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