Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/ftbase.cs @ 9280

Last change on this file since 9280 was 2645, checked in by mkommend, 15 years ago

extracted external libraries and adapted dependent plugins (ticket #837)
File size: 60.2 KB
Line 
1/*************************************************************************
2Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class ftbase
26    {
27        public struct ftplan
28        {
29            public int[] plan;
30            public double[] precomputed;
31            public double[] tmpbuf;
32            public double[] stackbuf;
33        };
34
35
36
37
38        public const int ftbaseplanentrysize = 8;
39        public const int ftbasecffttask = 0;
40        public const int ftbaserfhttask = 1;
41        public const int ftbaserffttask = 2;
42        public const int fftcooleytukeyplan = 0;
43        public const int fftbluesteinplan = 1;
44        public const int fftcodeletplan = 2;
45        public const int fhtcooleytukeyplan = 3;
46        public const int fhtcodeletplan = 4;
47        public const int fftrealcooleytukeyplan = 5;
48        public const int fftemptyplan = 6;
49        public const int fhtn2plan = 999;
50        public const int ftbaseupdatetw = 4;
51        public const int ftbasecodeletmax = 5;
52        public const int ftbasecodeletrecommended = 5;
53        public const double ftbaseinefficiencyfactor = 1.3;
54        public const int ftbasemaxsmoothfactor = 5;
55
56
57        /*************************************************************************
58        This subroutine generates FFT plan - a decomposition of a N-length FFT to
59        the more simpler operations. Plan consists of the root entry and the child
60        entries.
61
62        Subroutine parameters:
63            N               task size
64           
65        Output parameters:
66            Plan            plan
67
68          -- ALGLIB --
69             Copyright 01.05.2009 by Bochkanov Sergey
70        *************************************************************************/
71        public static void ftbasegeneratecomplexfftplan(int n,
72            ref ftplan plan)
73        {
74            int planarraysize = 0;
75            int plansize = 0;
76            int precomputedsize = 0;
77            int tmpmemsize = 0;
78            int stackmemsize = 0;
79            int stackptr = 0;
80
81            planarraysize = 1;
82            plansize = 0;
83            precomputedsize = 0;
84            stackmemsize = 0;
85            stackptr = 0;
86            tmpmemsize = 2*n;
87            plan.plan = new int[planarraysize];
88            ftbasegenerateplanrec(n, ftbasecffttask, ref plan, ref plansize, ref precomputedsize, ref planarraysize, ref tmpmemsize, ref stackmemsize, stackptr);
89            System.Diagnostics.Debug.Assert(stackptr==0, "Internal error in FTBaseGenerateComplexFFTPlan: stack ptr!");
90            plan.stackbuf = new double[Math.Max(stackmemsize, 1)];
91            plan.tmpbuf = new double[Math.Max(tmpmemsize, 1)];
92            plan.precomputed = new double[Math.Max(precomputedsize, 1)];
93            stackptr = 0;
94            ftbaseprecomputeplanrec(ref plan, 0, stackptr);
95            System.Diagnostics.Debug.Assert(stackptr==0, "Internal error in FTBaseGenerateComplexFFTPlan: stack ptr!");
96        }
97
98
99        /*************************************************************************
100        Generates real FFT plan
101        *************************************************************************/
102        public static void ftbasegeneraterealfftplan(int n,
103            ref ftplan plan)
104        {
105            int planarraysize = 0;
106            int plansize = 0;
107            int precomputedsize = 0;
108            int tmpmemsize = 0;
109            int stackmemsize = 0;
110            int stackptr = 0;
111
112            planarraysize = 1;
113            plansize = 0;
114            precomputedsize = 0;
115            stackmemsize = 0;
116            stackptr = 0;
117            tmpmemsize = 2*n;
118            plan.plan = new int[planarraysize];
119            ftbasegenerateplanrec(n, ftbaserffttask, ref plan, ref plansize, ref precomputedsize, ref planarraysize, ref tmpmemsize, ref stackmemsize, stackptr);
120            System.Diagnostics.Debug.Assert(stackptr==0, "Internal error in FTBaseGenerateRealFFTPlan: stack ptr!");
121            plan.stackbuf = new double[Math.Max(stackmemsize, 1)];
122            plan.tmpbuf = new double[Math.Max(tmpmemsize, 1)];
123            plan.precomputed = new double[Math.Max(precomputedsize, 1)];
124            stackptr = 0;
125            ftbaseprecomputeplanrec(ref plan, 0, stackptr);
126            System.Diagnostics.Debug.Assert(stackptr==0, "Internal error in FTBaseGenerateRealFFTPlan: stack ptr!");
127        }
128
129
130        /*************************************************************************
131        Generates real FHT plan
132        *************************************************************************/
133        public static void ftbasegeneraterealfhtplan(int n,
134            ref ftplan plan)
135        {
136            int planarraysize = 0;
137            int plansize = 0;
138            int precomputedsize = 0;
139            int tmpmemsize = 0;
140            int stackmemsize = 0;
141            int stackptr = 0;
142
143            planarraysize = 1;
144            plansize = 0;
145            precomputedsize = 0;
146            stackmemsize = 0;
147            stackptr = 0;
148            tmpmemsize = n;
149            plan.plan = new int[planarraysize];
150            ftbasegenerateplanrec(n, ftbaserfhttask, ref plan, ref plansize, ref precomputedsize, ref planarraysize, ref tmpmemsize, ref stackmemsize, stackptr);
151            System.Diagnostics.Debug.Assert(stackptr==0, "Internal error in FTBaseGenerateRealFHTPlan: stack ptr!");
152            plan.stackbuf = new double[Math.Max(stackmemsize, 1)];
153            plan.tmpbuf = new double[Math.Max(tmpmemsize, 1)];
154            plan.precomputed = new double[Math.Max(precomputedsize, 1)];
155            stackptr = 0;
156            ftbaseprecomputeplanrec(ref plan, 0, stackptr);
157            System.Diagnostics.Debug.Assert(stackptr==0, "Internal error in FTBaseGenerateRealFHTPlan: stack ptr!");
158        }
159
160
161        /*************************************************************************
162        This subroutine executes FFT/FHT plan.
163
164        If Plan is a:
165        * complex FFT plan  -   sizeof(A)=2*N,
166                                A contains interleaved real/imaginary values
167        * real FFT plan     -   sizeof(A)=2*N,
168                                A contains real values interleaved with zeros
169        * real FHT plan     -   sizeof(A)=2*N,
170                                A contains real values interleaved with zeros
171
172          -- ALGLIB --
173             Copyright 01.05.2009 by Bochkanov Sergey
174        *************************************************************************/
175        public static void ftbaseexecuteplan(ref double[] a,
176            int aoffset,
177            int n,
178            ref ftplan plan)
179        {
180            int stackptr = 0;
181
182            stackptr = 0;
183            ftbaseexecuteplanrec(ref a, aoffset, ref plan, 0, stackptr);
184        }
185
186
187        /*************************************************************************
188        Recurrent subroutine for the FTBaseExecutePlan
189
190        Parameters:
191            A           FFT'ed array
192            AOffset     offset of the FFT'ed part (distance is measured in doubles)
193
194          -- ALGLIB --
195             Copyright 01.05.2009 by Bochkanov Sergey
196        *************************************************************************/
197        public static void ftbaseexecuteplanrec(ref double[] a,
198            int aoffset,
199            ref ftplan plan,
200            int entryoffset,
201            int stackptr)
202        {
203            int i = 0;
204            int j = 0;
205            int k = 0;
206            int idx = 0;
207            int n1 = 0;
208            int n2 = 0;
209            int n = 0;
210            int m = 0;
211            int offs = 0;
212            int offs1 = 0;
213            int offs2 = 0;
214            int offsa = 0;
215            int offsb = 0;
216            int offsp = 0;
217            double hk = 0;
218            double hnk = 0;
219            double x = 0;
220            double y = 0;
221            double bx = 0;
222            double by = 0;
223            double v = 0;
224            double[] emptyarray = new double[0];
225            double a0x = 0;
226            double a0y = 0;
227            double a1x = 0;
228            double a1y = 0;
229            double a2x = 0;
230            double a2y = 0;
231            double a3x = 0;
232            double a3y = 0;
233            double a4x = 0;
234            double a4y = 0;
235            double v0 = 0;
236            double v1 = 0;
237            double v2 = 0;
238            double v3 = 0;
239            double t1x = 0;
240            double t1y = 0;
241            double t2x = 0;
242            double t2y = 0;
243            double t3x = 0;
244            double t3y = 0;
245            double t4x = 0;
246            double t4y = 0;
247            double t5x = 0;
248            double t5y = 0;
249            double m1x = 0;
250            double m1y = 0;
251            double m2x = 0;
252            double m2y = 0;
253            double m3x = 0;
254            double m3y = 0;
255            double m4x = 0;
256            double m4y = 0;
257            double m5x = 0;
258            double m5y = 0;
259            double s1x = 0;
260            double s1y = 0;
261            double s2x = 0;
262            double s2y = 0;
263            double s3x = 0;
264            double s3y = 0;
265            double s4x = 0;
266            double s4y = 0;
267            double s5x = 0;
268            double s5y = 0;
269            double c1 = 0;
270            double c2 = 0;
271            double c3 = 0;
272            double c4 = 0;
273            double c5 = 0;
274            double[] tmp = new double[0];
275            int i_ = 0;
276            int i1_ = 0;
277
278            if( plan.plan[entryoffset+3]==fftemptyplan )
279            {
280                return;
281            }
282            if( plan.plan[entryoffset+3]==fftcooleytukeyplan )
283            {
284               
285                //
286                // Cooley-Tukey plan
287                // * transposition
288                // * row-wise FFT
289                // * twiddle factors:
290                //   - TwBase is a basis twiddle factor for I=1, J=1
291                //   - TwRow is a twiddle factor for a second element in a row (J=1)
292                //   - Tw is a twiddle factor for a current element
293                // * transposition again
294                // * row-wise FFT again
295                //
296                n1 = plan.plan[entryoffset+1];
297                n2 = plan.plan[entryoffset+2];
298                internalcomplexlintranspose(ref a, n1, n2, aoffset, ref plan.tmpbuf);
299                for(i=0; i<=n2-1; i++)
300                {
301                    ftbaseexecuteplanrec(ref a, aoffset+i*n1*2, ref plan, plan.plan[entryoffset+5], stackptr);
302                }
303                ffttwcalc(ref a, aoffset, n1, n2);
304                internalcomplexlintranspose(ref a, n2, n1, aoffset, ref plan.tmpbuf);
305                for(i=0; i<=n1-1; i++)
306                {
307                    ftbaseexecuteplanrec(ref a, aoffset+i*n2*2, ref plan, plan.plan[entryoffset+6], stackptr);
308                }
309                internalcomplexlintranspose(ref a, n1, n2, aoffset, ref plan.tmpbuf);
310                return;
311            }
312            if( plan.plan[entryoffset+3]==fftrealcooleytukeyplan )
313            {
314               
315                //
316                // Cooley-Tukey plan
317                // * transposition
318                // * row-wise FFT
319                // * twiddle factors:
320                //   - TwBase is a basis twiddle factor for I=1, J=1
321                //   - TwRow is a twiddle factor for a second element in a row (J=1)
322                //   - Tw is a twiddle factor for a current element
323                // * transposition again
324                // * row-wise FFT again
325                //
326                n1 = plan.plan[entryoffset+1];
327                n2 = plan.plan[entryoffset+2];
328                internalcomplexlintranspose(ref a, n2, n1, aoffset, ref plan.tmpbuf);
329                for(i=0; i<=n1/2-1; i++)
330                {
331                   
332                    //
333                    // pack two adjacent smaller real FFT's together,
334                    // make one complex FFT,
335                    // unpack result
336                    //
337                    offs = aoffset+2*i*n2*2;
338                    for(k=0; k<=n2-1; k++)
339                    {
340                        a[offs+2*k+1] = a[offs+2*n2+2*k+0];
341                    }
342                    ftbaseexecuteplanrec(ref a, offs, ref plan, plan.plan[entryoffset+6], stackptr);
343                    plan.tmpbuf[0] = a[offs+0];
344                    plan.tmpbuf[1] = 0;
345                    plan.tmpbuf[2*n2+0] = a[offs+1];
346                    plan.tmpbuf[2*n2+1] = 0;
347                    for(k=1; k<=n2-1; k++)
348                    {
349                        offs1 = 2*k;
350                        offs2 = 2*n2+2*k;
351                        hk = a[offs+2*k+0];
352                        hnk = a[offs+2*(n2-k)+0];
353                        plan.tmpbuf[offs1+0] = +(0.5*(hk+hnk));
354                        plan.tmpbuf[offs2+1] = -(0.5*(hk-hnk));
355                        hk = a[offs+2*k+1];
356                        hnk = a[offs+2*(n2-k)+1];
357                        plan.tmpbuf[offs2+0] = +(0.5*(hk+hnk));
358                        plan.tmpbuf[offs1+1] = +(0.5*(hk-hnk));
359                    }
360                    i1_ = (0) - (offs);
361                    for(i_=offs; i_<=offs+2*n2*2-1;i_++)
362                    {
363                        a[i_] = plan.tmpbuf[i_+i1_];
364                    }
365                }
366                if( n1%2!=0 )
367                {
368                    ftbaseexecuteplanrec(ref a, aoffset+(n1-1)*n2*2, ref plan, plan.plan[entryoffset+6], stackptr);
369                }
370                ffttwcalc(ref a, aoffset, n2, n1);
371                internalcomplexlintranspose(ref a, n1, n2, aoffset, ref plan.tmpbuf);
372                for(i=0; i<=n2-1; i++)
373                {
374                    ftbaseexecuteplanrec(ref a, aoffset+i*n1*2, ref plan, plan.plan[entryoffset+5], stackptr);
375                }
376                internalcomplexlintranspose(ref a, n2, n1, aoffset, ref plan.tmpbuf);
377                return;
378            }
379            if( plan.plan[entryoffset+3]==fhtcooleytukeyplan )
380            {
381               
382                //
383                // Cooley-Tukey FHT plan:
384                // * transpose                    \
385                // * smaller FHT's                |
386                // * pre-process                  |
387                // * multiply by twiddle factors  | corresponds to multiplication by H1
388                // * post-process                 |
389                // * transpose again              /
390                // * multiply by H2 (smaller FHT's)
391                // * final transposition
392                //
393                // For more details see Vitezslav Vesely, "Fast algorithms
394                // of Fourier and Hartley transform and their implementation in MATLAB",
395                // page 31.
396                //
397                n1 = plan.plan[entryoffset+1];
398                n2 = plan.plan[entryoffset+2];
399                n = n1*n2;
400                internalreallintranspose(ref a, n1, n2, aoffset, ref plan.tmpbuf);
401                for(i=0; i<=n2-1; i++)
402                {
403                    ftbaseexecuteplanrec(ref a, aoffset+i*n1, ref plan, plan.plan[entryoffset+5], stackptr);
404                }
405                for(i=0; i<=n2-1; i++)
406                {
407                    for(j=0; j<=n1-1; j++)
408                    {
409                        offsa = aoffset+i*n1;
410                        hk = a[offsa+j];
411                        hnk = a[offsa+(n1-j)%n1];
412                        offs = 2*(i*n1+j);
413                        plan.tmpbuf[offs+0] = -(0.5*(hnk-hk));
414                        plan.tmpbuf[offs+1] = +(0.5*(hk+hnk));
415                    }
416                }
417                ffttwcalc(ref plan.tmpbuf, 0, n1, n2);
418                for(j=0; j<=n1-1; j++)
419                {
420                    a[aoffset+j] = plan.tmpbuf[2*j+0]+plan.tmpbuf[2*j+1];
421                }
422                if( n2%2==0 )
423                {
424                    offs = 2*(n2/2)*n1;
425                    offsa = aoffset+n2/2*n1;
426                    for(j=0; j<=n1-1; j++)
427                    {
428                        a[offsa+j] = plan.tmpbuf[offs+2*j+0]+plan.tmpbuf[offs+2*j+1];
429                    }
430                }
431                for(i=1; i<=(n2+1)/2-1; i++)
432                {
433                    offs = 2*i*n1;
434                    offs2 = 2*(n2-i)*n1;
435                    offsa = aoffset+i*n1;
436                    for(j=0; j<=n1-1; j++)
437                    {
438                        a[offsa+j] = plan.tmpbuf[offs+2*j+1]+plan.tmpbuf[offs2+2*j+0];
439                    }
440                    offsa = aoffset+(n2-i)*n1;
441                    for(j=0; j<=n1-1; j++)
442                    {
443                        a[offsa+j] = plan.tmpbuf[offs+2*j+0]+plan.tmpbuf[offs2+2*j+1];
444                    }
445                }
446                internalreallintranspose(ref a, n2, n1, aoffset, ref plan.tmpbuf);
447                for(i=0; i<=n1-1; i++)
448                {
449                    ftbaseexecuteplanrec(ref a, aoffset+i*n2, ref plan, plan.plan[entryoffset+6], stackptr);
450                }
451                internalreallintranspose(ref a, n1, n2, aoffset, ref plan.tmpbuf);
452                return;
453            }
454            if( plan.plan[entryoffset+3]==fhtn2plan )
455            {
456               
457                //
458                // Cooley-Tukey FHT plan
459                //
460                n1 = plan.plan[entryoffset+1];
461                n2 = plan.plan[entryoffset+2];
462                n = n1*n2;
463                reffht(ref a, n, aoffset);
464                return;
465            }
466            if( plan.plan[entryoffset+3]==fftcodeletplan )
467            {
468                n1 = plan.plan[entryoffset+1];
469                n2 = plan.plan[entryoffset+2];
470                n = n1*n2;
471                if( n==2 )
472                {
473                    a0x = a[aoffset+0];
474                    a0y = a[aoffset+1];
475                    a1x = a[aoffset+2];
476                    a1y = a[aoffset+3];
477                    v0 = a0x+a1x;
478                    v1 = a0y+a1y;
479                    v2 = a0x-a1x;
480                    v3 = a0y-a1y;
481                    a[aoffset+0] = v0;
482                    a[aoffset+1] = v1;
483                    a[aoffset+2] = v2;
484                    a[aoffset+3] = v3;
485                    return;
486                }
487                if( n==3 )
488                {
489                    offs = plan.plan[entryoffset+7];
490                    c1 = plan.precomputed[offs+0];
491                    c2 = plan.precomputed[offs+1];
492                    a0x = a[aoffset+0];
493                    a0y = a[aoffset+1];
494                    a1x = a[aoffset+2];
495                    a1y = a[aoffset+3];
496                    a2x = a[aoffset+4];
497                    a2y = a[aoffset+5];
498                    t1x = a1x+a2x;
499                    t1y = a1y+a2y;
500                    a0x = a0x+t1x;
501                    a0y = a0y+t1y;
502                    m1x = c1*t1x;
503                    m1y = c1*t1y;
504                    m2x = c2*(a1y-a2y);
505                    m2y = c2*(a2x-a1x);
506                    s1x = a0x+m1x;
507                    s1y = a0y+m1y;
508                    a1x = s1x+m2x;
509                    a1y = s1y+m2y;
510                    a2x = s1x-m2x;
511                    a2y = s1y-m2y;
512                    a[aoffset+0] = a0x;
513                    a[aoffset+1] = a0y;
514                    a[aoffset+2] = a1x;
515                    a[aoffset+3] = a1y;
516                    a[aoffset+4] = a2x;
517                    a[aoffset+5] = a2y;
518                    return;
519                }
520                if( n==4 )
521                {
522                    a0x = a[aoffset+0];
523                    a0y = a[aoffset+1];
524                    a1x = a[aoffset+2];
525                    a1y = a[aoffset+3];
526                    a2x = a[aoffset+4];
527                    a2y = a[aoffset+5];
528                    a3x = a[aoffset+6];
529                    a3y = a[aoffset+7];
530                    t1x = a0x+a2x;
531                    t1y = a0y+a2y;
532                    t2x = a1x+a3x;
533                    t2y = a1y+a3y;
534                    m2x = a0x-a2x;
535                    m2y = a0y-a2y;
536                    m3x = a1y-a3y;
537                    m3y = a3x-a1x;
538                    a[aoffset+0] = t1x+t2x;
539                    a[aoffset+1] = t1y+t2y;
540                    a[aoffset+4] = t1x-t2x;
541                    a[aoffset+5] = t1y-t2y;
542                    a[aoffset+2] = m2x+m3x;
543                    a[aoffset+3] = m2y+m3y;
544                    a[aoffset+6] = m2x-m3x;
545                    a[aoffset+7] = m2y-m3y;
546                    return;
547                }
548                if( n==5 )
549                {
550                    offs = plan.plan[entryoffset+7];
551                    c1 = plan.precomputed[offs+0];
552                    c2 = plan.precomputed[offs+1];
553                    c3 = plan.precomputed[offs+2];
554                    c4 = plan.precomputed[offs+3];
555                    c5 = plan.precomputed[offs+4];
556                    t1x = a[aoffset+2]+a[aoffset+8];
557                    t1y = a[aoffset+3]+a[aoffset+9];
558                    t2x = a[aoffset+4]+a[aoffset+6];
559                    t2y = a[aoffset+5]+a[aoffset+7];
560                    t3x = a[aoffset+2]-a[aoffset+8];
561                    t3y = a[aoffset+3]-a[aoffset+9];
562                    t4x = a[aoffset+6]-a[aoffset+4];
563                    t4y = a[aoffset+7]-a[aoffset+5];
564                    t5x = t1x+t2x;
565                    t5y = t1y+t2y;
566                    a[aoffset+0] = a[aoffset+0]+t5x;
567                    a[aoffset+1] = a[aoffset+1]+t5y;
568                    m1x = c1*t5x;
569                    m1y = c1*t5y;
570                    m2x = c2*(t1x-t2x);
571                    m2y = c2*(t1y-t2y);
572                    m3x = -(c3*(t3y+t4y));
573                    m3y = c3*(t3x+t4x);
574                    m4x = -(c4*t4y);
575                    m4y = c4*t4x;
576                    m5x = -(c5*t3y);
577                    m5y = c5*t3x;
578                    s3x = m3x-m4x;
579                    s3y = m3y-m4y;
580                    s5x = m3x+m5x;
581                    s5y = m3y+m5y;
582                    s1x = a[aoffset+0]+m1x;
583                    s1y = a[aoffset+1]+m1y;
584                    s2x = s1x+m2x;
585                    s2y = s1y+m2y;
586                    s4x = s1x-m2x;
587                    s4y = s1y-m2y;
588                    a[aoffset+2] = s2x+s3x;
589                    a[aoffset+3] = s2y+s3y;
590                    a[aoffset+4] = s4x+s5x;
591                    a[aoffset+5] = s4y+s5y;
592                    a[aoffset+6] = s4x-s5x;
593                    a[aoffset+7] = s4y-s5y;
594                    a[aoffset+8] = s2x-s3x;
595                    a[aoffset+9] = s2y-s3y;
596                    return;
597                }
598            }
599            if( plan.plan[entryoffset+3]==fhtcodeletplan )
600            {
601                n1 = plan.plan[entryoffset+1];
602                n2 = plan.plan[entryoffset+2];
603                n = n1*n2;
604                if( n==2 )
605                {
606                    a0x = a[aoffset+0];
607                    a1x = a[aoffset+1];
608                    a[aoffset+0] = a0x+a1x;
609                    a[aoffset+1] = a0x-a1x;
610                    return;
611                }
612                if( n==3 )
613                {
614                    offs = plan.plan[entryoffset+7];
615                    c1 = plan.precomputed[offs+0];
616                    c2 = plan.precomputed[offs+1];
617                    a0x = a[aoffset+0];
618                    a1x = a[aoffset+1];
619                    a2x = a[aoffset+2];
620                    t1x = a1x+a2x;
621                    a0x = a0x+t1x;
622                    m1x = c1*t1x;
623                    m2y = c2*(a2x-a1x);
624                    s1x = a0x+m1x;
625                    a[aoffset+0] = a0x;
626                    a[aoffset+1] = s1x-m2y;
627                    a[aoffset+2] = s1x+m2y;
628                    return;
629                }
630                if( n==4 )
631                {
632                    a0x = a[aoffset+0];
633                    a1x = a[aoffset+1];
634                    a2x = a[aoffset+2];
635                    a3x = a[aoffset+3];
636                    t1x = a0x+a2x;
637                    t2x = a1x+a3x;
638                    m2x = a0x-a2x;
639                    m3y = a3x-a1x;
640                    a[aoffset+0] = t1x+t2x;
641                    a[aoffset+1] = m2x-m3y;
642                    a[aoffset+2] = t1x-t2x;
643                    a[aoffset+3] = m2x+m3y;
644                    return;
645                }
646                if( n==5 )
647                {
648                    offs = plan.plan[entryoffset+7];
649                    c1 = plan.precomputed[offs+0];
650                    c2 = plan.precomputed[offs+1];
651                    c3 = plan.precomputed[offs+2];
652                    c4 = plan.precomputed[offs+3];
653                    c5 = plan.precomputed[offs+4];
654                    t1x = a[aoffset+1]+a[aoffset+4];
655                    t2x = a[aoffset+2]+a[aoffset+3];
656                    t3x = a[aoffset+1]-a[aoffset+4];
657                    t4x = a[aoffset+3]-a[aoffset+2];
658                    t5x = t1x+t2x;
659                    v0 = a[aoffset+0]+t5x;
660                    a[aoffset+0] = v0;
661                    m2x = c2*(t1x-t2x);
662                    m3y = c3*(t3x+t4x);
663                    s3y = m3y-c4*t4x;
664                    s5y = m3y+c5*t3x;
665                    s1x = v0+c1*t5x;
666                    s2x = s1x+m2x;
667                    s4x = s1x-m2x;
668                    a[aoffset+1] = s2x-s3y;
669                    a[aoffset+2] = s4x-s5y;
670                    a[aoffset+3] = s4x+s5y;
671                    a[aoffset+4] = s2x+s3y;
672                    return;
673                }
674            }
675            if( plan.plan[entryoffset+3]==fftbluesteinplan )
676            {
677               
678                //
679                // Bluestein plan:
680                // 1. multiply by precomputed coefficients
681                // 2. make convolution: forward FFT, multiplication by precomputed FFT
682                //    and backward FFT. backward FFT is represented as
683                //
684                //        invfft(x) = fft(x')'/M
685                //
686                //    for performance reasons reduction of inverse FFT to
687                //    forward FFT is merged with multiplication of FFT components
688                //    and last stage of Bluestein's transformation.
689                // 3. post-multiplication by Bluestein factors
690                //
691                n = plan.plan[entryoffset+1];
692                m = plan.plan[entryoffset+4];
693                offs = plan.plan[entryoffset+7];
694                for(i=stackptr+2*n; i<=stackptr+2*m-1; i++)
695                {
696                    plan.stackbuf[i] = 0;
697                }
698                offsp = offs+2*m;
699                offsa = aoffset;
700                offsb = stackptr;
701                for(i=0; i<=n-1; i++)
702                {
703                    bx = plan.precomputed[offsp+0];
704                    by = plan.precomputed[offsp+1];
705                    x = a[offsa+0];
706                    y = a[offsa+1];
707                    plan.stackbuf[offsb+0] = x*bx-y*-by;
708                    plan.stackbuf[offsb+1] = x*-by+y*bx;
709                    offsp = offsp+2;
710                    offsa = offsa+2;
711                    offsb = offsb+2;
712                }
713                ftbaseexecuteplanrec(ref plan.stackbuf, stackptr, ref plan, plan.plan[entryoffset+5], stackptr+2*2*m);
714                offsb = stackptr;
715                offsp = offs;
716                for(i=0; i<=m-1; i++)
717                {
718                    x = plan.stackbuf[offsb+0];
719                    y = plan.stackbuf[offsb+1];
720                    bx = plan.precomputed[offsp+0];
721                    by = plan.precomputed[offsp+1];
722                    plan.stackbuf[offsb+0] = x*bx-y*by;
723                    plan.stackbuf[offsb+1] = -(x*by+y*bx);
724                    offsb = offsb+2;
725                    offsp = offsp+2;
726                }
727                ftbaseexecuteplanrec(ref plan.stackbuf, stackptr, ref plan, plan.plan[entryoffset+5], stackptr+2*2*m);
728                offsb = stackptr;
729                offsp = offs+2*m;
730                offsa = aoffset;
731                for(i=0; i<=n-1; i++)
732                {
733                    x = +(plan.stackbuf[offsb+0]/m);
734                    y = -(plan.stackbuf[offsb+1]/m);
735                    bx = plan.precomputed[offsp+0];
736                    by = plan.precomputed[offsp+1];
737                    a[offsa+0] = x*bx-y*-by;
738                    a[offsa+1] = x*-by+y*bx;
739                    offsp = offsp+2;
740                    offsa = offsa+2;
741                    offsb = offsb+2;
742                }
743                return;
744            }
745        }
746
747
748        /*************************************************************************
749        Returns good factorization N=N1*N2.
750
751        Usually N1<=N2 (but not always - small N's may be exception).
752        if N1<>1 then N2<>1.
753
754        Factorization is chosen depending on task type and codelets we have.
755
756          -- ALGLIB --
757             Copyright 01.05.2009 by Bochkanov Sergey
758        *************************************************************************/
759        public static void ftbasefactorize(int n,
760            int tasktype,
761            ref int n1,
762            ref int n2)
763        {
764            int j = 0;
765
766            n1 = 0;
767            n2 = 0;
768           
769            //
770            // try to find good codelet
771            //
772            if( n1*n2!=n )
773            {
774                for(j=ftbasecodeletrecommended; j>=2; j--)
775                {
776                    if( n%j==0 )
777                    {
778                        n1 = j;
779                        n2 = n/j;
780                        break;
781                    }
782                }
783            }
784           
785            //
786            // try to factorize N
787            //
788            if( n1*n2!=n )
789            {
790                for(j=ftbasecodeletrecommended+1; j<=n-1; j++)
791                {
792                    if( n%j==0 )
793                    {
794                        n1 = j;
795                        n2 = n/j;
796                        break;
797                    }
798                }
799            }
800           
801            //
802            // looks like N is prime :(
803            //
804            if( n1*n2!=n )
805            {
806                n1 = 1;
807                n2 = n;
808            }
809           
810            //
811            // normalize
812            //
813            if( n2==1 & n1!=1 )
814            {
815                n2 = n1;
816                n1 = 1;
817            }
818        }
819
820
821        /*************************************************************************
822        Is number smooth?
823
824          -- ALGLIB --
825             Copyright 01.05.2009 by Bochkanov Sergey
826        *************************************************************************/
827        public static bool ftbaseissmooth(int n)
828        {
829            bool result = new bool();
830            int i = 0;
831
832            for(i=2; i<=ftbasemaxsmoothfactor; i++)
833            {
834                while( n%i==0 )
835                {
836                    n = n/i;
837                }
838            }
839            result = n==1;
840            return result;
841        }
842
843
844        /*************************************************************************
845        Returns smallest smooth (divisible only by 2, 3, 5) number that is greater
846        than or equal to max(N,2)
847
848          -- ALGLIB --
849             Copyright 01.05.2009 by Bochkanov Sergey
850        *************************************************************************/
851        public static int ftbasefindsmooth(int n)
852        {
853            int result = 0;
854            int best = 0;
855
856            best = 2;
857            while( best<n )
858            {
859                best = 2*best;
860            }
861            ftbasefindsmoothrec(n, 1, 2, ref best);
862            result = best;
863            return result;
864        }
865
866
867        /*************************************************************************
868        Returns  smallest  smooth  (divisible only by 2, 3, 5) even number that is
869        greater than or equal to max(N,2)
870
871          -- ALGLIB --
872             Copyright 01.05.2009 by Bochkanov Sergey
873        *************************************************************************/
874        public static int ftbasefindsmootheven(int n)
875        {
876            int result = 0;
877            int best = 0;
878
879            best = 2;
880            while( best<n )
881            {
882                best = 2*best;
883            }
884            ftbasefindsmoothrec(n, 2, 2, ref best);
885            result = best;
886            return result;
887        }
888
889
890        /*************************************************************************
891        Returns estimate of FLOP count for the FFT.
892
893        It is only an estimate based on operations count for the PERFECT FFT
894        and relative inefficiency of the algorithm actually used.
895
896        N should be power of 2, estimates are badly wrong for non-power-of-2 N's.
897
898          -- ALGLIB --
899             Copyright 01.05.2009 by Bochkanov Sergey
900        *************************************************************************/
901        public static double ftbasegetflopestimate(int n)
902        {
903            double result = 0;
904
905            result = ftbaseinefficiencyfactor*(4*n*Math.Log(n)/Math.Log(2)-6*n+8);
906            return result;
907        }
908
909
910        /*************************************************************************
911        Recurrent subroutine for the FFTGeneratePlan:
912
913        PARAMETERS:
914            N                   plan size
915            IsReal              whether input is real or not.
916                                subroutine MUST NOT ignore this flag because real
917                                inputs comes with non-initialized imaginary parts,
918                                so ignoring this flag will result in corrupted output
919            HalfOut             whether full output or only half of it from 0 to
920                                floor(N/2) is needed. This flag may be ignored if
921                                doing so will simplify calculations
922            Plan                plan array
923            PlanSize            size of used part (in integers)
924            PrecomputedSize     size of precomputed array allocated yet
925            PlanArraySize       plan array size (actual)
926            TmpMemSize          temporary memory required size
927            BluesteinMemSize    temporary memory required size
928
929          -- ALGLIB --
930             Copyright 01.05.2009 by Bochkanov Sergey
931        *************************************************************************/
932        private static void ftbasegenerateplanrec(int n,
933            int tasktype,
934            ref ftplan plan,
935            ref int plansize,
936            ref int precomputedsize,
937            ref int planarraysize,
938            ref int tmpmemsize,
939            ref int stackmemsize,
940            int stackptr)
941        {
942            int j = 0;
943            int k = 0;
944            int m = 0;
945            int n1 = 0;
946            int n2 = 0;
947            int esize = 0;
948            int entryoffset = 0;
949
950           
951            //
952            // prepare
953            //
954            if( plansize+ftbaseplanentrysize>planarraysize )
955            {
956                fftarrayresize(ref plan.plan, ref planarraysize, 8*planarraysize);
957            }
958            entryoffset = plansize;
959            esize = ftbaseplanentrysize;
960            plansize = plansize+esize;
961           
962            //
963            // if N=1, generate empty plan and exit
964            //
965            if( n==1 )
966            {
967                plan.plan[entryoffset+0] = esize;
968                plan.plan[entryoffset+1] = -1;
969                plan.plan[entryoffset+2] = -1;
970                plan.plan[entryoffset+3] = fftemptyplan;
971                plan.plan[entryoffset+4] = -1;
972                plan.plan[entryoffset+5] = -1;
973                plan.plan[entryoffset+6] = -1;
974                plan.plan[entryoffset+7] = -1;
975                return;
976            }
977           
978            //
979            // generate plans
980            //
981            ftbasefactorize(n, tasktype, ref n1, ref n2);
982            if( tasktype==ftbasecffttask | tasktype==ftbaserffttask )
983            {
984               
985                //
986                // complex FFT plans
987                //
988                if( n1!=1 )
989                {
990                   
991                    //
992                    // Cooley-Tukey plan (real or complex)
993                    //
994                    // Note that child plans are COMPLEX
995                    // (whether plan itself is complex or not).
996                    //
997                    tmpmemsize = Math.Max(tmpmemsize, 2*n1*n2);
998                    plan.plan[entryoffset+0] = esize;
999                    plan.plan[entryoffset+1] = n1;
1000                    plan.plan[entryoffset+2] = n2;
1001                    if( tasktype==ftbasecffttask )
1002                    {
1003                        plan.plan[entryoffset+3] = fftcooleytukeyplan;
1004                    }
1005                    else
1006                    {
1007                        plan.plan[entryoffset+3] = fftrealcooleytukeyplan;
1008                    }
1009                    plan.plan[entryoffset+4] = 0;
1010                    plan.plan[entryoffset+5] = plansize;
1011                    ftbasegenerateplanrec(n1, ftbasecffttask, ref plan, ref plansize, ref precomputedsize, ref planarraysize, ref tmpmemsize, ref stackmemsize, stackptr);
1012                    plan.plan[entryoffset+6] = plansize;
1013                    ftbasegenerateplanrec(n2, ftbasecffttask, ref plan, ref plansize, ref precomputedsize, ref planarraysize, ref tmpmemsize, ref stackmemsize, stackptr);
1014                    plan.plan[entryoffset+7] = -1;
1015                    return;
1016                }
1017                else
1018                {
1019                    if( n==2 | n==3 | n==4 | n==5 )
1020                    {
1021                       
1022                        //
1023                        // hard-coded plan
1024                        //
1025                        plan.plan[entryoffset+0] = esize;
1026                        plan.plan[entryoffset+1] = n1;
1027                        plan.plan[entryoffset+2] = n2;
1028                        plan.plan[entryoffset+3] = fftcodeletplan;
1029                        plan.plan[entryoffset+4] = 0;
1030                        plan.plan[entryoffset+5] = -1;
1031                        plan.plan[entryoffset+6] = -1;
1032                        plan.plan[entryoffset+7] = precomputedsize;
1033                        if( n==3 )
1034                        {
1035                            precomputedsize = precomputedsize+2;
1036                        }
1037                        if( n==5 )
1038                        {
1039                            precomputedsize = precomputedsize+5;
1040                        }
1041                        return;
1042                    }
1043                    else
1044                    {
1045                       
1046                        //
1047                        // Bluestein's plan
1048                        //
1049                        // Select such M that M>=2*N-1, M is composite, and M's
1050                        // factors are 2, 3, 5
1051                        //
1052                        k = 2*n2-1;
1053                        m = ftbasefindsmooth(k);
1054                        tmpmemsize = Math.Max(tmpmemsize, 2*m);
1055                        plan.plan[entryoffset+0] = esize;
1056                        plan.plan[entryoffset+1] = n2;
1057                        plan.plan[entryoffset+2] = -1;
1058                        plan.plan[entryoffset+3] = fftbluesteinplan;
1059                        plan.plan[entryoffset+4] = m;
1060                        plan.plan[entryoffset+5] = plansize;
1061                        stackptr = stackptr+2*2*m;
1062                        stackmemsize = Math.Max(stackmemsize, stackptr);
1063                        ftbasegenerateplanrec(m, ftbasecffttask, ref plan, ref plansize, ref precomputedsize, ref planarraysize, ref tmpmemsize, ref stackmemsize, stackptr);
1064                        stackptr = stackptr-2*2*m;
1065                        plan.plan[entryoffset+6] = -1;
1066                        plan.plan[entryoffset+7] = precomputedsize;
1067                        precomputedsize = precomputedsize+2*m+2*n;
1068                        return;
1069                    }
1070                }
1071            }
1072            if( tasktype==ftbaserfhttask )
1073            {
1074               
1075                //
1076                // real FHT plans
1077                //
1078                if( n1!=1 )
1079                {
1080                   
1081                    //
1082                    // Cooley-Tukey plan
1083                    //
1084                    //
1085                    tmpmemsize = Math.Max(tmpmemsize, 2*n1*n2);
1086                    plan.plan[entryoffset+0] = esize;
1087                    plan.plan[entryoffset+1] = n1;
1088                    plan.plan[entryoffset+2] = n2;
1089                    plan.plan[entryoffset+3] = fhtcooleytukeyplan;
1090                    plan.plan[entryoffset+4] = 0;
1091                    plan.plan[entryoffset+5] = plansize;
1092                    ftbasegenerateplanrec(n1, tasktype, ref plan, ref plansize, ref precomputedsize, ref planarraysize, ref tmpmemsize, ref stackmemsize, stackptr);
1093                    plan.plan[entryoffset+6] = plansize;
1094                    ftbasegenerateplanrec(n2, tasktype, ref plan, ref plansize, ref precomputedsize, ref planarraysize, ref tmpmemsize, ref stackmemsize, stackptr);
1095                    plan.plan[entryoffset+7] = -1;
1096                    return;
1097                }
1098                else
1099                {
1100                   
1101                    //
1102                    // N2 plan
1103                    //
1104                    plan.plan[entryoffset+0] = esize;
1105                    plan.plan[entryoffset+1] = n1;
1106                    plan.plan[entryoffset+2] = n2;
1107                    plan.plan[entryoffset+3] = fhtn2plan;
1108                    plan.plan[entryoffset+4] = 0;
1109                    plan.plan[entryoffset+5] = -1;
1110                    plan.plan[entryoffset+6] = -1;
1111                    plan.plan[entryoffset+7] = -1;
1112                    if( n==2 | n==3 | n==4 | n==5 )
1113                    {
1114                       
1115                        //
1116                        // hard-coded plan
1117                        //
1118                        plan.plan[entryoffset+0] = esize;
1119                        plan.plan[entryoffset+1] = n1;
1120                        plan.plan[entryoffset+2] = n2;
1121                        plan.plan[entryoffset+3] = fhtcodeletplan;
1122                        plan.plan[entryoffset+4] = 0;
1123                        plan.plan[entryoffset+5] = -1;
1124                        plan.plan[entryoffset+6] = -1;
1125                        plan.plan[entryoffset+7] = precomputedsize;
1126                        if( n==3 )
1127                        {
1128                            precomputedsize = precomputedsize+2;
1129                        }
1130                        if( n==5 )
1131                        {
1132                            precomputedsize = precomputedsize+5;
1133                        }
1134                        return;
1135                    }
1136                    return;
1137                }
1138            }
1139        }
1140
1141
1142        /*************************************************************************
1143        Recurrent subroutine for precomputing FFT plans
1144
1145          -- ALGLIB --
1146             Copyright 01.05.2009 by Bochkanov Sergey
1147        *************************************************************************/
1148        private static void ftbaseprecomputeplanrec(ref ftplan plan,
1149            int entryoffset,
1150            int stackptr)
1151        {
1152            int i = 0;
1153            int j = 0;
1154            int idx = 0;
1155            int n1 = 0;
1156            int n2 = 0;
1157            int n = 0;
1158            int m = 0;
1159            int offs = 0;
1160            double v = 0;
1161            double[] emptyarray = new double[0];
1162            double bx = 0;
1163            double by = 0;
1164
1165            if( plan.plan[entryoffset+3]==fftcooleytukeyplan | plan.plan[entryoffset+3]==fftrealcooleytukeyplan | plan.plan[entryoffset+3]==fhtcooleytukeyplan )
1166            {
1167                ftbaseprecomputeplanrec(ref plan, plan.plan[entryoffset+5], stackptr);
1168                ftbaseprecomputeplanrec(ref plan, plan.plan[entryoffset+6], stackptr);
1169                return;
1170            }
1171            if( plan.plan[entryoffset+3]==fftcodeletplan | plan.plan[entryoffset+3]==fhtcodeletplan )
1172            {
1173                n1 = plan.plan[entryoffset+1];
1174                n2 = plan.plan[entryoffset+2];
1175                n = n1*n2;
1176                if( n==3 )
1177                {
1178                    offs = plan.plan[entryoffset+7];
1179                    plan.precomputed[offs+0] = Math.Cos(2*Math.PI/3)-1;
1180                    plan.precomputed[offs+1] = Math.Sin(2*Math.PI/3);
1181                    return;
1182                }
1183                if( n==5 )
1184                {
1185                    offs = plan.plan[entryoffset+7];
1186                    v = 2*Math.PI/5;
1187                    plan.precomputed[offs+0] = (Math.Cos(v)+Math.Cos(2*v))/2-1;
1188                    plan.precomputed[offs+1] = (Math.Cos(v)-Math.Cos(2*v))/2;
1189                    plan.precomputed[offs+2] = -Math.Sin(v);
1190                    plan.precomputed[offs+3] = -(Math.Sin(v)+Math.Sin(2*v));
1191                    plan.precomputed[offs+4] = Math.Sin(v)-Math.Sin(2*v);
1192                    return;
1193                }
1194            }
1195            if( plan.plan[entryoffset+3]==fftbluesteinplan )
1196            {
1197                ftbaseprecomputeplanrec(ref plan, plan.plan[entryoffset+5], stackptr);
1198                n = plan.plan[entryoffset+1];
1199                m = plan.plan[entryoffset+4];
1200                offs = plan.plan[entryoffset+7];
1201                for(i=0; i<=2*m-1; i++)
1202                {
1203                    plan.precomputed[offs+i] = 0;
1204                }
1205                for(i=0; i<=n-1; i++)
1206                {
1207                    bx = Math.Cos(Math.PI*AP.Math.Sqr(i)/n);
1208                    by = Math.Sin(Math.PI*AP.Math.Sqr(i)/n);
1209                    plan.precomputed[offs+2*i+0] = bx;
1210                    plan.precomputed[offs+2*i+1] = by;
1211                    plan.precomputed[offs+2*m+2*i+0] = bx;
1212                    plan.precomputed[offs+2*m+2*i+1] = by;
1213                    if( i>0 )
1214                    {
1215                        plan.precomputed[offs+2*(m-i)+0] = bx;
1216                        plan.precomputed[offs+2*(m-i)+1] = by;
1217                    }
1218                }
1219                ftbaseexecuteplanrec(ref plan.precomputed, offs, ref plan, plan.plan[entryoffset+5], stackptr);
1220                return;
1221            }
1222        }
1223
1224
1225        /*************************************************************************
1226        Twiddle factors calculation
1227
1228          -- ALGLIB --
1229             Copyright 01.05.2009 by Bochkanov Sergey
1230        *************************************************************************/
1231        private static void ffttwcalc(ref double[] a,
1232            int aoffset,
1233            int n1,
1234            int n2)
1235        {
1236            int i = 0;
1237            int j = 0;
1238            int n = 0;
1239            int idx = 0;
1240            int offs = 0;
1241            double x = 0;
1242            double y = 0;
1243            double twxm1 = 0;
1244            double twy = 0;
1245            double twbasexm1 = 0;
1246            double twbasey = 0;
1247            double twrowxm1 = 0;
1248            double twrowy = 0;
1249            double tmpx = 0;
1250            double tmpy = 0;
1251            double v = 0;
1252
1253            n = n1*n2;
1254            v = -(2*Math.PI/n);
1255            twbasexm1 = -(2*AP.Math.Sqr(Math.Sin(0.5*v)));
1256            twbasey = Math.Sin(v);
1257            twrowxm1 = 0;
1258            twrowy = 0;
1259            for(i=0; i<=n2-1; i++)
1260            {
1261                twxm1 = 0;
1262                twy = 0;
1263                for(j=0; j<=n1-1; j++)
1264                {
1265                    idx = i*n1+j;
1266                    offs = aoffset+2*idx;
1267                    x = a[offs+0];
1268                    y = a[offs+1];
1269                    tmpx = x*twxm1-y*twy;
1270                    tmpy = x*twy+y*twxm1;
1271                    a[offs+0] = x+tmpx;
1272                    a[offs+1] = y+tmpy;
1273                   
1274                    //
1275                    // update Tw: Tw(new) = Tw(old)*TwRow
1276                    //
1277                    if( j<n1-1 )
1278                    {
1279                        if( j%ftbaseupdatetw==0 )
1280                        {
1281                            v = -(2*Math.PI*i*(j+1)/n);
1282                            twxm1 = -(2*AP.Math.Sqr(Math.Sin(0.5*v)));
1283                            twy = Math.Sin(v);
1284                        }
1285                        else
1286                        {
1287                            tmpx = twrowxm1+twxm1*twrowxm1-twy*twrowy;
1288                            tmpy = twrowy+twxm1*twrowy+twy*twrowxm1;
1289                            twxm1 = twxm1+tmpx;
1290                            twy = twy+tmpy;
1291                        }
1292                    }
1293                }
1294               
1295                //
1296                // update TwRow: TwRow(new) = TwRow(old)*TwBase
1297                //
1298                if( i<n2-1 )
1299                {
1300                    if( j%ftbaseupdatetw==0 )
1301                    {
1302                        v = -(2*Math.PI*(i+1)/n);
1303                        twrowxm1 = -(2*AP.Math.Sqr(Math.Sin(0.5*v)));
1304                        twrowy = Math.Sin(v);
1305                    }
1306                    else
1307                    {
1308                        tmpx = twbasexm1+twrowxm1*twbasexm1-twrowy*twbasey;
1309                        tmpy = twbasey+twrowxm1*twbasey+twrowy*twbasexm1;
1310                        twrowxm1 = twrowxm1+tmpx;
1311                        twrowy = twrowy+tmpy;
1312                    }
1313                }
1314            }
1315        }
1316
1317
1318        /*************************************************************************
1319        Linear transpose: transpose complex matrix stored in 1-dimensional array
1320
1321          -- ALGLIB --
1322             Copyright 01.05.2009 by Bochkanov Sergey
1323        *************************************************************************/
1324        private static void internalcomplexlintranspose(ref double[] a,
1325            int m,
1326            int n,
1327            int astart,
1328            ref double[] buf)
1329        {
1330            int i_ = 0;
1331            int i1_ = 0;
1332
1333            ffticltrec(ref a, astart, n, ref buf, 0, m, m, n);
1334            i1_ = (0) - (astart);
1335            for(i_=astart; i_<=astart+2*m*n-1;i_++)
1336            {
1337                a[i_] = buf[i_+i1_];
1338            }
1339        }
1340
1341
1342        /*************************************************************************
1343        Linear transpose: transpose real matrix stored in 1-dimensional array
1344
1345          -- ALGLIB --
1346             Copyright 01.05.2009 by Bochkanov Sergey
1347        *************************************************************************/
1348        private static void internalreallintranspose(ref double[] a,
1349            int m,
1350            int n,
1351            int astart,
1352            ref double[] buf)
1353        {
1354            int i_ = 0;
1355            int i1_ = 0;
1356
1357            fftirltrec(ref a, astart, n, ref buf, 0, m, m, n);
1358            i1_ = (0) - (astart);
1359            for(i_=astart; i_<=astart+m*n-1;i_++)
1360            {
1361                a[i_] = buf[i_+i1_];
1362            }
1363        }
1364
1365
1366        /*************************************************************************
1367        Recurrent subroutine for a InternalComplexLinTranspose
1368
1369        Write A^T to B, where:
1370        * A is m*n complex matrix stored in array A as pairs of real/image values,
1371          beginning from AStart position, with AStride stride
1372        * B is n*m complex matrix stored in array B as pairs of real/image values,
1373          beginning from BStart position, with BStride stride
1374        stride is measured in complex numbers, i.e. in real/image pairs.
1375
1376          -- ALGLIB --
1377             Copyright 01.05.2009 by Bochkanov Sergey
1378        *************************************************************************/
1379        private static void ffticltrec(ref double[] a,
1380            int astart,
1381            int astride,
1382            ref double[] b,
1383            int bstart,
1384            int bstride,
1385            int m,
1386            int n)
1387        {
1388            int i = 0;
1389            int j = 0;
1390            int idx1 = 0;
1391            int idx2 = 0;
1392            int m2 = 0;
1393            int m1 = 0;
1394            int n1 = 0;
1395
1396            if( m==0 | n==0 )
1397            {
1398                return;
1399            }
1400            if( Math.Max(m, n)<=8 )
1401            {
1402                m2 = 2*bstride;
1403                for(i=0; i<=m-1; i++)
1404                {
1405                    idx1 = bstart+2*i;
1406                    idx2 = astart+2*i*astride;
1407                    for(j=0; j<=n-1; j++)
1408                    {
1409                        b[idx1+0] = a[idx2+0];
1410                        b[idx1+1] = a[idx2+1];
1411                        idx1 = idx1+m2;
1412                        idx2 = idx2+2;
1413                    }
1414                }
1415                return;
1416            }
1417            if( n>m )
1418            {
1419               
1420                //
1421                // New partition:
1422                //
1423                // "A^T -> B" becomes "(A1 A2)^T -> ( B1 )
1424                //                                  ( B2 )
1425                //
1426                n1 = n/2;
1427                if( n-n1>=8 & n1%8!=0 )
1428                {
1429                    n1 = n1+(8-n1%8);
1430                }
1431                System.Diagnostics.Debug.Assert(n-n1>0);
1432                ffticltrec(ref a, astart, astride, ref b, bstart, bstride, m, n1);
1433                ffticltrec(ref a, astart+2*n1, astride, ref b, bstart+2*n1*bstride, bstride, m, n-n1);
1434            }
1435            else
1436            {
1437               
1438                //
1439                // New partition:
1440                //
1441                // "A^T -> B" becomes "( A1 )^T -> ( B1 B2 )
1442                //                     ( A2 )
1443                //
1444                m1 = m/2;
1445                if( m-m1>=8 & m1%8!=0 )
1446                {
1447                    m1 = m1+(8-m1%8);
1448                }
1449                System.Diagnostics.Debug.Assert(m-m1>0);
1450                ffticltrec(ref a, astart, astride, ref b, bstart, bstride, m1, n);
1451                ffticltrec(ref a, astart+2*m1*astride, astride, ref b, bstart+2*m1, bstride, m-m1, n);
1452            }
1453        }
1454
1455
1456        /*************************************************************************
1457        Recurrent subroutine for a InternalRealLinTranspose
1458
1459
1460          -- ALGLIB --
1461             Copyright 01.05.2009 by Bochkanov Sergey
1462        *************************************************************************/
1463        private static void fftirltrec(ref double[] a,
1464            int astart,
1465            int astride,
1466            ref double[] b,
1467            int bstart,
1468            int bstride,
1469            int m,
1470            int n)
1471        {
1472            int i = 0;
1473            int j = 0;
1474            int idx1 = 0;
1475            int idx2 = 0;
1476            int m1 = 0;
1477            int n1 = 0;
1478
1479            if( m==0 | n==0 )
1480            {
1481                return;
1482            }
1483            if( Math.Max(m, n)<=8 )
1484            {
1485                for(i=0; i<=m-1; i++)
1486                {
1487                    idx1 = bstart+i;
1488                    idx2 = astart+i*astride;
1489                    for(j=0; j<=n-1; j++)
1490                    {
1491                        b[idx1] = a[idx2];
1492                        idx1 = idx1+bstride;
1493                        idx2 = idx2+1;
1494                    }
1495                }
1496                return;
1497            }
1498            if( n>m )
1499            {
1500               
1501                //
1502                // New partition:
1503                //
1504                // "A^T -> B" becomes "(A1 A2)^T -> ( B1 )
1505                //                                  ( B2 )
1506                //
1507                n1 = n/2;
1508                if( n-n1>=8 & n1%8!=0 )
1509                {
1510                    n1 = n1+(8-n1%8);
1511                }
1512                System.Diagnostics.Debug.Assert(n-n1>0);
1513                fftirltrec(ref a, astart, astride, ref b, bstart, bstride, m, n1);
1514                fftirltrec(ref a, astart+n1, astride, ref b, bstart+n1*bstride, bstride, m, n-n1);
1515            }
1516            else
1517            {
1518               
1519                //
1520                // New partition:
1521                //
1522                // "A^T -> B" becomes "( A1 )^T -> ( B1 B2 )
1523                //                     ( A2 )
1524                //
1525                m1 = m/2;
1526                if( m-m1>=8 & m1%8!=0 )
1527                {
1528                    m1 = m1+(8-m1%8);
1529                }
1530                System.Diagnostics.Debug.Assert(m-m1>0);
1531                fftirltrec(ref a, astart, astride, ref b, bstart, bstride, m1, n);
1532                fftirltrec(ref a, astart+m1*astride, astride, ref b, bstart+m1, bstride, m-m1, n);
1533            }
1534        }
1535
1536
1537        /*************************************************************************
1538        recurrent subroutine for FFTFindSmoothRec
1539
1540          -- ALGLIB --
1541             Copyright 01.05.2009 by Bochkanov Sergey
1542        *************************************************************************/
1543        private static void ftbasefindsmoothrec(int n,
1544            int seed,
1545            int leastfactor,
1546            ref int best)
1547        {
1548            System.Diagnostics.Debug.Assert(ftbasemaxsmoothfactor<=5, "FTBaseFindSmoothRec: internal error!");
1549            if( seed>=n )
1550            {
1551                best = Math.Min(best, seed);
1552                return;
1553            }
1554            if( leastfactor<=2 )
1555            {
1556                ftbasefindsmoothrec(n, seed*2, 2, ref best);
1557            }
1558            if( leastfactor<=3 )
1559            {
1560                ftbasefindsmoothrec(n, seed*3, 3, ref best);
1561            }
1562            if( leastfactor<=5 )
1563            {
1564                ftbasefindsmoothrec(n, seed*5, 5, ref best);
1565            }
1566        }
1567
1568
1569        /*************************************************************************
1570        Internal subroutine: array resize
1571
1572          -- ALGLIB --
1573             Copyright 01.05.2009 by Bochkanov Sergey
1574        *************************************************************************/
1575        private static void fftarrayresize(ref int[] a,
1576            ref int asize,
1577            int newasize)
1578        {
1579            int[] tmp = new int[0];
1580            int i = 0;
1581
1582            tmp = new int[asize];
1583            for(i=0; i<=asize-1; i++)
1584            {
1585                tmp[i] = a[i];
1586            }
1587            a = new int[newasize];
1588            for(i=0; i<=asize-1; i++)
1589            {
1590                a[i] = tmp[i];
1591            }
1592            asize = newasize;
1593        }
1594
1595
1596        /*************************************************************************
1597        Reference FHT stub
1598        *************************************************************************/
1599        private static void reffht(ref double[] a,
1600            int n,
1601            int offs)
1602        {
1603            double[] buf = new double[0];
1604            int i = 0;
1605            int j = 0;
1606            double v = 0;
1607
1608            System.Diagnostics.Debug.Assert(n>0, "RefFHTR1D: incorrect N!");
1609            buf = new double[n];
1610            for(i=0; i<=n-1; i++)
1611            {
1612                v = 0;
1613                for(j=0; j<=n-1; j++)
1614                {
1615                    v = v+a[offs+j]*(Math.Cos(2*Math.PI*i*j/n)+Math.Sin(2*Math.PI*i*j/n));
1616                }
1617                buf[i] = v;
1618            }
1619            for(i=0; i<=n-1; i++)
1620            {
1621                a[offs+i] = buf[i];
1622            }
1623        }
1624    }
1625}
Note: See TracBrowser for help on using the repository browser.