001/* JOrbis
002 * Copyright (C) 2000 ymnk, JCraft,Inc.
003 *  
004 * Written by: 2000 ymnk<ymnk@jcaft.com>
005 *   
006 * Many thanks to 
007 *   Monty <monty@xiph.org> and 
008 *   The XIPHOPHORUS Company http://www.xiph.org/ .
009 * JOrbis has been based on their awesome works, Vorbis codec.
010 *   
011 * This program is free software; you can redistribute it and/or
012 * modify it under the terms of the GNU Library General Public License
013 * as published by the Free Software Foundation; either version 2 of
014 * the License, or (at your option) any later version.
015   
016 * This program is distributed in the hope that it will be useful,
017 * but WITHOUT ANY WARRANTY; without even the implied warranty of
018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
019 * GNU Library General Public License for more details.
020 * 
021 * You should have received a copy of the GNU Library General Public
022 * License along with this program; if not, write to the Free Software
023 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
024 */
025
026package com.jcraft.jorbis;
027
028class Drft{
029  int n;
030  float[] trigcache;    
031  int[] splitcache;
032
033  void backward(float[] data){
034    //System.err.println("Drft.backward");
035    if(n==1)return;
036    drftb1(n,data,trigcache,trigcache,n,splitcache);
037  }
038
039  void init(int n){
040    //System.err.println("Drft.init");
041    this.n=n;
042    trigcache=new float[3*n];
043    splitcache=new int[32];
044    fdrffti(n, trigcache, splitcache);
045  }
046
047  void clear(){
048    //System.err.println("Drft.clear");
049    if(trigcache!=null)trigcache=null;
050    if(splitcache!=null)splitcache=null;
051//    memset(l,0,sizeof(drft_lookup));
052  }
053
054  static int[] ntryh = { 4,2,3,5 };
055  static float tpi = 6.28318530717958647692528676655900577f;
056  static float hsqt2 = .70710678118654752440084436210485f;
057  static float taui = .86602540378443864676372317075293618f;
058  static float taur = -.5f;
059  static float sqrt2 = 1.4142135623730950488016887242097f;
060
061  static void drfti1(int n, float[] wa, int index, int[] ifac){
062    float arg,argh,argld,fi;
063    int ntry=0,i,j=-1;
064    int k1, l1, l2, ib;
065    int ld, ii, ip, is, nq, nr;
066    int ido, ipm, nfm1;
067    int nl=n;
068    int nf=0;
069
070    int state=101;
071
072    loop: while(true){
073System.err.print("");
074      switch(state){
075      case 101:
076        j++;
077        if (j < 4)
078          ntry=ntryh[j];
079        else
080          ntry+=2;
081      case 104:
082        nq=nl/ntry;
083        nr=nl-ntry*nq;
084        if(nr!=0){
085          state=101;
086          break;
087        }
088        nf++;
089        ifac[nf+1]=ntry;
090        nl=nq;
091        if(ntry!=2){
092          state=107;
093          break;
094        }
095        if(nf==1){
096          state=107;
097          break;
098        }
099
100        for(i=1;i<nf;i++){
101          ib=nf-i+1;
102          ifac[ib+1]=ifac[ib];
103        }
104        ifac[2] = 2;
105      case 107:
106        if(nl!=1){
107          state=104;
108          break;
109        }
110        ifac[0]=n;
111        ifac[1]=nf;
112        argh=tpi/n;
113        is=0;
114        nfm1=nf-1;
115        l1=1;
116
117        if(nfm1==0)return;
118
119        for (k1=0;k1<nfm1;k1++){
120          ip=ifac[k1+2];
121          ld=0;
122          l2=l1*ip;
123          ido=n/l2;
124          ipm=ip-1;
125
126          for (j=0;j<ipm;j++){
127            ld+=l1;
128            i=is;
129            argld=(float)ld*argh;
130            fi=0.f;
131            for (ii=2;ii<ido;ii+=2){
132              fi+=1.f;
133              arg=fi*argld;
134              wa[index+i++]=(float)Math.cos(arg);
135              wa[index+i++]=(float)Math.sin(arg);
136            }
137            is+=ido;
138          }
139          l1=l2;
140        }
141        break loop;
142      }
143    }
144  }
145
146  static void fdrffti(int n, float[] wsave, int[] ifac){
147//System.err.println("fdrffti: n="+n);
148    if(n == 1) return;
149    drfti1(n, wsave, n, ifac);
150  }
151
152  static void dradf2(int ido,int l1,float[] cc, float[] ch, float[] wa1, int index){
153    int i,k;
154    float ti2,tr2;
155    int t0,t1,t2,t3,t4,t5,t6;
156
157    t1=0;
158    t0=(t2=l1*ido);
159    t3=ido<<1;
160    for(k=0;k<l1;k++){
161      ch[t1<<1]=cc[t1]+cc[t2];
162      ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
163      t1+=ido;
164      t2+=ido;
165    }
166    
167    if(ido<2)return;
168
169    if(ido!=2){
170      t1=0;
171      t2=t0;
172      for(k=0;k<l1;k++){
173        t3=t2;
174        t4=(t1<<1)+(ido<<1);
175        t5=t1;
176        t6=t1+t1;
177        for(i=2;i<ido;i+=2){
178          t3+=2;
179          t4-=2;
180          t5+=2;
181          t6+=2;
182          tr2=wa1[index+i-2]*cc[t3-1]+wa1[index+i-1]*cc[t3];
183          ti2=wa1[index+i-2]*cc[t3]-wa1[index+i-1]*cc[t3-1];
184          ch[t6]=cc[t5]+ti2;
185          ch[t4]=ti2-cc[t5];
186          ch[t6-1]=cc[t5-1]+tr2;
187          ch[t4-1]=cc[t5-1]-tr2;
188        }
189        t1+=ido;
190        t2+=ido;
191      }
192      if(ido%2==1)return;
193    }
194    
195    t3=(t2=(t1=ido)-1);
196    t2+=t0;
197    for(k=0;k<l1;k++){
198      ch[t1]=-cc[t2];
199      ch[t1-1]=cc[t3];
200      t1+=ido<<1;
201      t2+=ido;
202      t3+=ido;
203    }
204  }
205
206  static void dradf4(int ido,int l1,float[] cc, float[] ch, 
207                     float[] wa1, int index1,
208                     float[] wa2, int index2,
209                     float[] wa3, int index3){
210    int i,k,t0,t1,t2,t3,t4,t5,t6;
211    float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
212    t0=l1*ido;
213  
214    t1=t0;
215    t4=t1<<1;
216    t2=t1+(t1<<1);
217    t3=0;
218
219    for(k=0;k<l1;k++){
220      tr1=cc[t1]+cc[t2];
221      tr2=cc[t3]+cc[t4];
222
223      ch[t5=t3<<2]=tr1+tr2;
224      ch[(ido<<2)+t5-1]=tr2-tr1;
225      ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
226      ch[t5]=cc[t2]-cc[t1];
227      
228      t1+=ido;
229      t2+=ido;
230      t3+=ido;
231      t4+=ido;
232    }
233    if(ido<2)return;
234
235    if(ido!=2){
236      t1=0;
237      for(k=0;k<l1;k++){
238        t2=t1;
239        t4=t1<<2;
240        t5=(t6=ido<<1)+t4;
241        for(i=2;i<ido;i+=2){
242          t3=(t2+=2);
243          t4+=2;
244          t5-=2;
245
246          t3+=t0;
247          cr2=wa1[index1+i-2]*cc[t3-1]+wa1[index1+i-1]*cc[t3];
248          ci2=wa1[index1+i-2]*cc[t3]-wa1[index1+i-1]*cc[t3-1];
249          t3+=t0;
250          cr3=wa2[index2+i-2]*cc[t3-1]+wa2[index2+i-1]*cc[t3];
251          ci3=wa2[index2+i-2]*cc[t3]-wa2[index2+i-1]*cc[t3-1];
252          t3+=t0;
253          cr4=wa3[index3+i-2]*cc[t3-1]+wa3[index3+i-1]*cc[t3];
254          ci4=wa3[index3+i-2]*cc[t3]-wa3[index3+i-1]*cc[t3-1];
255
256          tr1=cr2+cr4;
257          tr4=cr4-cr2;
258          ti1=ci2+ci4;
259          ti4=ci2-ci4;
260
261          ti2=cc[t2]+ci3;
262          ti3=cc[t2]-ci3;
263          tr2=cc[t2-1]+cr3;
264          tr3=cc[t2-1]-cr3;
265          
266          ch[t4-1]=tr1+tr2;
267          ch[t4]=ti1+ti2;
268
269          ch[t5-1]=tr3-ti4;
270          ch[t5]=tr4-ti3;
271
272          ch[t4+t6-1]=ti4+tr3;
273          ch[t4+t6]=tr4+ti3;
274
275          ch[t5+t6-1]=tr2-tr1;
276          ch[t5+t6]=ti1-ti2;
277        }
278        t1+=ido;
279      }
280      if((ido&1)!=0)return;
281    }
282
283    t2=(t1=t0+ido-1)+(t0<<1);
284    t3=ido<<2;
285    t4=ido;
286    t5=ido<<1;
287    t6=ido;
288
289    for(k=0;k<l1;k++){
290      ti1=-hsqt2*(cc[t1]+cc[t2]);
291      tr1=hsqt2*(cc[t1]-cc[t2]);
292
293      ch[t4-1]=tr1+cc[t6-1];
294      ch[t4+t5-1]=cc[t6-1]-tr1;
295
296      ch[t4]=ti1-cc[t1+t0];
297      ch[t4+t5]=ti1+cc[t1+t0];
298      
299      t1+=ido;
300      t2+=ido;
301      t4+=t3;
302      t6+=ido;
303    }
304  }
305
306  static void dradfg(int ido,int ip,int l1,int idl1,float[] cc,float[] c1,
307                     float[] c2, float[] ch, float[] ch2, float[] wa, int index){
308    int idij,ipph,i,j,k,l,ic,ik,is;
309    int t0,t1,t2=0,t3,t4,t5,t6,t7,t8,t9,t10;
310    float dc2,ai1,ai2,ar1,ar2,ds2;
311    int nbd;
312    float dcp=0,arg,dsp=0,ar1h,ar2h;
313    int idp2,ipp2;
314  
315    arg=tpi/(float)ip;
316    dcp=(float)Math.cos(arg);
317    dsp=(float)Math.sin(arg);
318    ipph=(ip+1)>>1;
319    ipp2=ip;
320    idp2=ido;
321    nbd=(ido-1)>>1;
322    t0=l1*ido;
323    t10=ip*ido;
324
325    int state=100;
326    loop: while(true){
327      switch(state){
328      case 101:
329        if(ido==1){
330          state=119;
331          break;
332        }
333        for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
334
335        t1=0;
336        for(j=1;j<ip;j++){
337          t1+=t0;
338          t2=t1;
339          for(k=0;k<l1;k++){
340            ch[t2]=c1[t2];
341            t2+=ido;
342          }
343        }
344
345        is=-ido;
346        t1=0;
347        if(nbd>l1){
348          for(j=1;j<ip;j++){
349            t1+=t0;
350            is+=ido;
351            t2= -ido+t1;
352            for(k=0;k<l1;k++){
353              idij=is-1;
354              t2+=ido;
355              t3=t2;
356              for(i=2;i<ido;i+=2){
357                idij+=2;
358                t3+=2;
359                ch[t3-1]=wa[index+idij-1]*c1[t3-1]+wa[index+idij]*c1[t3];
360                ch[t3]=wa[index+idij-1]*c1[t3]-wa[index+idij]*c1[t3-1];
361              }
362            }
363          }
364        }
365        else{
366
367          for(j=1;j<ip;j++){
368            is+=ido;
369            idij=is-1;
370            t1+=t0;
371            t2=t1;
372            for(i=2;i<ido;i+=2){
373              idij+=2;
374              t2+=2;
375              t3=t2;
376              for(k=0;k<l1;k++){
377                ch[t3-1]=wa[index+idij-1]*c1[t3-1]+wa[index+idij]*c1[t3];
378                ch[t3]=wa[index+idij-1]*c1[t3]-wa[index+idij]*c1[t3-1];
379                t3+=ido;
380              }
381            }
382          }
383        }
384
385        t1=0;
386        t2=ipp2*t0;
387        if(nbd<l1){
388          for(j=1;j<ipph;j++){
389            t1+=t0;
390            t2-=t0;
391            t3=t1;
392            t4=t2;
393            for(i=2;i<ido;i+=2){
394              t3+=2;
395              t4+=2;
396              t5=t3-ido;
397              t6=t4-ido;
398              for(k=0;k<l1;k++){
399                t5+=ido;
400                t6+=ido;
401                c1[t5-1]=ch[t5-1]+ch[t6-1];
402                c1[t6-1]=ch[t5]-ch[t6];
403                c1[t5]=ch[t5]+ch[t6];
404                c1[t6]=ch[t6-1]-ch[t5-1];
405              }
406            }
407          }
408        }
409        else{
410          for(j=1;j<ipph;j++){
411            t1+=t0;
412            t2-=t0;
413            t3=t1;
414            t4=t2;
415            for(k=0;k<l1;k++){
416              t5=t3;
417              t6=t4;
418              for(i=2;i<ido;i+=2){
419                t5+=2;
420                t6+=2;
421                c1[t5-1]=ch[t5-1]+ch[t6-1];
422                c1[t6-1]=ch[t5]-ch[t6];
423                c1[t5]=ch[t5]+ch[t6];
424                c1[t6]=ch[t6-1]-ch[t5-1];
425              }
426              t3+=ido;
427              t4+=ido;
428            }
429          }
430        }
431      case 119:
432        for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
433
434        t1=0;
435        t2=ipp2*idl1;
436        for(j=1;j<ipph;j++){
437          t1+=t0;
438          t2-=t0;
439          t3=t1-ido;
440          t4=t2-ido;
441          for(k=0;k<l1;k++){
442            t3+=ido;
443            t4+=ido;
444            c1[t3]=ch[t3]+ch[t4];
445            c1[t4]=ch[t4]-ch[t3];
446          }
447        }
448
449        ar1=1.f;
450        ai1=0.f;
451        t1=0;
452        t2=ipp2*idl1;
453        t3=(ip-1)*idl1;
454        for(l=1;l<ipph;l++){
455          t1+=idl1;
456          t2-=idl1;
457          ar1h=dcp*ar1-dsp*ai1;
458          ai1=dcp*ai1+dsp*ar1;
459          ar1=ar1h;
460          t4=t1;
461          t5=t2;
462          t6=t3;
463          t7=idl1;
464
465          for(ik=0;ik<idl1;ik++){
466            ch2[t4++]=c2[ik]+ar1*c2[t7++];
467            ch2[t5++]=ai1*c2[t6++];
468          }
469
470          dc2=ar1;
471          ds2=ai1;
472          ar2=ar1;
473          ai2=ai1;
474
475          t4=idl1;
476          t5=(ipp2-1)*idl1;
477          for(j=2;j<ipph;j++){
478            t4+=idl1;
479            t5-=idl1;
480
481            ar2h=dc2*ar2-ds2*ai2;
482            ai2=dc2*ai2+ds2*ar2;
483            ar2=ar2h;
484
485            t6=t1;
486            t7=t2;
487            t8=t4;
488            t9=t5;
489            for(ik=0;ik<idl1;ik++){
490              ch2[t6++]+=ar2*c2[t8++];
491              ch2[t7++]+=ai2*c2[t9++];
492            }
493          }
494        }
495        t1=0;
496        for(j=1;j<ipph;j++){
497          t1+=idl1;
498          t2=t1;
499          for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
500        }
501
502        if(ido<l1){
503          state=132;
504          break;
505        }
506
507        t1=0;
508        t2=0;
509        for(k=0;k<l1;k++){
510          t3=t1;
511          t4=t2;
512          for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
513          t1+=ido;
514          t2+=t10;
515        }
516        state=135;
517        break;
518
519      case 132:
520        for(i=0;i<ido;i++){
521          t1=i;
522          t2=i;
523          for(k=0;k<l1;k++){
524            cc[t2]=ch[t1];
525            t1+=ido;
526            t2+=t10;
527          }
528        }
529      case 135:
530        t1=0;
531        t2=ido<<1;
532        t3=0;
533        t4=ipp2*t0;
534        for(j=1;j<ipph;j++){
535          t1+=t2;
536          t3+=t0;
537          t4-=t0;
538
539          t5=t1;
540          t6=t3;
541          t7=t4;
542
543          for(k=0;k<l1;k++){
544            cc[t5-1]=ch[t6];
545            cc[t5]=ch[t7];
546            t5+=t10;
547            t6+=ido;
548            t7+=ido;
549          }
550        }
551
552        if(ido==1)return;
553        if(nbd<l1){
554          state=141;
555          break;
556        }
557
558        t1=-ido;
559        t3=0;
560        t4=0;
561        t5=ipp2*t0;
562        for(j=1;j<ipph;j++){
563          t1+=t2;
564          t3+=t2;
565          t4+=t0;
566          t5-=t0;
567          t6=t1;
568          t7=t3;
569          t8=t4;
570          t9=t5;
571          for(k=0;k<l1;k++){
572            for(i=2;i<ido;i+=2){
573              ic=idp2-i;
574              cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
575              cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
576              cc[i+t7]=ch[i+t8]+ch[i+t9];
577              cc[ic+t6]=ch[i+t9]-ch[i+t8];
578            }
579            t6+=t10;
580            t7+=t10;
581            t8+=ido;
582            t9+=ido;
583          }
584        }
585        return;
586      case 141:
587        t1=-ido;
588        t3=0;
589        t4=0;
590        t5=ipp2*t0;
591        for(j=1;j<ipph;j++){
592          t1+=t2;
593          t3+=t2;
594          t4+=t0;
595          t5-=t0;
596          for(i=2;i<ido;i+=2){
597            t6=idp2+t1-i;
598            t7=i+t3;
599            t8=i+t4;
600            t9=i+t5;
601            for(k=0;k<l1;k++){
602              cc[t7-1]=ch[t8-1]+ch[t9-1];
603              cc[t6-1]=ch[t8-1]-ch[t9-1];
604              cc[t7]=ch[t8]+ch[t9];
605              cc[t6]=ch[t9]-ch[t8];
606              t6+=t10;
607              t7+=t10;
608              t8+=ido;
609              t9+=ido;
610            }
611          }
612        }
613        break loop;
614      }
615    }
616  }
617
618  static void drftf1(int n,float[] c, float[] ch, float[] wa, int[] ifac){
619    int i,k1,l1,l2;
620    int na,kh,nf;
621    int ip,iw,ido,idl1,ix2,ix3;
622
623    nf=ifac[1];
624    na=1;
625    l2=n;
626    iw=n;
627
628    for(k1=0;k1<nf;k1++){
629      kh=nf-k1;
630      ip=ifac[kh+1];
631      l1=l2/ip;
632      ido=n/l2;
633      idl1=ido*l1;
634      iw-=(ip-1)*ido;
635      na=1-na;
636
637      int state=100;
638      loop: while(true){
639        switch(state){
640        case 100:
641          if(ip!=4){
642            state=102;
643            break;
644          }
645
646          ix2=iw+ido;
647          ix3=ix2+ido;
648          if(na!=0)
649            dradf4(ido,l1,ch,c,wa,iw-1,wa,ix2-1,wa,ix3-1);
650          else
651            dradf4(ido,l1,c,ch,wa,iw-1,wa,ix2-1,wa,ix3-1);
652          state=110;
653          break;
654        case 102:
655          if(ip!=2){
656            state=104;
657            break;
658          }
659          if(na!=0){
660            state=103;
661            break;
662          }
663          dradf2(ido,l1,c,ch,wa, iw-1);
664          state=110;
665          break;
666        case 103:
667          dradf2(ido,l1,ch,c,wa, iw-1);
668        case 104:
669          if(ido==1)na=1-na;
670          if(na!=0){
671            state=109;
672            break;
673          }
674          dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa,iw-1);
675          na=1;
676          state=110;
677          break;
678        case 109:
679          dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa,iw-1);
680          na=0;
681        case 110:
682          l2=l1;
683          break loop;
684        }
685      }
686    }
687    if(na==1)return;
688    for(i=0;i<n;i++)c[i]=ch[i];
689  }
690
691  static void dradb2(int ido,int l1,float[] cc,float[] ch,float[] wa1, int index){
692    int i,k,t0,t1,t2,t3,t4,t5,t6;
693    float ti2,tr2;
694
695    t0=l1*ido;
696    
697    t1=0;
698    t2=0;
699    t3=(ido<<1)-1;
700    for(k=0;k<l1;k++){
701      ch[t1]=cc[t2]+cc[t3+t2];
702      ch[t1+t0]=cc[t2]-cc[t3+t2];
703      t2=(t1+=ido)<<1;
704    }
705
706    if(ido<2)return;
707    if(ido!=2){
708      t1=0;
709      t2=0;
710      for(k=0;k<l1;k++){
711        t3=t1;
712        t5=(t4=t2)+(ido<<1);
713        t6=t0+t1;
714        for(i=2;i<ido;i+=2){
715          t3+=2;
716          t4+=2;
717          t5-=2;
718          t6+=2;
719          ch[t3-1]=cc[t4-1]+cc[t5-1];
720          tr2=cc[t4-1]-cc[t5-1];
721          ch[t3]=cc[t4]-cc[t5];
722          ti2=cc[t4]+cc[t5];
723          ch[t6-1]=wa1[index+i-2]*tr2-wa1[index+i-1]*ti2;
724          ch[t6]=wa1[index+i-2]*ti2+wa1[index+i-1]*tr2;
725        }
726        t2=(t1+=ido)<<1;
727      }
728      if((ido%2)==1)return;
729    }
730
731    t1=ido-1;
732    t2=ido-1;
733    for(k=0;k<l1;k++){
734      ch[t1]=cc[t2]+cc[t2];
735      ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
736      t1+=ido;
737      t2+=ido<<1;
738    }
739  }
740
741  static void dradb3(int ido,int l1,float[] cc,float[] ch,
742                     float[] wa1, int index1,
743                     float[] wa2, int index2){
744    int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
745    float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
746    t0=l1*ido;
747
748    t1=0;
749    t2=t0<<1;
750    t3=ido<<1;
751    t4=ido+(ido<<1);
752    t5=0;
753    for(k=0;k<l1;k++){
754      tr2=cc[t3-1]+cc[t3-1];
755      cr2=cc[t5]+(taur*tr2);
756      ch[t1]=cc[t5]+tr2;
757      ci3=taui*(cc[t3]+cc[t3]);
758      ch[t1+t0]=cr2-ci3;
759      ch[t1+t2]=cr2+ci3;
760      t1+=ido;
761      t3+=t4;
762      t5+=t4;
763    }
764
765    if(ido==1)return;
766
767    t1=0;
768    t3=ido<<1;
769    for(k=0;k<l1;k++){
770      t7=t1+(t1<<1);
771      t6=(t5=t7+t3);
772      t8=t1;
773      t10=(t9=t1+t0)+t0;
774
775      for(i=2;i<ido;i+=2){
776        t5+=2;
777        t6-=2;
778        t7+=2;
779        t8+=2;
780        t9+=2;
781        t10+=2;
782        tr2=cc[t5-1]+cc[t6-1];
783        cr2=cc[t7-1]+(taur*tr2);
784        ch[t8-1]=cc[t7-1]+tr2;
785        ti2=cc[t5]-cc[t6];
786        ci2=cc[t7]+(taur*ti2);
787        ch[t8]=cc[t7]+ti2;
788        cr3=taui*(cc[t5-1]-cc[t6-1]);
789        ci3=taui*(cc[t5]+cc[t6]);
790        dr2=cr2-ci3;
791        dr3=cr2+ci3;
792        di2=ci2+cr3;
793        di3=ci2-cr3;
794        ch[t9-1]=wa1[index1+i-2]*dr2-wa1[index1+i-1]*di2;
795        ch[t9]=wa1[index1+i-2]*di2+wa1[index1+i-1]*dr2;
796        ch[t10-1]=wa2[index2+i-2]*dr3-wa2[index2+i-1]*di3;
797        ch[t10]=wa2[index2+i-2]*di3+wa2[index2+i-1]*dr3;
798      }
799      t1+=ido;
800    }
801  }
802
803  static void dradb4(int ido,int l1,float[] cc,float[] ch,
804                     float[] wa1, int index1,
805                     float[] wa2, int index2,
806                     float[] wa3, int index3){
807    int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
808    float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
809    t0=l1*ido;
810  
811    t1=0;
812    t2=ido<<2;
813    t3=0;
814    t6=ido<<1;
815    for(k=0;k<l1;k++){
816      t4=t3+t6;
817      t5=t1;
818      tr3=cc[t4-1]+cc[t4-1];
819      tr4=cc[t4]+cc[t4]; 
820      tr1=cc[t3]-cc[(t4+=t6)-1];
821      tr2=cc[t3]+cc[t4-1];
822      ch[t5]=tr2+tr3;
823      ch[t5+=t0]=tr1-tr4;
824      ch[t5+=t0]=tr2-tr3;
825      ch[t5+=t0]=tr1+tr4;
826      t1+=ido;
827      t3+=t2;
828    }
829
830    if(ido<2)return;
831    if(ido!=2){
832      t1=0;
833      for(k=0;k<l1;k++){
834        t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
835        t7=t1;
836        for(i=2;i<ido;i+=2){
837          t2+=2;
838          t3+=2;
839          t4-=2;
840          t5-=2;
841          t7+=2;
842          ti1=cc[t2]+cc[t5];
843          ti2=cc[t2]-cc[t5];
844          ti3=cc[t3]-cc[t4];
845          tr4=cc[t3]+cc[t4];
846          tr1=cc[t2-1]-cc[t5-1];
847          tr2=cc[t2-1]+cc[t5-1];
848          ti4=cc[t3-1]-cc[t4-1];
849          tr3=cc[t3-1]+cc[t4-1];
850          ch[t7-1]=tr2+tr3;
851          cr3=tr2-tr3;
852          ch[t7]=ti2+ti3;
853          ci3=ti2-ti3;
854          cr2=tr1-tr4;
855          cr4=tr1+tr4;
856          ci2=ti1+ti4;
857          ci4=ti1-ti4;
858
859          ch[(t8=t7+t0)-1]=wa1[index1+i-2]*cr2-wa1[index1+i-1]*ci2;
860          ch[t8]=wa1[index1+i-2]*ci2+wa1[index1+i-1]*cr2;
861          ch[(t8+=t0)-1]=wa2[index2+i-2]*cr3-wa2[index2+i-1]*ci3;
862          ch[t8]=wa2[index2+i-2]*ci3+wa2[index2+i-1]*cr3;
863          ch[(t8+=t0)-1]=wa3[index3+i-2]*cr4-wa3[index3+i-1]*ci4;
864          ch[t8]=wa3[index3+i-2]*ci4+wa3[index3+i-1]*cr4;
865        }
866        t1+=ido;
867      }
868      if(ido%2 == 1)return;
869    }
870
871    t1=ido;
872    t2=ido<<2;
873    t3=ido-1;
874    t4=ido+(ido<<1);
875    for(k=0;k<l1;k++){
876      t5=t3;
877      ti1=cc[t1]+cc[t4];
878      ti2=cc[t4]-cc[t1];
879      tr1=cc[t1-1]-cc[t4-1];
880      tr2=cc[t1-1]+cc[t4-1];
881      ch[t5]=tr2+tr2;
882      ch[t5+=t0]=sqrt2*(tr1-ti1);
883      ch[t5+=t0]=ti2+ti2;
884      ch[t5+=t0]=-sqrt2*(tr1+ti1);
885      
886      t3+=ido;
887      t1+=t2;
888      t4+=t2;
889    }
890  }
891
892  static void dradbg(int ido,int ip,int l1,int idl1,float[] cc,float[] c1,
893                     float[] c2,float[] ch,float[] ch2,float[] wa, int index ){
894
895    int idij,ipph=0,i,j,k,l,ik,is,t0=0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10=0,
896      t11,t12;
897    float dc2,ai1,ai2,ar1,ar2,ds2;
898    int nbd=0;
899    float dcp=0,arg,dsp=0,ar1h,ar2h;
900    int ipp2=0;
901
902    int state=100;
903
904    loop: while(true){
905      switch(state){
906      case 100:
907        t10=ip*ido;
908        t0=l1*ido;
909        arg=tpi/(float)ip;
910        dcp=(float)Math.cos(arg);
911        dsp=(float)Math.sin(arg);
912        nbd=(ido-1)>>>1;
913        ipp2=ip;
914        ipph=(ip+1)>>>1;
915        if(ido<l1){
916          state=103;
917          break;
918        }
919        t1=0;
920        t2=0;
921        for(k=0;k<l1;k++){
922          t3=t1;
923          t4=t2;
924          for(i=0;i<ido;i++){
925            ch[t3]=cc[t4];
926            t3++;
927            t4++;
928          }
929          t1+=ido;
930          t2+=t10;
931        }
932        state=106;
933        break;
934      case 103:
935        t1=0;
936        for(i=0;i<ido;i++){
937          t2=t1;
938          t3=t1;
939          for(k=0;k<l1;k++){
940            ch[t2]=cc[t3];
941            t2+=ido;
942            t3+=t10;
943          }
944          t1++;
945        }
946      case 106:
947        t1=0;
948        t2=ipp2*t0;
949        t7=(t5=ido<<1);
950        for(j=1;j<ipph;j++){
951          t1+=t0;
952          t2-=t0;
953          t3=t1;
954          t4=t2;
955          t6=t5;
956          for(k=0;k<l1;k++){
957            ch[t3]=cc[t6-1]+cc[t6-1];
958            ch[t4]=cc[t6]+cc[t6];
959            t3+=ido;
960            t4+=ido;
961            t6+=t10;
962          }
963          t5+=t7;
964        }
965        if (ido == 1){
966          state=116;
967          break;
968        }
969        if(nbd<l1){
970          state=112;
971          break;
972        }
973
974        t1=0;
975        t2=ipp2*t0;
976        t7=0;
977        for(j=1;j<ipph;j++){
978          t1+=t0;
979          t2-=t0;
980          t3=t1;
981          t4=t2;
982
983          t7+=(ido<<1);
984          t8=t7;
985          for(k=0;k<l1;k++){
986            t5=t3;
987            t6=t4;
988            t9=t8;
989            t11=t8;
990            for(i=2;i<ido;i+=2){
991              t5+=2;
992              t6+=2;
993              t9+=2;
994              t11-=2;
995              ch[t5-1]=cc[t9-1]+cc[t11-1];
996              ch[t6-1]=cc[t9-1]-cc[t11-1];
997              ch[t5]=cc[t9]-cc[t11];
998              ch[t6]=cc[t9]+cc[t11];
999            }
1000            t3+=ido;
1001            t4+=ido;
1002            t8+=t10;
1003          }
1004        }
1005        state=116;
1006        break;
1007      case 112:
1008        t1=0;
1009        t2=ipp2*t0;
1010        t7=0;
1011        for(j=1;j<ipph;j++){
1012          t1+=t0;
1013          t2-=t0;
1014          t3=t1;
1015          t4=t2;
1016          t7+=(ido<<1);
1017          t8=t7;
1018          t9=t7;
1019          for(i=2;i<ido;i+=2){
1020            t3+=2;
1021            t4+=2;
1022            t8+=2;
1023            t9-=2;
1024            t5=t3;
1025            t6=t4;
1026            t11=t8;
1027            t12=t9;
1028            for(k=0;k<l1;k++){
1029              ch[t5-1]=cc[t11-1]+cc[t12-1];
1030              ch[t6-1]=cc[t11-1]-cc[t12-1];
1031              ch[t5]=cc[t11]-cc[t12];
1032              ch[t6]=cc[t11]+cc[t12];
1033              t5+=ido;
1034              t6+=ido;
1035              t11+=t10;
1036              t12+=t10;
1037            }
1038          }
1039        }
1040      case 116:
1041        ar1=1.f;
1042        ai1=0.f;
1043        t1=0;
1044        t9=(t2=ipp2*idl1);
1045        t3=(ip-1)*idl1;
1046        for(l=1;l<ipph;l++){
1047          t1+=idl1;
1048          t2-=idl1;
1049          
1050          ar1h=dcp*ar1-dsp*ai1;
1051          ai1=dcp*ai1+dsp*ar1;
1052          ar1=ar1h;
1053          t4=t1;
1054          t5=t2;
1055          t6=0;
1056          t7=idl1;
1057          t8=t3;
1058          for(ik=0;ik<idl1;ik++){
1059            c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
1060            c2[t5++]=ai1*ch2[t8++];
1061          }
1062          dc2=ar1;
1063          ds2=ai1;
1064          ar2=ar1;
1065          ai2=ai1;
1066
1067          t6=idl1;
1068          t7=t9-idl1;
1069          for(j=2;j<ipph;j++){
1070            t6+=idl1;
1071            t7-=idl1;
1072            ar2h=dc2*ar2-ds2*ai2;
1073            ai2=dc2*ai2+ds2*ar2;
1074            ar2=ar2h;
1075            t4=t1;
1076            t5=t2;
1077            t11=t6;
1078            t12=t7;
1079            for(ik=0;ik<idl1;ik++){
1080              c2[t4++]+=ar2*ch2[t11++];
1081              c2[t5++]+=ai2*ch2[t12++];
1082            }
1083          }
1084        }
1085
1086        t1=0;
1087        for(j=1;j<ipph;j++){
1088          t1+=idl1;
1089          t2=t1;
1090          for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
1091        }
1092
1093        t1=0;
1094        t2=ipp2*t0;
1095        for(j=1;j<ipph;j++){
1096          t1+=t0;
1097          t2-=t0;
1098          t3=t1;
1099          t4=t2;
1100          for(k=0;k<l1;k++){
1101            ch[t3]=c1[t3]-c1[t4];
1102            ch[t4]=c1[t3]+c1[t4];
1103            t3+=ido;
1104            t4+=ido;
1105          }
1106        }
1107
1108        if(ido==1){
1109          state=132;
1110          break;
1111        }
1112        if(nbd<l1){
1113          state=128;
1114          break;
1115        }
1116
1117        t1=0;
1118        t2=ipp2*t0;
1119        for(j=1;j<ipph;j++){
1120          t1+=t0;
1121          t2-=t0;
1122          t3=t1;
1123          t4=t2;
1124          for(k=0;k<l1;k++){
1125            t5=t3;
1126            t6=t4;
1127            for(i=2;i<ido;i+=2){
1128              t5+=2;
1129              t6+=2;
1130              ch[t5-1]=c1[t5-1]-c1[t6];
1131              ch[t6-1]=c1[t5-1]+c1[t6];
1132              ch[t5]=c1[t5]+c1[t6-1];
1133              ch[t6]=c1[t5]-c1[t6-1];
1134            }
1135            t3+=ido;
1136            t4+=ido;
1137          }
1138        }
1139        state=132;
1140        break;
1141      case 128:
1142        t1=0;
1143        t2=ipp2*t0;
1144        for(j=1;j<ipph;j++){
1145          t1+=t0;
1146          t2-=t0;
1147          t3=t1;
1148          t4=t2;
1149          for(i=2;i<ido;i+=2){
1150            t3+=2;
1151            t4+=2;
1152            t5=t3;
1153            t6=t4;
1154            for(k=0;k<l1;k++){
1155              ch[t5-1]=c1[t5-1]-c1[t6];
1156              ch[t6-1]=c1[t5-1]+c1[t6];
1157              ch[t5]=c1[t5]+c1[t6-1];
1158              ch[t6]=c1[t5]-c1[t6-1];
1159              t5+=ido;
1160              t6+=ido;
1161            }
1162          }
1163        }
1164      case 132:
1165        if(ido==1)return;
1166
1167        for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
1168
1169        t1=0;
1170        for(j=1;j<ip;j++){
1171          t2=(t1+=t0);
1172          for(k=0;k<l1;k++){
1173            c1[t2]=ch[t2];
1174            t2+=ido;
1175          }
1176        }
1177
1178        if(nbd>l1){
1179          state=139;
1180          break;
1181        }
1182
1183        is= -ido-1;
1184        t1=0;
1185        for(j=1;j<ip;j++){
1186          is+=ido;
1187          t1+=t0;
1188          idij=is;
1189          t2=t1;
1190          for(i=2;i<ido;i+=2){
1191            t2+=2;
1192            idij+=2;
1193            t3=t2;
1194            for(k=0;k<l1;k++){
1195              c1[t3-1]=wa[index+idij-1]*ch[t3-1]-wa[index+idij]*ch[t3];
1196              c1[t3]=wa[index+idij-1]*ch[t3]+wa[index+idij]*ch[t3-1];
1197              t3+=ido;
1198            }
1199          }
1200        }
1201        return;
1202
1203      case 139:
1204        is= -ido-1;
1205        t1=0;
1206        for(j=1;j<ip;j++){
1207          is+=ido;
1208          t1+=t0;
1209          t2=t1;
1210          for(k=0;k<l1;k++){
1211            idij=is;
1212            t3=t2;
1213            for(i=2;i<ido;i+=2){
1214              idij+=2;
1215              t3+=2;
1216              c1[t3-1]=wa[index+idij-1]*ch[t3-1]-wa[index+idij]*ch[t3];
1217              c1[t3]=wa[index+idij-1]*ch[t3]+wa[index+idij]*ch[t3-1];
1218            }
1219            t2+=ido;
1220          }
1221        }
1222        break loop; 
1223      }
1224    }
1225  }
1226
1227  static void drftb1(int n, float[] c, float[] ch, float[] wa, int index, int[] ifac){
1228    int i,k1,l1,l2=0;
1229    int na;
1230    int nf,ip=0,iw,ix2,ix3,ido=0,idl1=0;
1231
1232    nf=ifac[1];
1233    na=0;
1234    l1=1;
1235    iw=1;
1236
1237    for(k1=0;k1<nf;k1++){
1238      int state=100;
1239      loop: while(true){
1240        switch(state){
1241        case 100:
1242          ip=ifac[k1 + 2];
1243          l2=ip*l1;
1244          ido=n/l2;
1245          idl1=ido*l1;
1246          if(ip!=4){
1247            state=103;
1248            break;
1249          }
1250          ix2=iw+ido;
1251          ix3=ix2+ido;
1252
1253          if(na!=0)
1254            dradb4(ido,l1,ch,c,wa,index+iw-1,wa,index+ix2-1,wa,index+ix3-1);
1255          else
1256            dradb4(ido,l1,c,ch,wa,index+iw-1,wa,index+ix2-1,wa,index+ix3-1);
1257          na=1-na;
1258          state=115;
1259          break;
1260        case 103:
1261          if(ip!=2){
1262            state=106;
1263            break;
1264          }
1265
1266          if(na!=0)
1267            dradb2(ido,l1,ch,c,wa,index+iw-1);
1268          else
1269            dradb2(ido,l1,c,ch,wa,index+iw-1);
1270          na=1-na;
1271          state=115;
1272          break;
1273
1274        case 106:
1275          if(ip!=3){
1276            state=109;
1277            break;
1278          }
1279
1280          ix2=iw+ido;
1281          if(na!=0)
1282            dradb3(ido,l1,ch,c,wa,index+iw-1,wa,index+ix2-1);
1283          else
1284            dradb3(ido,l1,c,ch,wa,index+iw-1,wa,index+ix2-1);
1285          na=1-na;
1286          state=115;
1287          break;
1288        case 109:
1289          //    The radix five case can be translated later.....
1290          //    if(ip!=5)goto L112;
1291          //
1292          //ix2=iw+ido;
1293          //ix3=ix2+ido;
1294          //ix4=ix3+ido;
1295          //if(na!=0)
1296          //  dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
1297          //else
1298          //  dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
1299          //na=1-na;
1300          //state=115; 
1301          //break;
1302          if(na!=0)
1303            dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa,index+iw-1);
1304          else
1305            dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa,index+iw-1);
1306          if(ido==1)na=1-na;
1307
1308        case 115:
1309          l1=l2;
1310          iw+=(ip-1)*ido;
1311          break loop;
1312        }
1313      }
1314    }
1315    if(na==0)return;
1316    for(i=0;i<n;i++)c[i]=ch[i];
1317  }
1318}