1 /*

     2 * Portions Copyright (c) 2009 Nokia Corporation and/or its subsidiary(-ies).

     3 * All rights reserved.

     4 * This component and the accompanying materials are made available

     5 * under the terms of "Eclipse Public License v1.0"

     6 * which accompanies this distribution, and is available

     7 * at the URL "http://www.eclipse.org/legal/epl-v10.html".

     8 *

     9 * Initial Contributors:

    10 * Nokia Corporation - initial contribution.

    11 *

    12 * Contributors:

    13 *

    14 * Description: 

    15 * The original of this file was released into the public domain, see below for details

    16 */

    17 

    18 

    19 /* Notes from RFB: 

    20 

    21    Looks like the user-level routines are:

    22 

    23    Real FFT

    24 

    25    void __ogg_fdrffti(int n, double *wsave, int *ifac)

    26    void __ogg_fdrfftf(int n,double *r,double *wsave,int *ifac)

    27    void __ogg_fdrfftb(int n, double *r, double *wsave, int *ifac)

    28 

    29      __ogg_fdrffti == initialization

    30      __ogg_fdrfftf == forward transform

    31      __ogg_fdrfftb == backward transform

    32 

    33      Parameters are

    34      n == length of sequence

    35      r == sequence to be transformed (input)

    36        == transformed sequence (output)

    37      wsave == work array of length 2n (allocated by caller)

    38      ifac == work array of length 15 (allocated by caller)

    39 

    40    Cosine quarter-wave FFT

    41 

    42    void __ogg_fdcosqi(int n, double *wsave, int *ifac)

    43    void __ogg_fdcosqf(int n,double *x,double *wsave,int *ifac)

    44    void __ogg_fdcosqb(int n,double *x,double *wsave,int *ifac)

    45 */

    46 

    47 /********************************************************************

    48  *                                                                  *

    49  * THIS FILE IS PART OF THE OggSQUISH SOFTWARE CODEC SOURCE CODE.   *

    50  *                                                                  *

    51  ********************************************************************

    52 

    53   file: fft.c

    54   function: Fast discrete Fourier and cosine transforms and inverses

    55   author: Monty <xiphmont@mit.edu>

    56   modifications by: Monty

    57   last modification date: Jul 1 1996

    58 

    59  ********************************************************************/

    60 

    61 /* These Fourier routines were originally based on the Fourier

    62    routines of the same names from the NETLIB bihar and fftpack

    63    fortran libraries developed by Paul N. Swarztrauber at the National

    64    Center for Atmospheric Research in Boulder, CO USA.  They have been

    65    reimplemented in C and optimized in a few ways for OggSquish. */

    66 

    67 /* As the original fortran libraries are public domain, the C Fourier

    68    routines in this file are hereby released to the public domain as

    69    well.  The C routines here produce output exactly equivalent to the

    70    original fortran routines.  Of particular interest are the facts

    71    that (like the original fortran), these routines can work on

    72    arbitrary length vectors that need not be powers of two in

    73    length. */

    74 

    75 #include "openc.h"

    76 #define STIN static

    77 

    78 static void drfti1(int n, double *wa, int *ifac){

    79   static int ntryh[4] = { 4,2,3,5 };

    80   static double tpi = 6.28318530717958647692528676655900577;

    81   double arg,argh,argld,fi;

    82   int ntry=0,i,j=-1;

    83   int k1, l1, l2, ib;

    84   int ld, ii, ip, is, nq, nr;

    85   int ido, ipm, nfm1;

    86   int nl=n;

    87   int nf=0;

    88 

    89  L101:

    90   j++;

    91   if (j < 4)

    92     ntry=ntryh[j];

    93   else

    94     ntry+=2;

    95 

    96  L104:

    97   nq=nl/ntry;

    98   nr=nl-ntry*nq;

    99   if (nr!=0) goto L101;

   100 

   101   nf++;

   102   ifac[nf+1]=ntry;

   103   nl=nq;

   104   if(ntry!=2)goto L107;

   105   if(nf==1)goto L107;

   106 

   107   for (i=1;i<nf;i++){

   108     ib=nf-i+1;

   109     ifac[ib+1]=ifac[ib];

   110   }

   111   ifac[2] = 2;

   112 

   113  L107:

   114   if(nl!=1)goto L104;

   115   ifac[0]=n;

   116   ifac[1]=nf;

   117   argh=tpi/n;

   118   is=0;

   119   nfm1=nf-1;

   120   l1=1;

   121 

   122   if(nfm1==0)return;

   123 

   124   for (k1=0;k1<nfm1;k1++){

   125     ip=ifac[k1+2];

   126     ld=0;

   127     l2=l1*ip;

   128     ido=n/l2;

   129     ipm=ip-1;

   130 

   131     for (j=0;j<ipm;j++){

   132       ld+=l1;

   133       i=is;

   134       argld=(double)ld*argh;

   135       fi=0.;

   136       for (ii=2;ii<ido;ii+=2){

   137 	fi+=1.;

   138 	arg=fi*argld;

   139 	wa[i++]=cos(arg);

   140 	wa[i++]=sin(arg);

   141       }

   142       is+=ido;

   143     }

   144     l1=l2;

   145   }

   146 }

   147 

   148 void __ogg_fdrffti(int n, double *wsave, int *ifac){

   149 

   150   if (n == 1) return;

   151   drfti1(n, wsave+n, ifac);

   152 }

   153 

   154 void __ogg_fdcosqi(int n, double *wsave, int *ifac){

   155   static double pih = 1.57079632679489661923132169163975;

   156   static int k;

   157   static double fk, dt;

   158 

   159   dt=pih/n;

   160   fk=0.;

   161   for(k=0;k<n;k++){

   162     fk+=1.;

   163     wsave[k] = cos(fk*dt);

   164   }

   165 

   166   __ogg_fdrffti(n, wsave+n,ifac);

   167 }

   168 

   169 STIN void dradf2(int ido,int l1,double *cc,double *ch,double *wa1){

   170   int i,k;

   171   double ti2,tr2;

   172   int t0,t1,t2,t3,t4,t5,t6;

   173   

   174   t1=0;

   175   t0=(t2=l1*ido);

   176   t3=ido<<1;

   177   for(k=0;k<l1;k++){

   178     ch[t1<<1]=cc[t1]+cc[t2];

   179     ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];

   180     t1+=ido;

   181     t2+=ido;

   182   }

   183 

   184   if(ido<2)return;

   185   if(ido==2)goto L105;

   186 

   187   t1=0;

   188   t2=t0;

   189   for(k=0;k<l1;k++){

   190     t3=t2;

   191     t4=(t1<<1)+(ido<<1);

   192     t5=t1;

   193     t6=t1+t1;

   194     for(i=2;i<ido;i+=2){

   195       t3+=2;

   196       t4-=2;

   197       t5+=2;

   198       t6+=2;

   199       tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];

   200       ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];

   201       ch[t6]=cc[t5]+ti2;

   202       ch[t4]=ti2-cc[t5];

   203       ch[t6-1]=cc[t5-1]+tr2;

   204       ch[t4-1]=cc[t5-1]-tr2;

   205     }

   206     t1+=ido;

   207     t2+=ido;

   208   }

   209 

   210   if(ido%2==1)return;

   211 

   212  L105:

   213   t3=(t2=(t1=ido)-1);

   214   t2+=t0;

   215   for(k=0;k<l1;k++){

   216     ch[t1]=-cc[t2];

   217     ch[t1-1]=cc[t3];

   218     t1+=ido<<1;

   219     t2+=ido;

   220     t3+=ido;

   221   }

   222 }

   223 

   224 STIN void dradf4(int ido,int l1,double *cc,double *ch,double *wa1,

   225 	    double *wa2,double *wa3){

   226   static double hsqt2 = .70710678118654752440084436210485;

   227   int i,k,t0,t1,t2,t3,t4,t5,t6;

   228   double ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;

   229   t0=l1*ido;

   230   

   231   t1=t0;

   232   t4=t1<<1;

   233   t2=t1+(t1<<1);

   234   t3=0;

   235 

   236   for(k=0;k<l1;k++){

   237     tr1=cc[t1]+cc[t2];

   238     tr2=cc[t3]+cc[t4];

   239     ch[t5=t3<<2]=tr1+tr2;

   240     ch[(ido<<2)+t5-1]=tr2-tr1;

   241     ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];

   242     ch[t5]=cc[t2]-cc[t1];

   243 

   244     t1+=ido;

   245     t2+=ido;

   246     t3+=ido;

   247     t4+=ido;

   248   }

   249 

   250   if(ido<2)return;

   251   if(ido==2)goto L105;

   252 

   253   t1=0;

   254   for(k=0;k<l1;k++){

   255     t2=t1;

   256     t4=t1<<2;

   257     t5=(t6=ido<<1)+t4;

   258     for(i=2;i<ido;i+=2){

   259       t3=(t2+=2);

   260       t4+=2;

   261       t5-=2;

   262 

   263       t3+=t0;

   264       cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];

   265       ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];

   266       t3+=t0;

   267       cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];

   268       ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];

   269       t3+=t0;

   270       cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];

   271       ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];

   272 

   273       tr1=cr2+cr4;

   274       tr4=cr4-cr2;

   275       ti1=ci2+ci4;

   276       ti4=ci2-ci4;

   277       ti2=cc[t2]+ci3;

   278       ti3=cc[t2]-ci3;

   279       tr2=cc[t2-1]+cr3;

   280       tr3=cc[t2-1]-cr3;

   281 

   282       

   283       ch[t4-1]=tr1+tr2;

   284       ch[t4]=ti1+ti2;

   285 

   286       ch[t5-1]=tr3-ti4;

   287       ch[t5]=tr4-ti3;

   288 

   289       ch[t4+t6-1]=ti4+tr3;

   290       ch[t4+t6]=tr4+ti3;

   291 

   292       ch[t5+t6-1]=tr2-tr1;

   293       ch[t5+t6]=ti1-ti2;

   294     }

   295     t1+=ido;

   296   }

   297   if(ido%2==1)return;

   298 

   299  L105:

   300   

   301   t2=(t1=t0+ido-1)+(t0<<1);

   302   t3=ido<<2;

   303   t4=ido;

   304   t5=ido<<1;

   305   t6=ido;

   306 

   307   for(k=0;k<l1;k++){

   308     ti1=-hsqt2*(cc[t1]+cc[t2]);

   309     tr1=hsqt2*(cc[t1]-cc[t2]);

   310     ch[t4-1]=tr1+cc[t6-1];

   311     ch[t4+t5-1]=cc[t6-1]-tr1;

   312     ch[t4]=ti1-cc[t1+t0];

   313     ch[t4+t5]=ti1+cc[t1+t0];

   314     t1+=ido;

   315     t2+=ido;

   316     t4+=t3;

   317     t6+=ido;

   318   }

   319 }

   320 

   321 STIN void dradfg(int ido,int ip,int l1,int idl1,double *cc,double *c1,

   322 			  double *c2,double *ch,double *ch2,double *wa){

   323 

   324   static double tpi=6.28318530717958647692528676655900577;

   325   int idij,ipph,i,j,k,l,ic,ik,is;

   326   int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;

   327   double dc2,ai1,ai2,ar1,ar2,ds2;

   328   int nbd;

   329   double dcp,arg,dsp,ar1h,ar2h;

   330   int idp2,ipp2;

   331   

   332   arg=tpi/(double)ip;

   333   dcp=cos(arg);

   334   dsp=sin(arg);

   335   ipph=(ip+1)>>1;

   336   ipp2=ip;

   337   idp2=ido;

   338   nbd=(ido-1)>>1;

   339   t0=l1*ido;

   340   t10=ip*ido;

   341 

   342   if(ido==1)goto L119;

   343   for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];

   344 

   345   t1=0;

   346   for(j=1;j<ip;j++){

   347     t1+=t0;

   348     t2=t1;

   349     for(k=0;k<l1;k++){

   350       ch[t2]=c1[t2];

   351       t2+=ido;

   352     }

   353   }

   354 

   355   is=-ido;

   356   t1=0;

   357   if(nbd>l1){

   358     for(j=1;j<ip;j++){

   359       t1+=t0;

   360       is+=ido;

   361       t2= -ido+t1;

   362       for(k=0;k<l1;k++){

   363 	idij=is-1;

   364 	t2+=ido;

   365 	t3=t2;

   366 	for(i=2;i<ido;i+=2){

   367 	  idij+=2;

   368 	  t3+=2;

   369 	  ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];

   370 	  ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];

   371 	}

   372       }

   373     }

   374   }else{

   375 

   376     for(j=1;j<ip;j++){

   377       is+=ido;

   378       idij=is-1;

   379       t1+=t0;

   380       t2=t1;

   381       for(i=2;i<ido;i+=2){

   382 	idij+=2;

   383 	t2+=2;

   384 	t3=t2;

   385 	for(k=0;k<l1;k++){

   386 	  ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];

   387 	  ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];

   388 	  t3+=ido;

   389 	}

   390       }

   391     }

   392   }

   393 

   394   t1=0;

   395   t2=ipp2*t0;

   396   if(nbd<l1){

   397     for(j=1;j<ipph;j++){

   398       t1+=t0;

   399       t2-=t0;

   400       t3=t1;

   401       t4=t2;

   402       for(i=2;i<ido;i+=2){

   403 	t3+=2;

   404 	t4+=2;

   405 	t5=t3-ido;

   406 	t6=t4-ido;

   407 	for(k=0;k<l1;k++){

   408 	  t5+=ido;

   409 	  t6+=ido;

   410 	  c1[t5-1]=ch[t5-1]+ch[t6-1];

   411 	  c1[t6-1]=ch[t5]-ch[t6];

   412 	  c1[t5]=ch[t5]+ch[t6];

   413 	  c1[t6]=ch[t6-1]-ch[t5-1];

   414 	}

   415       }

   416     }

   417   }else{

   418     for(j=1;j<ipph;j++){

   419       t1+=t0;

   420       t2-=t0;

   421       t3=t1;

   422       t4=t2;

   423       for(k=0;k<l1;k++){

   424 	t5=t3;

   425 	t6=t4;

   426 	for(i=2;i<ido;i+=2){

   427 	  t5+=2;

   428 	  t6+=2;

   429 	  c1[t5-1]=ch[t5-1]+ch[t6-1];

   430 	  c1[t6-1]=ch[t5]-ch[t6];

   431 	  c1[t5]=ch[t5]+ch[t6];

   432 	  c1[t6]=ch[t6-1]-ch[t5-1];

   433 	}

   434 	t3+=ido;

   435 	t4+=ido;

   436       }

   437     }

   438   }

   439 

   440 L119:

   441   for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];

   442 

   443   t1=0;

   444   t2=ipp2*idl1;

   445   for(j=1;j<ipph;j++){

   446     t1+=t0;

   447     t2-=t0;

   448     t3=t1-ido;

   449     t4=t2-ido;

   450     for(k=0;k<l1;k++){

   451       t3+=ido;

   452       t4+=ido;

   453       c1[t3]=ch[t3]+ch[t4];

   454       c1[t4]=ch[t4]-ch[t3];

   455     }

   456   }

   457 

   458   ar1=1.;

   459   ai1=0.;

   460   t1=0;

   461   t2=ipp2*idl1;

   462   t3=(ip-1)*idl1;

   463   for(l=1;l<ipph;l++){

   464     t1+=idl1;

   465     t2-=idl1;

   466     ar1h=dcp*ar1-dsp*ai1;

   467     ai1=dcp*ai1+dsp*ar1;

   468     ar1=ar1h;

   469     t4=t1;

   470     t5=t2;

   471     t6=t3;

   472     t7=idl1;

   473 

   474     for(ik=0;ik<idl1;ik++){

   475       ch2[t4++]=c2[ik]+ar1*c2[t7++];

   476       ch2[t5++]=ai1*c2[t6++];

   477     }

   478 

   479     dc2=ar1;

   480     ds2=ai1;

   481     ar2=ar1;

   482     ai2=ai1;

   483 

   484     t4=idl1;

   485     t5=(ipp2-1)*idl1;

   486     for(j=2;j<ipph;j++){

   487       t4+=idl1;

   488       t5-=idl1;

   489 

   490       ar2h=dc2*ar2-ds2*ai2;

   491       ai2=dc2*ai2+ds2*ar2;

   492       ar2=ar2h;

   493 

   494       t6=t1;

   495       t7=t2;

   496       t8=t4;

   497       t9=t5;

   498       for(ik=0;ik<idl1;ik++){

   499 	ch2[t6++]+=ar2*c2[t8++];

   500 	ch2[t7++]+=ai2*c2[t9++];

   501       }

   502     }

   503   }

   504 

   505   t1=0;

   506   for(j=1;j<ipph;j++){

   507     t1+=idl1;

   508     t2=t1;

   509     for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];

   510   }

   511 

   512   if(ido<l1)goto L132;

   513 

   514   t1=0;

   515   t2=0;

   516   for(k=0;k<l1;k++){

   517     t3=t1;

   518     t4=t2;

   519     for(i=0;i<ido;i++)cc[t4++]=ch[t3++];

   520     t1+=ido;

   521     t2+=t10;

   522   }

   523 

   524   goto L135;

   525 

   526  L132:

   527   for(i=0;i<ido;i++){

   528     t1=i;

   529     t2=i;

   530     for(k=0;k<l1;k++){

   531       cc[t2]=ch[t1];

   532       t1+=ido;

   533       t2+=t10;

   534     }

   535   }

   536 

   537  L135:

   538   t1=0;

   539   t2=ido<<1;

   540   t3=0;

   541   t4=ipp2*t0;

   542   for(j=1;j<ipph;j++){

   543 

   544     t1+=t2;

   545     t3+=t0;

   546     t4-=t0;

   547 

   548     t5=t1;

   549     t6=t3;

   550     t7=t4;

   551 

   552     for(k=0;k<l1;k++){

   553       cc[t5-1]=ch[t6];

   554       cc[t5]=ch[t7];

   555       t5+=t10;

   556       t6+=ido;

   557       t7+=ido;

   558     }

   559   }

   560 

   561   if(ido==1)return;

   562   if(nbd<l1)goto L141;

   563 

   564   t1=-ido;

   565   t3=0;

   566   t4=0;

   567   t5=ipp2*t0;

   568   for(j=1;j<ipph;j++){

   569     t1+=t2;

   570     t3+=t2;

   571     t4+=t0;

   572     t5-=t0;

   573     t6=t1;

   574     t7=t3;

   575     t8=t4;

   576     t9=t5;

   577     for(k=0;k<l1;k++){

   578       for(i=2;i<ido;i+=2){

   579 	ic=idp2-i;

   580 	cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];

   581 	cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];

   582 	cc[i+t7]=ch[i+t8]+ch[i+t9];

   583 	cc[ic+t6]=ch[i+t9]-ch[i+t8];

   584       }

   585       t6+=t10;

   586       t7+=t10;

   587       t8+=ido;

   588       t9+=ido;

   589     }

   590   }

   591   return;

   592 

   593  L141:

   594 

   595   t1=-ido;

   596   t3=0;

   597   t4=0;

   598   t5=ipp2*t0;

   599   for(j=1;j<ipph;j++){

   600     t1+=t2;

   601     t3+=t2;

   602     t4+=t0;

   603     t5-=t0;

   604     for(i=2;i<ido;i+=2){

   605       t6=idp2+t1-i;

   606       t7=i+t3;

   607       t8=i+t4;

   608       t9=i+t5;

   609       for(k=0;k<l1;k++){

   610 	cc[t7-1]=ch[t8-1]+ch[t9-1];

   611 	cc[t6-1]=ch[t8-1]-ch[t9-1];

   612 	cc[t7]=ch[t8]+ch[t9];

   613 	cc[t6]=ch[t9]-ch[t8];

   614 	t6+=t10;

   615 	t7+=t10;

   616 	t8+=ido;

   617 	t9+=ido;

   618       }

   619     }

   620   }

   621 }

   622 

   623 STIN void drftf1(int n,double *c,double *ch,double *wa,int *ifac){

   624   int i,k1,l1,l2;

   625   int na,kh,nf;

   626   int ip,iw,ido,idl1,ix2,ix3;

   627 

   628   nf=ifac[1];

   629   na=1;

   630   l2=n;

   631   iw=n;

   632 

   633   for(k1=0;k1<nf;k1++){

   634     kh=nf-k1;

   635     ip=ifac[kh+1];

   636     l1=l2/ip;

   637     ido=n/l2;

   638     idl1=ido*l1;

   639     iw-=(ip-1)*ido;

   640     na=1-na;

   641 

   642     if(ip!=4)goto L102;

   643 

   644     ix2=iw+ido;

   645     ix3=ix2+ido;

   646     if(na!=0)

   647       dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);

   648     else

   649       dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);

   650     goto L110;

   651 

   652  L102:

   653     if(ip!=2)goto L104;

   654     if(na!=0)goto L103;

   655 

   656     dradf2(ido,l1,c,ch,wa+iw-1);

   657     goto L110;

   658 

   659   L103:

   660     dradf2(ido,l1,ch,c,wa+iw-1);

   661     goto L110;

   662 

   663   L104:

   664     if(ido==1)na=1-na;

   665     if(na!=0)goto L109;

   666 

   667     dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);

   668     na=1;

   669     goto L110;

   670 

   671   L109:

   672     dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);

   673     na=0;

   674 

   675   L110:

   676     l2=l1;

   677   }

   678 

   679   if(na==1)return;

   680 

   681   for(i=0;i<n;i++)c[i]=ch[i];

   682 }

   683 

   684 void __ogg_fdrfftf(int n,double *r,double *wsave,int *ifac){

   685   if(n==1)return;

   686   drftf1(n,r,wsave,wsave+n,ifac);

   687 }

   688 

   689 STIN void dcsqf1(int n,double *x,double *w,double *xh,int *ifac){

   690   int modn,i,k,kc;

   691   int np2,ns2;

   692   double xim1;

   693 

   694   ns2=(n+1)>>1;

   695   np2=n;

   696 

   697   kc=np2;

   698   for(k=1;k<ns2;k++){

   699     kc--;

   700     xh[k]=x[k]+x[kc];

   701     xh[kc]=x[k]-x[kc];

   702   }

   703 

   704   modn=n%2;

   705   if(modn==0)xh[ns2]=x[ns2]+x[ns2];

   706 

   707   for(k=1;k<ns2;k++){

   708     kc=np2-k;

   709     x[k]=w[k-1]*xh[kc]+w[kc-1]*xh[k];

   710     x[kc]=w[k-1]*xh[k]-w[kc-1]*xh[kc];

   711   }

   712 

   713   if(modn==0)x[ns2]=w[ns2-1]*xh[ns2];

   714 

   715   __ogg_fdrfftf(n,x,xh,ifac);

   716 

   717   for(i=2;i<n;i+=2){

   718     xim1=x[i-1]-x[i];

   719     x[i]=x[i-1]+x[i];

   720     x[i-1]=xim1;

   721   }

   722 }

   723 

   724 void __ogg_fdcosqf(int n,double *x,double *wsave,int *ifac){

   725     static double sqrt2=1.4142135623730950488016887242097;

   726     double tsqx;

   727 

   728   switch(n){

   729   case 0:case 1:

   730     return;

   731   case 2:

   732     tsqx=sqrt2*x[1];

   733     x[1]=x[0]-tsqx;

   734     x[0]+=tsqx;

   735     return;

   736   default:

   737     dcsqf1(n,x,wsave,wsave+n,ifac);

   738     return;

   739   }

   740 }

   741 

   742 STIN void dradb2(int ido,int l1,double *cc,double *ch,double *wa1){

   743   int i,k,t0,t1,t2,t3,t4,t5,t6;

   744   double ti2,tr2;

   745 

   746   t0=l1*ido;

   747   

   748   t1=0;

   749   t2=0;

   750   t3=(ido<<1)-1;

   751   for(k=0;k<l1;k++){

   752     ch[t1]=cc[t2]+cc[t3+t2];

   753     ch[t1+t0]=cc[t2]-cc[t3+t2];

   754     t2=(t1+=ido)<<1;

   755   }

   756 

   757   if(ido<2)return;

   758   if(ido==2)goto L105;

   759 

   760   t1=0;

   761   t2=0;

   762   for(k=0;k<l1;k++){

   763     t3=t1;

   764     t5=(t4=t2)+(ido<<1);

   765     t6=t0+t1;

   766     for(i=2;i<ido;i+=2){

   767       t3+=2;

   768       t4+=2;

   769       t5-=2;

   770       t6+=2;

   771       ch[t3-1]=cc[t4-1]+cc[t5-1];

   772       tr2=cc[t4-1]-cc[t5-1];

   773       ch[t3]=cc[t4]-cc[t5];

   774       ti2=cc[t4]+cc[t5];

   775       ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;

   776       ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;

   777     }

   778     t2=(t1+=ido)<<1;

   779   }

   780 

   781   if(ido%2==1)return;

   782 

   783 L105:

   784   t1=ido-1;

   785   t2=ido-1;

   786   for(k=0;k<l1;k++){

   787     ch[t1]=cc[t2]+cc[t2];

   788     ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);

   789     t1+=ido;

   790     t2+=ido<<1;

   791   }

   792 }

   793 

   794 STIN void dradb3(int ido,int l1,double *cc,double *ch,double *wa1,

   795 			  double *wa2){

   796   static double taur = -.5;

   797   static double taui = .86602540378443864676372317075293618;

   798   int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;

   799   double ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;

   800   t0=l1*ido;

   801 

   802   t1=0;

   803   t2=t0<<1;

   804   t3=ido<<1;

   805   t4=ido+(ido<<1);

   806   t5=0;

   807   for(k=0;k<l1;k++){

   808     tr2=cc[t3-1]+cc[t3-1];

   809     cr2=cc[t5]+(taur*tr2);

   810     ch[t1]=cc[t5]+tr2;

   811     ci3=taui*(cc[t3]+cc[t3]);

   812     ch[t1+t0]=cr2-ci3;

   813     ch[t1+t2]=cr2+ci3;

   814     t1+=ido;

   815     t3+=t4;

   816     t5+=t4;

   817   }

   818 

   819   if(ido==1)return;

   820 

   821   t1=0;

   822   t3=ido<<1;

   823   for(k=0;k<l1;k++){

   824     t7=t1+(t1<<1);

   825     t6=(t5=t7+t3);

   826     t8=t1;

   827     t10=(t9=t1+t0)+t0;

   828 

   829     for(i=2;i<ido;i+=2){

   830       t5+=2;

   831       t6-=2;

   832       t7+=2;

   833       t8+=2;

   834       t9+=2;

   835       t10+=2;

   836       tr2=cc[t5-1]+cc[t6-1];

   837       cr2=cc[t7-1]+(taur*tr2);

   838       ch[t8-1]=cc[t7-1]+tr2;

   839       ti2=cc[t5]-cc[t6];

   840       ci2=cc[t7]+(taur*ti2);

   841       ch[t8]=cc[t7]+ti2;

   842       cr3=taui*(cc[t5-1]-cc[t6-1]);

   843       ci3=taui*(cc[t5]+cc[t6]);

   844       dr2=cr2-ci3;

   845       dr3=cr2+ci3;

   846       di2=ci2+cr3;

   847       di3=ci2-cr3;

   848       ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;

   849       ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;

   850       ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;

   851       ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;

   852     }

   853     t1+=ido;

   854   }

   855 }

   856 

   857 STIN void dradb4(int ido,int l1,double *cc,double *ch,double *wa1,

   858 			  double *wa2,double *wa3){

   859   static double sqrt2=1.4142135623730950488016887242097;

   860   int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;

   861   double ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;

   862   t0=l1*ido;

   863   

   864   t1=0;

   865   t2=ido<<2;

   866   t3=0;

   867   t6=ido<<1;

   868   for(k=0;k<l1;k++){

   869     t4=t3+t6;

   870     t5=t1;

   871     tr3=cc[t4-1]+cc[t4-1];

   872     tr4=cc[t4]+cc[t4]; 

   873     tr1=cc[t3]-cc[(t4+=t6)-1];

   874     tr2=cc[t3]+cc[t4-1];

   875     ch[t5]=tr2+tr3;

   876     ch[t5+=t0]=tr1-tr4;

   877     ch[t5+=t0]=tr2-tr3;

   878     ch[t5+=t0]=tr1+tr4;

   879     t1+=ido;

   880     t3+=t2;

   881   }

   882 

   883   if(ido<2)return;

   884   if(ido==2)goto L105;

   885 

   886   t1=0;

   887   for(k=0;k<l1;k++){

   888     t5=(t4=(t3=(t2=t1<<2)+t6))+t6;

   889     t7=t1;

   890     for(i=2;i<ido;i+=2){

   891       t2+=2;

   892       t3+=2;

   893       t4-=2;

   894       t5-=2;

   895       t7+=2;

   896       ti1=cc[t2]+cc[t5];

   897       ti2=cc[t2]-cc[t5];

   898       ti3=cc[t3]-cc[t4];

   899       tr4=cc[t3]+cc[t4];

   900       tr1=cc[t2-1]-cc[t5-1];

   901       tr2=cc[t2-1]+cc[t5-1];

   902       ti4=cc[t3-1]-cc[t4-1];

   903       tr3=cc[t3-1]+cc[t4-1];

   904       ch[t7-1]=tr2+tr3;

   905       cr3=tr2-tr3;

   906       ch[t7]=ti2+ti3;

   907       ci3=ti2-ti3;

   908       cr2=tr1-tr4;

   909       cr4=tr1+tr4;

   910       ci2=ti1+ti4;

   911       ci4=ti1-ti4;

   912 

   913       ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;

   914       ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;

   915       ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;

   916       ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;

   917       ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;

   918       ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;

   919     }

   920     t1+=ido;

   921   }

   922 

   923   if(ido%2 == 1)return;

   924 

   925  L105:

   926 

   927   t1=ido;

   928   t2=ido<<2;

   929   t3=ido-1;

   930   t4=ido+(ido<<1);

   931   for(k=0;k<l1;k++){

   932     t5=t3;

   933     ti1=cc[t1]+cc[t4];

   934     ti2=cc[t4]-cc[t1];

   935     tr1=cc[t1-1]-cc[t4-1];

   936     tr2=cc[t1-1]+cc[t4-1];

   937     ch[t5]=tr2+tr2;

   938     ch[t5+=t0]=sqrt2*(tr1-ti1);

   939     ch[t5+=t0]=ti2+ti2;

   940     ch[t5+=t0]=-sqrt2*(tr1+ti1);

   941 

   942     t3+=ido;

   943     t1+=t2;

   944     t4+=t2;

   945   }

   946 }

   947 

   948 STIN void dradbg(int ido,int ip,int l1,int idl1,double *cc,double *c1,

   949 	    double *c2,double *ch,double *ch2,double *wa){

   950   static double tpi=6.28318530717958647692528676655900577;

   951   int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,

   952       t11,t12;

   953   double dc2,ai1,ai2,ar1,ar2,ds2;

   954   int nbd;

   955   double dcp,arg,dsp,ar1h,ar2h;

   956   int ipp2;

   957 

   958   t10=ip*ido;

   959   t0=l1*ido;

   960   arg=tpi/(double)ip;

   961   dcp=cos(arg);

   962   dsp=sin(arg);

   963   nbd=(ido-1)>>1;

   964   ipp2=ip;

   965   ipph=(ip+1)>>1;

   966   if(ido<l1)goto L103;

   967   

   968   t1=0;

   969   t2=0;

   970   for(k=0;k<l1;k++){

   971     t3=t1;

   972     t4=t2;

   973     for(i=0;i<ido;i++){

   974       ch[t3]=cc[t4];

   975       t3++;

   976       t4++;

   977     }

   978     t1+=ido;

   979     t2+=t10;

   980   }

   981   goto L106;

   982 

   983  L103:

   984   t1=0;

   985   for(i=0;i<ido;i++){

   986     t2=t1;

   987     t3=t1;

   988     for(k=0;k<l1;k++){

   989       ch[t2]=cc[t3];

   990       t2+=ido;

   991       t3+=t10;

   992     }

   993     t1++;

   994   }

   995 

   996  L106:

   997   t1=0;

   998   t2=ipp2*t0;

   999   t7=(t5=ido<<1);

  1000   for(j=1;j<ipph;j++){

  1001     t1+=t0;

  1002     t2-=t0;

  1003     t3=t1;

  1004     t4=t2;

  1005     t6=t5;

  1006     for(k=0;k<l1;k++){

  1007       ch[t3]=cc[t6-1]+cc[t6-1];

  1008       ch[t4]=cc[t6]+cc[t6];

  1009       t3+=ido;

  1010       t4+=ido;

  1011       t6+=t10;

  1012     }

  1013     t5+=t7;

  1014   }

  1015 

  1016   if (ido == 1)goto L116;

  1017   if(nbd<l1)goto L112;

  1018 

  1019   t1=0;

  1020   t2=ipp2*t0;

  1021   t7=0;

  1022   for(j=1;j<ipph;j++){

  1023     t1+=t0;

  1024     t2-=t0;

  1025     t3=t1;

  1026     t4=t2;

  1027 

  1028     t7+=(ido<<1);

  1029     t8=t7;

  1030     for(k=0;k<l1;k++){

  1031       t5=t3;

  1032       t6=t4;

  1033       t9=t8;

  1034       t11=t8;

  1035       for(i=2;i<ido;i+=2){

  1036 	t5+=2;

  1037 	t6+=2;

  1038 	t9+=2;

  1039 	t11-=2;

  1040 	ch[t5-1]=cc[t9-1]+cc[t11-1];

  1041 	ch[t6-1]=cc[t9-1]-cc[t11-1];

  1042 	ch[t5]=cc[t9]-cc[t11];

  1043 	ch[t6]=cc[t9]+cc[t11];

  1044       }

  1045       t3+=ido;

  1046       t4+=ido;

  1047       t8+=t10;

  1048     }

  1049   }

  1050   goto L116;

  1051 

  1052  L112:

  1053   t1=0;

  1054   t2=ipp2*t0;

  1055   t7=0;

  1056   for(j=1;j<ipph;j++){

  1057     t1+=t0;

  1058     t2-=t0;

  1059     t3=t1;

  1060     t4=t2;

  1061     t7+=(ido<<1);

  1062     t8=t7;

  1063     t9=t7;

  1064     for(i=2;i<ido;i+=2){

  1065       t3+=2;

  1066       t4+=2;

  1067       t8+=2;

  1068       t9-=2;

  1069       t5=t3;

  1070       t6=t4;

  1071       t11=t8;

  1072       t12=t9;

  1073       for(k=0;k<l1;k++){

  1074 	ch[t5-1]=cc[t11-1]+cc[t12-1];

  1075 	ch[t6-1]=cc[t11-1]-cc[t12-1];

  1076 	ch[t5]=cc[t11]-cc[t12];

  1077 	ch[t6]=cc[t11]+cc[t12];

  1078 	t5+=ido;

  1079 	t6+=ido;

  1080 	t11+=t10;

  1081 	t12+=t10;

  1082       }

  1083     }

  1084   }

  1085 

  1086 L116:

  1087   ar1=1.;

  1088   ai1=0.;

  1089   t1=0;

  1090   t9=(t2=ipp2*idl1);

  1091   t3=(ip-1)*idl1;

  1092   for(l=1;l<ipph;l++){

  1093     t1+=idl1;

  1094     t2-=idl1;

  1095 

  1096     ar1h=dcp*ar1-dsp*ai1;

  1097     ai1=dcp*ai1+dsp*ar1;

  1098     ar1=ar1h;

  1099     t4=t1;

  1100     t5=t2;

  1101     t6=0;

  1102     t7=idl1;

  1103     t8=t3;

  1104     for(ik=0;ik<idl1;ik++){

  1105       c2[t4++]=ch2[t6++]+ar1*ch2[t7++];

  1106       c2[t5++]=ai1*ch2[t8++];

  1107     }

  1108     dc2=ar1;

  1109     ds2=ai1;

  1110     ar2=ar1;

  1111     ai2=ai1;

  1112 

  1113     t6=idl1;

  1114     t7=t9-idl1;

  1115     for(j=2;j<ipph;j++){

  1116       t6+=idl1;

  1117       t7-=idl1;

  1118       ar2h=dc2*ar2-ds2*ai2;

  1119       ai2=dc2*ai2+ds2*ar2;

  1120       ar2=ar2h;

  1121       t4=t1;

  1122       t5=t2;

  1123       t11=t6;

  1124       t12=t7;

  1125       for(ik=0;ik<idl1;ik++){

  1126 	c2[t4++]+=ar2*ch2[t11++];

  1127 	c2[t5++]+=ai2*ch2[t12++];

  1128       }

  1129     }

  1130   }

  1131 

  1132   t1=0;

  1133   for(j=1;j<ipph;j++){

  1134     t1+=idl1;