/* Copyright (c) 2007-2008 CSIRO

   Copyright (c) 2007-2009 Xiph.Org Foundation

   Copyright (c) 2008-2009 Gregory Maxwell

   Written by Jean-Marc Valin and Gregory Maxwell */

/*

   Redistribution and use in source and binary forms, with or without

   modification, are permitted provided that the following conditions

   are met:

   - Redistributions of source code must retain the above copyright

   notice, this list of conditions and the following disclaimer.

   - Redistributions in binary form must reproduce the above copyright

   notice, this list of conditions and the following disclaimer in the

   documentation and/or other materials provided with the distribution.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR

   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER

   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,

   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,

   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR

   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF

   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING

   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS

   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/

#ifdef HAVE_CONFIG_H

#include "config.h"

#endif

#include <math.h>

#include "bands.h"

#include "modes.h"

#include "vq.h"

#include "cwrs.h"

#include "stack_alloc.h"

#include "os_support.h"

#include "mathops.h"

#include "rate.h"

#include "quant_bands.h"

#include "pitch.h"

int hysteresis_decision(opus_val16 val, const opus_val16 *thresholds, const opus_val16 *hysteresis, int N, int prev)

{

   int i;

   for (i=0;i<N;i++)

   {

      if (val < thresholds[i])

         break;

   }

   if (i>prev && val < thresholds[prev]+hysteresis[prev])

      i=prev;

   if (i<prev && val > thresholds[prev-1]-hysteresis[prev-1])

      i=prev;

   return i;

}

opus_uint32 celt_lcg_rand(opus_uint32 seed)

{

   return 1664525 * seed + 1013904223;

}

/* This is a cos() approximation designed to be bit-exact on any platform. Bit exactness

   with this approximation is important because it has an impact on the bit allocation */

opus_int16 bitexact_cos(opus_int16 x)

{

   opus_int32 tmp;

   opus_int16 x2;

   tmp = (4096+((opus_int32)(x)*(x)))>>13;

   celt_sig_assert(tmp<=32767);

   x2 = tmp;

   x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));

   celt_sig_assert(x2<=32766);

   return 1+x2;

}

int bitexact_log2tan(int isin,int icos)

{

   int lc;

   int ls;

   lc=EC_ILOG(icos);

   ls=EC_ILOG(isin);

   icos<<=15-lc;

   isin<<=15-ls;

   return (ls-lc)*(1<<11)

         +FRAC_MUL16(isin, FRAC_MUL16(isin, -2597) + 7932)

         -FRAC_MUL16(icos, FRAC_MUL16(icos, -2597) + 7932);

}

#ifdef FIXED_POINT

/* Compute the amplitude (sqrt energy) in each of the bands */

void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int LM, int arch)

{

   int i, c, N;

   const opus_int16 *eBands = m->eBands;

   (void)arch;

   N = m->shortMdctSize<<LM;

   c=0; do {

      for (i=0;i<end;i++)

      {

         int j;

         opus_val32 maxval=0;

         opus_val32 sum = 0;

         maxval = celt_maxabs32(&X[c*N+(eBands[i]<<LM)], (eBands[i+1]-eBands[i])<<LM);

         if (maxval > 0)

         {

            int shift = celt_ilog2(maxval) - 14 + (((m->logN[i]>>BITRES)+LM+1)>>1);

            j=eBands[i]<<LM;

            if (shift>0)

            {

               do {

                  sum = MAC16_16(sum, EXTRACT16(SHR32(X[j+c*N],shift)),

                        EXTRACT16(SHR32(X[j+c*N],shift)));

               } while (++j<eBands[i+1]<<LM);

            } else {

               do {

                  sum = MAC16_16(sum, EXTRACT16(SHL32(X[j+c*N],-shift)),

                        EXTRACT16(SHL32(X[j+c*N],-shift)));

               } while (++j<eBands[i+1]<<LM);

            }

            /* We're adding one here to ensure the normalized band isn't larger than unity norm */

            bandE[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);

         } else {

            bandE[i+c*m->nbEBands] = EPSILON;

         }

         /*printf ("%f ", bandE[i+c*m->nbEBands]);*/

      }

   } while (++c<C);

   /*printf ("\n");*/

}

/* Normalise each band such that the energy is one. */

void normalise_bands(const CELTMode *m, const celt_sig * OPUS_RESTRICT freq, celt_norm * OPUS_RESTRICT X, const celt_ener *bandE, int end, int C, int M)

{

   int i, c, N;

   const opus_int16 *eBands = m->eBands;

   N = M*m->shortMdctSize;

   c=0; do {

      i=0; do {

         opus_val16 g;

         int j,shift;

         opus_val16 E;

         shift = celt_zlog2(bandE[i+c*m->nbEBands])-13;

         E = VSHR32(bandE[i+c*m->nbEBands], shift);

         g = EXTRACT16(celt_rcp(SHL32(E,3)));

         j=M*eBands[i]; do {

            X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);

         } while (++j<M*eBands[i+1]);

      } while (++i<end);

   } while (++c<C);

}

#else /* FIXED_POINT */

/* Compute the amplitude (sqrt energy) in each of the bands */

void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int LM, int arch)

{

   int i, c, N;

   const opus_int16 *eBands = m->eBands;

   N = m->shortMdctSize<<LM;

   c=0; do {

      for (i=0;i<end;i++)

      {

         opus_val32 sum;

         sum = 1e-27f + celt_inner_prod(&X[c*N+(eBands[i]<<LM)], &X[c*N+(eBands[i]<<LM)], (eBands[i+1]-eBands[i])<<LM, arch);

         bandE[i+c*m->nbEBands] = celt_sqrt(sum);

         /*printf ("%f ", bandE[i+c*m->nbEBands]);*/

      }

   } while (++c<C);

   /*printf ("\n");*/

}

/* Normalise each band such that the energy is one. */

void normalise_bands(const CELTMode *m, const celt_sig * OPUS_RESTRICT freq, celt_norm * OPUS_RESTRICT X, const celt_ener *bandE, int end, int C, int M)

{

   int i, c, N;

   const opus_int16 *eBands = m->eBands;

   N = M*m->shortMdctSize;

   c=0; do {

      for (i=0;i<end;i++)

      {

         int j;

         opus_val16 g = 1.f/(1e-27f+bandE[i+c*m->nbEBands]);

         for (j=M*eBands[i];j<M*eBands[i+1];j++)

            X[j+c*N] = freq[j+c*N]*g;

      }

   } while (++c<C);

}

#endif /* FIXED_POINT */

/* De-normalise the energy to produce the synthesis from the unit-energy bands */

void denormalise_bands(const CELTMode *m, const celt_norm * OPUS_RESTRICT X,

      celt_sig * OPUS_RESTRICT freq, const opus_val16 *bandLogE, int start,

      int end, int M, int downsample, int silence)

{

   int i, N;

   int bound;

   celt_sig * OPUS_RESTRICT f;

   const celt_norm * OPUS_RESTRICT x;

   const opus_int16 *eBands = m->eBands;

   N = M*m->shortMdctSize;

   bound = M*eBands[end];

   if (downsample!=1)

      bound = IMIN(bound, N/downsample);

   if (silence)

   {

      bound = 0;

      start = end = 0;

   }

   f = freq;

   x = X+M*eBands[start];

   for (i=0;i<M*eBands[start];i++)

      *f++ = 0;

   for (i=start;i<end;i++)

   {

      int j, band_end;

      opus_val16 g;

      opus_val16 lg;

#ifdef FIXED_POINT

      int shift;

#endif

      j=M*eBands[i];

      band_end = M*eBands[i+1];

      lg = SATURATE16(ADD32(bandLogE[i], SHL32((opus_val32)eMeans[i],6)));

#ifndef FIXED_POINT

      g = celt_exp2(MIN32(32.f, lg));

#else

      /* Handle the integer part of the log energy */

      shift = 16-(lg>>DB_SHIFT);

      if (shift>31)

      {

         shift=0;

         g=0;

      } else {

         /* Handle the fractional part. */

         g = celt_exp2_frac(lg&((1<<DB_SHIFT)-1));

      }

      /* Handle extreme gains with negative shift. */

      if (shift<0)

      {

         /* For shift <= -2 and g > 16384 we'd be likely to overflow, so we're

            capping the gain here, which is equivalent to a cap of 18 on lg.

            This shouldn't trigger unless the bitstream is already corrupted. */

         if (shift <= -2)

         {

            g = 16384;

            shift = -2;

         }

         do {

            *f++ = SHL32(MULT16_16(*x++, g), -shift);

         } while (++j<band_end);

      } else

#endif

         /* Be careful of the fixed-point "else" just above when changing this code */

         do {

            *f++ = SHR32(MULT16_16(*x++, g), shift);

         } while (++j<band_end);

   }

   celt_assert(start <= end);

   OPUS_CLEAR(&freq[bound], N-bound);

}

/* This prevents energy collapse for transients with multiple short MDCTs */

void anti_collapse(const CELTMode *m, celt_norm *X_, unsigned char *collapse_masks, int LM, int C, int size,

      int start, int end, const opus_val16 *logE, const opus_val16 *prev1logE,

      const opus_val16 *prev2logE, const int *pulses, opus_uint32 seed, int arch)

{

   int c, i, j, k;

   for (i=start;i<end;i++)

   {

      int N0;

      opus_val16 thresh, sqrt_1;

      int depth;

#ifdef FIXED_POINT

      int shift;

      opus_val32 thresh32;

#endif

      N0 = m->eBands[i+1]-m->eBands[i];

      /* depth in 1/8 bits */

      celt_sig_assert(pulses[i]>=0);

      depth = celt_udiv(1+pulses[i], (m->eBands[i+1]-m->eBands[i]))>>LM;

#ifdef FIXED_POINT

      thresh32 = SHR32(celt_exp2(-SHL16(depth, 10-BITRES)),1);

      thresh = MULT16_32_Q15(QCONST16(0.5f, 15), MIN32(32767,thresh32));

      {

         opus_val32 t;

         t = N0<<LM;

         shift = celt_ilog2(t)>>1;

         t = SHL32(t, (7-shift)<<1);

         sqrt_1 = celt_rsqrt_norm(t);

      }

#else

      thresh = .5f*celt_exp2(-.125f*depth);

      sqrt_1 = celt_rsqrt(N0<<LM);

#endif

      c=0; do

      {

         celt_norm *X;

         opus_val16 prev1;

         opus_val16 prev2;

         opus_val32 Ediff;

         opus_val16 r;

         int renormalize=0;

         prev1 = prev1logE[c*m->nbEBands+i];

         prev2 = prev2logE[c*m->nbEBands+i];

         if (C==1)

         {

            prev1 = MAX16(prev1,prev1logE[m->nbEBands+i]);

            prev2 = MAX16(prev2,prev2logE[m->nbEBands+i]);

         }

         Ediff = EXTEND32(logE[c*m->nbEBands+i])-EXTEND32(MIN16(prev1,prev2));

         Ediff = MAX32(0, Ediff);

#ifdef FIXED_POINT

         if (Ediff < 16384)

         {

            opus_val32 r32 = SHR32(celt_exp2(-EXTRACT16(Ediff)),1);

            r = 2*MIN16(16383,r32);

         } else {

            r = 0;

         }

         if (LM==3)

            r = MULT16_16_Q14(23170, MIN32(23169, r));

         r = SHR16(MIN16(thresh, r),1);

         r = SHR32(MULT16_16_Q15(sqrt_1, r),shift);

#else

         /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because

            short blocks don't have the same energy as long */

         r = 2.f*celt_exp2(-Ediff);

         if (LM==3)

            r *= 1.41421356f;

         r = MIN16(thresh, r);

         r = r*sqrt_1;

#endif

         X = X_+c*size+(m->eBands[i]<<LM);

         for (k=0;k<1<<LM;k++)

         {

            /* Detect collapse */

            if (!(collapse_masks[i*C+c]&1<<k))

            {

               /* Fill with noise */

               for (j=0;j<N0;j++)

               {

                  seed = celt_lcg_rand(seed);

                  X[(j<<LM)+k] = (seed&0x8000 ? r : -r);

               }

               renormalize = 1;

            }

         }

         /* We just added some energy, so we need to renormalise */

         if (renormalize)

            renormalise_vector(X, N0<<LM, Q15ONE, arch);

      } while (++c<C);

   }

}

/* Compute the weights to use for optimizing normalized distortion across

   channels. We use the amplitude to weight square distortion, which means

   that we use the square root of the value we would have been using if we

   wanted to minimize the MSE in the non-normalized domain. This roughly

   corresponds to some quick-and-dirty perceptual experiments I ran to

   measure inter-aural masking (there doesn't seem to be any published data

   on the topic). */

static void compute_channel_weights(celt_ener Ex, celt_ener Ey, opus_val16 w[2])

{

   celt_ener minE;

#ifdef FIXED_POINT

   int shift;

#endif

   minE = MIN32(Ex, Ey);

   /* Adjustment to make the weights a bit more conservative. */

   Ex = ADD32(Ex, minE/3);

   Ey = ADD32(Ey, minE/3);

#ifdef FIXED_POINT

   shift = celt_ilog2(EPSILON+MAX32(Ex, Ey))-14;

#endif

   w[0] = VSHR32(Ex, shift);

   w[1] = VSHR32(Ey, shift);

}

static void intensity_stereo(const CELTMode *m, celt_norm * OPUS_RESTRICT X, const celt_norm * OPUS_RESTRICT Y, const celt_ener *bandE, int bandID, int N)

{

   int i = bandID;

   int j;

   opus_val16 a1, a2;

   opus_val16 left, right;

   opus_val16 norm;

#ifdef FIXED_POINT

   int shift = celt_zlog2(MAX32(bandE[i], bandE[i+m->nbEBands]))-13;

#endif

   left = VSHR32(bandE[i],shift);

   right = VSHR32(bandE[i+m->nbEBands],shift);

   norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));

   a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);

   a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);

   for (j=0;j<N;j++)

   {

      celt_norm r, l;

      l = X[j];

      r = Y[j];

      X[j] = EXTRACT16(SHR32(MAC16_16(MULT16_16(a1, l), a2, r), 14));

      /* Side is not encoded, no need to calculate */

   }

}

static void stereo_split(celt_norm * OPUS_RESTRICT X, celt_norm * OPUS_RESTRICT Y, int N)

{

   int j;

   for (j=0;j<N;j++)

   {

      opus_val32 r, l;

      l = MULT16_16(QCONST16(.70710678f, 15), X[j]);

      r = MULT16_16(QCONST16(.70710678f, 15), Y[j]);

      X[j] = EXTRACT16(SHR32(ADD32(l, r), 15));

      Y[j] = EXTRACT16(SHR32(SUB32(r, l), 15));

   }

}

static void stereo_merge(celt_norm * OPUS_RESTRICT X, celt_norm * OPUS_RESTRICT Y, opus_val16 mid, int N, int arch)

{

   int j;

   opus_val32 xp=0, side=0;

   opus_val32 El, Er;

   opus_val16 mid2;

#ifdef FIXED_POINT

   int kl, kr;

#endif

   opus_val32 t, lgain, rgain;

   /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */

   dual_inner_prod(Y, X, Y, N, &xp, &side, arch);

   /* Compensating for the mid normalization */

   xp = MULT16_32_Q15(mid, xp);

   /* mid and side are in Q15, not Q14 like X and Y */

   mid2 = SHR16(mid, 1);

   El = MULT16_16(mid2, mid2) + side - 2*xp;

   Er = MULT16_16(mid2, mid2) + side + 2*xp;

   if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28))

   {

      OPUS_COPY(Y, X, N);

      return;

   }

#ifdef FIXED_POINT

   kl = celt_ilog2(El)>>1;

   kr = celt_ilog2(Er)>>1;

#endif

   t = VSHR32(El, (kl-7)<<1);

   lgain = celt_rsqrt_norm(t);

   t = VSHR32(Er, (kr-7)<<1);

   rgain = celt_rsqrt_norm(t);

#ifdef FIXED_POINT

   if (kl < 7)

      kl = 7;

   if (kr < 7)

      kr = 7;

#endif

   for (j=0;j<N;j++)

   {

      celt_norm r, l;

      /* Apply mid scaling (side is already scaled) */

      l = MULT16_16_P15(mid, X[j]);

      r = Y[j];

      X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));

      Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));

   }

}

/* Decide whether we should spread the pulses in the current frame */

int spreading_decision(const CELTMode *m, const celt_norm *X, int *average,

      int last_decision, int *hf_average, int *tapset_decision, int update_hf,

      int end, int C, int M, const int *spread_weight)

{

   int i, c, N0;

   int sum = 0, nbBands=0;

   const opus_int16 * OPUS_RESTRICT eBands = m->eBands;

   int decision;

   int hf_sum=0;

   celt_assert(end>0);

   N0 = M*m->shortMdctSize;

   if (M*(eBands[end]-eBands[end-1]) <= 8)

      return SPREAD_NONE;

   c=0; do {

      for (i=0;i<end;i++)

      {

         int j, N, tmp=0;

         int tcount[3] = {0,0,0};

         const celt_norm * OPUS_RESTRICT x = X+M*eBands[i]+c*N0;

         N = M*(eBands[i+1]-eBands[i]);

         if (N<=8)

            continue;

         /* Compute rough CDF of |x[j]| */

         for (j=0;j<N;j++)

         {

            opus_val32 x2N; /* Q13 */

            x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);

            if (x2N < QCONST16(0.25f,13))

               tcount[0]++;

            if (x2N < QCONST16(0.0625f,13))

               tcount[1]++;

            if (x2N < QCONST16(0.015625f,13))

               tcount[2]++;

         }

         /* Only include four last bands (8 kHz and up) */

         if (i>m->nbEBands-4)

            hf_sum += celt_udiv(32*(tcount[1]+tcount[0]), N);

         tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);

         sum += tmp*spread_weight[i];

         nbBands+=spread_weight[i];

      }

   } while (++c<C);

   if (update_hf)

   {

      if (hf_sum)

         hf_sum = celt_udiv(hf_sum, C*(4-m->nbEBands+end));

      *hf_average = (*hf_average+hf_sum)>>1;

      hf_sum = *hf_average;

      if (*tapset_decision==2)

         hf_sum += 4;

      else if (*tapset_decision==0)

         hf_sum -= 4;

      if (hf_sum > 22)

         *tapset_decision=2;

      else if (hf_sum > 18)

         *tapset_decision=1;

      else

         *tapset_decision=0;

   }

   /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/

   celt_assert(nbBands>0); /* end has to be non-zero */

   celt_assert(sum>=0);

   sum = celt_udiv((opus_int32)sum<<8, nbBands);

   /* Recursive averaging */

   sum = (sum+*average)>>1;

   *average = sum;

   /* Hysteresis */

   sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2;

   if (sum < 80)

   {

      decision = SPREAD_AGGRESSIVE;

   } else if (sum < 256)

   {

      decision = SPREAD_NORMAL;

   } else if (sum < 384)

   {

      decision = SPREAD_LIGHT;

   } else {

      decision = SPREAD_NONE;

   }

#ifdef FUZZING

   decision = rand()&0x3;

   *tapset_decision=rand()%3;

#endif

   return decision;

}

/* Indexing table for converting from natural Hadamard to ordery Hadamard

   This is essentially a bit-reversed Gray, on top of which we've added

   an inversion of the order because we want the DC at the end rather than

   the beginning. The lines are for N=2, 4, 8, 16 */

static const int ordery_table[] = {

       1,  0,

       3,  0,  2,  1,

       7,  0,  4,  3,  6,  1,  5,  2,

      15,  0,  8,  7, 12,  3, 11,  4, 14,  1,  9,  6, 13,  2, 10,  5,

};

static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)

{

   int i,j;

   VARDECL(celt_norm, tmp);

   int N;

   SAVE_STACK;

   N = N0*stride;

   ALLOC(tmp, N, celt_norm);

   celt_assert(stride>0);

   if (hadamard)

   {

      const int *ordery = ordery_table+stride-2;

      for (i=0;i<stride;i++)

      {

         for (j=0;j<N0;j++)

            tmp[ordery[i]*N0+j] = X[j*stride+i];

      }

   } else {

      for (i=0;i<stride;i++)

         for (j=0;j<N0;j++)

            tmp[i*N0+j] = X[j*stride+i];

   }

   OPUS_COPY(X, tmp, N);

   RESTORE_STACK;

}

static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)

{

   int i,j;

   VARDECL(celt_norm, tmp);

   int N;

   SAVE_STACK;

   N = N0*stride;

   ALLOC(tmp, N, celt_norm);

   if (hadamard)

   {

      const int *ordery = ordery_table+stride-2;

      for (i=0;i<stride;i++)

         for (j=0;j<N0;j++)

            tmp[j*stride+i] = X[ordery[i]*N0+j];

   } else {

      for (i=0;i<stride;i++)

         for (j=0;j<N0;j++)

            tmp[j*stride+i] = X[i*N0+j];

   }

   OPUS_COPY(X, tmp, N);

   RESTORE_STACK;

}

void haar1(celt_norm *X, int N0, int stride)

{

   int i, j;

   N0 >>= 1;

   for (i=0;i<stride;i++)

      for (j=0;j<N0;j++)

      {

         opus_val32 tmp1, tmp2;

         tmp1 = MULT16_16(QCONST16(.70710678f,15), X[stride*2*j+i]);

         tmp2 = MULT16_16(QCONST16(.70710678f,15), X[stride*(2*j+1)+i]);

         X[stride*2*j+i] = EXTRACT16(PSHR32(ADD32(tmp1, tmp2), 15));

         X[stride*(2*j+1)+i] = EXTRACT16(PSHR32(SUB32(tmp1, tmp2), 15));

      }

}

static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo)

{

   static const opus_int16 exp2_table8[8] =

      {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};

   int qn, qb;

   int N2 = 2*N-1;

   if (stereo && N==2)

      N2--;

   /* The upper limit ensures that in a stereo split with itheta==16384, we'll

       always have enough bits left over to code at least one pulse in the

       side; otherwise it would collapse, since it doesn't get folded. */

   qb = celt_sudiv(b+N2*offset, N2);

   qb = IMIN(b-pulse_cap-(4<<BITRES), qb);

   qb = IMIN(8<<BITRES, qb);

   if (qb<(1<<BITRES>>1)) {

      qn = 1;

   } else {

      qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));

      qn = (qn+1)>>1<<1;

   }

   celt_assert(qn <= 256);

   return qn;

}

struct band_ctx {

   int encode;

   int resynth;

   const CELTMode *m;

   int i;

   int intensity;

   int spread;

   int tf_change;

   ec_ctx *ec;

   opus_int32 remaining_bits;

   const celt_ener *bandE;

   opus_uint32 seed;

   int arch;

   int theta_round;

   int disable_inv;

   int avoid_split_noise;

};

struct split_ctx {

   int inv;

   int imid;

   int iside;

   int delta;

   int itheta;

   int qalloc;

};

static void compute_theta(struct band_ctx *ctx, struct split_ctx *sctx,

      celt_norm *X, celt_norm *Y, int N, int *b, int B, int B0,

      int LM,

      int stereo, int *fill)

{

   int qn;

   int itheta=0;

   int delta;

   int imid, iside;

   int qalloc;

   int pulse_cap;

   int offset;

   opus_int32 tell;

   int inv=0;

   int encode;

   const CELTMode *m;

   int i;

   int intensity;

   ec_ctx *ec;

   const celt_ener *bandE;

   encode = ctx->encode;

   m = ctx->m;

   i = ctx->i;

   intensity = ctx->intensity;

   ec = ctx->ec;

   bandE = ctx->bandE;

   /* Decide on the resolution to give to the split parameter theta */

   pulse_cap = m->logN[i]+LM*(1<<BITRES);

   offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);

   qn = compute_qn(N, *b, offset, pulse_cap, stereo);

   if (stereo && i>=intensity)

      qn = 1;

   if (encode)

   {

      /* theta is the atan() of the ratio between the (normalized)

         side and mid. With just that parameter, we can re-scale both

         mid and side because we know that 1) they have unit norm and

         2) they are orthogonal. */

      itheta = stereo_itheta(X, Y, stereo, N, ctx->arch);

   }

   tell = ec_tell_frac(ec);

   if (qn!=1)

   {

      if (encode)

      {

         if (!stereo || ctx->theta_round == 0)

         {

            itheta = (itheta*(opus_int32)qn+8192)>>14;

            if (!stereo && ctx->avoid_split_noise && itheta > 0 && itheta < qn)

            {

               /* Check if the selected value of theta will cause the bit allocation

                  to inject noise on one side. If so, make sure the energy of that side

                  is zero. */

               int unquantized = celt_udiv((opus_int32)itheta*16384, qn);

               imid = bitexact_cos((opus_int16)unquantized);

               iside = bitexact_cos((opus_int16)(16384-unquantized));

               delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid));

               if (delta > *b)

                  itheta = qn;

               else if (delta < -*b)

                  itheta = 0;

            }

         } else {

            int down;

            /* Bias quantization towards itheta=0 and itheta=16384. */

            int bias = itheta > 8192 ? 32767/qn : -32767/qn;

            down = IMIN(qn-1, IMAX(0, (itheta*(opus_int32)qn + bias)>>14));

            if (ctx->theta_round < 0)

               itheta = down;

            else

               itheta = down+1;

         }

      }

      /* Entropy coding of the angle. We use a uniform pdf for the

         time split, a step for stereo, and a triangular one for the rest. */

      if (stereo && N>2)

      {

         int p0 = 3;

         int x = itheta;

         int x0 = qn/2;

         int ft = p0*(x0+1) + x0;

         /* Use a probability of p0 up to itheta=8192 and then use 1 after */

         if (encode)

         {

            ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);

         } else {

            int fs;

            fs=ec_decode(ec,ft);

            if (fs<(x0+1)*p0)

               x=fs/p0;

            else

               x=x0+1+(fs-(x0+1)*p0);

            ec_dec_update(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);

            itheta = x;

         }

      } else if (B0>1 || stereo) {

         /* Uniform pdf */

         if (encode)

            ec_enc_uint(ec, itheta, qn+1);

         else

            itheta = ec_dec_uint(ec, qn+1);

      } else {

         int fs=1, ft;

         ft = ((qn>>1)+1)*((qn>>1)+1);

         if (encode)

         {

            int fl;

            fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;

            fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :

             ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);

            ec_encode(ec, fl, fl+fs, ft);

         } else {

            /* Triangular pdf */

            int fl=0;

            int fm;

            fm = ec_decode(ec, ft);

            if (fm < ((qn>>1)*((qn>>1) + 1)>>1))

            {

               itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;

               fs = itheta + 1;

               fl = itheta*(itheta + 1)>>1;

            }

            else

            {

               itheta = (2*(qn + 1)

                - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;

               fs = qn + 1 - itheta;

               fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);

            }

            ec_dec_update(ec, fl, fl+fs, ft);

         }

      }

      celt_assert(itheta>=0);

      itheta = celt_udiv((opus_int32)itheta*16384, qn);

      if (encode && stereo)

      {

         if (itheta==0)

            intensity_stereo(m, X, Y, bandE, i, N);

         else

            stereo_split(X, Y, N);

      }

      /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.

               Let's do that at higher complexity */

   } else if (stereo) {

      if (encode)

      {

         inv = itheta > 8192 && !ctx->disable_inv;

         if (inv)

         {

            int j;

            for (j=0;j<N;j++)

               Y[j] = -Y[j];

         }

         intensity_stereo(m, X, Y, bandE, i, N);

      }

      if (*b>2<<BITRES && ctx->remaining_bits > 2<<BITRES)

      {

         if (encode)

            ec_enc_bit_logp(ec, inv, 2);

         else

            inv = ec_dec_bit_logp(ec, 2);

      } else

         inv = 0;

      /* inv flag override to avoid problems with downmixing. */

      if (ctx->disable_inv)

         inv = 0;

      itheta = 0;

   }

   qalloc = ec_tell_frac(ec) - tell;

   *b -= qalloc;

   if (itheta == 0)

   {

      imid = 32767;

      iside = 0;

      *fill &= (1<<B)-1;

      delta = -16384;

   } else if (itheta == 16384)

   {

      imid = 0;

      iside = 32767;

      *fill &= ((1<<B)-1)<<B;

      delta = 16384;

   } else {

      imid = bitexact_cos((opus_int16)itheta);