jorbis/src/com/jcraft/jorbis/Drft.java

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