Update copyright headers
[qt:qt.git] / demos / mobile / guitartuner / src / fftpack.c
1 /****************************************************************************
2 **
3 ** Copyright (C) 2015 The Qt Company Ltd.
4 ** Contact: http://www.qt.io/licensing/
5 **
6 ** This file is part of the QtDeclarative module of the Qt Toolkit.
7 **
8 ** $QT_BEGIN_LICENSE:BSD$
9 ** You may use this file under the terms of the BSD license as follows:
10 **
11 ** "Redistribution and use in source and binary forms, with or without
12 ** modification, are permitted provided that the following conditions are
13 ** met:
14 **   * Redistributions of source code must retain the above copyright
15 **     notice, this list of conditions and the following disclaimer.
16 **   * Redistributions in binary form must reproduce the above copyright
17 **     notice, this list of conditions and the following disclaimer in
18 **     the documentation and/or other materials provided with the
19 **     distribution.
20 **   * Neither the name of The Qt Company Ltd nor the names of its
21 **     contributors may be used to endorse or promote products derived
22 **     from this software without specific prior written permission.
23 **
24 **
25 ** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
26 ** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
27 ** LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
28 ** A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
29 ** OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
30 ** SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
31 ** LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
32 ** DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
33 ** THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
34 ** (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
35 ** OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE."
36 **
37 ** $QT_END_LICENSE$
38 **
39 ****************************************************************************/
40
41 /********************************************************************
42
43   The routines in this file are from http://www.netlib.org/fftpack/.
44   According to the comments in the original file (which are provided below)
45   and according to the wikipedia article about the FFTPACK[1],
46   they are released as public domain.
47
48   For information about the author of the original, Fortran routines,
49   and the author of the C translation these routines,
50   see the comments below.
51
52   [1] http://en.wikipedia.org/wiki/FFTPACK, referenced 2010-12-21.
53
54  ********************************************************************/
55
56 /********************************************************************
57
58   file: fft.c
59   function: Fast discrete Fourier and cosine transforms and inverses
60   author: Monty <xiphmont@mit.edu>
61   modifications by: Monty
62   last modification date: Jul 1 1996
63
64  ********************************************************************/
65
66 /* These Fourier routines were originally based on the Fourier
67    routines of the same names from the NETLIB bihar and fftpack
68    fortran libraries developed by Paul N. Swarztrauber at the National
69    Center for Atmospheric Research in Boulder, CO USA.  They have been
70    reimplemented in C and optimized in a few ways for OggSquish. */
71
72 /* As the original fortran libraries are public domain, the C Fourier
73    routines in this file are hereby released to the public domain as
74    well.  The C routines here produce output exactly equivalent to the
75    original fortran routines.  Of particular interest are the facts
76    that (like the original fortran), these routines can work on
77    arbitrary length vectors that need not be powers of two in
78    length. */
79
80 #include <math.h>
81
82 __STATIC void drfti1(int n, float *wa, int *ifac){
83   static int ntryh[4] = { 4,2,3,5 };
84   static float tpi = 6.28318530717958647692528676655900577;
85   float arg,argh,argld,fi;
86   int ntry=0,i,j=-1;
87   int k1, l1, l2, ib;
88   int ld, ii, ip, is, nq, nr;
89   int ido, ipm, nfm1;
90   int nl=n;
91   int nf=0;
92
93  L101:
94   j++;
95   if (j < 4)
96     ntry=ntryh[j];
97   else
98     ntry+=2;
99
100  L104:
101   nq=nl/ntry;
102   nr=nl-ntry*nq;
103   if (nr!=0) goto L101;
104
105   nf++;
106   ifac[nf+1]=ntry;
107   nl=nq;
108   if (ntry!=2) goto L107;
109   if (nf==1) goto L107;
110
111   for (i=1;i<nf;i++){
112     ib=nf-i+1;
113     ifac[ib+1]=ifac[ib];
114   }
115   ifac[2] = 2;
116
117  L107:
118   if (nl!=1) goto L104;
119   ifac[0]=n;
120   ifac[1]=nf;
121   argh=tpi/n;
122   is=0;
123   nfm1=nf-1;
124   l1=1;
125
126   if (nfm1==0) return;
127
128   for (k1=0;k1<nfm1;k1++){
129     ip=ifac[k1+2];
130     ld=0;
131     l2=l1*ip;
132     ido=n/l2;
133     ipm=ip-1;
134
135     for (j=0;j<ipm;j++){
136       ld+=l1;
137       i=is;
138       argld=(float)ld*argh;
139       fi=0.;
140       for (ii=2;ii<ido;ii+=2){
141     fi+=1.;
142     arg=fi*argld;
143     wa[i++]=cos(arg);
144     wa[i++]=sin(arg);
145       }
146       is+=ido;
147     }
148     l1=l2;
149   }
150 }
151
152 void __ogg_fdrffti(int n, float *wsave, int *ifac){
153
154   if (n == 1) return;
155   drfti1(n, wsave+n, ifac);
156 }
157
158 void __ogg_fdcosqi(int n, float *wsave, int *ifac){
159   static float pih = 1.57079632679489661923132169163975;
160   static int k;
161   static float fk, dt;
162
163   dt=pih/n;
164   fk=0.;
165   for (k=0;k<n;k++){
166     fk+=1.;
167     wsave[k] = cos(fk*dt);
168   }
169
170   __ogg_fdrffti(n, wsave+n,ifac);
171 }
172
173 STIN void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){
174   int i,k;
175   float ti2,tr2;
176   int t0,t1,t2,t3,t4,t5,t6;
177
178   t1=0;
179   t0=(t2=l1*ido);
180   t3=ido<<1;
181   for (k=0;k<l1;k++){
182     ch[t1<<1]=cc[t1]+cc[t2];
183     ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
184     t1+=ido;
185     t2+=ido;
186   }
187
188   if (ido<2) return;
189   if (ido==2) goto L105;
190
191   t1=0;
192   t2=t0;
193   for (k=0;k<l1;k++){
194     t3=t2;
195     t4=(t1<<1)+(ido<<1);
196     t5=t1;
197     t6=t1+t1;
198     for (i=2;i<ido;i+=2){
199       t3+=2;
200       t4-=2;
201       t5+=2;
202       t6+=2;
203       tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
204       ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
205       ch[t6]=cc[t5]+ti2;
206       ch[t4]=ti2-cc[t5];
207       ch[t6-1]=cc[t5-1]+tr2;
208       ch[t4-1]=cc[t5-1]-tr2;
209     }
210     t1+=ido;
211     t2+=ido;
212   }
213
214   if (ido%2==1) return;
215
216  L105:
217   t3=(t2=(t1=ido)-1);
218   t2+=t0;
219   for (k=0;k<l1;k++){
220     ch[t1]=-cc[t2];
221     ch[t1-1]=cc[t3];
222     t1+=ido<<1;
223     t2+=ido;
224     t3+=ido;
225   }
226 }
227
228 STIN void dradf4(int ido,int l1,float *cc,float *ch,float *wa1,
229         float *wa2,float *wa3){
230   static float hsqt2 = .70710678118654752440084436210485;
231   int i,k,t0,t1,t2,t3,t4,t5,t6;
232   float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
233   t0=l1*ido;
234
235   t1=t0;
236   t4=t1<<1;
237   t2=t1+(t1<<1);
238   t3=0;
239
240   for (k=0;k<l1;k++){
241     tr1=cc[t1]+cc[t2];
242     tr2=cc[t3]+cc[t4];
243     ch[t5=t3<<2]=tr1+tr2;
244     ch[(ido<<2)+t5-1]=tr2-tr1;
245     ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
246     ch[t5]=cc[t2]-cc[t1];
247
248     t1+=ido;
249     t2+=ido;
250     t3+=ido;
251     t4+=ido;
252   }
253
254   if (ido<2) return;
255   if (ido==2) goto L105;
256
257   t1=0;
258   for (k=0;k<l1;k++){
259     t2=t1;
260     t4=t1<<2;
261     t5=(t6=ido<<1)+t4;
262     for (i=2;i<ido;i+=2){
263       t3=(t2+=2);
264       t4+=2;
265       t5-=2;
266
267       t3+=t0;
268       cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
269       ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
270       t3+=t0;
271       cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];
272       ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];
273       t3+=t0;
274       cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];
275       ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];
276
277       tr1=cr2+cr4;
278       tr4=cr4-cr2;
279       ti1=ci2+ci4;
280       ti4=ci2-ci4;
281       ti2=cc[t2]+ci3;
282       ti3=cc[t2]-ci3;
283       tr2=cc[t2-1]+cr3;
284       tr3=cc[t2-1]-cr3;
285
286
287       ch[t4-1]=tr1+tr2;
288       ch[t4]=ti1+ti2;
289
290       ch[t5-1]=tr3-ti4;
291       ch[t5]=tr4-ti3;
292
293       ch[t4+t6-1]=ti4+tr3;
294       ch[t4+t6]=tr4+ti3;
295
296       ch[t5+t6-1]=tr2-tr1;
297       ch[t5+t6]=ti1-ti2;
298     }
299     t1+=ido;
300   }
301   if (ido%2==1) return;
302
303  L105:
304
305   t2=(t1=t0+ido-1)+(t0<<1);
306   t3=ido<<2;
307   t4=ido;
308   t5=ido<<1;
309   t6=ido;
310
311   for (k=0;k<l1;k++){
312     ti1=-hsqt2*(cc[t1]+cc[t2]);
313     tr1=hsqt2*(cc[t1]-cc[t2]);
314     ch[t4-1]=tr1+cc[t6-1];
315     ch[t4+t5-1]=cc[t6-1]-tr1;
316     ch[t4]=ti1-cc[t1+t0];
317     ch[t4+t5]=ti1+cc[t1+t0];
318     t1+=ido;
319     t2+=ido;
320     t4+=t3;
321     t6+=ido;
322   }
323 }
324
325 STIN void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
326               float *c2,float *ch,float *ch2,float *wa){
327
328   static float tpi=6.28318530717958647692528676655900577;
329   int idij,ipph,i,j,k,l,ic,ik,is;
330   int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
331   float dc2,ai1,ai2,ar1,ar2,ds2;
332   int nbd;
333   float dcp,arg,dsp,ar1h,ar2h;
334   int idp2,ipp2;
335
336   arg=tpi/(float)ip;
337   dcp=cos(arg);
338   dsp=sin(arg);
339   ipph=(ip+1)>>1;
340   ipp2=ip;
341   idp2=ido;
342   nbd=(ido-1)>>1;
343   t0=l1*ido;
344   t10=ip*ido;
345
346   if (ido==1) goto L119;
347   for (ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
348
349   t1=0;
350   for (j=1;j<ip;j++){
351     t1+=t0;
352     t2=t1;
353     for (k=0;k<l1;k++){
354       ch[t2]=c1[t2];
355       t2+=ido;
356     }
357   }
358
359   is=-ido;
360   t1=0;
361   if (nbd>l1){
362     for (j=1;j<ip;j++){
363       t1+=t0;
364       is+=ido;
365       t2= -ido+t1;
366       for (k=0;k<l1;k++){
367     idij=is-1;
368     t2+=ido;
369     t3=t2;
370     for (i=2;i<ido;i+=2){
371       idij+=2;
372       t3+=2;
373       ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
374       ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
375     }
376       }
377     }
378   }else{
379
380     for (j=1;j<ip;j++){
381       is+=ido;
382       idij=is-1;
383       t1+=t0;
384       t2=t1;
385       for (i=2;i<ido;i+=2){
386     idij+=2;
387     t2+=2;
388     t3=t2;
389     for (k=0;k<l1;k++){
390       ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
391       ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
392       t3+=ido;
393     }
394       }
395     }
396   }
397
398   t1=0;
399   t2=ipp2*t0;
400   if (nbd<l1){
401     for (j=1;j<ipph;j++){
402       t1+=t0;
403       t2-=t0;
404       t3=t1;
405       t4=t2;
406       for (i=2;i<ido;i+=2){
407     t3+=2;
408     t4+=2;
409     t5=t3-ido;
410     t6=t4-ido;
411     for (k=0;k<l1;k++){
412       t5+=ido;
413       t6+=ido;
414       c1[t5-1]=ch[t5-1]+ch[t6-1];
415       c1[t6-1]=ch[t5]-ch[t6];
416       c1[t5]=ch[t5]+ch[t6];
417       c1[t6]=ch[t6-1]-ch[t5-1];
418     }
419       }
420     }
421   }else{
422     for (j=1;j<ipph;j++){
423       t1+=t0;
424       t2-=t0;
425       t3=t1;
426       t4=t2;
427       for (k=0;k<l1;k++){
428     t5=t3;
429     t6=t4;
430     for (i=2;i<ido;i+=2){
431       t5+=2;
432       t6+=2;
433       c1[t5-1]=ch[t5-1]+ch[t6-1];
434       c1[t6-1]=ch[t5]-ch[t6];
435       c1[t5]=ch[t5]+ch[t6];
436       c1[t6]=ch[t6-1]-ch[t5-1];
437     }
438     t3+=ido;
439     t4+=ido;
440       }
441     }
442   }
443
444 L119:
445   for (ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
446
447   t1=0;
448   t2=ipp2*idl1;
449   for (j=1;j<ipph;j++){
450     t1+=t0;
451     t2-=t0;
452     t3=t1-ido;
453     t4=t2-ido;
454     for (k=0;k<l1;k++){
455       t3+=ido;
456       t4+=ido;
457       c1[t3]=ch[t3]+ch[t4];
458       c1[t4]=ch[t4]-ch[t3];
459     }
460   }
461
462   ar1=1.;
463   ai1=0.;
464   t1=0;
465   t2=ipp2*idl1;
466   t3=(ip-1)*idl1;
467   for (l=1;l<ipph;l++){
468     t1+=idl1;
469     t2-=idl1;
470     ar1h=dcp*ar1-dsp*ai1;
471     ai1=dcp*ai1+dsp*ar1;
472     ar1=ar1h;
473     t4=t1;
474     t5=t2;
475     t6=t3;
476     t7=idl1;
477
478     for (ik=0;ik<idl1;ik++){
479       ch2[t4++]=c2[ik]+ar1*c2[t7++];
480       ch2[t5++]=ai1*c2[t6++];
481     }
482
483     dc2=ar1;
484     ds2=ai1;
485     ar2=ar1;
486     ai2=ai1;
487
488     t4=idl1;
489     t5=(ipp2-1)*idl1;
490     for (j=2;j<ipph;j++){
491       t4+=idl1;
492       t5-=idl1;
493
494       ar2h=dc2*ar2-ds2*ai2;
495       ai2=dc2*ai2+ds2*ar2;
496       ar2=ar2h;
497
498       t6=t1;
499       t7=t2;
500       t8=t4;
501       t9=t5;
502       for (ik=0;ik<idl1;ik++){
503     ch2[t6++]+=ar2*c2[t8++];
504     ch2[t7++]+=ai2*c2[t9++];
505       }
506     }
507   }
508
509   t1=0;
510   for (j=1;j<ipph;j++){
511     t1+=idl1;
512     t2=t1;
513     for (ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
514   }
515
516   if (ido<l1) goto L132;
517
518   t1=0;
519   t2=0;
520   for (k=0;k<l1;k++){
521     t3=t1;
522     t4=t2;
523     for (i=0;i<ido;i++)cc[t4++]=ch[t3++];
524     t1+=ido;
525     t2+=t10;
526   }
527
528   goto L135;
529
530  L132:
531   for (i=0;i<ido;i++){
532     t1=i;
533     t2=i;
534     for (k=0;k<l1;k++){
535       cc[t2]=ch[t1];
536       t1+=ido;
537       t2+=t10;
538     }
539   }
540
541  L135:
542   t1=0;
543   t2=ido<<1;
544   t3=0;
545   t4=ipp2*t0;
546   for (j=1;j<ipph;j++){
547
548     t1+=t2;
549     t3+=t0;
550     t4-=t0;
551
552     t5=t1;
553     t6=t3;
554     t7=t4;
555
556     for (k=0;k<l1;k++){
557       cc[t5-1]=ch[t6];
558       cc[t5]=ch[t7];
559       t5+=t10;
560       t6+=ido;
561       t7+=ido;
562     }
563   }
564
565   if (ido==1) return;
566   if (nbd<l1) goto L141;
567
568   t1=-ido;
569   t3=0;
570   t4=0;
571   t5=ipp2*t0;
572   for (j=1;j<ipph;j++){
573     t1+=t2;
574     t3+=t2;
575     t4+=t0;
576     t5-=t0;
577     t6=t1;
578     t7=t3;
579     t8=t4;
580     t9=t5;
581     for (k=0;k<l1;k++){
582       for (i=2;i<ido;i+=2){
583     ic=idp2-i;
584     cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
585     cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
586     cc[i+t7]=ch[i+t8]+ch[i+t9];
587     cc[ic+t6]=ch[i+t9]-ch[i+t8];
588       }
589       t6+=t10;
590       t7+=t10;
591       t8+=ido;
592       t9+=ido;
593     }
594   }
595   return;
596
597  L141:
598
599   t1=-ido;
600   t3=0;
601   t4=0;
602   t5=ipp2*t0;
603   for (j=1;j<ipph;j++){
604     t1+=t2;
605     t3+=t2;
606     t4+=t0;
607     t5-=t0;
608     for (i=2;i<ido;i+=2){
609       t6=idp2+t1-i;
610       t7=i+t3;
611       t8=i+t4;
612       t9=i+t5;
613       for (k=0;k<l1;k++){
614     cc[t7-1]=ch[t8-1]+ch[t9-1];
615     cc[t6-1]=ch[t8-1]-ch[t9-1];
616     cc[t7]=ch[t8]+ch[t9];
617     cc[t6]=ch[t9]-ch[t8];
618     t6+=t10;
619     t7+=t10;
620     t8+=ido;
621     t9+=ido;
622       }
623     }
624   }
625 }
626
627 STIN void drftf1(int n,float *c,float *ch,float *wa,int *ifac){
628   int i,k1,l1,l2;
629   int na,kh,nf;
630   int ip,iw,ido,idl1,ix2,ix3;
631
632   nf=ifac[1];
633   na=1;
634   l2=n;
635   iw=n;
636
637   for (k1=0;k1<nf;k1++){
638     kh=nf-k1;
639     ip=ifac[kh+1];
640     l1=l2/ip;
641     ido=n/l2;
642     idl1=ido*l1;
643     iw-=(ip-1)*ido;
644     na=1-na;
645
646     if (ip!=4) goto L102;
647
648     ix2=iw+ido;
649     ix3=ix2+ido;
650     if (na!=0)
651       dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
652     else
653       dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
654     goto L110;
655
656  L102:
657     if (ip!=2) goto L104;
658     if (na!=0) goto L103;
659
660     dradf2(ido,l1,c,ch,wa+iw-1);
661     goto L110;
662
663   L103:
664     dradf2(ido,l1,ch,c,wa+iw-1);
665     goto L110;
666
667   L104:
668     if (ido==1)na=1-na;
669     if (na!=0) goto L109;
670
671     dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
672     na=1;
673     goto L110;
674
675   L109:
676     dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
677     na=0;
678
679   L110:
680     l2=l1;
681   }
682
683   if (na==1) return;
684
685   for (i=0;i<n;i++)c[i]=ch[i];
686 }
687
688 void __ogg_fdrfftf(int n,float *r,float *wsave,int *ifac){
689   if (n==1) return;
690   drftf1(n,r,wsave,wsave+n,ifac);
691 }
692
693 STIN void dcsqf1(int n,float *x,float *w,float *xh,int *ifac){
694   int modn,i,k,kc;
695   int np2,ns2;
696   float xim1;
697
698   ns2=(n+1)>>1;
699   np2=n;
700
701   kc=np2;
702   for (k=1;k<ns2;k++){
703     kc--;
704     xh[k]=x[k]+x[kc];
705     xh[kc]=x[k]-x[kc];
706   }
707
708   modn=n%2;
709   if (modn==0)xh[ns2]=x[ns2]+x[ns2];
710
711   for (k=1;k<ns2;k++){
712     kc=np2-k;
713     x[k]=w[k-1]*xh[kc]+w[kc-1]*xh[k];
714     x[kc]=w[k-1]*xh[k]-w[kc-1]*xh[kc];
715   }
716
717   if (modn==0)x[ns2]=w[ns2-1]*xh[ns2];
718
719   __ogg_fdrfftf(n,x,xh,ifac);
720
721   for (i=2;i<n;i+=2){
722     xim1=x[i-1]-x[i];
723     x[i]=x[i-1]+x[i];
724     x[i-1]=xim1;
725   }
726 }
727
728 void __ogg_fdcosqf(int n,float *x,float *wsave,int *ifac){
729     static float sqrt2=1.4142135623730950488016887242097;
730     float tsqx;
731
732   switch (n){
733   case 0:case 1:
734     return;
735   case 2:
736     tsqx=sqrt2*x[1];
737     x[1]=x[0]-tsqx;
738     x[0]+=tsqx;
739     return;
740   default:
741     dcsqf1(n,x,wsave,wsave+n,ifac);
742     return;
743   }
744 }
745
746 STIN void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){
747   int i,k,t0,t1,t2,t3,t4,t5,t6;
748   float ti2,tr2;
749
750   t0=l1*ido;
751
752   t1=0;
753   t2=0;
754   t3=(ido<<1)-1;
755   for (k=0;k<l1;k++){
756     ch[t1]=cc[t2]+cc[t3+t2];
757     ch[t1+t0]=cc[t2]-cc[t3+t2];
758     t2=(t1+=ido)<<1;
759   }
760
761   if (ido<2) return;
762   if (ido==2) goto L105;
763
764   t1=0;
765   t2=0;
766   for (k=0;k<l1;k++){
767     t3=t1;
768     t5=(t4=t2)+(ido<<1);
769     t6=t0+t1;
770     for (i=2;i<ido;i+=2){
771       t3+=2;
772       t4+=2;
773       t5-=2;
774       t6+=2;
775       ch[t3-1]=cc[t4-1]+cc[t5-1];
776       tr2=cc[t4-1]-cc[t5-1];
777       ch[t3]=cc[t4]-cc[t5];
778       ti2=cc[t4]+cc[t5];
779       ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
780       ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
781     }
782     t2=(t1+=ido)<<1;
783   }
784
785   if (ido%2==1) return;
786
787 L105:
788   t1=ido-1;
789   t2=ido-1;
790   for (k=0;k<l1;k++){
791     ch[t1]=cc[t2]+cc[t2];
792     ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
793     t1+=ido;
794     t2+=ido<<1;
795   }
796 }
797
798 STIN void dradb3(int ido,int l1,float *cc,float *ch,float *wa1,
799               float *wa2){
800   static float taur = -.5;
801   static float taui = .86602540378443864676372317075293618;
802   int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
803   float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
804   t0=l1*ido;
805
806   t1=0;
807   t2=t0<<1;
808   t3=ido<<1;
809   t4=ido+(ido<<1);
810   t5=0;
811   for (k=0;k<l1;k++){
812     tr2=cc[t3-1]+cc[t3-1];
813     cr2=cc[t5]+(taur*tr2);
814     ch[t1]=cc[t5]+tr2;
815     ci3=taui*(cc[t3]+cc[t3]);
816     ch[t1+t0]=cr2-ci3;
817     ch[t1+t2]=cr2+ci3;
818     t1+=ido;
819     t3+=t4;
820     t5+=t4;
821   }
822
823   if (ido==1) return;
824
825   t1=0;
826   t3=ido<<1;
827   for (k=0;k<l1;k++){
828     t7=t1+(t1<<1);
829     t6=(t5=t7+t3);
830     t8=t1;
831     t10=(t9=t1+t0)+t0;
832
833     for (i=2;i<ido;i+=2){
834       t5+=2;
835       t6-=2;
836       t7+=2;
837       t8+=2;
838       t9+=2;
839       t10+=2;
840       tr2=cc[t5-1]+cc[t6-1];
841       cr2=cc[t7-1]+(taur*tr2);
842       ch[t8-1]=cc[t7-1]+tr2;
843       ti2=cc[t5]-cc[t6];
844       ci2=cc[t7]+(taur*ti2);
845       ch[t8]=cc[t7]+ti2;
846       cr3=taui*(cc[t5-1]-cc[t6-1]);
847       ci3=taui*(cc[t5]+cc[t6]);
848       dr2=cr2-ci3;
849       dr3=cr2+ci3;
850       di2=ci2+cr3;
851       di3=ci2-cr3;
852       ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;
853       ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;
854       ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;
855       ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;
856     }
857     t1+=ido;
858   }
859 }
860
861 STIN void dradb4(int ido,int l1,float *cc,float *ch,float *wa1,
862               float *wa2,float *wa3){
863   static float sqrt2=1.4142135623730950488016887242097;
864   int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
865   float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
866   t0=l1*ido;
867
868   t1=0;
869   t2=ido<<2;
870   t3=0;
871   t6=ido<<1;
872   for (k=0;k<l1;k++){
873     t4=t3+t6;
874     t5=t1;
875     tr3=cc[t4-1]+cc[t4-1];
876     tr4=cc[t4]+cc[t4];
877     tr1=cc[t3]-cc[(t4+=t6)-1];
878     tr2=cc[t3]+cc[t4-1];
879     ch[t5]=tr2+tr3;
880     ch[t5+=t0]=tr1-tr4;
881     ch[t5+=t0]=tr2-tr3;
882     ch[t5+=t0]=tr1+tr4;
883     t1+=ido;
884     t3+=t2;
885   }
886
887   if (ido<2) return;
888   if (ido==2) goto L105;
889
890   t1=0;
891   for (k=0;k<l1;k++){
892     t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
893     t7=t1;
894     for (i=2;i<ido;i+=2){
895       t2+=2;
896       t3+=2;
897       t4-=2;
898       t5-=2;
899       t7+=2;
900       ti1=cc[t2]+cc[t5];
901       ti2=cc[t2]-cc[t5];
902       ti3=cc[t3]-cc[t4];
903       tr4=cc[t3]+cc[t4];
904       tr1=cc[t2-1]-cc[t5-1];
905       tr2=cc[t2-1]+cc[t5-1];
906       ti4=cc[t3-1]-cc[t4-1];
907       tr3=cc[t3-1]+cc[t4-1];
908       ch[t7-1]=tr2+tr3;
909       cr3=tr2-tr3;
910       ch[t7]=ti2+ti3;
911       ci3=ti2-ti3;
912       cr2=tr1-tr4;
913       cr4=tr1+tr4;
914       ci2=ti1+ti4;
915       ci4=ti1-ti4;
916
917       ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
918       ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
919       ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
920       ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
921       ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
922       ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
923     }
924     t1+=ido;
925   }
926
927   if (ido%2 == 1) return;
928
929  L105:
930
931   t1=ido;
932   t2=ido<<2;
933   t3=ido-1;
934   t4=ido+(ido<<1);
935   for (k=0;k<l1;k++){
936     t5=t3;
937     ti1=cc[t1]+cc[t4];
938     ti2=cc[t4]-cc[t1];
939     tr1=cc[t1-1]-cc[t4-1];
940     tr2=cc[t1-1]+cc[t4-1];
941     ch[t5]=tr2+tr2;
942     ch[t5+=t0]=sqrt2*(tr1-ti1);
943     ch[t5+=t0]=ti2+ti2;
944     ch[t5+=t0]=-sqrt2*(tr1+ti1);
945
946     t3+=ido;
947     t1+=t2;
948     t4+=t2;
949   }
950 }
951
952 STIN void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
953         float *c2,float *ch,float *ch2,float *wa){
954   static float tpi=6.28318530717958647692528676655900577;
955   int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
956       t11,t12;
957   float dc2,ai1,ai2,ar1,ar2,ds2;
958   int nbd;
959   float dcp,arg,dsp,ar1h,ar2h;
960   int ipp2;
961
962   t10=ip*ido;
963   t0=l1*ido;
964   arg=tpi/(float)ip;
965   dcp=cos(arg);
966   dsp=sin(arg);
967   nbd=(ido-1)>>1;
968   ipp2=ip;
969   ipph=(ip+1)>>1;
970   if (ido<l1) goto L103;
971
972   t1=0;
973   t2=0;
974   for (k=0;k<l1;k++){
975     t3=t1;
976     t4=t2;
977     for (i=0;i<ido;i++){
978       ch[t3]=cc[t4];
979       t3++;
980       t4++;
981     }
982     t1+=ido;
983     t2+=t10;
984   }
985   goto L106;
986
987  L103:
988   t1=0;
989   for (i=0;i<ido;i++){
990     t2=t1;
991     t3=t1;
992     for (k=0;k<l1;k++){
993       ch[t2]=cc[t3];
994       t2+=ido;
995       t3+=t10;
996     }
997     t1++;
998   }
999
1000  L106:
1001   t1=0;
1002   t2=ipp2*t0;
1003   t7=(t5=ido<<1);
1004   for (j=1;j<ipph;j++){
1005     t1+=t0;
1006     t2-=t0;
1007     t3=t1;
1008     t4=t2;
1009     t6=t5;
1010     for (k=0;k<l1;k++){
1011       ch[t3]=cc[t6-1]+cc[t6-1];
1012       ch[t4]=cc[t6]+cc[t6];
1013       t3+=ido;
1014       t4+=ido;
1015       t6+=t10;
1016     }
1017     t5+=t7;
1018   }
1019
1020   if (ido == 1) goto L116;
1021   if (nbd<l1) goto L112;
1022
1023   t1=0;
1024   t2=ipp2*t0;
1025   t7=0;
1026   for (j=1;j<ipph;j++){
1027     t1+=t0;
1028     t2-=t0;
1029     t3=t1;
1030     t4=t2;
1031
1032     t7+=(ido<<1);
1033     t8=t7;
1034     for (k=0;k<l1;k++){
1035       t5=t3;
1036       t6=t4;
1037       t9=t8;
1038       t11=t8;
1039       for (i=2;i<ido;i+=2){
1040     t5+=2;
1041     t6+=2;
1042     t9+=2;
1043     t11-=2;
1044     ch[t5-1]=cc[t9-1]+cc[t11-1];
1045     ch[t6-1]=cc[t9-1]-cc[t11-1];
1046     ch[t5]=cc[t9]-cc[t11];
1047     ch[t6]=cc[t9]+cc[t11];
1048       }
1049       t3+=ido;
1050       t4+=ido;
1051       t8+=t10;
1052     }
1053   }
1054   goto L116;
1055
1056  L112:
1057   t1=0;
1058   t2=ipp2*t0;
1059   t7=0;
1060   for (j=1;j<ipph;j++){
1061     t1+=t0;
1062     t2-=t0;
1063     t3=t1;
1064     t4=t2;
1065     t7+=(ido<<1);
1066     t8=t7;
1067     t9=t7;
1068     for (i=2;i<ido;i+=2){
1069       t3+=2;
1070       t4+=2;
1071       t8+=2;
1072       t9-=2;
1073       t5=t3;
1074       t6=t4;
1075       t11=t8;
1076       t12=t9;
1077       for (k=0;k<l1;k++){
1078     ch[t5-1]=cc[t11-1]+cc[t12-1];
1079     ch[t6-1]=cc[t11-1]-cc[t12-1];
1080     ch[t5]=cc[t11]-cc[t12];
1081     ch[t6]=cc[t11]+cc[t12];
1082     t5+=ido;
1083     t6+=ido;
1084     t11+=t10;
1085     t12+=t10;
1086       }
1087     }
1088   }
1089
1090 L116:
1091   ar1=1.;
1092   ai1=0.;
1093   t1=0;
1094   t9=(t2=ipp2*idl1);
1095   t3=(ip-1)*idl1;
1096   for (l=1;l<ipph;l++){
1097     t1+=idl1;
1098     t2-=idl1;
1099
1100     ar1h=dcp*ar1-dsp*ai1;
1101     ai1=dcp*ai1+dsp*ar1;
1102     ar1=ar1h;
1103     t4=t1;
1104     t5=t2;
1105     t6=0;
1106     t7=idl1;
1107     t8=t3;
1108     for (ik=0;ik<idl1;ik++){
1109       c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
1110       c2[t5++]=ai1*ch2[t8++];
1111     }
1112     dc2=ar1;
1113     ds2=ai1;
1114     ar2=ar1;
1115     ai2=ai1;
1116
1117     t6=idl1;
1118     t7=t9-idl1;
1119     for (j=2;j<ipph;j++){
1120       t6+=idl1;
1121       t7-=idl1;
1122       ar2h=dc2*ar2-ds2*ai2;
1123       ai2=dc2*ai2+ds2*ar2;
1124       ar2=ar2h;
1125       t4=t1;
1126       t5=t2;
1127       t11=t6;
1128       t12=t7;
1129       for (ik=0;ik<idl1;ik++){
1130     c2[t4++]+=ar2*ch2[t11++];
1131     c2[t5++]+=ai2*ch2[t12++];
1132       }
1133     }
1134   }
1135
1136   t1=0;
1137   for (j=1;j<ipph;j++){
1138     t1+=idl1;
1139     t2=t1;
1140     for (ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
1141   }
1142
1143   t1=0;
1144   t2=ipp2*t0;
1145   for (j=1;j<ipph;j++){
1146     t1+=t0;
1147     t2-=t0;
1148     t3=t1;
1149     t4=t2;
1150     for (k=0;k<l1;k++){
1151       ch[t3]=c1[t3]-c1[t4];
1152       ch[t4]=c1[t3]+c1[t4];
1153       t3+=ido;
1154       t4+=ido;
1155     }
1156   }
1157
1158   if (ido==1) goto L132;
1159   if (nbd<l1) goto L128;
1160
1161   t1=0;
1162   t2=ipp2*t0;
1163   for (j=1;j<ipph;j++){
1164     t1+=t0;
1165     t2-=t0;
1166     t3=t1;
1167     t4=t2;
1168     for (k=0;k<l1;k++){
1169       t5=t3;
1170       t6=t4;
1171       for (i=2;i<ido;i+=2){
1172     t5+=2;
1173     t6+=2;
1174     ch[t5-1]=c1[t5-1]-c1[t6];
1175     ch[t6-1]=c1[t5-1]+c1[t6];
1176     ch[t5]=c1[t5]+c1[t6-1];
1177     ch[t6]=c1[t5]-c1[t6-1];
1178       }
1179       t3+=ido;
1180       t4+=ido;
1181     }
1182   }
1183   goto L132;
1184
1185  L128:
1186   t1=0;
1187   t2=ipp2*t0;
1188   for (j=1;j<ipph;j++){
1189     t1+=t0;
1190     t2-=t0;
1191     t3=t1;
1192     t4=t2;
1193     for (i=2;i<ido;i+=2){
1194       t3+=2;
1195       t4+=2;
1196       t5=t3;
1197       t6=t4;
1198       for (k=0;k<l1;k++){
1199     ch[t5-1]=c1[t5-1]-c1[t6];
1200     ch[t6-1]=c1[t5-1]+c1[t6];
1201     ch[t5]=c1[t5]+c1[t6-1];
1202     ch[t6]=c1[t5]-c1[t6-1];
1203     t5+=ido;
1204     t6+=ido;
1205       }
1206     }
1207   }
1208
1209 L132:
1210   if (ido==1) return;
1211
1212   for (ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
1213
1214   t1=0;
1215   for (j=1;j<ip;j++){
1216     t2=(t1+=t0);
1217     for (k=0;k<l1;k++){
1218       c1[t2]=ch[t2];
1219       t2+=ido;
1220     }
1221   }
1222
1223   if (nbd>l1) goto L139;
1224
1225   is= -ido-1;
1226   t1=0;
1227   for (j=1;j<ip;j++){
1228     is+=ido;
1229     t1+=t0;
1230     idij=is;
1231     t2=t1;
1232     for (i=2;i<ido;i+=2){
1233       t2+=2;
1234       idij+=2;
1235       t3=t2;
1236       for (k=0;k<l1;k++){
1237     c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
1238     c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
1239     t3+=ido;
1240       }
1241     }
1242   }
1243   return;
1244
1245  L139:
1246   is= -ido-1;
1247   t1=0;
1248   for (j=1;j<ip;j++){
1249     is+=ido;
1250     t1+=t0;
1251     t2=t1;
1252     for (k=0;k<l1;k++){
1253       idij=is;
1254       t3=t2;
1255       for (i=2;i<ido;i+=2){
1256     idij+=2;
1257     t3+=2;
1258     c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
1259     c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
1260       }
1261       t2+=ido;
1262     }
1263   }
1264 }
1265
1266 STIN void drftb1(int n, float *c, float *ch, float *wa, int *ifac){
1267   int i,k1,l1,l2;
1268   int na;
1269   int nf,ip,iw,ix2,ix3,ido,idl1;
1270
1271   nf=ifac[1];
1272   na=0;
1273   l1=1;
1274   iw=1;
1275
1276   for (k1=0;k1<nf;k1++){
1277     ip=ifac[k1 + 2];
1278     l2=ip*l1;
1279     ido=n/l2;
1280     idl1=ido*l1;
1281     if (ip!=4) goto L103;
1282     ix2=iw+ido;
1283     ix3=ix2+ido;
1284
1285     if (na!=0)
1286       dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
1287     else
1288       dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
1289     na=1-na;
1290     goto L115;
1291
1292   L103:
1293     if (ip!=2) goto L106;
1294
1295     if (na!=0)
1296       dradb2(ido,l1,ch,c,wa+iw-1);
1297     else
1298       dradb2(ido,l1,c,ch,wa+iw-1);
1299     na=1-na;
1300     goto L115;
1301
1302   L106:
1303     if (ip!=3) goto L109;
1304
1305     ix2=iw+ido;
1306     if (na!=0)
1307       dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);
1308     else
1309       dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);
1310     na=1-na;
1311     goto L115;
1312
1313   L109:
1314 /*    The radix five case can be translated later..... */
1315 /*    if (ip!=5) goto L112;
1316
1317     ix2=iw+ido;
1318     ix3=ix2+ido;
1319     ix4=ix3+ido;
1320     if (na!=0)
1321       dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
1322     else
1323       dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
1324     na=1-na;
1325     goto L115;
1326
1327   L112:*/
1328     if (na!=0)
1329       dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
1330     else
1331       dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
1332     if (ido==1)na=1-na;
1333
1334   L115:
1335     l1=l2;
1336     iw+=(ip-1)*ido;
1337   }
1338
1339   if (na==0) return;
1340
1341   for (i=0;i<n;i++)c[i]=ch[i];
1342 }
1343
1344 void __ogg_fdrfftb(int n, float *r, float *wsave, int *ifac){
1345   if (n == 1) return;
1346   drftb1(n, r, wsave, wsave+n, ifac);
1347 }
1348
1349 STIN void dcsqb1(int n,float *x,float *w,float *xh,int *ifac){
1350   int modn,i,k,kc;
1351   int np2,ns2;
1352   float xim1;
1353
1354   ns2=(n+1)>>1;
1355   np2=n;
1356
1357   for (i=2;i<n;i+=2){
1358     xim1=x[i-1]+x[i];
1359     x[i]-=x[i-1];
1360     x[i-1]=xim1;
1361   }
1362
1363   x[0]+=x[0];
1364   modn=n%2;
1365   if (modn==0)x[n-1]+=x[n-1];
1366
1367   __ogg_fdrfftb(n,x,xh,ifac);
1368
1369   kc=np2;
1370   for (k=1;k<ns2;k++){
1371     kc--;
1372     xh[k]=w[k-1]*x[kc]+w[kc-1]*x[k];
1373     xh[kc]=w[k-1]*x[k]-w[kc-1]*x[kc];
1374   }
1375
1376   if (modn==0)x[ns2]=w[ns2-1]*(x[ns2]+x[ns2]);
1377
1378   kc=np2;
1379   for (k=1;k<ns2;k++){
1380     kc--;
1381     x[k]=xh[k]+xh[kc];
1382     x[kc]=xh[k]-xh[kc];
1383   }
1384   x[0]+=x[0];
1385 }
1386
1387 void __ogg_fdcosqb(int n,float *x,float *wsave,int *ifac){
1388   static float tsqrt2 = 2.8284271247461900976033774484194;
1389   float x1;
1390
1391   if (n<2){
1392     x[0]*=4;
1393     return;
1394   }
1395   if (n==2){
1396     x1=(x[0]+x[1])*4;
1397     x[1]=tsqrt2*(x[0]-x[1]);
1398     x[0]=x1;
1399     return;
1400   }
1401
1402   dcsqb1(n,x,wsave,wsave+n,ifac);
1403 }
1404
1405
1406