tor-browser

The Tor Browser
git clone https://git.dasho.dev/tor-browser.git
Log | Files | Refs | README | LICENSE

vorbis_smallft.c (22185B)

      1 /********************************************************************
      2 *                                                                  *
      3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
      4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
      5 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
      6 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
      7 *                                                                  *
      8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
      9 * by the Xiph.Org Foundation https://xiph.org/                     *
     10 *                                                                  *
     11 ********************************************************************
     12 
     13 function: *unnormalized* fft transform
     14 
     15 ********************************************************************/
     16 
     17 /* FFT implementation from OggSquish, minus cosine transforms,
     18 * minus all but radix 2/4 case.  In Vorbis we only need this
     19 * cut-down version.
     20 *
     21 * To do more than just power-of-two sized vectors, see the full
     22 * version I wrote for NetLib.
     23 *
     24 * Note that the packing is a little strange; rather than the FFT r/i
     25 * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1,
     26 * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
     27 * FORTRAN version
     28 */
     29 
     30 #include <stdlib.h>
     31 #include <string.h>
     32 #include <math.h>
     33 #include "smallft.h"
     34 #include "os.h"
     35 #include "misc.h"
     36 
     37 static void drfti1(int n, float *wa, int *ifac){
     38  static int ntryh[4] = { 4,2,3,5 };
     39  static float tpi = 6.28318530717958648f;
     40  float arg,argh,argld,fi;
     41  int ntry=0,i,j=-1;
     42  int k1, l1, l2, ib;
     43  int ld, ii, ip, is, nq, nr;
     44  int ido, ipm, nfm1;
     45  int nl=n;
     46  int nf=0;
     47 
     48 L101:
     49  j++;
     50  if (j < 4)
     51    ntry=ntryh[j];
     52  else
     53    ntry+=2;
     54 
     55 L104:
     56  nq=nl/ntry;
     57  nr=nl-ntry*nq;
     58  if (nr!=0) goto L101;
     59 
     60  nf++;
     61  ifac[nf+1]=ntry;
     62  nl=nq;
     63  if(ntry!=2)goto L107;
     64  if(nf==1)goto L107;
     65 
     66  for (i=1;i<nf;i++){
     67    ib=nf-i+1;
     68    ifac[ib+1]=ifac[ib];
     69  }
     70  ifac[2] = 2;
     71 
     72 L107:
     73  if(nl!=1)goto L104;
     74  ifac[0]=n;
     75  ifac[1]=nf;
     76  argh=tpi/n;
     77  is=0;
     78  nfm1=nf-1;
     79  l1=1;
     80 
     81  if(nfm1==0)return;
     82 
     83  for (k1=0;k1<nfm1;k1++){
     84    ip=ifac[k1+2];
     85    ld=0;
     86    l2=l1*ip;
     87    ido=n/l2;
     88    ipm=ip-1;
     89 
     90    for (j=0;j<ipm;j++){
     91      ld+=l1;
     92      i=is;
     93      argld=(float)ld*argh;
     94      fi=0.f;
     95      for (ii=2;ii<ido;ii+=2){
     96        fi+=1.f;
     97        arg=fi*argld;
     98        wa[i++]=cos(arg);
     99        wa[i++]=sin(arg);
    100      }
    101      is+=ido;
    102    }
    103    l1=l2;
    104  }
    105 }
    106 
    107 static void fdrffti(int n, float *wsave, int *ifac){
    108 
    109  if (n == 1) return;
    110  drfti1(n, wsave+n, ifac);
    111 }
    112 
    113 static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){
    114  int i,k;
    115  float ti2,tr2;
    116  int t0,t1,t2,t3,t4,t5,t6;
    117 
    118  t1=0;
    119  t0=(t2=l1*ido);
    120  t3=ido<<1;
    121  for(k=0;k<l1;k++){
    122    ch[t1<<1]=cc[t1]+cc[t2];
    123    ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
    124    t1+=ido;
    125    t2+=ido;
    126  }
    127 
    128  if(ido<2)return;
    129  if(ido==2)goto L105;
    130 
    131  t1=0;
    132  t2=t0;
    133  for(k=0;k<l1;k++){
    134    t3=t2;
    135    t4=(t1<<1)+(ido<<1);
    136    t5=t1;
    137    t6=t1+t1;
    138    for(i=2;i<ido;i+=2){
    139      t3+=2;
    140      t4-=2;
    141      t5+=2;
    142      t6+=2;
    143      tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
    144      ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
    145      ch[t6]=cc[t5]+ti2;
    146      ch[t4]=ti2-cc[t5];
    147      ch[t6-1]=cc[t5-1]+tr2;
    148      ch[t4-1]=cc[t5-1]-tr2;
    149    }
    150    t1+=ido;
    151    t2+=ido;
    152  }
    153 
    154  if(ido%2==1)return;
    155 
    156 L105:
    157  t3=(t2=(t1=ido)-1);
    158  t2+=t0;
    159  for(k=0;k<l1;k++){
    160    ch[t1]=-cc[t2];
    161    ch[t1-1]=cc[t3];
    162    t1+=ido<<1;
    163    t2+=ido;
    164    t3+=ido;
    165  }
    166 }
    167 
    168 static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1,
    169            float *wa2,float *wa3){
    170  static float hsqt2 = .70710678118654752f;
    171  int i,k,t0,t1,t2,t3,t4,t5,t6;
    172  float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
    173  t0=l1*ido;
    174 
    175  t1=t0;
    176  t4=t1<<1;
    177  t2=t1+(t1<<1);
    178  t3=0;
    179 
    180  for(k=0;k<l1;k++){
    181    tr1=cc[t1]+cc[t2];
    182    tr2=cc[t3]+cc[t4];
    183 
    184    ch[t5=t3<<2]=tr1+tr2;
    185    ch[(ido<<2)+t5-1]=tr2-tr1;
    186    ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
    187    ch[t5]=cc[t2]-cc[t1];
    188 
    189    t1+=ido;
    190    t2+=ido;
    191    t3+=ido;
    192    t4+=ido;
    193  }
    194 
    195  if(ido<2)return;
    196  if(ido==2)goto L105;
    197 
    198 
    199  t1=0;
    200  for(k=0;k<l1;k++){
    201    t2=t1;
    202    t4=t1<<2;
    203    t5=(t6=ido<<1)+t4;
    204    for(i=2;i<ido;i+=2){
    205      t3=(t2+=2);
    206      t4+=2;
    207      t5-=2;
    208 
    209      t3+=t0;
    210      cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
    211      ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
    212      t3+=t0;
    213      cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];
    214      ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];
    215      t3+=t0;
    216      cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];
    217      ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];
    218 
    219      tr1=cr2+cr4;
    220      tr4=cr4-cr2;
    221      ti1=ci2+ci4;
    222      ti4=ci2-ci4;
    223 
    224      ti2=cc[t2]+ci3;
    225      ti3=cc[t2]-ci3;
    226      tr2=cc[t2-1]+cr3;
    227      tr3=cc[t2-1]-cr3;
    228 
    229      ch[t4-1]=tr1+tr2;
    230      ch[t4]=ti1+ti2;
    231 
    232      ch[t5-1]=tr3-ti4;
    233      ch[t5]=tr4-ti3;
    234 
    235      ch[t4+t6-1]=ti4+tr3;
    236      ch[t4+t6]=tr4+ti3;
    237 
    238      ch[t5+t6-1]=tr2-tr1;
    239      ch[t5+t6]=ti1-ti2;
    240    }
    241    t1+=ido;
    242  }
    243  if(ido&1)return;
    244 
    245 L105:
    246 
    247  t2=(t1=t0+ido-1)+(t0<<1);
    248  t3=ido<<2;
    249  t4=ido;
    250  t5=ido<<1;
    251  t6=ido;
    252 
    253  for(k=0;k<l1;k++){
    254    ti1=-hsqt2*(cc[t1]+cc[t2]);
    255    tr1=hsqt2*(cc[t1]-cc[t2]);
    256 
    257    ch[t4-1]=tr1+cc[t6-1];
    258    ch[t4+t5-1]=cc[t6-1]-tr1;
    259 
    260    ch[t4]=ti1-cc[t1+t0];
    261    ch[t4+t5]=ti1+cc[t1+t0];
    262 
    263    t1+=ido;
    264    t2+=ido;
    265    t4+=t3;
    266    t6+=ido;
    267  }
    268 }
    269 
    270 static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
    271                          float *c2,float *ch,float *ch2,float *wa){
    272 
    273  static float tpi=6.283185307179586f;
    274  int idij,ipph,i,j,k,l,ic,ik,is;
    275  int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
    276  float dc2,ai1,ai2,ar1,ar2,ds2;
    277  int nbd;
    278  float dcp,arg,dsp,ar1h,ar2h;
    279  int idp2,ipp2;
    280 
    281  arg=tpi/(float)ip;
    282  dcp=cos(arg);
    283  dsp=sin(arg);
    284  ipph=(ip+1)>>1;
    285  ipp2=ip;
    286  idp2=ido;
    287  nbd=(ido-1)>>1;
    288  t0=l1*ido;
    289  t10=ip*ido;
    290 
    291  if(ido==1)goto L119;
    292  for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
    293 
    294  t1=0;
    295  for(j=1;j<ip;j++){
    296    t1+=t0;
    297    t2=t1;
    298    for(k=0;k<l1;k++){
    299      ch[t2]=c1[t2];
    300      t2+=ido;
    301    }
    302  }
    303 
    304  is=-ido;
    305  t1=0;
    306  if(nbd>l1){
    307    for(j=1;j<ip;j++){
    308      t1+=t0;
    309      is+=ido;
    310      t2= -ido+t1;
    311      for(k=0;k<l1;k++){
    312        idij=is-1;
    313        t2+=ido;
    314        t3=t2;
    315        for(i=2;i<ido;i+=2){
    316          idij+=2;
    317          t3+=2;
    318          ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
    319          ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
    320        }
    321      }
    322    }
    323  }else{
    324 
    325    for(j=1;j<ip;j++){
    326      is+=ido;
    327      idij=is-1;
    328      t1+=t0;
    329      t2=t1;
    330      for(i=2;i<ido;i+=2){
    331        idij+=2;
    332        t2+=2;
    333        t3=t2;
    334        for(k=0;k<l1;k++){
    335          ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
    336          ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
    337          t3+=ido;
    338        }
    339      }
    340    }
    341  }
    342 
    343  t1=0;
    344  t2=ipp2*t0;
    345  if(nbd<l1){
    346    for(j=1;j<ipph;j++){
    347      t1+=t0;
    348      t2-=t0;
    349      t3=t1;
    350      t4=t2;
    351      for(i=2;i<ido;i+=2){
    352        t3+=2;
    353        t4+=2;
    354        t5=t3-ido;
    355        t6=t4-ido;
    356        for(k=0;k<l1;k++){
    357          t5+=ido;
    358          t6+=ido;
    359          c1[t5-1]=ch[t5-1]+ch[t6-1];
    360          c1[t6-1]=ch[t5]-ch[t6];
    361          c1[t5]=ch[t5]+ch[t6];
    362          c1[t6]=ch[t6-1]-ch[t5-1];
    363        }
    364      }
    365    }
    366  }else{
    367    for(j=1;j<ipph;j++){
    368      t1+=t0;
    369      t2-=t0;
    370      t3=t1;
    371      t4=t2;
    372      for(k=0;k<l1;k++){
    373        t5=t3;
    374        t6=t4;
    375        for(i=2;i<ido;i+=2){
    376          t5+=2;
    377          t6+=2;
    378          c1[t5-1]=ch[t5-1]+ch[t6-1];
    379          c1[t6-1]=ch[t5]-ch[t6];
    380          c1[t5]=ch[t5]+ch[t6];
    381          c1[t6]=ch[t6-1]-ch[t5-1];
    382        }
    383        t3+=ido;
    384        t4+=ido;
    385      }
    386    }
    387  }
    388 
    389 L119:
    390  for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
    391 
    392  t1=0;
    393  t2=ipp2*idl1;
    394  for(j=1;j<ipph;j++){
    395    t1+=t0;
    396    t2-=t0;
    397    t3=t1-ido;
    398    t4=t2-ido;
    399    for(k=0;k<l1;k++){
    400      t3+=ido;
    401      t4+=ido;
    402      c1[t3]=ch[t3]+ch[t4];
    403      c1[t4]=ch[t4]-ch[t3];
    404    }
    405  }
    406 
    407  ar1=1.f;
    408  ai1=0.f;
    409  t1=0;
    410  t2=ipp2*idl1;
    411  t3=(ip-1)*idl1;
    412  for(l=1;l<ipph;l++){
    413    t1+=idl1;
    414    t2-=idl1;
    415    ar1h=dcp*ar1-dsp*ai1;
    416    ai1=dcp*ai1+dsp*ar1;
    417    ar1=ar1h;
    418    t4=t1;
    419    t5=t2;
    420    t6=t3;
    421    t7=idl1;
    422 
    423    for(ik=0;ik<idl1;ik++){
    424      ch2[t4++]=c2[ik]+ar1*c2[t7++];
    425      ch2[t5++]=ai1*c2[t6++];
    426    }
    427 
    428    dc2=ar1;
    429    ds2=ai1;
    430    ar2=ar1;
    431    ai2=ai1;
    432 
    433    t4=idl1;
    434    t5=(ipp2-1)*idl1;
    435    for(j=2;j<ipph;j++){
    436      t4+=idl1;
    437      t5-=idl1;
    438 
    439      ar2h=dc2*ar2-ds2*ai2;
    440      ai2=dc2*ai2+ds2*ar2;
    441      ar2=ar2h;
    442 
    443      t6=t1;
    444      t7=t2;
    445      t8=t4;
    446      t9=t5;
    447      for(ik=0;ik<idl1;ik++){
    448        ch2[t6++]+=ar2*c2[t8++];
    449        ch2[t7++]+=ai2*c2[t9++];
    450      }
    451    }
    452  }
    453 
    454  t1=0;
    455  for(j=1;j<ipph;j++){
    456    t1+=idl1;
    457    t2=t1;
    458    for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
    459  }
    460 
    461  if(ido<l1)goto L132;
    462 
    463  t1=0;
    464  t2=0;
    465  for(k=0;k<l1;k++){
    466    t3=t1;
    467    t4=t2;
    468    for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
    469    t1+=ido;
    470    t2+=t10;
    471  }
    472 
    473  goto L135;
    474 
    475 L132:
    476  for(i=0;i<ido;i++){
    477    t1=i;
    478    t2=i;
    479    for(k=0;k<l1;k++){
    480      cc[t2]=ch[t1];
    481      t1+=ido;
    482      t2+=t10;
    483    }
    484  }
    485 
    486 L135:
    487  t1=0;
    488  t2=ido<<1;
    489  t3=0;
    490  t4=ipp2*t0;
    491  for(j=1;j<ipph;j++){
    492 
    493    t1+=t2;
    494    t3+=t0;
    495    t4-=t0;
    496 
    497    t5=t1;
    498    t6=t3;
    499    t7=t4;
    500 
    501    for(k=0;k<l1;k++){
    502      cc[t5-1]=ch[t6];
    503      cc[t5]=ch[t7];
    504      t5+=t10;
    505      t6+=ido;
    506      t7+=ido;
    507    }
    508  }
    509 
    510  if(ido==1)return;
    511  if(nbd<l1)goto L141;
    512 
    513  t1=-ido;
    514  t3=0;
    515  t4=0;
    516  t5=ipp2*t0;
    517  for(j=1;j<ipph;j++){
    518    t1+=t2;
    519    t3+=t2;
    520    t4+=t0;
    521    t5-=t0;
    522    t6=t1;
    523    t7=t3;
    524    t8=t4;
    525    t9=t5;
    526    for(k=0;k<l1;k++){
    527      for(i=2;i<ido;i+=2){
    528        ic=idp2-i;
    529        cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
    530        cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
    531        cc[i+t7]=ch[i+t8]+ch[i+t9];
    532        cc[ic+t6]=ch[i+t9]-ch[i+t8];
    533      }
    534      t6+=t10;
    535      t7+=t10;
    536      t8+=ido;
    537      t9+=ido;
    538    }
    539  }
    540  return;
    541 
    542 L141:
    543 
    544  t1=-ido;
    545  t3=0;
    546  t4=0;
    547  t5=ipp2*t0;
    548  for(j=1;j<ipph;j++){
    549    t1+=t2;
    550    t3+=t2;
    551    t4+=t0;
    552    t5-=t0;
    553    for(i=2;i<ido;i+=2){
    554      t6=idp2+t1-i;
    555      t7=i+t3;
    556      t8=i+t4;
    557      t9=i+t5;
    558      for(k=0;k<l1;k++){
    559        cc[t7-1]=ch[t8-1]+ch[t9-1];
    560        cc[t6-1]=ch[t8-1]-ch[t9-1];
    561        cc[t7]=ch[t8]+ch[t9];
    562        cc[t6]=ch[t9]-ch[t8];
    563        t6+=t10;
    564        t7+=t10;
    565        t8+=ido;
    566        t9+=ido;
    567      }
    568    }
    569  }
    570 }
    571 
    572 static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){
    573  int i,k1,l1,l2;
    574  int na,kh,nf;
    575  int ip,iw,ido,idl1,ix2,ix3;
    576 
    577  nf=ifac[1];
    578  na=1;
    579  l2=n;
    580  iw=n;
    581 
    582  for(k1=0;k1<nf;k1++){
    583    kh=nf-k1;
    584    ip=ifac[kh+1];
    585    l1=l2/ip;
    586    ido=n/l2;
    587    idl1=ido*l1;
    588    iw-=(ip-1)*ido;
    589    na=1-na;
    590 
    591    if(ip!=4)goto L102;
    592 
    593    ix2=iw+ido;
    594    ix3=ix2+ido;
    595    if(na!=0)
    596      dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
    597    else
    598      dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
    599    goto L110;
    600 
    601 L102:
    602    if(ip!=2)goto L104;
    603    if(na!=0)goto L103;
    604 
    605    dradf2(ido,l1,c,ch,wa+iw-1);
    606    goto L110;
    607 
    608  L103:
    609    dradf2(ido,l1,ch,c,wa+iw-1);
    610    goto L110;
    611 
    612  L104:
    613    if(ido==1)na=1-na;
    614    if(na!=0)goto L109;
    615 
    616    dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
    617    na=1;
    618    goto L110;
    619 
    620  L109:
    621    dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
    622    na=0;
    623 
    624  L110:
    625    l2=l1;
    626  }
    627 
    628  if(na==1)return;
    629 
    630  for(i=0;i<n;i++)c[i]=ch[i];
    631 }
    632 
    633 static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){
    634  int i,k,t0,t1,t2,t3,t4,t5,t6;
    635  float ti2,tr2;
    636 
    637  t0=l1*ido;
    638 
    639  t1=0;
    640  t2=0;
    641  t3=(ido<<1)-1;
    642  for(k=0;k<l1;k++){
    643    ch[t1]=cc[t2]+cc[t3+t2];
    644    ch[t1+t0]=cc[t2]-cc[t3+t2];
    645    t2=(t1+=ido)<<1;
    646  }
    647 
    648  if(ido<2)return;
    649  if(ido==2)goto L105;
    650 
    651  t1=0;
    652  t2=0;
    653  for(k=0;k<l1;k++){
    654    t3=t1;
    655    t5=(t4=t2)+(ido<<1);
    656    t6=t0+t1;
    657    for(i=2;i<ido;i+=2){
    658      t3+=2;
    659      t4+=2;
    660      t5-=2;
    661      t6+=2;
    662      ch[t3-1]=cc[t4-1]+cc[t5-1];
    663      tr2=cc[t4-1]-cc[t5-1];
    664      ch[t3]=cc[t4]-cc[t5];
    665      ti2=cc[t4]+cc[t5];
    666      ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
    667      ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
    668    }
    669    t2=(t1+=ido)<<1;
    670  }
    671 
    672  if(ido%2==1)return;
    673 
    674 L105:
    675  t1=ido-1;
    676  t2=ido-1;
    677  for(k=0;k<l1;k++){
    678    ch[t1]=cc[t2]+cc[t2];
    679    ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
    680    t1+=ido;
    681    t2+=ido<<1;
    682  }
    683 }
    684 
    685 static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1,
    686                          float *wa2){
    687  static float taur = -.5f;
    688  static float taui = .8660254037844386f;
    689  int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
    690  float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
    691  t0=l1*ido;
    692 
    693  t1=0;
    694  t2=t0<<1;
    695  t3=ido<<1;
    696  t4=ido+(ido<<1);
    697  t5=0;
    698  for(k=0;k<l1;k++){
    699    tr2=cc[t3-1]+cc[t3-1];
    700    cr2=cc[t5]+(taur*tr2);
    701    ch[t1]=cc[t5]+tr2;
    702    ci3=taui*(cc[t3]+cc[t3]);
    703    ch[t1+t0]=cr2-ci3;
    704    ch[t1+t2]=cr2+ci3;
    705    t1+=ido;
    706    t3+=t4;
    707    t5+=t4;
    708  }
    709 
    710  if(ido==1)return;
    711 
    712  t1=0;
    713  t3=ido<<1;
    714  for(k=0;k<l1;k++){
    715    t7=t1+(t1<<1);
    716    t6=(t5=t7+t3);
    717    t8=t1;
    718    t10=(t9=t1+t0)+t0;
    719 
    720    for(i=2;i<ido;i+=2){
    721      t5+=2;
    722      t6-=2;
    723      t7+=2;
    724      t8+=2;
    725      t9+=2;
    726      t10+=2;
    727      tr2=cc[t5-1]+cc[t6-1];
    728      cr2=cc[t7-1]+(taur*tr2);
    729      ch[t8-1]=cc[t7-1]+tr2;
    730      ti2=cc[t5]-cc[t6];
    731      ci2=cc[t7]+(taur*ti2);
    732      ch[t8]=cc[t7]+ti2;
    733      cr3=taui*(cc[t5-1]-cc[t6-1]);
    734      ci3=taui*(cc[t5]+cc[t6]);
    735      dr2=cr2-ci3;
    736      dr3=cr2+ci3;
    737      di2=ci2+cr3;
    738      di3=ci2-cr3;
    739      ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;
    740      ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;
    741      ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;
    742      ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;
    743    }
    744    t1+=ido;
    745  }
    746 }
    747 
    748 static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1,
    749                          float *wa2,float *wa3){
    750  static float sqrt2=1.414213562373095f;
    751  int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
    752  float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
    753  t0=l1*ido;
    754 
    755  t1=0;
    756  t2=ido<<2;
    757  t3=0;
    758  t6=ido<<1;
    759  for(k=0;k<l1;k++){
    760    t4=t3+t6;
    761    t5=t1;
    762    tr3=cc[t4-1]+cc[t4-1];
    763    tr4=cc[t4]+cc[t4];
    764    tr1=cc[t3]-cc[(t4+=t6)-1];
    765    tr2=cc[t3]+cc[t4-1];
    766    ch[t5]=tr2+tr3;
    767    ch[t5+=t0]=tr1-tr4;
    768    ch[t5+=t0]=tr2-tr3;
    769    ch[t5+=t0]=tr1+tr4;
    770    t1+=ido;
    771    t3+=t2;
    772  }
    773 
    774  if(ido<2)return;
    775  if(ido==2)goto L105;
    776 
    777  t1=0;
    778  for(k=0;k<l1;k++){
    779    t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
    780    t7=t1;
    781    for(i=2;i<ido;i+=2){
    782      t2+=2;
    783      t3+=2;
    784      t4-=2;
    785      t5-=2;
    786      t7+=2;
    787      ti1=cc[t2]+cc[t5];
    788      ti2=cc[t2]-cc[t5];
    789      ti3=cc[t3]-cc[t4];
    790      tr4=cc[t3]+cc[t4];
    791      tr1=cc[t2-1]-cc[t5-1];
    792      tr2=cc[t2-1]+cc[t5-1];
    793      ti4=cc[t3-1]-cc[t4-1];
    794      tr3=cc[t3-1]+cc[t4-1];
    795      ch[t7-1]=tr2+tr3;
    796      cr3=tr2-tr3;
    797      ch[t7]=ti2+ti3;
    798      ci3=ti2-ti3;
    799      cr2=tr1-tr4;
    800      cr4=tr1+tr4;
    801      ci2=ti1+ti4;
    802      ci4=ti1-ti4;
    803 
    804      ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
    805      ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
    806      ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
    807      ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
    808      ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
    809      ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
    810    }
    811    t1+=ido;
    812  }
    813 
    814  if(ido%2 == 1)return;
    815 
    816 L105:
    817 
    818  t1=ido;
    819  t2=ido<<2;
    820  t3=ido-1;
    821  t4=ido+(ido<<1);
    822  for(k=0;k<l1;k++){
    823    t5=t3;
    824    ti1=cc[t1]+cc[t4];
    825    ti2=cc[t4]-cc[t1];
    826    tr1=cc[t1-1]-cc[t4-1];
    827    tr2=cc[t1-1]+cc[t4-1];
    828    ch[t5]=tr2+tr2;
    829    ch[t5+=t0]=sqrt2*(tr1-ti1);
    830    ch[t5+=t0]=ti2+ti2;
    831    ch[t5+=t0]=-sqrt2*(tr1+ti1);
    832 
    833    t3+=ido;
    834    t1+=t2;
    835    t4+=t2;
    836  }
    837 }
    838 
    839 static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
    840            float *c2,float *ch,float *ch2,float *wa){
    841  static float tpi=6.283185307179586f;
    842  int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
    843      t11,t12;
    844  float dc2,ai1,ai2,ar1,ar2,ds2;
    845  int nbd;
    846  float dcp,arg,dsp,ar1h,ar2h;
    847  int ipp2;
    848 
    849  t10=ip*ido;
    850  t0=l1*ido;
    851  arg=tpi/(float)ip;
    852  dcp=cos(arg);
    853  dsp=sin(arg);
    854  nbd=(ido-1)>>1;
    855  ipp2=ip;
    856  ipph=(ip+1)>>1;
    857  if(ido<l1)goto L103;
    858 
    859  t1=0;
    860  t2=0;
    861  for(k=0;k<l1;k++){
    862    t3=t1;
    863    t4=t2;
    864    for(i=0;i<ido;i++){
    865      ch[t3]=cc[t4];
    866      t3++;
    867      t4++;
    868    }
    869    t1+=ido;
    870    t2+=t10;
    871  }
    872  goto L106;
    873 
    874 L103:
    875  t1=0;
    876  for(i=0;i<ido;i++){
    877    t2=t1;
    878    t3=t1;
    879    for(k=0;k<l1;k++){
    880      ch[t2]=cc[t3];
    881      t2+=ido;
    882      t3+=t10;
    883    }
    884    t1++;
    885  }
    886 
    887 L106:
    888  t1=0;
    889  t2=ipp2*t0;
    890  t7=(t5=ido<<1);
    891  for(j=1;j<ipph;j++){
    892    t1+=t0;
    893    t2-=t0;
    894    t3=t1;
    895    t4=t2;
    896    t6=t5;
    897    for(k=0;k<l1;k++){
    898      ch[t3]=cc[t6-1]+cc[t6-1];
    899      ch[t4]=cc[t6]+cc[t6];
    900      t3+=ido;
    901      t4+=ido;
    902      t6+=t10;
    903    }
    904    t5+=t7;
    905  }
    906 
    907  if (ido == 1)goto L116;
    908  if(nbd<l1)goto L112;
    909 
    910  t1=0;
    911  t2=ipp2*t0;
    912  t7=0;
    913  for(j=1;j<ipph;j++){
    914    t1+=t0;
    915    t2-=t0;
    916    t3=t1;
    917    t4=t2;
    918 
    919    t7+=(ido<<1);
    920    t8=t7;
    921    for(k=0;k<l1;k++){
    922      t5=t3;
    923      t6=t4;
    924      t9=t8;
    925      t11=t8;
    926      for(i=2;i<ido;i+=2){
    927        t5+=2;
    928        t6+=2;
    929        t9+=2;
    930        t11-=2;
    931        ch[t5-1]=cc[t9-1]+cc[t11-1];
    932        ch[t6-1]=cc[t9-1]-cc[t11-1];
    933        ch[t5]=cc[t9]-cc[t11];
    934        ch[t6]=cc[t9]+cc[t11];
    935      }
    936      t3+=ido;
    937      t4+=ido;
    938      t8+=t10;
    939    }
    940  }
    941  goto L116;
    942 
    943 L112:
    944  t1=0;
    945  t2=ipp2*t0;
    946  t7=0;
    947  for(j=1;j<ipph;j++){
    948    t1+=t0;
    949    t2-=t0;
    950    t3=t1;
    951    t4=t2;
    952    t7+=(ido<<1);
    953    t8=t7;
    954    t9=t7;
    955    for(i=2;i<ido;i+=2){
    956      t3+=2;
    957      t4+=2;
    958      t8+=2;
    959      t9-=2;
    960      t5=t3;
    961      t6=t4;
    962      t11=t8;
    963      t12=t9;
    964      for(k=0;k<l1;k++){
    965        ch[t5-1]=cc[t11-1]+cc[t12-1];
    966        ch[t6-1]=cc[t11-1]-cc[t12-1];
    967        ch[t5]=cc[t11]-cc[t12];
    968        ch[t6]=cc[t11]+cc[t12];
    969        t5+=ido;
    970        t6+=ido;
    971        t11+=t10;
    972        t12+=t10;
    973      }
    974    }
    975  }
    976 
    977 L116:
    978  ar1=1.f;
    979  ai1=0.f;
    980  t1=0;
    981  t9=(t2=ipp2*idl1);
    982  t3=(ip-1)*idl1;
    983  for(l=1;l<ipph;l++){
    984    t1+=idl1;
    985    t2-=idl1;
    986 
    987    ar1h=dcp*ar1-dsp*ai1;
    988    ai1=dcp*ai1+dsp*ar1;
    989    ar1=ar1h;
    990    t4=t1;
    991    t5=t2;
    992    t6=0;
    993    t7=idl1;
    994    t8=t3;
    995    for(ik=0;ik<idl1;ik++){
    996      c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
    997      c2[t5++]=ai1*ch2[t8++];
    998    }
    999    dc2=ar1;
   1000    ds2=ai1;
   1001    ar2=ar1;
   1002    ai2=ai1;
   1003 
   1004    t6=idl1;
   1005    t7=t9-idl1;
   1006    for(j=2;j<ipph;j++){
   1007      t6+=idl1;
   1008      t7-=idl1;
   1009      ar2h=dc2*ar2-ds2*ai2;
   1010      ai2=dc2*ai2+ds2*ar2;
   1011      ar2=ar2h;
   1012      t4=t1;
   1013      t5=t2;
   1014      t11=t6;
   1015      t12=t7;
   1016      for(ik=0;ik<idl1;ik++){
   1017        c2[t4++]+=ar2*ch2[t11++];
   1018        c2[t5++]+=ai2*ch2[t12++];
   1019      }
   1020    }
   1021  }
   1022 
   1023  t1=0;
   1024  for(j=1;j<ipph;j++){
   1025    t1+=idl1;
   1026    t2=t1;
   1027    for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
   1028  }
   1029 
   1030  t1=0;
   1031  t2=ipp2*t0;
   1032  for(j=1;j<ipph;j++){
   1033    t1+=t0;
   1034    t2-=t0;
   1035    t3=t1;
   1036    t4=t2;
   1037    for(k=0;k<l1;k++){
   1038      ch[t3]=c1[t3]-c1[t4];
   1039      ch[t4]=c1[t3]+c1[t4];
   1040      t3+=ido;
   1041      t4+=ido;
   1042    }
   1043  }
   1044 
   1045  if(ido==1)goto L132;
   1046  if(nbd<l1)goto L128;
   1047 
   1048  t1=0;
   1049  t2=ipp2*t0;
   1050  for(j=1;j<ipph;j++){
   1051    t1+=t0;
   1052    t2-=t0;
   1053    t3=t1;
   1054    t4=t2;
   1055    for(k=0;k<l1;k++){
   1056      t5=t3;
   1057      t6=t4;
   1058      for(i=2;i<ido;i+=2){
   1059        t5+=2;
   1060        t6+=2;
   1061        ch[t5-1]=c1[t5-1]-c1[t6];
   1062        ch[t6-1]=c1[t5-1]+c1[t6];
   1063        ch[t5]=c1[t5]+c1[t6-1];
   1064        ch[t6]=c1[t5]-c1[t6-1];
   1065      }
   1066      t3+=ido;
   1067      t4+=ido;
   1068    }
   1069  }
   1070  goto L132;
   1071 
   1072 L128:
   1073  t1=0;
   1074  t2=ipp2*t0;
   1075  for(j=1;j<ipph;j++){
   1076    t1+=t0;
   1077    t2-=t0;
   1078    t3=t1;
   1079    t4=t2;
   1080    for(i=2;i<ido;i+=2){
   1081      t3+=2;
   1082      t4+=2;
   1083      t5=t3;
   1084      t6=t4;
   1085      for(k=0;k<l1;k++){
   1086        ch[t5-1]=c1[t5-1]-c1[t6];
   1087        ch[t6-1]=c1[t5-1]+c1[t6];
   1088        ch[t5]=c1[t5]+c1[t6-1];
   1089        ch[t6]=c1[t5]-c1[t6-1];
   1090        t5+=ido;
   1091        t6+=ido;
   1092      }
   1093    }
   1094  }
   1095 
   1096 L132:
   1097  if(ido==1)return;
   1098 
   1099  for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
   1100 
   1101  t1=0;
   1102  for(j=1;j<ip;j++){
   1103    t2=(t1+=t0);
   1104    for(k=0;k<l1;k++){
   1105      c1[t2]=ch[t2];
   1106      t2+=ido;
   1107    }
   1108  }
   1109 
   1110  if(nbd>l1)goto L139;
   1111 
   1112  is= -ido-1;
   1113  t1=0;
   1114  for(j=1;j<ip;j++){
   1115    is+=ido;
   1116    t1+=t0;
   1117    idij=is;
   1118    t2=t1;
   1119    for(i=2;i<ido;i+=2){
   1120      t2+=2;
   1121      idij+=2;
   1122      t3=t2;
   1123      for(k=0;k<l1;k++){
   1124        c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
   1125        c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
   1126        t3+=ido;
   1127      }
   1128    }
   1129  }
   1130  return;
   1131 
   1132 L139:
   1133  is= -ido-1;
   1134  t1=0;
   1135  for(j=1;j<ip;j++){
   1136    is+=ido;
   1137    t1+=t0;
   1138    t2=t1;
   1139    for(k=0;k<l1;k++){
   1140      idij=is;
   1141      t3=t2;
   1142      for(i=2;i<ido;i+=2){
   1143        idij+=2;
   1144        t3+=2;
   1145        c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
   1146        c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
   1147      }
   1148      t2+=ido;
   1149    }
   1150  }
   1151 }
   1152 
   1153 static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){
   1154  int i,k1,l1,l2;
   1155  int na;
   1156  int nf,ip,iw,ix2,ix3,ido,idl1;
   1157 
   1158  nf=ifac[1];
   1159  na=0;
   1160  l1=1;
   1161  iw=1;
   1162 
   1163  for(k1=0;k1<nf;k1++){
   1164    ip=ifac[k1 + 2];
   1165    l2=ip*l1;
   1166    ido=n/l2;
   1167    idl1=ido*l1;
   1168    if(ip!=4)goto L103;
   1169    ix2=iw+ido;
   1170    ix3=ix2+ido;
   1171 
   1172    if(na!=0)
   1173      dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
   1174    else
   1175      dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
   1176    na=1-na;
   1177    goto L115;
   1178 
   1179  L103:
   1180    if(ip!=2)goto L106;
   1181 
   1182    if(na!=0)
   1183      dradb2(ido,l1,ch,c,wa+iw-1);
   1184    else
   1185      dradb2(ido,l1,c,ch,wa+iw-1);
   1186    na=1-na;
   1187    goto L115;
   1188 
   1189  L106:
   1190    if(ip!=3)goto L109;
   1191 
   1192    ix2=iw+ido;
   1193    if(na!=0)
   1194      dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);
   1195    else
   1196      dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);
   1197    na=1-na;
   1198    goto L115;
   1199 
   1200  L109:
   1201 /*    The radix five case can be translated later..... */
   1202 /*    if(ip!=5)goto L112;
   1203 
   1204    ix2=iw+ido;
   1205    ix3=ix2+ido;
   1206    ix4=ix3+ido;
   1207    if(na!=0)
   1208      dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
   1209    else
   1210      dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
   1211    na=1-na;
   1212    goto L115;
   1213 
   1214  L112:*/
   1215    if(na!=0)
   1216      dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
   1217    else
   1218      dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
   1219    if(ido==1)na=1-na;
   1220 
   1221  L115:
   1222    l1=l2;
   1223    iw+=(ip-1)*ido;
   1224  }
   1225 
   1226  if(na==0)return;
   1227 
   1228  for(i=0;i<n;i++)c[i]=ch[i];
   1229 }
   1230 
   1231 void drft_forward(drft_lookup *l,float *data){
   1232  if(l->n==1)return;
   1233  drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
   1234 }
   1235 
   1236 void drft_backward(drft_lookup *l,float *data){
   1237  if (l->n==1)return;
   1238  drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
   1239 }
   1240 
   1241 void drft_init(drft_lookup *l,int n){
   1242  l->n=n;
   1243  l->trigcache=_ogg_calloc(3*n,sizeof(*l->trigcache));
   1244  l->splitcache=_ogg_calloc(32,sizeof(*l->splitcache));
   1245  fdrffti(n, l->trigcache, l->splitcache);
   1246 }
   1247 
   1248 void drft_clear(drft_lookup *l){
   1249  if(l){
   1250    if(l->trigcache)_ogg_free(l->trigcache);
   1251    if(l->splitcache)_ogg_free(l->splitcache);
   1252    memset(l,0,sizeof(*l));
   1253  }
   1254 }