00001
00002
00003
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099 #define WANT_STREAM
00100
00101 #define WANT_MATH
00102
00103 #include "newmatap.h"
00104
00105 #ifdef use_namespace
00106 namespace NEWMAT {
00107 #endif
00108
00109 #ifdef DO_REPORT
00110 #define REPORT { static ExeCounter ExeCount(__LINE__,20); ++ExeCount; }
00111 #else
00112 #define REPORT {}
00113 #endif
00114
00115 inline Real square(Real x) { return x*x; }
00116 inline int square(int x) { return x*x; }
00117
00118 static void GR_1D_FS (int PTS, int N_SYM, int N_UN_SYM,
00119 const SimpleIntArray& SYM, int P_SYM, const SimpleIntArray& UN_SYM,
00120 Real* X, Real* Y);
00121 static void GR_1D_FT (int N, int N_FACTOR, const SimpleIntArray& FACTOR,
00122 Real* X, Real* Y);
00123 static void R_P_FTK (int N, int M, int P, Real* X, Real* Y);
00124 static void R_2_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1);
00125 static void R_3_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1,
00126 Real* X2, Real* Y2);
00127 static void R_4_FTK (int N, int M,
00128 Real* X0, Real* Y0, Real* X1, Real* Y1,
00129 Real* X2, Real* Y2, Real* X3, Real* Y3);
00130 static void R_5_FTK (int N, int M,
00131 Real* X0, Real* Y0, Real* X1, Real* Y1, Real* X2, Real* Y2,
00132 Real* X3, Real* Y3, Real* X4, Real* Y4);
00133 static void R_8_FTK (int N, int M,
00134 Real* X0, Real* Y0, Real* X1, Real* Y1,
00135 Real* X2, Real* Y2, Real* X3, Real* Y3,
00136 Real* X4, Real* Y4, Real* X5, Real* Y5,
00137 Real* X6, Real* Y6, Real* X7, Real* Y7);
00138 static void R_16_FTK (int N, int M,
00139 Real* X0, Real* Y0, Real* X1, Real* Y1,
00140 Real* X2, Real* Y2, Real* X3, Real* Y3,
00141 Real* X4, Real* Y4, Real* X5, Real* Y5,
00142 Real* X6, Real* Y6, Real* X7, Real* Y7,
00143 Real* X8, Real* Y8, Real* X9, Real* Y9,
00144 Real* X10, Real* Y10, Real* X11, Real* Y11,
00145 Real* X12, Real* Y12, Real* X13, Real* Y13,
00146 Real* X14, Real* Y14, Real* X15, Real* Y15);
00147 static int BitReverse(int x, int prod, int n, const SimpleIntArray& f);
00148
00149
00150 bool FFT_Controller::ar_1d_ft (int PTS, Real* X, Real *Y)
00151 {
00152
00153
00154 REPORT
00155
00156 int F,J,N,NF,P,PMAX,P_SYM,P_TWO,Q,R,TWO_GRP;
00157
00158
00159
00160 const int NP = 16, NQ = 10;
00161 SimpleIntArray PP(NP), QQ(NQ);
00162
00163 TWO_GRP=16; PMAX=19;
00164
00165
00166
00167
00168
00169
00170 if (PTS<=1) return true;
00171 N=PTS; P_SYM=1; F=2; P=0; Q=0;
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187 while (N > 1)
00188 {
00189 bool fail = true;
00190 for (J=F; J<=PMAX; J++)
00191 if (N % J == 0) { fail = false; F=J; break; }
00192 if (fail || P >= NP || Q >= NQ) return false;
00193 N /= F;
00194 if (N % F != 0) QQ[Q++] = F;
00195 else { N /= F; PP[P++] = F; P_SYM *= F; }
00196 }
00197
00198 R = (Q == 0) ? 0 : 1;
00199
00200 NF = 2*P + Q;
00201 SimpleIntArray FACTOR(NF + 1), SYM(2*P + R);
00202 FACTOR[NF] = 0;
00203
00204
00205 for (J=0; J<P; J++)
00206 { SYM[J]=FACTOR[J]=PP[P-1-J]; FACTOR[P+Q+J]=SYM[P+R+J]=PP[J]; }
00207
00208 if (Q>0)
00209 {
00210 REPORT
00211 for (J=0; J<Q; J++) FACTOR[P+J]=QQ[J];
00212 SYM[P]=PTS/square(P_SYM);
00213 }
00214
00215
00216 P_TWO = 1;
00217 for (J=0; J < NF; J++)
00218 {
00219 if (FACTOR[J]!=2) continue;
00220 P_TWO=P_TWO*2; FACTOR[J]=1;
00221 if (P_TWO<TWO_GRP && FACTOR[J+1]==2) continue;
00222 FACTOR[J]=P_TWO; P_TWO=1;
00223 }
00224
00225 if (P==0) R=0;
00226 if (Q<=1) Q=0;
00227
00228
00229 GR_1D_FT(PTS,NF,FACTOR,X,Y);
00230 GR_1D_FS(PTS,2*P+R,Q,SYM,P_SYM,QQ,X,Y);
00231
00232 return true;
00233
00234 }
00235
00236 static void GR_1D_FS (int PTS, int N_SYM, int N_UN_SYM,
00237 const SimpleIntArray& SYM, int P_SYM, const SimpleIntArray& UN_SYM,
00238 Real* X, Real* Y)
00239 {
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249 REPORT
00250
00251 Real T;
00252 int JJ,KK,P_UN_SYM;
00253
00254
00255
00256
00257 if (N_SYM > 0)
00258 {
00259 REPORT
00260 SimpleIntArray U(N_SYM);
00261 for(MultiRadixCounter MRC(N_SYM, SYM, U); !MRC.Finish(); ++MRC)
00262 {
00263 if (MRC.Swap())
00264 {
00265 int P = MRC.Reverse(); int JJ = MRC.Counter(); Real T;
00266 T=X[JJ]; X[JJ]=X[P]; X[P]=T; T=Y[JJ]; Y[JJ]=Y[P]; Y[P]=T;
00267 }
00268 }
00269 }
00270
00271 int J,JL,K,L,M,MS;
00272
00273
00274
00275
00276
00277
00278
00279 if (N_UN_SYM==0) { REPORT return; }
00280 P_UN_SYM=PTS/square(P_SYM); JL=(P_UN_SYM-3)*P_SYM; MS=P_UN_SYM*P_SYM;
00281
00282 for (J = P_SYM; J<=JL; J+=P_SYM)
00283 {
00284 K=J;
00285 do K = P_SYM * BitReverse(K / P_SYM, P_UN_SYM, N_UN_SYM, UN_SYM);
00286 while (K<J);
00287
00288 if (K!=J)
00289 {
00290 REPORT
00291 for (L=0; L<P_SYM; L++) for (M=L; M<PTS; M+=MS)
00292 {
00293 JJ=M+J; KK=M+K;
00294 T=X[JJ]; X[JJ]=X[KK]; X[KK]=T; T=Y[JJ]; Y[JJ]=Y[KK]; Y[KK]=T;
00295 }
00296 }
00297 }
00298
00299 return;
00300 }
00301
00302 static void GR_1D_FT (int N, int N_FACTOR, const SimpleIntArray& FACTOR,
00303 Real* X, Real* Y)
00304 {
00305
00306
00307 REPORT
00308
00309 int M = N;
00310
00311 for (int i = 0; i < N_FACTOR; i++)
00312 {
00313 int P = FACTOR[i]; M /= P;
00314
00315 switch(P)
00316 {
00317 case 1: REPORT break;
00318 case 2: REPORT R_2_FTK (N,M,X,Y,X+M,Y+M); break;
00319 case 3: REPORT R_3_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M); break;
00320 case 4: REPORT R_4_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,X+3*M,Y+3*M); break;
00321 case 5:
00322 REPORT
00323 R_5_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,X+3*M,Y+3*M,X+4*M,Y+4*M);
00324 break;
00325 case 8:
00326 REPORT
00327 R_8_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,
00328 X+3*M,Y+3*M,X+4*M,Y+4*M,X+5*M,Y+5*M,
00329 X+6*M,Y+6*M,X+7*M,Y+7*M);
00330 break;
00331 case 16:
00332 REPORT
00333 R_16_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,
00334 X+3*M,Y+3*M,X+4*M,Y+4*M,X+5*M,Y+5*M,
00335 X+6*M,Y+6*M,X+7*M,Y+7*M,X+8*M,Y+8*M,
00336 X+9*M,Y+9*M,X+10*M,Y+10*M,X+11*M,Y+11*M,
00337 X+12*M,Y+12*M,X+13*M,Y+13*M,X+14*M,Y+14*M,
00338 X+15*M,Y+15*M);
00339 break;
00340 default: REPORT R_P_FTK (N,M,P,X,Y); break;
00341 }
00342 }
00343
00344 }
00345
00346 static void R_P_FTK (int N, int M, int P, Real* X, Real* Y)
00347
00348
00349 {
00350 REPORT
00351 bool NO_FOLD,ZERO;
00352 Real ANGLE,IS,IU,RS,RU,T,TWOPI,XT,YT;
00353 int J,JJ,K0,K,M_OVER_2,MP,PM,PP,U,V;
00354
00355 Real AA [9][9], BB [9][9];
00356 Real A [18], B [18], C [18], S [18];
00357 Real IA [9], IB [9], RA [9], RB [9];
00358
00359 TWOPI=8.0*atan(1.0);
00360 M_OVER_2=M/2+1; MP=M*P; PP=P/2; PM=P-1;
00361
00362 for (U=0; U<PP; U++)
00363 {
00364 ANGLE=TWOPI*Real(U+1)/Real(P);
00365 JJ=P-U-2;
00366 A[U]=cos(ANGLE); B[U]=sin(ANGLE);
00367 A[JJ]=A[U]; B[JJ]= -B[U];
00368 }
00369
00370 for (U=1; U<=PP; U++)
00371 {
00372 for (V=1; V<=PP; V++)
00373 { JJ=U*V-U*V/P*P; AA[V-1][U-1]=A[JJ-1]; BB[V-1][U-1]=B[JJ-1]; }
00374 }
00375
00376 for (J=0; J<M_OVER_2; J++)
00377 {
00378 NO_FOLD = (J==0 || 2*J==M);
00379 K0=J;
00380 ANGLE=TWOPI*Real(J)/Real(MP); ZERO=ANGLE==0.0;
00381 C[0]=cos(ANGLE); S[0]=sin(ANGLE);
00382 for (U=1; U<PM; U++)
00383 {
00384 C[U]=C[U-1]*C[0]-S[U-1]*S[0];
00385 S[U]=S[U-1]*C[0]+C[U-1]*S[0];
00386 }
00387 goto L700;
00388 L500:
00389 REPORT
00390 if (NO_FOLD) { REPORT goto L1500; }
00391 REPORT
00392 NO_FOLD=true; K0=M-J;
00393 for (U=0; U<PM; U++)
00394 { T=C[U]*A[U]+S[U]*B[U]; S[U]= -S[U]*A[U]+C[U]*B[U]; C[U]=T; }
00395 L700:
00396 REPORT
00397 for (K=K0; K<N; K+=MP)
00398 {
00399 XT=X[K]; YT=Y[K];
00400 for (U=1; U<=PP; U++)
00401 {
00402 RA[U-1]=XT; IA[U-1]=YT;
00403 RB[U-1]=0.0; IB[U-1]=0.0;
00404 }
00405 for (U=1; U<=PP; U++)
00406 {
00407 JJ=P-U;
00408 RS=X[K+M*U]+X[K+M*JJ]; IS=Y[K+M*U]+Y[K+M*JJ];
00409 RU=X[K+M*U]-X[K+M*JJ]; IU=Y[K+M*U]-Y[K+M*JJ];
00410 XT=XT+RS; YT=YT+IS;
00411 for (V=0; V<PP; V++)
00412 {
00413 RA[V]=RA[V]+RS*AA[V][U-1]; IA[V]=IA[V]+IS*AA[V][U-1];
00414 RB[V]=RB[V]+RU*BB[V][U-1]; IB[V]=IB[V]+IU*BB[V][U-1];
00415 }
00416 }
00417 X[K]=XT; Y[K]=YT;
00418 for (U=1; U<=PP; U++)
00419 {
00420 if (!ZERO)
00421 {
00422 REPORT
00423 XT=RA[U-1]+IB[U-1]; YT=IA[U-1]-RB[U-1];
00424 X[K+M*U]=XT*C[U-1]+YT*S[U-1]; Y[K+M*U]=YT*C[U-1]-XT*S[U-1];
00425 JJ=P-U;
00426 XT=RA[U-1]-IB[U-1]; YT=IA[U-1]+RB[U-1];
00427 X[K+M*JJ]=XT*C[JJ-1]+YT*S[JJ-1];
00428 Y[K+M*JJ]=YT*C[JJ-1]-XT*S[JJ-1];
00429 }
00430 else
00431 {
00432 REPORT
00433 X[K+M*U]=RA[U-1]+IB[U-1]; Y[K+M*U]=IA[U-1]-RB[U-1];
00434 JJ=P-U;
00435 X[K+M*JJ]=RA[U-1]-IB[U-1]; Y[K+M*JJ]=IA[U-1]+RB[U-1];
00436 }
00437 }
00438 }
00439 goto L500;
00440 L1500: ;
00441 }
00442 return;
00443 }
00444
00445 static void R_2_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1)
00446
00447 {
00448 REPORT
00449 bool NO_FOLD,ZERO;
00450 int J,K,K0,M2,M_OVER_2;
00451 Real ANGLE,C,IS,IU,RS,RU,S,TWOPI;
00452
00453 M2=M*2; M_OVER_2=M/2+1;
00454 TWOPI=8.0*atan(1.0);
00455
00456 for (J=0; J<M_OVER_2; J++)
00457 {
00458 NO_FOLD = (J==0 || 2*J==M);
00459 K0=J;
00460 ANGLE=TWOPI*Real(J)/Real(M2); ZERO=ANGLE==0.0;
00461 C=cos(ANGLE); S=sin(ANGLE);
00462 goto L200;
00463 L100:
00464 REPORT
00465 if (NO_FOLD) { REPORT goto L600; }
00466 REPORT
00467 NO_FOLD=true; K0=M-J; C= -C;
00468 L200:
00469 REPORT
00470 for (K=K0; K<N; K+=M2)
00471 {
00472 RS=X0[K]+X1[K]; IS=Y0[K]+Y1[K];
00473 RU=X0[K]-X1[K]; IU=Y0[K]-Y1[K];
00474 X0[K]=RS; Y0[K]=IS;
00475 if (!ZERO) { X1[K]=RU*C+IU*S; Y1[K]=IU*C-RU*S; }
00476 else { X1[K]=RU; Y1[K]=IU; }
00477 }
00478 goto L100;
00479 L600: ;
00480 }
00481
00482 return;
00483 }
00484
00485 static void R_3_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1,
00486 Real* X2, Real* Y2)
00487
00488 {
00489 REPORT
00490 bool NO_FOLD,ZERO;
00491 int J,K,K0,M3,M_OVER_2;
00492 Real ANGLE,A,B,C1,C2,S1,S2,T,TWOPI;
00493 Real I0,I1,I2,IA,IB,IS,R0,R1,R2,RA,RB,RS;
00494
00495 M3=M*3; M_OVER_2=M/2+1; TWOPI=8.0*atan(1.0);
00496 A=cos(TWOPI/3.0); B=sin(TWOPI/3.0);
00497
00498 for (J=0; J<M_OVER_2; J++)
00499 {
00500 NO_FOLD = (J==0 || 2*J==M);
00501 K0=J;
00502 ANGLE=TWOPI*Real(J)/Real(M3); ZERO=ANGLE==0.0;
00503 C1=cos(ANGLE); S1=sin(ANGLE);
00504 C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
00505 goto L200;
00506 L100:
00507 REPORT
00508 if (NO_FOLD) { REPORT goto L600; }
00509 REPORT
00510 NO_FOLD=true; K0=M-J;
00511 T=C1*A+S1*B; S1=C1*B-S1*A; C1=T;
00512 T=C2*A-S2*B; S2= -C2*B-S2*A; C2=T;
00513 L200:
00514 REPORT
00515 for (K=K0; K<N; K+=M3)
00516 {
00517 R0 = X0[K]; I0 = Y0[K];
00518 RS=X1[K]+X2[K]; IS=Y1[K]+Y2[K];
00519 X0[K]=R0+RS; Y0[K]=I0+IS;
00520 RA=R0+RS*A; IA=I0+IS*A;
00521 RB=(X1[K]-X2[K])*B; IB=(Y1[K]-Y2[K])*B;
00522 if (!ZERO)
00523 {
00524 REPORT
00525 R1=RA+IB; I1=IA-RB; R2=RA-IB; I2=IA+RB;
00526 X1[K]=R1*C1+I1*S1; Y1[K]=I1*C1-R1*S1;
00527 X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
00528 }
00529 else { REPORT X1[K]=RA+IB; Y1[K]=IA-RB; X2[K]=RA-IB; Y2[K]=IA+RB; }
00530 }
00531 goto L100;
00532 L600: ;
00533 }
00534
00535 return;
00536 }
00537
00538 static void R_4_FTK (int N, int M,
00539 Real* X0, Real* Y0, Real* X1, Real* Y1,
00540 Real* X2, Real* Y2, Real* X3, Real* Y3)
00541
00542 {
00543 REPORT
00544 bool NO_FOLD,ZERO;
00545 int J,K,K0,M4,M_OVER_2;
00546 Real ANGLE,C1,C2,C3,S1,S2,S3,T,TWOPI;
00547 Real I1,I2,I3,IS0,IS1,IU0,IU1,R1,R2,R3,RS0,RS1,RU0,RU1;
00548
00549 M4=M*4; M_OVER_2=M/2+1;
00550 TWOPI=8.0*atan(1.0);
00551
00552 for (J=0; J<M_OVER_2; J++)
00553 {
00554 NO_FOLD = (J==0 || 2*J==M);
00555 K0=J;
00556 ANGLE=TWOPI*Real(J)/Real(M4); ZERO=ANGLE==0.0;
00557 C1=cos(ANGLE); S1=sin(ANGLE);
00558 C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
00559 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00560 goto L200;
00561 L100:
00562 REPORT
00563 if (NO_FOLD) { REPORT goto L600; }
00564 REPORT
00565 NO_FOLD=true; K0=M-J;
00566 T=C1; C1=S1; S1=T;
00567 C2= -C2;
00568 T=C3; C3= -S3; S3= -T;
00569 L200:
00570 REPORT
00571 for (K=K0; K<N; K+=M4)
00572 {
00573 RS0=X0[K]+X2[K]; IS0=Y0[K]+Y2[K];
00574 RU0=X0[K]-X2[K]; IU0=Y0[K]-Y2[K];
00575 RS1=X1[K]+X3[K]; IS1=Y1[K]+Y3[K];
00576 RU1=X1[K]-X3[K]; IU1=Y1[K]-Y3[K];
00577 X0[K]=RS0+RS1; Y0[K]=IS0+IS1;
00578 if (!ZERO)
00579 {
00580 REPORT
00581 R1=RU0+IU1; I1=IU0-RU1;
00582 R2=RS0-RS1; I2=IS0-IS1;
00583 R3=RU0-IU1; I3=IU0+RU1;
00584 X2[K]=R1*C1+I1*S1; Y2[K]=I1*C1-R1*S1;
00585 X1[K]=R2*C2+I2*S2; Y1[K]=I2*C2-R2*S2;
00586 X3[K]=R3*C3+I3*S3; Y3[K]=I3*C3-R3*S3;
00587 }
00588 else
00589 {
00590 REPORT
00591 X2[K]=RU0+IU1; Y2[K]=IU0-RU1;
00592 X1[K]=RS0-RS1; Y1[K]=IS0-IS1;
00593 X3[K]=RU0-IU1; Y3[K]=IU0+RU1;
00594 }
00595 }
00596 goto L100;
00597 L600: ;
00598 }
00599
00600 return;
00601 }
00602
00603 static void R_5_FTK (int N, int M,
00604 Real* X0, Real* Y0, Real* X1, Real* Y1, Real* X2, Real* Y2,
00605 Real* X3, Real* Y3, Real* X4, Real* Y4)
00606
00607
00608 {
00609 REPORT
00610 bool NO_FOLD,ZERO;
00611 int J,K,K0,M5,M_OVER_2;
00612 Real ANGLE,A1,A2,B1,B2,C1,C2,C3,C4,S1,S2,S3,S4,T,TWOPI;
00613 Real R0,R1,R2,R3,R4,RA1,RA2,RB1,RB2,RS1,RS2,RU1,RU2;
00614 Real I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2;
00615
00616 M5=M*5; M_OVER_2=M/2+1;
00617 TWOPI=8.0*atan(1.0);
00618 A1=cos(TWOPI/5.0); B1=sin(TWOPI/5.0);
00619 A2=cos(2.0*TWOPI/5.0); B2=sin(2.0*TWOPI/5.0);
00620
00621 for (J=0; J<M_OVER_2; J++)
00622 {
00623 NO_FOLD = (J==0 || 2*J==M);
00624 K0=J;
00625 ANGLE=TWOPI*Real(J)/Real(M5); ZERO=ANGLE==0.0;
00626 C1=cos(ANGLE); S1=sin(ANGLE);
00627 C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
00628 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00629 C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
00630 goto L200;
00631 L100:
00632 REPORT
00633 if (NO_FOLD) { REPORT goto L600; }
00634 REPORT
00635 NO_FOLD=true; K0=M-J;
00636 T=C1*A1+S1*B1; S1=C1*B1-S1*A1; C1=T;
00637 T=C2*A2+S2*B2; S2=C2*B2-S2*A2; C2=T;
00638 T=C3*A2-S3*B2; S3= -C3*B2-S3*A2; C3=T;
00639 T=C4*A1-S4*B1; S4= -C4*B1-S4*A1; C4=T;
00640 L200:
00641 REPORT
00642 for (K=K0; K<N; K+=M5)
00643 {
00644 R0=X0[K]; I0=Y0[K];
00645 RS1=X1[K]+X4[K]; IS1=Y1[K]+Y4[K];
00646 RU1=X1[K]-X4[K]; IU1=Y1[K]-Y4[K];
00647 RS2=X2[K]+X3[K]; IS2=Y2[K]+Y3[K];
00648 RU2=X2[K]-X3[K]; IU2=Y2[K]-Y3[K];
00649 X0[K]=R0+RS1+RS2; Y0[K]=I0+IS1+IS2;
00650 RA1=R0+RS1*A1+RS2*A2; IA1=I0+IS1*A1+IS2*A2;
00651 RA2=R0+RS1*A2+RS2*A1; IA2=I0+IS1*A2+IS2*A1;
00652 RB1=RU1*B1+RU2*B2; IB1=IU1*B1+IU2*B2;
00653 RB2=RU1*B2-RU2*B1; IB2=IU1*B2-IU2*B1;
00654 if (!ZERO)
00655 {
00656 REPORT
00657 R1=RA1+IB1; I1=IA1-RB1;
00658 R2=RA2+IB2; I2=IA2-RB2;
00659 R3=RA2-IB2; I3=IA2+RB2;
00660 R4=RA1-IB1; I4=IA1+RB1;
00661 X1[K]=R1*C1+I1*S1; Y1[K]=I1*C1-R1*S1;
00662 X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
00663 X3[K]=R3*C3+I3*S3; Y3[K]=I3*C3-R3*S3;
00664 X4[K]=R4*C4+I4*S4; Y4[K]=I4*C4-R4*S4;
00665 }
00666 else
00667 {
00668 REPORT
00669 X1[K]=RA1+IB1; Y1[K]=IA1-RB1;
00670 X2[K]=RA2+IB2; Y2[K]=IA2-RB2;
00671 X3[K]=RA2-IB2; Y3[K]=IA2+RB2;
00672 X4[K]=RA1-IB1; Y4[K]=IA1+RB1;
00673 }
00674 }
00675 goto L100;
00676 L600: ;
00677 }
00678
00679 return;
00680 }
00681
00682 static void R_8_FTK (int N, int M,
00683 Real* X0, Real* Y0, Real* X1, Real* Y1,
00684 Real* X2, Real* Y2, Real* X3, Real* Y3,
00685 Real* X4, Real* Y4, Real* X5, Real* Y5,
00686 Real* X6, Real* Y6, Real* X7, Real* Y7)
00687
00688 {
00689 REPORT
00690 bool NO_FOLD,ZERO;
00691 int J,K,K0,M8,M_OVER_2;
00692 Real ANGLE,C1,C2,C3,C4,C5,C6,C7,E,S1,S2,S3,S4,S5,S6,S7,T,TWOPI;
00693 Real R1,R2,R3,R4,R5,R6,R7,RS0,RS1,RS2,RS3,RU0,RU1,RU2,RU3;
00694 Real I1,I2,I3,I4,I5,I6,I7,IS0,IS1,IS2,IS3,IU0,IU1,IU2,IU3;
00695 Real RSS0,RSS1,RSU0,RSU1,RUS0,RUS1,RUU0,RUU1;
00696 Real ISS0,ISS1,ISU0,ISU1,IUS0,IUS1,IUU0,IUU1;
00697
00698 M8=M*8; M_OVER_2=M/2+1;
00699 TWOPI=8.0*atan(1.0); E=cos(TWOPI/8.0);
00700
00701 for (J=0;J<M_OVER_2;J++)
00702 {
00703 NO_FOLD= (J==0 || 2*J==M);
00704 K0=J;
00705 ANGLE=TWOPI*Real(J)/Real(M8); ZERO=ANGLE==0.0;
00706 C1=cos(ANGLE); S1=sin(ANGLE);
00707 C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
00708 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00709 C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
00710 C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
00711 C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
00712 C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
00713 goto L200;
00714 L100:
00715 REPORT
00716 if (NO_FOLD) { REPORT goto L600; }
00717 REPORT
00718 NO_FOLD=true; K0=M-J;
00719 T=(C1+S1)*E; S1=(C1-S1)*E; C1=T;
00720 T=S2; S2=C2; C2=T;
00721 T=(-C3+S3)*E; S3=(C3+S3)*E; C3=T;
00722 C4= -C4;
00723 T= -(C5+S5)*E; S5=(-C5+S5)*E; C5=T;
00724 T= -S6; S6= -C6; C6=T;
00725 T=(C7-S7)*E; S7= -(C7+S7)*E; C7=T;
00726 L200:
00727 REPORT
00728 for (K=K0; K<N; K+=M8)
00729 {
00730 RS0=X0[K]+X4[K]; IS0=Y0[K]+Y4[K];
00731 RU0=X0[K]-X4[K]; IU0=Y0[K]-Y4[K];
00732 RS1=X1[K]+X5[K]; IS1=Y1[K]+Y5[K];
00733 RU1=X1[K]-X5[K]; IU1=Y1[K]-Y5[K];
00734 RS2=X2[K]+X6[K]; IS2=Y2[K]+Y6[K];
00735 RU2=X2[K]-X6[K]; IU2=Y2[K]-Y6[K];
00736 RS3=X3[K]+X7[K]; IS3=Y3[K]+Y7[K];
00737 RU3=X3[K]-X7[K]; IU3=Y3[K]-Y7[K];
00738 RSS0=RS0+RS2; ISS0=IS0+IS2;
00739 RSU0=RS0-RS2; ISU0=IS0-IS2;
00740 RSS1=RS1+RS3; ISS1=IS1+IS3;
00741 RSU1=RS1-RS3; ISU1=IS1-IS3;
00742 RUS0=RU0-IU2; IUS0=IU0+RU2;
00743 RUU0=RU0+IU2; IUU0=IU0-RU2;
00744 RUS1=RU1-IU3; IUS1=IU1+RU3;
00745 RUU1=RU1+IU3; IUU1=IU1-RU3;
00746 T=(RUS1+IUS1)*E; IUS1=(IUS1-RUS1)*E; RUS1=T;
00747 T=(RUU1+IUU1)*E; IUU1=(IUU1-RUU1)*E; RUU1=T;
00748 X0[K]=RSS0+RSS1; Y0[K]=ISS0+ISS1;
00749 if (!ZERO)
00750 {
00751 REPORT
00752 R1=RUU0+RUU1; I1=IUU0+IUU1;
00753 R2=RSU0+ISU1; I2=ISU0-RSU1;
00754 R3=RUS0+IUS1; I3=IUS0-RUS1;
00755 R4=RSS0-RSS1; I4=ISS0-ISS1;
00756 R5=RUU0-RUU1; I5=IUU0-IUU1;
00757 R6=RSU0-ISU1; I6=ISU0+RSU1;
00758 R7=RUS0-IUS1; I7=IUS0+RUS1;
00759 X4[K]=R1*C1+I1*S1; Y4[K]=I1*C1-R1*S1;
00760 X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
00761 X6[K]=R3*C3+I3*S3; Y6[K]=I3*C3-R3*S3;
00762 X1[K]=R4*C4+I4*S4; Y1[K]=I4*C4-R4*S4;
00763 X5[K]=R5*C5+I5*S5; Y5[K]=I5*C5-R5*S5;
00764 X3[K]=R6*C6+I6*S6; Y3[K]=I6*C6-R6*S6;
00765 X7[K]=R7*C7+I7*S7; Y7[K]=I7*C7-R7*S7;
00766 }
00767 else
00768 {
00769 REPORT
00770 X4[K]=RUU0+RUU1; Y4[K]=IUU0+IUU1;
00771 X2[K]=RSU0+ISU1; Y2[K]=ISU0-RSU1;
00772 X6[K]=RUS0+IUS1; Y6[K]=IUS0-RUS1;
00773 X1[K]=RSS0-RSS1; Y1[K]=ISS0-ISS1;
00774 X5[K]=RUU0-RUU1; Y5[K]=IUU0-IUU1;
00775 X3[K]=RSU0-ISU1; Y3[K]=ISU0+RSU1;
00776 X7[K]=RUS0-IUS1; Y7[K]=IUS0+RUS1;
00777 }
00778 }
00779 goto L100;
00780 L600: ;
00781 }
00782
00783 return;
00784 }
00785
00786 static void R_16_FTK (int N, int M,
00787 Real* X0, Real* Y0, Real* X1, Real* Y1,
00788 Real* X2, Real* Y2, Real* X3, Real* Y3,
00789 Real* X4, Real* Y4, Real* X5, Real* Y5,
00790 Real* X6, Real* Y6, Real* X7, Real* Y7,
00791 Real* X8, Real* Y8, Real* X9, Real* Y9,
00792 Real* X10, Real* Y10, Real* X11, Real* Y11,
00793 Real* X12, Real* Y12, Real* X13, Real* Y13,
00794 Real* X14, Real* Y14, Real* X15, Real* Y15)
00795
00796 {
00797 REPORT
00798 bool NO_FOLD,ZERO;
00799 int J,K,K0,M16,M_OVER_2;
00800 Real ANGLE,EI1,ER1,E2,EI3,ER3,EI5,ER5,T,TWOPI;
00801 Real RS0,RS1,RS2,RS3,RS4,RS5,RS6,RS7;
00802 Real IS0,IS1,IS2,IS3,IS4,IS5,IS6,IS7;
00803 Real RU0,RU1,RU2,RU3,RU4,RU5,RU6,RU7;
00804 Real IU0,IU1,IU2,IU3,IU4,IU5,IU6,IU7;
00805 Real RUS0,RUS1,RUS2,RUS3,RUU0,RUU1,RUU2,RUU3;
00806 Real ISS0,ISS1,ISS2,ISS3,ISU0,ISU1,ISU2,ISU3;
00807 Real RSS0,RSS1,RSS2,RSS3,RSU0,RSU1,RSU2,RSU3;
00808 Real IUS0,IUS1,IUS2,IUS3,IUU0,IUU1,IUU2,IUU3;
00809 Real RSSS0,RSSS1,RSSU0,RSSU1,RSUS0,RSUS1,RSUU0,RSUU1;
00810 Real ISSS0,ISSS1,ISSU0,ISSU1,ISUS0,ISUS1,ISUU0,ISUU1;
00811 Real RUSS0,RUSS1,RUSU0,RUSU1,RUUS0,RUUS1,RUUU0,RUUU1;
00812 Real IUSS0,IUSS1,IUSU0,IUSU1,IUUS0,IUUS1,IUUU0,IUUU1;
00813 Real R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15;
00814 Real I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14,I15;
00815 Real C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15;
00816 Real S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,S15;
00817
00818 M16=M*16; M_OVER_2=M/2+1;
00819 TWOPI=8.0*atan(1.0);
00820 ER1=cos(TWOPI/16.0); EI1=sin(TWOPI/16.0);
00821 E2=cos(TWOPI/8.0);
00822 ER3=cos(3.0*TWOPI/16.0); EI3=sin(3.0*TWOPI/16.0);
00823 ER5=cos(5.0*TWOPI/16.0); EI5=sin(5.0*TWOPI/16.0);
00824
00825 for (J=0; J<M_OVER_2; J++)
00826 {
00827 NO_FOLD = (J==0 || 2*J==M);
00828 K0=J;
00829 ANGLE=TWOPI*Real(J)/Real(M16);
00830 ZERO=ANGLE==0.0;
00831 C1=cos(ANGLE); S1=sin(ANGLE);
00832 C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
00833 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00834 C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
00835 C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
00836 C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
00837 C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
00838 C8=C4*C4-S4*S4; S8=C4*S4+S4*C4;
00839 C9=C8*C1-S8*S1; S9=S8*C1+C8*S1;
00840 C10=C8*C2-S8*S2; S10=S8*C2+C8*S2;
00841 C11=C8*C3-S8*S3; S11=S8*C3+C8*S3;
00842 C12=C8*C4-S8*S4; S12=S8*C4+C8*S4;
00843 C13=C8*C5-S8*S5; S13=S8*C5+C8*S5;
00844 C14=C8*C6-S8*S6; S14=S8*C6+C8*S6;
00845 C15=C8*C7-S8*S7; S15=S8*C7+C8*S7;
00846 goto L200;
00847 L100:
00848 REPORT
00849 if (NO_FOLD) { REPORT goto L600; }
00850 REPORT
00851 NO_FOLD=true; K0=M-J;
00852 T=C1*ER1+S1*EI1; S1= -S1*ER1+C1*EI1; C1=T;
00853 T=(C2+S2)*E2; S2=(C2-S2)*E2; C2=T;
00854 T=C3*ER3+S3*EI3; S3= -S3*ER3+C3*EI3; C3=T;
00855 T=S4; S4=C4; C4=T;
00856 T=S5*ER1-C5*EI1; S5=C5*ER1+S5*EI1; C5=T;
00857 T=(-C6+S6)*E2; S6=(C6+S6)*E2; C6=T;
00858 T=S7*ER3-C7*EI3; S7=C7*ER3+S7*EI3; C7=T;
00859 C8= -C8;
00860 T= -(C9*ER1+S9*EI1); S9=S9*ER1-C9*EI1; C9=T;
00861 T= -(C10+S10)*E2; S10=(-C10+S10)*E2; C10=T;
00862 T= -(C11*ER3+S11*EI3); S11=S11*ER3-C11*EI3; C11=T;
00863 T= -S12; S12= -C12; C12=T;
00864 T= -S13*ER1+C13*EI1; S13= -(C13*ER1+S13*EI1); C13=T;
00865 T=(C14-S14)*E2; S14= -(C14+S14)*E2; C14=T;
00866 T= -S15*ER3+C15*EI3; S15= -(C15*ER3+S15*EI3); C15=T;
00867 L200:
00868 REPORT
00869 for (K=K0; K<N; K+=M16)
00870 {
00871 RS0=X0[K]+X8[K]; IS0=Y0[K]+Y8[K];
00872 RU0=X0[K]-X8[K]; IU0=Y0[K]-Y8[K];
00873 RS1=X1[K]+X9[K]; IS1=Y1[K]+Y9[K];
00874 RU1=X1[K]-X9[K]; IU1=Y1[K]-Y9[K];
00875 RS2=X2[K]+X10[K]; IS2=Y2[K]+Y10[K];
00876 RU2=X2[K]-X10[K]; IU2=Y2[K]-Y10[K];
00877 RS3=X3[K]+X11[K]; IS3=Y3[K]+Y11[K];
00878 RU3=X3[K]-X11[K]; IU3=Y3[K]-Y11[K];
00879 RS4=X4[K]+X12[K]; IS4=Y4[K]+Y12[K];
00880 RU4=X4[K]-X12[K]; IU4=Y4[K]-Y12[K];
00881 RS5=X5[K]+X13[K]; IS5=Y5[K]+Y13[K];
00882 RU5=X5[K]-X13[K]; IU5=Y5[K]-Y13[K];
00883 RS6=X6[K]+X14[K]; IS6=Y6[K]+Y14[K];
00884 RU6=X6[K]-X14[K]; IU6=Y6[K]-Y14[K];
00885 RS7=X7[K]+X15[K]; IS7=Y7[K]+Y15[K];
00886 RU7=X7[K]-X15[K]; IU7=Y7[K]-Y15[K];
00887 RSS0=RS0+RS4; ISS0=IS0+IS4;
00888 RSS1=RS1+RS5; ISS1=IS1+IS5;
00889 RSS2=RS2+RS6; ISS2=IS2+IS6;
00890 RSS3=RS3+RS7; ISS3=IS3+IS7;
00891 RSU0=RS0-RS4; ISU0=IS0-IS4;
00892 RSU1=RS1-RS5; ISU1=IS1-IS5;
00893 RSU2=RS2-RS6; ISU2=IS2-IS6;
00894 RSU3=RS3-RS7; ISU3=IS3-IS7;
00895 RUS0=RU0-IU4; IUS0=IU0+RU4;
00896 RUS1=RU1-IU5; IUS1=IU1+RU5;
00897 RUS2=RU2-IU6; IUS2=IU2+RU6;
00898 RUS3=RU3-IU7; IUS3=IU3+RU7;
00899 RUU0=RU0+IU4; IUU0=IU0-RU4;
00900 RUU1=RU1+IU5; IUU1=IU1-RU5;
00901 RUU2=RU2+IU6; IUU2=IU2-RU6;
00902 RUU3=RU3+IU7; IUU3=IU3-RU7;
00903 T=(RSU1+ISU1)*E2; ISU1=(ISU1-RSU1)*E2; RSU1=T;
00904 T=(RSU3+ISU3)*E2; ISU3=(ISU3-RSU3)*E2; RSU3=T;
00905 T=RUS1*ER3+IUS1*EI3; IUS1=IUS1*ER3-RUS1*EI3; RUS1=T;
00906 T=(RUS2+IUS2)*E2; IUS2=(IUS2-RUS2)*E2; RUS2=T;
00907 T=RUS3*ER5+IUS3*EI5; IUS3=IUS3*ER5-RUS3*EI5; RUS3=T;
00908 T=RUU1*ER1+IUU1*EI1; IUU1=IUU1*ER1-RUU1*EI1; RUU1=T;
00909 T=(RUU2+IUU2)*E2; IUU2=(IUU2-RUU2)*E2; RUU2=T;
00910 T=RUU3*ER3+IUU3*EI3; IUU3=IUU3*ER3-RUU3*EI3; RUU3=T;
00911 RSSS0=RSS0+RSS2; ISSS0=ISS0+ISS2;
00912 RSSS1=RSS1+RSS3; ISSS1=ISS1+ISS3;
00913 RSSU0=RSS0-RSS2; ISSU0=ISS0-ISS2;
00914 RSSU1=RSS1-RSS3; ISSU1=ISS1-ISS3;
00915 RSUS0=RSU0-ISU2; ISUS0=ISU0+RSU2;
00916 RSUS1=RSU1-ISU3; ISUS1=ISU1+RSU3;
00917 RSUU0=RSU0+ISU2; ISUU0=ISU0-RSU2;
00918 RSUU1=RSU1+ISU3; ISUU1=ISU1-RSU3;
00919 RUSS0=RUS0-IUS2; IUSS0=IUS0+RUS2;
00920 RUSS1=RUS1-IUS3; IUSS1=IUS1+RUS3;
00921 RUSU0=RUS0+IUS2; IUSU0=IUS0-RUS2;
00922 RUSU1=RUS1+IUS3; IUSU1=IUS1-RUS3;
00923 RUUS0=RUU0+RUU2; IUUS0=IUU0+IUU2;
00924 RUUS1=RUU1+RUU3; IUUS1=IUU1+IUU3;
00925 RUUU0=RUU0-RUU2; IUUU0=IUU0-IUU2;
00926 RUUU1=RUU1-RUU3; IUUU1=IUU1-IUU3;
00927 X0[K]=RSSS0+RSSS1; Y0[K]=ISSS0+ISSS1;
00928 if (!ZERO)
00929 {
00930 REPORT
00931 R1=RUUS0+RUUS1; I1=IUUS0+IUUS1;
00932 R2=RSUU0+RSUU1; I2=ISUU0+ISUU1;
00933 R3=RUSU0+RUSU1; I3=IUSU0+IUSU1;
00934 R4=RSSU0+ISSU1; I4=ISSU0-RSSU1;
00935 R5=RUUU0+IUUU1; I5=IUUU0-RUUU1;
00936 R6=RSUS0+ISUS1; I6=ISUS0-RSUS1;
00937 R7=RUSS0+IUSS1; I7=IUSS0-RUSS1;
00938 R8=RSSS0-RSSS1; I8=ISSS0-ISSS1;
00939 R9=RUUS0-RUUS1; I9=IUUS0-IUUS1;
00940 R10=RSUU0-RSUU1; I10=ISUU0-ISUU1;
00941 R11=RUSU0-RUSU1; I11=IUSU0-IUSU1;
00942 R12=RSSU0-ISSU1; I12=ISSU0+RSSU1;
00943 R13=RUUU0-IUUU1; I13=IUUU0+RUUU1;
00944 R14=RSUS0-ISUS1; I14=ISUS0+RSUS1;
00945 R15=RUSS0-IUSS1; I15=IUSS0+RUSS1;
00946 X8[K]=R1*C1+I1*S1; Y8[K]=I1*C1-R1*S1;
00947 X4[K]=R2*C2+I2*S2; Y4[K]=I2*C2-R2*S2;
00948 X12[K]=R3*C3+I3*S3; Y12[K]=I3*C3-R3*S3;
00949 X2[K]=R4*C4+I4*S4; Y2[K]=I4*C4-R4*S4;
00950 X10[K]=R5*C5+I5*S5; Y10[K]=I5*C5-R5*S5;
00951 X6[K]=R6*C6+I6*S6; Y6[K]=I6*C6-R6*S6;
00952 X14[K]=R7*C7+I7*S7; Y14[K]=I7*C7-R7*S7;
00953 X1[K]=R8*C8+I8*S8; Y1[K]=I8*C8-R8*S8;
00954 X9[K]=R9*C9+I9*S9; Y9[K]=I9*C9-R9*S9;
00955 X5[K]=R10*C10+I10*S10; Y5[K]=I10*C10-R10*S10;
00956 X13[K]=R11*C11+I11*S11; Y13[K]=I11*C11-R11*S11;
00957 X3[K]=R12*C12+I12*S12; Y3[K]=I12*C12-R12*S12;
00958 X11[K]=R13*C13+I13*S13; Y11[K]=I13*C13-R13*S13;
00959 X7[K]=R14*C14+I14*S14; Y7[K]=I14*C14-R14*S14;
00960 X15[K]=R15*C15+I15*S15; Y15[K]=I15*C15-R15*S15;
00961 }
00962 else
00963 {
00964 REPORT
00965 X8[K]=RUUS0+RUUS1; Y8[K]=IUUS0+IUUS1;
00966 X4[K]=RSUU0+RSUU1; Y4[K]=ISUU0+ISUU1;
00967 X12[K]=RUSU0+RUSU1; Y12[K]=IUSU0+IUSU1;
00968 X2[K]=RSSU0+ISSU1; Y2[K]=ISSU0-RSSU1;
00969 X10[K]=RUUU0+IUUU1; Y10[K]=IUUU0-RUUU1;
00970 X6[K]=RSUS0+ISUS1; Y6[K]=ISUS0-RSUS1;
00971 X14[K]=RUSS0+IUSS1; Y14[K]=IUSS0-RUSS1;
00972 X1[K]=RSSS0-RSSS1; Y1[K]=ISSS0-ISSS1;
00973 X9[K]=RUUS0-RUUS1; Y9[K]=IUUS0-IUUS1;
00974 X5[K]=RSUU0-RSUU1; Y5[K]=ISUU0-ISUU1;
00975 X13[K]=RUSU0-RUSU1; Y13[K]=IUSU0-IUSU1;
00976 X3[K]=RSSU0-ISSU1; Y3[K]=ISSU0+RSSU1;
00977 X11[K]=RUUU0-IUUU1; Y11[K]=IUUU0+RUUU1;
00978 X7[K]=RSUS0-ISUS1; Y7[K]=ISUS0+RSUS1;
00979 X15[K]=RUSS0-IUSS1; Y15[K]=IUSS0+RUSS1;
00980 }
00981 }
00982 goto L100;
00983 L600: ;
00984 }
00985
00986 return;
00987 }
00988
00989
00990
00991
00992 bool FFT_Controller::CanFactor(int PTS)
00993 {
00994 REPORT
00995 const int NP = 16, NQ = 10, PMAX=19;
00996
00997 if (PTS<=1) { REPORT return true; }
00998
00999 int N = PTS, F = 2, P = 0, Q = 0;
01000
01001 while (N > 1)
01002 {
01003 bool fail = true;
01004 for (int J = F; J <= PMAX; J++)
01005 if (N % J == 0) { fail = false; F=J; break; }
01006 if (fail || P >= NP || Q >= NQ) { REPORT return false; }
01007 N /= F;
01008 if (N % F != 0) Q++; else { N /= F; P++; }
01009 }
01010
01011 return true;
01012
01013 }
01014
01015 bool FFT_Controller::OnlyOldFFT;
01016
01017
01018
01019 MultiRadixCounter::MultiRadixCounter(int nx, const SimpleIntArray& rx,
01020 SimpleIntArray& vx)
01021 : Radix(rx), Value(vx), n(nx), reverse(0),
01022 product(1), counter(0), finish(false)
01023 {
01024 REPORT for (int k = 0; k < n; k++) { Value[k] = 0; product *= Radix[k]; }
01025 }
01026
01027 void MultiRadixCounter::operator++()
01028 {
01029 REPORT
01030 counter++; int p = product;
01031 for (int k = 0; k < n; k++)
01032 {
01033 Value[k]++; int p1 = p / Radix[k]; reverse += p1;
01034 if (Value[k] == Radix[k]) { REPORT Value[k] = 0; reverse -= p; p = p1; }
01035 else { REPORT return; }
01036 }
01037 finish = true;
01038 }
01039
01040
01041 static int BitReverse(int x, int prod, int n, const SimpleIntArray& f)
01042 {
01043
01044
01045
01046
01047
01048 REPORT
01049 const int* d = f.Data() + n; int sum = 0; int q = 1;
01050 while (n--)
01051 {
01052 prod /= *(--d);
01053 int c = x / prod; x-= c * prod;
01054 sum += q * c; q *= *d;
01055 }
01056 return sum;
01057 }
01058
01059
01060 #ifdef use_namespace
01061 }
01062 #endif
01063
01065