Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/ftbase.cs @ 3021

Last change on this file since 3021 was 2806, checked in by gkronber, 15 years ago

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