1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/libs/vorbis/smallft.c	Sun Nov 01 00:36:56 2015 +0200
     1.3 @@ -0,0 +1,1255 @@
     1.4 +/********************************************************************
     1.5 + *                                                                  *
     1.6 + * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
     1.7 + * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
     1.8 + * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
     1.9 + * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
    1.10 + *                                                                  *
    1.11 + * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
    1.12 + * by the Xiph.Org Foundation http://www.xiph.org/                  *
    1.13 + *                                                                  *
    1.14 + ********************************************************************
    1.15 +
    1.16 + function: *unnormalized* fft transform
    1.17 + last mod: $Id: smallft.c 16227 2009-07-08 06:58:46Z xiphmont $
    1.18 +
    1.19 + ********************************************************************/
    1.20 +
    1.21 +/* FFT implementation from OggSquish, minus cosine transforms,
    1.22 + * minus all but radix 2/4 case.  In Vorbis we only need this
    1.23 + * cut-down version.
    1.24 + *
    1.25 + * To do more than just power-of-two sized vectors, see the full
    1.26 + * version I wrote for NetLib.
    1.27 + *
    1.28 + * Note that the packing is a little strange; rather than the FFT r/i
    1.29 + * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1,
    1.30 + * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
    1.31 + * FORTRAN version
    1.32 + */
    1.33 +
    1.34 +#include <stdlib.h>
    1.35 +#include <string.h>
    1.36 +#include <math.h>
    1.37 +#include "smallft.h"
    1.38 +#include "os.h"
    1.39 +#include "misc.h"
    1.40 +
    1.41 +static void drfti1(int n, float *wa, int *ifac){
    1.42 +  static int ntryh[4] = { 4,2,3,5 };
    1.43 +  static float tpi = 6.28318530717958648f;
    1.44 +  float arg,argh,argld,fi;
    1.45 +  int ntry=0,i,j=-1;
    1.46 +  int k1, l1, l2, ib;
    1.47 +  int ld, ii, ip, is, nq, nr;
    1.48 +  int ido, ipm, nfm1;
    1.49 +  int nl=n;
    1.50 +  int nf=0;
    1.51 +
    1.52 + L101:
    1.53 +  j++;
    1.54 +  if (j < 4)
    1.55 +    ntry=ntryh[j];
    1.56 +  else
    1.57 +    ntry+=2;
    1.58 +
    1.59 + L104:
    1.60 +  nq=nl/ntry;
    1.61 +  nr=nl-ntry*nq;
    1.62 +  if (nr!=0) goto L101;
    1.63 +
    1.64 +  nf++;
    1.65 +  ifac[nf+1]=ntry;
    1.66 +  nl=nq;
    1.67 +  if(ntry!=2)goto L107;
    1.68 +  if(nf==1)goto L107;
    1.69 +
    1.70 +  for (i=1;i<nf;i++){
    1.71 +    ib=nf-i+1;
    1.72 +    ifac[ib+1]=ifac[ib];
    1.73 +  }
    1.74 +  ifac[2] = 2;
    1.75 +
    1.76 + L107:
    1.77 +  if(nl!=1)goto L104;
    1.78 +  ifac[0]=n;
    1.79 +  ifac[1]=nf;
    1.80 +  argh=tpi/n;
    1.81 +  is=0;
    1.82 +  nfm1=nf-1;
    1.83 +  l1=1;
    1.84 +
    1.85 +  if(nfm1==0)return;
    1.86 +
    1.87 +  for (k1=0;k1<nfm1;k1++){
    1.88 +    ip=ifac[k1+2];
    1.89 +    ld=0;
    1.90 +    l2=l1*ip;
    1.91 +    ido=n/l2;
    1.92 +    ipm=ip-1;
    1.93 +
    1.94 +    for (j=0;j<ipm;j++){
    1.95 +      ld+=l1;
    1.96 +      i=is;
    1.97 +      argld=(float)ld*argh;
    1.98 +      fi=0.f;
    1.99 +      for (ii=2;ii<ido;ii+=2){
   1.100 +        fi+=1.f;
   1.101 +        arg=fi*argld;
   1.102 +        wa[i++]=cos(arg);
   1.103 +        wa[i++]=sin(arg);
   1.104 +      }
   1.105 +      is+=ido;
   1.106 +    }
   1.107 +    l1=l2;
   1.108 +  }
   1.109 +}
   1.110 +
   1.111 +static void fdrffti(int n, float *wsave, int *ifac){
   1.112 +
   1.113 +  if (n == 1) return;
   1.114 +  drfti1(n, wsave+n, ifac);
   1.115 +}
   1.116 +
   1.117 +static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){
   1.118 +  int i,k;
   1.119 +  float ti2,tr2;
   1.120 +  int t0,t1,t2,t3,t4,t5,t6;
   1.121 +
   1.122 +  t1=0;
   1.123 +  t0=(t2=l1*ido);
   1.124 +  t3=ido<<1;
   1.125 +  for(k=0;k<l1;k++){
   1.126 +    ch[t1<<1]=cc[t1]+cc[t2];
   1.127 +    ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
   1.128 +    t1+=ido;
   1.129 +    t2+=ido;
   1.130 +  }
   1.131 +
   1.132 +  if(ido<2)return;
   1.133 +  if(ido==2)goto L105;
   1.134 +
   1.135 +  t1=0;
   1.136 +  t2=t0;
   1.137 +  for(k=0;k<l1;k++){
   1.138 +    t3=t2;
   1.139 +    t4=(t1<<1)+(ido<<1);
   1.140 +    t5=t1;
   1.141 +    t6=t1+t1;
   1.142 +    for(i=2;i<ido;i+=2){
   1.143 +      t3+=2;
   1.144 +      t4-=2;
   1.145 +      t5+=2;
   1.146 +      t6+=2;
   1.147 +      tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
   1.148 +      ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
   1.149 +      ch[t6]=cc[t5]+ti2;
   1.150 +      ch[t4]=ti2-cc[t5];
   1.151 +      ch[t6-1]=cc[t5-1]+tr2;
   1.152 +      ch[t4-1]=cc[t5-1]-tr2;
   1.153 +    }
   1.154 +    t1+=ido;
   1.155 +    t2+=ido;
   1.156 +  }
   1.157 +
   1.158 +  if(ido%2==1)return;
   1.159 +
   1.160 + L105:
   1.161 +  t3=(t2=(t1=ido)-1);
   1.162 +  t2+=t0;
   1.163 +  for(k=0;k<l1;k++){
   1.164 +    ch[t1]=-cc[t2];
   1.165 +    ch[t1-1]=cc[t3];
   1.166 +    t1+=ido<<1;
   1.167 +    t2+=ido;
   1.168 +    t3+=ido;
   1.169 +  }
   1.170 +}
   1.171 +
   1.172 +static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1,
   1.173 +            float *wa2,float *wa3){
   1.174 +  static float hsqt2 = .70710678118654752f;
   1.175 +  int i,k,t0,t1,t2,t3,t4,t5,t6;
   1.176 +  float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
   1.177 +  t0=l1*ido;
   1.178 +
   1.179 +  t1=t0;
   1.180 +  t4=t1<<1;
   1.181 +  t2=t1+(t1<<1);
   1.182 +  t3=0;
   1.183 +
   1.184 +  for(k=0;k<l1;k++){
   1.185 +    tr1=cc[t1]+cc[t2];
   1.186 +    tr2=cc[t3]+cc[t4];
   1.187 +
   1.188 +    ch[t5=t3<<2]=tr1+tr2;
   1.189 +    ch[(ido<<2)+t5-1]=tr2-tr1;
   1.190 +    ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
   1.191 +    ch[t5]=cc[t2]-cc[t1];
   1.192 +
   1.193 +    t1+=ido;
   1.194 +    t2+=ido;
   1.195 +    t3+=ido;
   1.196 +    t4+=ido;
   1.197 +  }
   1.198 +
   1.199 +  if(ido<2)return;
   1.200 +  if(ido==2)goto L105;
   1.201 +
   1.202 +
   1.203 +  t1=0;
   1.204 +  for(k=0;k<l1;k++){
   1.205 +    t2=t1;
   1.206 +    t4=t1<<2;
   1.207 +    t5=(t6=ido<<1)+t4;
   1.208 +    for(i=2;i<ido;i+=2){
   1.209 +      t3=(t2+=2);
   1.210 +      t4+=2;
   1.211 +      t5-=2;
   1.212 +
   1.213 +      t3+=t0;
   1.214 +      cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
   1.215 +      ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
   1.216 +      t3+=t0;
   1.217 +      cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];
   1.218 +      ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];
   1.219 +      t3+=t0;
   1.220 +      cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];
   1.221 +      ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];
   1.222 +
   1.223 +      tr1=cr2+cr4;
   1.224 +      tr4=cr4-cr2;
   1.225 +      ti1=ci2+ci4;
   1.226 +      ti4=ci2-ci4;
   1.227 +
   1.228 +      ti2=cc[t2]+ci3;
   1.229 +      ti3=cc[t2]-ci3;
   1.230 +      tr2=cc[t2-1]+cr3;
   1.231 +      tr3=cc[t2-1]-cr3;
   1.232 +
   1.233 +      ch[t4-1]=tr1+tr2;
   1.234 +      ch[t4]=ti1+ti2;
   1.235 +
   1.236 +      ch[t5-1]=tr3-ti4;
   1.237 +      ch[t5]=tr4-ti3;
   1.238 +
   1.239 +      ch[t4+t6-1]=ti4+tr3;
   1.240 +      ch[t4+t6]=tr4+ti3;
   1.241 +
   1.242 +      ch[t5+t6-1]=tr2-tr1;
   1.243 +      ch[t5+t6]=ti1-ti2;
   1.244 +    }
   1.245 +    t1+=ido;
   1.246 +  }
   1.247 +  if(ido&1)return;
   1.248 +
   1.249 + L105:
   1.250 +
   1.251 +  t2=(t1=t0+ido-1)+(t0<<1);
   1.252 +  t3=ido<<2;
   1.253 +  t4=ido;
   1.254 +  t5=ido<<1;
   1.255 +  t6=ido;
   1.256 +
   1.257 +  for(k=0;k<l1;k++){
   1.258 +    ti1=-hsqt2*(cc[t1]+cc[t2]);
   1.259 +    tr1=hsqt2*(cc[t1]-cc[t2]);
   1.260 +
   1.261 +    ch[t4-1]=tr1+cc[t6-1];
   1.262 +    ch[t4+t5-1]=cc[t6-1]-tr1;
   1.263 +
   1.264 +    ch[t4]=ti1-cc[t1+t0];
   1.265 +    ch[t4+t5]=ti1+cc[t1+t0];
   1.266 +
   1.267 +    t1+=ido;
   1.268 +    t2+=ido;
   1.269 +    t4+=t3;
   1.270 +    t6+=ido;
   1.271 +  }
   1.272 +}
   1.273 +
   1.274 +static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
   1.275 +                          float *c2,float *ch,float *ch2,float *wa){
   1.276 +
   1.277 +  static float tpi=6.283185307179586f;
   1.278 +  int idij,ipph,i,j,k,l,ic,ik,is;
   1.279 +  int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
   1.280 +  float dc2,ai1,ai2,ar1,ar2,ds2;
   1.281 +  int nbd;
   1.282 +  float dcp,arg,dsp,ar1h,ar2h;
   1.283 +  int idp2,ipp2;
   1.284 +
   1.285 +  arg=tpi/(float)ip;
   1.286 +  dcp=cos(arg);
   1.287 +  dsp=sin(arg);
   1.288 +  ipph=(ip+1)>>1;
   1.289 +  ipp2=ip;
   1.290 +  idp2=ido;
   1.291 +  nbd=(ido-1)>>1;
   1.292 +  t0=l1*ido;
   1.293 +  t10=ip*ido;
   1.294 +
   1.295 +  if(ido==1)goto L119;
   1.296 +  for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
   1.297 +
   1.298 +  t1=0;
   1.299 +  for(j=1;j<ip;j++){
   1.300 +    t1+=t0;
   1.301 +    t2=t1;
   1.302 +    for(k=0;k<l1;k++){
   1.303 +      ch[t2]=c1[t2];
   1.304 +      t2+=ido;
   1.305 +    }
   1.306 +  }
   1.307 +
   1.308 +  is=-ido;
   1.309 +  t1=0;
   1.310 +  if(nbd>l1){
   1.311 +    for(j=1;j<ip;j++){
   1.312 +      t1+=t0;
   1.313 +      is+=ido;
   1.314 +      t2= -ido+t1;
   1.315 +      for(k=0;k<l1;k++){
   1.316 +        idij=is-1;
   1.317 +        t2+=ido;
   1.318 +        t3=t2;
   1.319 +        for(i=2;i<ido;i+=2){
   1.320 +          idij+=2;
   1.321 +          t3+=2;
   1.322 +          ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
   1.323 +          ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
   1.324 +        }
   1.325 +      }
   1.326 +    }
   1.327 +  }else{
   1.328 +
   1.329 +    for(j=1;j<ip;j++){
   1.330 +      is+=ido;
   1.331 +      idij=is-1;
   1.332 +      t1+=t0;
   1.333 +      t2=t1;
   1.334 +      for(i=2;i<ido;i+=2){
   1.335 +        idij+=2;
   1.336 +        t2+=2;
   1.337 +        t3=t2;
   1.338 +        for(k=0;k<l1;k++){
   1.339 +          ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
   1.340 +          ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
   1.341 +          t3+=ido;
   1.342 +        }
   1.343 +      }
   1.344 +    }
   1.345 +  }
   1.346 +
   1.347 +  t1=0;
   1.348 +  t2=ipp2*t0;
   1.349 +  if(nbd<l1){
   1.350 +    for(j=1;j<ipph;j++){
   1.351 +      t1+=t0;
   1.352 +      t2-=t0;
   1.353 +      t3=t1;
   1.354 +      t4=t2;
   1.355 +      for(i=2;i<ido;i+=2){
   1.356 +        t3+=2;
   1.357 +        t4+=2;
   1.358 +        t5=t3-ido;
   1.359 +        t6=t4-ido;
   1.360 +        for(k=0;k<l1;k++){
   1.361 +          t5+=ido;
   1.362 +          t6+=ido;
   1.363 +          c1[t5-1]=ch[t5-1]+ch[t6-1];
   1.364 +          c1[t6-1]=ch[t5]-ch[t6];
   1.365 +          c1[t5]=ch[t5]+ch[t6];
   1.366 +          c1[t6]=ch[t6-1]-ch[t5-1];
   1.367 +        }
   1.368 +      }
   1.369 +    }
   1.370 +  }else{
   1.371 +    for(j=1;j<ipph;j++){
   1.372 +      t1+=t0;
   1.373 +      t2-=t0;
   1.374 +      t3=t1;
   1.375 +      t4=t2;
   1.376 +      for(k=0;k<l1;k++){
   1.377 +        t5=t3;
   1.378 +        t6=t4;
   1.379 +        for(i=2;i<ido;i+=2){
   1.380 +          t5+=2;
   1.381 +          t6+=2;
   1.382 +          c1[t5-1]=ch[t5-1]+ch[t6-1];
   1.383 +          c1[t6-1]=ch[t5]-ch[t6];
   1.384 +          c1[t5]=ch[t5]+ch[t6];
   1.385 +          c1[t6]=ch[t6-1]-ch[t5-1];
   1.386 +        }
   1.387 +        t3+=ido;
   1.388 +        t4+=ido;
   1.389 +      }
   1.390 +    }
   1.391 +  }
   1.392 +
   1.393 +L119:
   1.394 +  for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
   1.395 +
   1.396 +  t1=0;
   1.397 +  t2=ipp2*idl1;
   1.398 +  for(j=1;j<ipph;j++){
   1.399 +    t1+=t0;
   1.400 +    t2-=t0;
   1.401 +    t3=t1-ido;
   1.402 +    t4=t2-ido;
   1.403 +    for(k=0;k<l1;k++){
   1.404 +      t3+=ido;
   1.405 +      t4+=ido;
   1.406 +      c1[t3]=ch[t3]+ch[t4];
   1.407 +      c1[t4]=ch[t4]-ch[t3];
   1.408 +    }
   1.409 +  }
   1.410 +
   1.411 +  ar1=1.f;
   1.412 +  ai1=0.f;
   1.413 +  t1=0;
   1.414 +  t2=ipp2*idl1;
   1.415 +  t3=(ip-1)*idl1;
   1.416 +  for(l=1;l<ipph;l++){
   1.417 +    t1+=idl1;
   1.418 +    t2-=idl1;
   1.419 +    ar1h=dcp*ar1-dsp*ai1;
   1.420 +    ai1=dcp*ai1+dsp*ar1;
   1.421 +    ar1=ar1h;
   1.422 +    t4=t1;
   1.423 +    t5=t2;
   1.424 +    t6=t3;
   1.425 +    t7=idl1;
   1.426 +
   1.427 +    for(ik=0;ik<idl1;ik++){
   1.428 +      ch2[t4++]=c2[ik]+ar1*c2[t7++];
   1.429 +      ch2[t5++]=ai1*c2[t6++];
   1.430 +    }
   1.431 +
   1.432 +    dc2=ar1;
   1.433 +    ds2=ai1;
   1.434 +    ar2=ar1;
   1.435 +    ai2=ai1;
   1.436 +
   1.437 +    t4=idl1;
   1.438 +    t5=(ipp2-1)*idl1;
   1.439 +    for(j=2;j<ipph;j++){
   1.440 +      t4+=idl1;
   1.441 +      t5-=idl1;
   1.442 +
   1.443 +      ar2h=dc2*ar2-ds2*ai2;
   1.444 +      ai2=dc2*ai2+ds2*ar2;
   1.445 +      ar2=ar2h;
   1.446 +
   1.447 +      t6=t1;
   1.448 +      t7=t2;
   1.449 +      t8=t4;
   1.450 +      t9=t5;
   1.451 +      for(ik=0;ik<idl1;ik++){
   1.452 +        ch2[t6++]+=ar2*c2[t8++];
   1.453 +        ch2[t7++]+=ai2*c2[t9++];
   1.454 +      }
   1.455 +    }
   1.456 +  }
   1.457 +
   1.458 +  t1=0;
   1.459 +  for(j=1;j<ipph;j++){
   1.460 +    t1+=idl1;
   1.461 +    t2=t1;
   1.462 +    for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
   1.463 +  }
   1.464 +
   1.465 +  if(ido<l1)goto L132;
   1.466 +
   1.467 +  t1=0;
   1.468 +  t2=0;
   1.469 +  for(k=0;k<l1;k++){
   1.470 +    t3=t1;
   1.471 +    t4=t2;
   1.472 +    for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
   1.473 +    t1+=ido;
   1.474 +    t2+=t10;
   1.475 +  }
   1.476 +
   1.477 +  goto L135;
   1.478 +
   1.479 + L132:
   1.480 +  for(i=0;i<ido;i++){
   1.481 +    t1=i;
   1.482 +    t2=i;
   1.483 +    for(k=0;k<l1;k++){
   1.484 +      cc[t2]=ch[t1];
   1.485 +      t1+=ido;
   1.486 +      t2+=t10;
   1.487 +    }
   1.488 +  }
   1.489 +
   1.490 + L135:
   1.491 +  t1=0;
   1.492 +  t2=ido<<1;
   1.493 +  t3=0;
   1.494 +  t4=ipp2*t0;
   1.495 +  for(j=1;j<ipph;j++){
   1.496 +
   1.497 +    t1+=t2;
   1.498 +    t3+=t0;
   1.499 +    t4-=t0;
   1.500 +
   1.501 +    t5=t1;
   1.502 +    t6=t3;
   1.503 +    t7=t4;
   1.504 +
   1.505 +    for(k=0;k<l1;k++){
   1.506 +      cc[t5-1]=ch[t6];
   1.507 +      cc[t5]=ch[t7];
   1.508 +      t5+=t10;
   1.509 +      t6+=ido;
   1.510 +      t7+=ido;
   1.511 +    }
   1.512 +  }
   1.513 +
   1.514 +  if(ido==1)return;
   1.515 +  if(nbd<l1)goto L141;
   1.516 +
   1.517 +  t1=-ido;
   1.518 +  t3=0;
   1.519 +  t4=0;
   1.520 +  t5=ipp2*t0;
   1.521 +  for(j=1;j<ipph;j++){
   1.522 +    t1+=t2;
   1.523 +    t3+=t2;
   1.524 +    t4+=t0;
   1.525 +    t5-=t0;
   1.526 +    t6=t1;
   1.527 +    t7=t3;
   1.528 +    t8=t4;
   1.529 +    t9=t5;
   1.530 +    for(k=0;k<l1;k++){
   1.531 +      for(i=2;i<ido;i+=2){
   1.532 +        ic=idp2-i;
   1.533 +        cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
   1.534 +        cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
   1.535 +        cc[i+t7]=ch[i+t8]+ch[i+t9];
   1.536 +        cc[ic+t6]=ch[i+t9]-ch[i+t8];
   1.537 +      }
   1.538 +      t6+=t10;
   1.539 +      t7+=t10;
   1.540 +      t8+=ido;
   1.541 +      t9+=ido;
   1.542 +    }
   1.543 +  }
   1.544 +  return;
   1.545 +
   1.546 + L141:
   1.547 +
   1.548 +  t1=-ido;
   1.549 +  t3=0;
   1.550 +  t4=0;
   1.551 +  t5=ipp2*t0;
   1.552 +  for(j=1;j<ipph;j++){
   1.553 +    t1+=t2;
   1.554 +    t3+=t2;
   1.555 +    t4+=t0;
   1.556 +    t5-=t0;
   1.557 +    for(i=2;i<ido;i+=2){
   1.558 +      t6=idp2+t1-i;
   1.559 +      t7=i+t3;
   1.560 +      t8=i+t4;
   1.561 +      t9=i+t5;
   1.562 +      for(k=0;k<l1;k++){
   1.563 +        cc[t7-1]=ch[t8-1]+ch[t9-1];
   1.564 +        cc[t6-1]=ch[t8-1]-ch[t9-1];
   1.565 +        cc[t7]=ch[t8]+ch[t9];
   1.566 +        cc[t6]=ch[t9]-ch[t8];
   1.567 +        t6+=t10;
   1.568 +        t7+=t10;
   1.569 +        t8+=ido;
   1.570 +        t9+=ido;
   1.571 +      }
   1.572 +    }
   1.573 +  }
   1.574 +}
   1.575 +
   1.576 +static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){
   1.577 +  int i,k1,l1,l2;
   1.578 +  int na,kh,nf;
   1.579 +  int ip,iw,ido,idl1,ix2,ix3;
   1.580 +
   1.581 +  nf=ifac[1];
   1.582 +  na=1;
   1.583 +  l2=n;
   1.584 +  iw=n;
   1.585 +
   1.586 +  for(k1=0;k1<nf;k1++){
   1.587 +    kh=nf-k1;
   1.588 +    ip=ifac[kh+1];
   1.589 +    l1=l2/ip;
   1.590 +    ido=n/l2;
   1.591 +    idl1=ido*l1;
   1.592 +    iw-=(ip-1)*ido;
   1.593 +    na=1-na;
   1.594 +
   1.595 +    if(ip!=4)goto L102;
   1.596 +
   1.597 +    ix2=iw+ido;
   1.598 +    ix3=ix2+ido;
   1.599 +    if(na!=0)
   1.600 +      dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
   1.601 +    else
   1.602 +      dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
   1.603 +    goto L110;
   1.604 +
   1.605 + L102:
   1.606 +    if(ip!=2)goto L104;
   1.607 +    if(na!=0)goto L103;
   1.608 +
   1.609 +    dradf2(ido,l1,c,ch,wa+iw-1);
   1.610 +    goto L110;
   1.611 +
   1.612 +  L103:
   1.613 +    dradf2(ido,l1,ch,c,wa+iw-1);
   1.614 +    goto L110;
   1.615 +
   1.616 +  L104:
   1.617 +    if(ido==1)na=1-na;
   1.618 +    if(na!=0)goto L109;
   1.619 +
   1.620 +    dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
   1.621 +    na=1;
   1.622 +    goto L110;
   1.623 +
   1.624 +  L109:
   1.625 +    dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
   1.626 +    na=0;
   1.627 +
   1.628 +  L110:
   1.629 +    l2=l1;
   1.630 +  }
   1.631 +
   1.632 +  if(na==1)return;
   1.633 +
   1.634 +  for(i=0;i<n;i++)c[i]=ch[i];
   1.635 +}
   1.636 +
   1.637 +static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){
   1.638 +  int i,k,t0,t1,t2,t3,t4,t5,t6;
   1.639 +  float ti2,tr2;
   1.640 +
   1.641 +  t0=l1*ido;
   1.642 +
   1.643 +  t1=0;
   1.644 +  t2=0;
   1.645 +  t3=(ido<<1)-1;
   1.646 +  for(k=0;k<l1;k++){
   1.647 +    ch[t1]=cc[t2]+cc[t3+t2];
   1.648 +    ch[t1+t0]=cc[t2]-cc[t3+t2];
   1.649 +    t2=(t1+=ido)<<1;
   1.650 +  }
   1.651 +
   1.652 +  if(ido<2)return;
   1.653 +  if(ido==2)goto L105;
   1.654 +
   1.655 +  t1=0;
   1.656 +  t2=0;
   1.657 +  for(k=0;k<l1;k++){
   1.658 +    t3=t1;
   1.659 +    t5=(t4=t2)+(ido<<1);
   1.660 +    t6=t0+t1;
   1.661 +    for(i=2;i<ido;i+=2){
   1.662 +      t3+=2;
   1.663 +      t4+=2;
   1.664 +      t5-=2;
   1.665 +      t6+=2;
   1.666 +      ch[t3-1]=cc[t4-1]+cc[t5-1];
   1.667 +      tr2=cc[t4-1]-cc[t5-1];
   1.668 +      ch[t3]=cc[t4]-cc[t5];
   1.669 +      ti2=cc[t4]+cc[t5];
   1.670 +      ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
   1.671 +      ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
   1.672 +    }
   1.673 +    t2=(t1+=ido)<<1;
   1.674 +  }
   1.675 +
   1.676 +  if(ido%2==1)return;
   1.677 +
   1.678 +L105:
   1.679 +  t1=ido-1;
   1.680 +  t2=ido-1;
   1.681 +  for(k=0;k<l1;k++){
   1.682 +    ch[t1]=cc[t2]+cc[t2];
   1.683 +    ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
   1.684 +    t1+=ido;
   1.685 +    t2+=ido<<1;
   1.686 +  }
   1.687 +}
   1.688 +
   1.689 +static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1,
   1.690 +                          float *wa2){
   1.691 +  static float taur = -.5f;
   1.692 +  static float taui = .8660254037844386f;
   1.693 +  int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
   1.694 +  float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
   1.695 +  t0=l1*ido;
   1.696 +
   1.697 +  t1=0;
   1.698 +  t2=t0<<1;
   1.699 +  t3=ido<<1;
   1.700 +  t4=ido+(ido<<1);
   1.701 +  t5=0;
   1.702 +  for(k=0;k<l1;k++){
   1.703 +    tr2=cc[t3-1]+cc[t3-1];
   1.704 +    cr2=cc[t5]+(taur*tr2);
   1.705 +    ch[t1]=cc[t5]+tr2;
   1.706 +    ci3=taui*(cc[t3]+cc[t3]);
   1.707 +    ch[t1+t0]=cr2-ci3;
   1.708 +    ch[t1+t2]=cr2+ci3;
   1.709 +    t1+=ido;
   1.710 +    t3+=t4;
   1.711 +    t5+=t4;
   1.712 +  }
   1.713 +
   1.714 +  if(ido==1)return;
   1.715 +
   1.716 +  t1=0;
   1.717 +  t3=ido<<1;
   1.718 +  for(k=0;k<l1;k++){
   1.719 +    t7=t1+(t1<<1);
   1.720 +    t6=(t5=t7+t3);
   1.721 +    t8=t1;
   1.722 +    t10=(t9=t1+t0)+t0;
   1.723 +
   1.724 +    for(i=2;i<ido;i+=2){
   1.725 +      t5+=2;
   1.726 +      t6-=2;
   1.727 +      t7+=2;
   1.728 +      t8+=2;
   1.729 +      t9+=2;
   1.730 +      t10+=2;
   1.731 +      tr2=cc[t5-1]+cc[t6-1];
   1.732 +      cr2=cc[t7-1]+(taur*tr2);
   1.733 +      ch[t8-1]=cc[t7-1]+tr2;
   1.734 +      ti2=cc[t5]-cc[t6];
   1.735 +      ci2=cc[t7]+(taur*ti2);
   1.736 +      ch[t8]=cc[t7]+ti2;
   1.737 +      cr3=taui*(cc[t5-1]-cc[t6-1]);
   1.738 +      ci3=taui*(cc[t5]+cc[t6]);
   1.739 +      dr2=cr2-ci3;
   1.740 +      dr3=cr2+ci3;
   1.741 +      di2=ci2+cr3;
   1.742 +      di3=ci2-cr3;
   1.743 +      ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;
   1.744 +      ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;
   1.745 +      ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;
   1.746 +      ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;
   1.747 +    }
   1.748 +    t1+=ido;
   1.749 +  }
   1.750 +}
   1.751 +
   1.752 +static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1,
   1.753 +                          float *wa2,float *wa3){
   1.754 +  static float sqrt2=1.414213562373095f;
   1.755 +  int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
   1.756 +  float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
   1.757 +  t0=l1*ido;
   1.758 +
   1.759 +  t1=0;
   1.760 +  t2=ido<<2;
   1.761 +  t3=0;
   1.762 +  t6=ido<<1;
   1.763 +  for(k=0;k<l1;k++){
   1.764 +    t4=t3+t6;
   1.765 +    t5=t1;
   1.766 +    tr3=cc[t4-1]+cc[t4-1];
   1.767 +    tr4=cc[t4]+cc[t4];
   1.768 +    tr1=cc[t3]-cc[(t4+=t6)-1];
   1.769 +    tr2=cc[t3]+cc[t4-1];
   1.770 +    ch[t5]=tr2+tr3;
   1.771 +    ch[t5+=t0]=tr1-tr4;
   1.772 +    ch[t5+=t0]=tr2-tr3;
   1.773 +    ch[t5+=t0]=tr1+tr4;
   1.774 +    t1+=ido;
   1.775 +    t3+=t2;
   1.776 +  }
   1.777 +
   1.778 +  if(ido<2)return;
   1.779 +  if(ido==2)goto L105;
   1.780 +
   1.781 +  t1=0;
   1.782 +  for(k=0;k<l1;k++){
   1.783 +    t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
   1.784 +    t7=t1;
   1.785 +    for(i=2;i<ido;i+=2){
   1.786 +      t2+=2;
   1.787 +      t3+=2;
   1.788 +      t4-=2;
   1.789 +      t5-=2;
   1.790 +      t7+=2;
   1.791 +      ti1=cc[t2]+cc[t5];
   1.792 +      ti2=cc[t2]-cc[t5];
   1.793 +      ti3=cc[t3]-cc[t4];
   1.794 +      tr4=cc[t3]+cc[t4];
   1.795 +      tr1=cc[t2-1]-cc[t5-1];
   1.796 +      tr2=cc[t2-1]+cc[t5-1];
   1.797 +      ti4=cc[t3-1]-cc[t4-1];
   1.798 +      tr3=cc[t3-1]+cc[t4-1];
   1.799 +      ch[t7-1]=tr2+tr3;
   1.800 +      cr3=tr2-tr3;
   1.801 +      ch[t7]=ti2+ti3;
   1.802 +      ci3=ti2-ti3;
   1.803 +      cr2=tr1-tr4;
   1.804 +      cr4=tr1+tr4;
   1.805 +      ci2=ti1+ti4;
   1.806 +      ci4=ti1-ti4;
   1.807 +
   1.808 +      ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
   1.809 +      ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
   1.810 +      ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
   1.811 +      ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
   1.812 +      ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
   1.813 +      ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
   1.814 +    }
   1.815 +    t1+=ido;
   1.816 +  }
   1.817 +
   1.818 +  if(ido%2 == 1)return;
   1.819 +
   1.820 + L105:
   1.821 +
   1.822 +  t1=ido;
   1.823 +  t2=ido<<2;
   1.824 +  t3=ido-1;
   1.825 +  t4=ido+(ido<<1);
   1.826 +  for(k=0;k<l1;k++){
   1.827 +    t5=t3;
   1.828 +    ti1=cc[t1]+cc[t4];
   1.829 +    ti2=cc[t4]-cc[t1];
   1.830 +    tr1=cc[t1-1]-cc[t4-1];
   1.831 +    tr2=cc[t1-1]+cc[t4-1];
   1.832 +    ch[t5]=tr2+tr2;
   1.833 +    ch[t5+=t0]=sqrt2*(tr1-ti1);
   1.834 +    ch[t5+=t0]=ti2+ti2;
   1.835 +    ch[t5+=t0]=-sqrt2*(tr1+ti1);
   1.836 +
   1.837 +    t3+=ido;
   1.838 +    t1+=t2;
   1.839 +    t4+=t2;
   1.840 +  }
   1.841 +}
   1.842 +
   1.843 +static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
   1.844 +            float *c2,float *ch,float *ch2,float *wa){
   1.845 +  static float tpi=6.283185307179586f;
   1.846 +  int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
   1.847 +      t11,t12;
   1.848 +  float dc2,ai1,ai2,ar1,ar2,ds2;
   1.849 +  int nbd;
   1.850 +  float dcp,arg,dsp,ar1h,ar2h;
   1.851 +  int ipp2;
   1.852 +
   1.853 +  t10=ip*ido;
   1.854 +  t0=l1*ido;
   1.855 +  arg=tpi/(float)ip;
   1.856 +  dcp=cos(arg);
   1.857 +  dsp=sin(arg);
   1.858 +  nbd=(ido-1)>>1;
   1.859 +  ipp2=ip;
   1.860 +  ipph=(ip+1)>>1;
   1.861 +  if(ido<l1)goto L103;
   1.862 +
   1.863 +  t1=0;
   1.864 +  t2=0;
   1.865 +  for(k=0;k<l1;k++){
   1.866 +    t3=t1;
   1.867 +    t4=t2;
   1.868 +    for(i=0;i<ido;i++){
   1.869 +      ch[t3]=cc[t4];
   1.870 +      t3++;
   1.871 +      t4++;
   1.872 +    }
   1.873 +    t1+=ido;
   1.874 +    t2+=t10;
   1.875 +  }
   1.876 +  goto L106;
   1.877 +
   1.878 + L103:
   1.879 +  t1=0;
   1.880 +  for(i=0;i<ido;i++){
   1.881 +    t2=t1;
   1.882 +    t3=t1;
   1.883 +    for(k=0;k<l1;k++){
   1.884 +      ch[t2]=cc[t3];
   1.885 +      t2+=ido;
   1.886 +      t3+=t10;
   1.887 +    }
   1.888 +    t1++;
   1.889 +  }
   1.890 +
   1.891 + L106:
   1.892 +  t1=0;
   1.893 +  t2=ipp2*t0;
   1.894 +  t7=(t5=ido<<1);
   1.895 +  for(j=1;j<ipph;j++){
   1.896 +    t1+=t0;
   1.897 +    t2-=t0;
   1.898 +    t3=t1;
   1.899 +    t4=t2;
   1.900 +    t6=t5;
   1.901 +    for(k=0;k<l1;k++){
   1.902 +      ch[t3]=cc[t6-1]+cc[t6-1];
   1.903 +      ch[t4]=cc[t6]+cc[t6];
   1.904 +      t3+=ido;
   1.905 +      t4+=ido;
   1.906 +      t6+=t10;
   1.907 +    }
   1.908 +    t5+=t7;
   1.909 +  }
   1.910 +
   1.911 +  if (ido == 1)goto L116;
   1.912 +  if(nbd<l1)goto L112;
   1.913 +
   1.914 +  t1=0;
   1.915 +  t2=ipp2*t0;
   1.916 +  t7=0;
   1.917 +  for(j=1;j<ipph;j++){
   1.918 +    t1+=t0;
   1.919 +    t2-=t0;
   1.920 +    t3=t1;
   1.921 +    t4=t2;
   1.922 +
   1.923 +    t7+=(ido<<1);
   1.924 +    t8=t7;
   1.925 +    for(k=0;k<l1;k++){
   1.926 +      t5=t3;
   1.927 +      t6=t4;
   1.928 +      t9=t8;
   1.929 +      t11=t8;
   1.930 +      for(i=2;i<ido;i+=2){
   1.931 +        t5+=2;
   1.932 +        t6+=2;
   1.933 +        t9+=2;
   1.934 +        t11-=2;
   1.935 +        ch[t5-1]=cc[t9-1]+cc[t11-1];
   1.936 +        ch[t6-1]=cc[t9-1]-cc[t11-1];
   1.937 +        ch[t5]=cc[t9]-cc[t11];
   1.938 +        ch[t6]=cc[t9]+cc[t11];
   1.939 +      }
   1.940 +      t3+=ido;
   1.941 +      t4+=ido;
   1.942 +      t8+=t10;
   1.943 +    }
   1.944 +  }
   1.945 +  goto L116;
   1.946 +
   1.947 + L112:
   1.948 +  t1=0;
   1.949 +  t2=ipp2*t0;
   1.950 +  t7=0;
   1.951 +  for(j=1;j<ipph;j++){
   1.952 +    t1+=t0;
   1.953 +    t2-=t0;
   1.954 +    t3=t1;
   1.955 +    t4=t2;
   1.956 +    t7+=(ido<<1);
   1.957 +    t8=t7;
   1.958 +    t9=t7;
   1.959 +    for(i=2;i<ido;i+=2){
   1.960 +      t3+=2;
   1.961 +      t4+=2;
   1.962 +      t8+=2;
   1.963 +      t9-=2;
   1.964 +      t5=t3;
   1.965 +      t6=t4;
   1.966 +      t11=t8;
   1.967 +      t12=t9;
   1.968 +      for(k=0;k<l1;k++){
   1.969 +        ch[t5-1]=cc[t11-1]+cc[t12-1];
   1.970 +        ch[t6-1]=cc[t11-1]-cc[t12-1];
   1.971 +        ch[t5]=cc[t11]-cc[t12];
   1.972 +        ch[t6]=cc[t11]+cc[t12];
   1.973 +        t5+=ido;
   1.974 +        t6+=ido;
   1.975 +        t11+=t10;
   1.976 +        t12+=t10;
   1.977 +      }
   1.978 +    }
   1.979 +  }
   1.980 +
   1.981 +L116:
   1.982 +  ar1=1.f;
   1.983 +  ai1=0.f;
   1.984 +  t1=0;
   1.985 +  t9=(t2=ipp2*idl1);
   1.986 +  t3=(ip-1)*idl1;
   1.987 +  for(l=1;l<ipph;l++){
   1.988 +    t1+=idl1;
   1.989 +    t2-=idl1;
   1.990 +
   1.991 +    ar1h=dcp*ar1-dsp*ai1;
   1.992 +    ai1=dcp*ai1+dsp*ar1;
   1.993 +    ar1=ar1h;
   1.994 +    t4=t1;
   1.995 +    t5=t2;
   1.996 +    t6=0;
   1.997 +    t7=idl1;
   1.998 +    t8=t3;
   1.999 +    for(ik=0;ik<idl1;ik++){
  1.1000 +      c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
  1.1001 +      c2[t5++]=ai1*ch2[t8++];
  1.1002 +    }
  1.1003 +    dc2=ar1;
  1.1004 +    ds2=ai1;
  1.1005 +    ar2=ar1;
  1.1006 +    ai2=ai1;
  1.1007 +
  1.1008 +    t6=idl1;
  1.1009 +    t7=t9-idl1;
  1.1010 +    for(j=2;j<ipph;j++){
  1.1011 +      t6+=idl1;
  1.1012 +      t7-=idl1;
  1.1013 +      ar2h=dc2*ar2-ds2*ai2;
  1.1014 +      ai2=dc2*ai2+ds2*ar2;
  1.1015 +      ar2=ar2h;
  1.1016 +      t4=t1;
  1.1017 +      t5=t2;
  1.1018 +      t11=t6;
  1.1019 +      t12=t7;
  1.1020 +      for(ik=0;ik<idl1;ik++){
  1.1021 +        c2[t4++]+=ar2*ch2[t11++];
  1.1022 +        c2[t5++]+=ai2*ch2[t12++];
  1.1023 +      }
  1.1024 +    }
  1.1025 +  }
  1.1026 +
  1.1027 +  t1=0;
  1.1028 +  for(j=1;j<ipph;j++){
  1.1029 +    t1+=idl1;
  1.1030 +    t2=t1;
  1.1031 +    for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
  1.1032 +  }
  1.1033 +
  1.1034 +  t1=0;
  1.1035 +  t2=ipp2*t0;
  1.1036 +  for(j=1;j<ipph;j++){
  1.1037 +    t1+=t0;
  1.1038 +    t2-=t0;
  1.1039 +    t3=t1;
  1.1040 +    t4=t2;
  1.1041 +    for(k=0;k<l1;k++){
  1.1042 +      ch[t3]=c1[t3]-c1[t4];
  1.1043 +      ch[t4]=c1[t3]+c1[t4];
  1.1044 +      t3+=ido;
  1.1045 +      t4+=ido;
  1.1046 +    }
  1.1047 +  }
  1.1048 +
  1.1049 +  if(ido==1)goto L132;
  1.1050 +  if(nbd<l1)goto L128;
  1.1051 +
  1.1052 +  t1=0;
  1.1053 +  t2=ipp2*t0;
  1.1054 +  for(j=1;j<ipph;j++){
  1.1055 +    t1+=t0;
  1.1056 +    t2-=t0;
  1.1057 +    t3=t1;
  1.1058 +    t4=t2;
  1.1059 +    for(k=0;k<l1;k++){
  1.1060 +      t5=t3;
  1.1061 +      t6=t4;
  1.1062 +      for(i=2;i<ido;i+=2){
  1.1063 +        t5+=2;
  1.1064 +        t6+=2;
  1.1065 +        ch[t5-1]=c1[t5-1]-c1[t6];
  1.1066 +        ch[t6-1]=c1[t5-1]+c1[t6];
  1.1067 +        ch[t5]=c1[t5]+c1[t6-1];
  1.1068 +        ch[t6]=c1[t5]-c1[t6-1];
  1.1069 +      }
  1.1070 +      t3+=ido;
  1.1071 +      t4+=ido;
  1.1072 +    }
  1.1073 +  }
  1.1074 +  goto L132;
  1.1075 +
  1.1076 + L128:
  1.1077 +  t1=0;
  1.1078 +  t2=ipp2*t0;
  1.1079 +  for(j=1;j<ipph;j++){
  1.1080 +    t1+=t0;
  1.1081 +    t2-=t0;
  1.1082 +    t3=t1;
  1.1083 +    t4=t2;
  1.1084 +    for(i=2;i<ido;i+=2){
  1.1085 +      t3+=2;
  1.1086 +      t4+=2;
  1.1087 +      t5=t3;
  1.1088 +      t6=t4;
  1.1089 +      for(k=0;k<l1;k++){
  1.1090 +        ch[t5-1]=c1[t5-1]-c1[t6];
  1.1091 +        ch[t6-1]=c1[t5-1]+c1[t6];
  1.1092 +        ch[t5]=c1[t5]+c1[t6-1];
  1.1093 +        ch[t6]=c1[t5]-c1[t6-1];
  1.1094 +        t5+=ido;
  1.1095 +        t6+=ido;
  1.1096 +      }
  1.1097 +    }
  1.1098 +  }
  1.1099 +
  1.1100 +L132:
  1.1101 +  if(ido==1)return;
  1.1102 +
  1.1103 +  for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
  1.1104 +
  1.1105 +  t1=0;
  1.1106 +  for(j=1;j<ip;j++){
  1.1107 +    t2=(t1+=t0);
  1.1108 +    for(k=0;k<l1;k++){
  1.1109 +      c1[t2]=ch[t2];
  1.1110 +      t2+=ido;
  1.1111 +    }
  1.1112 +  }
  1.1113 +
  1.1114 +  if(nbd>l1)goto L139;
  1.1115 +
  1.1116 +  is= -ido-1;
  1.1117 +  t1=0;
  1.1118 +  for(j=1;j<ip;j++){
  1.1119 +    is+=ido;
  1.1120 +    t1+=t0;
  1.1121 +    idij=is;
  1.1122 +    t2=t1;
  1.1123 +    for(i=2;i<ido;i+=2){
  1.1124 +      t2+=2;
  1.1125 +      idij+=2;
  1.1126 +      t3=t2;
  1.1127 +      for(k=0;k<l1;k++){
  1.1128 +        c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
  1.1129 +        c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
  1.1130 +        t3+=ido;
  1.1131 +      }
  1.1132 +    }
  1.1133 +  }
  1.1134 +  return;
  1.1135 +
  1.1136 + L139:
  1.1137 +  is= -ido-1;
  1.1138 +  t1=0;
  1.1139 +  for(j=1;j<ip;j++){
  1.1140 +    is+=ido;
  1.1141 +    t1+=t0;
  1.1142 +    t2=t1;
  1.1143 +    for(k=0;k<l1;k++){
  1.1144 +      idij=is;
  1.1145 +      t3=t2;
  1.1146 +      for(i=2;i<ido;i+=2){
  1.1147 +        idij+=2;
  1.1148 +        t3+=2;
  1.1149 +        c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
  1.1150 +        c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
  1.1151 +      }
  1.1152 +      t2+=ido;
  1.1153 +    }
  1.1154 +  }
  1.1155 +}
  1.1156 +
  1.1157 +static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){
  1.1158 +  int i,k1,l1,l2;
  1.1159 +  int na;
  1.1160 +  int nf,ip,iw,ix2,ix3,ido,idl1;
  1.1161 +
  1.1162 +  nf=ifac[1];
  1.1163 +  na=0;
  1.1164 +  l1=1;
  1.1165 +  iw=1;
  1.1166 +
  1.1167 +  for(k1=0;k1<nf;k1++){
  1.1168 +    ip=ifac[k1 + 2];
  1.1169 +    l2=ip*l1;
  1.1170 +    ido=n/l2;
  1.1171 +    idl1=ido*l1;
  1.1172 +    if(ip!=4)goto L103;
  1.1173 +    ix2=iw+ido;
  1.1174 +    ix3=ix2+ido;
  1.1175 +
  1.1176 +    if(na!=0)
  1.1177 +      dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
  1.1178 +    else
  1.1179 +      dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
  1.1180 +    na=1-na;
  1.1181 +    goto L115;
  1.1182 +
  1.1183 +  L103:
  1.1184 +    if(ip!=2)goto L106;
  1.1185 +
  1.1186 +    if(na!=0)
  1.1187 +      dradb2(ido,l1,ch,c,wa+iw-1);
  1.1188 +    else
  1.1189 +      dradb2(ido,l1,c,ch,wa+iw-1);
  1.1190 +    na=1-na;
  1.1191 +    goto L115;
  1.1192 +
  1.1193 +  L106:
  1.1194 +    if(ip!=3)goto L109;
  1.1195 +
  1.1196 +    ix2=iw+ido;
  1.1197 +    if(na!=0)
  1.1198 +      dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);
  1.1199 +    else
  1.1200 +      dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);
  1.1201 +    na=1-na;
  1.1202 +    goto L115;
  1.1203 +
  1.1204 +  L109:
  1.1205 +/*    The radix five case can be translated later..... */
  1.1206 +/*    if(ip!=5)goto L112;
  1.1207 +
  1.1208 +    ix2=iw+ido;
  1.1209 +    ix3=ix2+ido;
  1.1210 +    ix4=ix3+ido;
  1.1211 +    if(na!=0)
  1.1212 +      dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
  1.1213 +    else
  1.1214 +      dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
  1.1215 +    na=1-na;
  1.1216 +    goto L115;
  1.1217 +
  1.1218 +  L112:*/
  1.1219 +    if(na!=0)
  1.1220 +      dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
  1.1221 +    else
  1.1222 +      dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
  1.1223 +    if(ido==1)na=1-na;
  1.1224 +
  1.1225 +  L115:
  1.1226 +    l1=l2;
  1.1227 +    iw+=(ip-1)*ido;
  1.1228 +  }
  1.1229 +
  1.1230 +  if(na==0)return;
  1.1231 +
  1.1232 +  for(i=0;i<n;i++)c[i]=ch[i];
  1.1233 +}
  1.1234 +
  1.1235 +void drft_forward(drft_lookup *l,float *data){
  1.1236 +  if(l->n==1)return;
  1.1237 +  drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
  1.1238 +}
  1.1239 +
  1.1240 +void drft_backward(drft_lookup *l,float *data){
  1.1241 +  if (l->n==1)return;
  1.1242 +  drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
  1.1243 +}
  1.1244 +
  1.1245 +void drft_init(drft_lookup *l,int n){
  1.1246 +  l->n=n;
  1.1247 +  l->trigcache=_ogg_calloc(3*n,sizeof(*l->trigcache));
  1.1248 +  l->splitcache=_ogg_calloc(32,sizeof(*l->splitcache));
  1.1249 +  fdrffti(n, l->trigcache, l->splitcache);
  1.1250 +}
  1.1251 +
  1.1252 +void drft_clear(drft_lookup *l){
  1.1253 +  if(l){
  1.1254 +    if(l->trigcache)_ogg_free(l->trigcache);
  1.1255 +    if(l->splitcache)_ogg_free(l->splitcache);
  1.1256 +    memset(l,0,sizeof(*l));
  1.1257 +  }
  1.1258 +}