+#include "lookup.h"
+#include "scales.h"
+
+/* three possible LSP to f curve functions; the exact computation
+ (float), a lookup based float implementation, and an integer
+ implementation. The float lookup is likely the optimal choice on
+ any machine with an FPU. The integer implementation is *not* fixed
+ point (due to the need for a large dynamic range and thus a
+ separately tracked exponent) and thus much more complex than the
+ relatively simple float implementations. It's mostly for future
+ work on a fully fixed point implementation for processors like the
+ ARM family. */
+
+/* define either of these (preferably FLOAT_LOOKUP) to have faster
+ but less precise implementation. */
+#undef FLOAT_LOOKUP
+#undef INT_LOOKUP
+
+#ifdef FLOAT_LOOKUP
+#include "lookup.c" /* catch this in the build system; we #include for
+ compilers (like gcc) that can't inline across
+ modules */
+
+/* side effect: changes *lsp to cosines of lsp */
+void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
+ float amp,float ampoffset){
+ int i;
+ float wdel=M_PI/ln;
+ vorbis_fpu_control fpu;
+
+ vorbis_fpu_setround(&fpu);
+ for(i=0;i<m;i++)lsp[i]=vorbis_coslook(lsp[i]);
+
+ i=0;
+ while(i<n){
+ int k=map[i];
+ int qexp;
+ float p=.7071067812f;
+ float q=.7071067812f;
+ float w=vorbis_coslook(wdel*k);
+ float *ftmp=lsp;
+ int c=m>>1;
+
+ while(c--){
+ q*=ftmp[0]-w;
+ p*=ftmp[1]-w;
+ ftmp+=2;
+ }
+
+ if(m&1){
+ /* odd order filter; slightly assymetric */
+ /* the last coefficient */
+ q*=ftmp[0]-w;
+ q*=q;
+ p*=p*(1.f-w*w);
+ }else{
+ /* even order filter; still symmetric */
+ q*=q*(1.f+w);
+ p*=p*(1.f-w);
+ }
+
+ q=frexp(p+q,&qexp);
+ q=vorbis_fromdBlook(amp*
+ vorbis_invsqlook(q)*
+ vorbis_invsq2explook(qexp+m)-
+ ampoffset);
+
+ do{
+ curve[i++]*=q;
+ }while(map[i]==k);
+ }
+ vorbis_fpu_restore(fpu);
+}
+
+#else
+
+#ifdef INT_LOOKUP
+#include "lookup.c" /* catch this in the build system; we #include for
+ compilers (like gcc) that can't inline across
+ modules */
+
+static const int MLOOP_1[64]={
+ 0,10,11,11, 12,12,12,12, 13,13,13,13, 13,13,13,13,
+ 14,14,14,14, 14,14,14,14, 14,14,14,14, 14,14,14,14,
+ 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
+ 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
+};
+
+static const int MLOOP_2[64]={
+ 0,4,5,5, 6,6,6,6, 7,7,7,7, 7,7,7,7,
+ 8,8,8,8, 8,8,8,8, 8,8,8,8, 8,8,8,8,
+ 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
+ 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
+};
+
+static const int MLOOP_3[8]={0,1,2,2,3,3,3,3};
+
+
+/* side effect: changes *lsp to cosines of lsp */
+void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
+ float amp,float ampoffset){
+
+ /* 0 <= m < 256 */
+
+ /* set up for using all int later */
+ int i;
+ int ampoffseti=rint(ampoffset*4096.f);
+ int ampi=rint(amp*16.f);
+ long *ilsp=alloca(m*sizeof(*ilsp));
+ for(i=0;i<m;i++)ilsp[i]=vorbis_coslook_i(lsp[i]/M_PI*65536.f+.5f);
+
+ i=0;
+ while(i<n){
+ int j,k=map[i];
+ unsigned long pi=46341; /* 2**-.5 in 0.16 */
+ unsigned long qi=46341;
+ int qexp=0,shift;
+ long wi=vorbis_coslook_i(k*65536/ln);
+
+ qi*=labs(ilsp[0]-wi);
+ pi*=labs(ilsp[1]-wi);
+
+ for(j=3;j<m;j+=2){
+ if(!(shift=MLOOP_1[(pi|qi)>>25]))
+ if(!(shift=MLOOP_2[(pi|qi)>>19]))
+ shift=MLOOP_3[(pi|qi)>>16];
+ qi=(qi>>shift)*labs(ilsp[j-1]-wi);
+ pi=(pi>>shift)*labs(ilsp[j]-wi);
+ qexp+=shift;
+ }
+ if(!(shift=MLOOP_1[(pi|qi)>>25]))
+ if(!(shift=MLOOP_2[(pi|qi)>>19]))
+ shift=MLOOP_3[(pi|qi)>>16];
+
+ /* pi,qi normalized collectively, both tracked using qexp */
+
+ if(m&1){
+ /* odd order filter; slightly assymetric */
+ /* the last coefficient */
+ qi=(qi>>shift)*labs(ilsp[j-1]-wi);
+ pi=(pi>>shift)<<14;
+ qexp+=shift;
+
+ if(!(shift=MLOOP_1[(pi|qi)>>25]))
+ if(!(shift=MLOOP_2[(pi|qi)>>19]))
+ shift=MLOOP_3[(pi|qi)>>16];
+
+ pi>>=shift;
+ qi>>=shift;
+ qexp+=shift-14*((m+1)>>1);
+
+ pi=((pi*pi)>>16);
+ qi=((qi*qi)>>16);
+ qexp=qexp*2+m;
+
+ pi*=(1<<14)-((wi*wi)>>14);
+ qi+=pi>>14;
+
+ }else{
+ /* even order filter; still symmetric */
+
+ /* p*=p(1-w), q*=q(1+w), let normalization drift because it isn't
+ worth tracking step by step */
+
+ pi>>=shift;
+ qi>>=shift;
+ qexp+=shift-7*m;
+
+ pi=((pi*pi)>>16);
+ qi=((qi*qi)>>16);
+ qexp=qexp*2+m;
+
+ pi*=(1<<14)-wi;
+ qi*=(1<<14)+wi;
+ qi=(qi+pi)>>14;
+
+ }
+
+
+ /* we've let the normalization drift because it wasn't important;
+ however, for the lookup, things must be normalized again. We
+ need at most one right shift or a number of left shifts */
+
+ if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
+ qi>>=1; qexp++;
+ }else
+ while(qi && !(qi&0x8000)){ /* checks for 0.0xxxxxxxxxxxxxxx or less*/
+ qi<<=1; qexp--;
+ }
+
+ amp=vorbis_fromdBlook_i(ampi* /* n.4 */
+ vorbis_invsqlook_i(qi,qexp)-
+ /* m.8, m+n<=8 */
+ ampoffseti); /* 8.12[0] */
+
+ curve[i]*=amp;
+ while(map[++i]==k)curve[i]*=amp;
+ }
+}
+
+#else
+
+/* old, nonoptimized but simple version for any poor sap who needs to
+ figure out what the hell this code does, or wants the other
+ fraction of a dB precision */