110 #define REPORT { static ExeCounter ExeCount(__LINE__,20); ++ExeCount; } 118 static void GR_1D_FS (
int PTS,
int N_SYM,
int N_UN_SYM,
127 static void R_4_FTK (
int N,
int M,
130 static void R_5_FTK (
int N,
int M,
133 static void R_8_FTK (
int N,
int M,
156 int F,J,N,NF,P,PMAX,P_SYM,P_TWO,Q,R,TWO_GRP;
160 const int NP = 16, NQ = 10;
170 if (PTS<=1)
return true;
171 N=PTS; P_SYM=1; F=2; P=0; Q=0;
190 for (J=F; J<=PMAX; J++)
191 if (N % J == 0) { fail =
false; F=J;
break; }
192 if (fail || P >= NP || Q >= NQ)
return false;
194 if (N % F != 0) QQ[Q++] = F;
195 else { N /= F; PP[P++] = F; P_SYM *= F; }
198 R = (Q == 0) ? 0 : 1;
206 { SYM[J]=FACTOR[J]=PP[P-1-J]; FACTOR[P+Q+J]=SYM[P+R+J]=PP[J]; }
211 for (J=0; J<Q; J++) FACTOR[P+J]=QQ[J];
217 for (J=0; J < NF; J++)
219 if (FACTOR[J]!=2)
continue;
220 P_TWO=P_TWO*2; FACTOR[J]=1;
221 if (P_TWO<TWO_GRP && FACTOR[J+1]==2)
continue;
222 FACTOR[J]=P_TWO; P_TWO=1;
230 GR_1D_FS(PTS,2*P+R,Q,SYM,P_SYM,QQ,X,Y);
236 static void GR_1D_FS (
int PTS,
int N_SYM,
int N_UN_SYM,
265 int P = MRC.Reverse();
int JJ = MRC.Counter();
Real T;
266 T=X[JJ]; X[JJ]=X[P]; X[P]=T; T=Y[JJ]; Y[JJ]=Y[P]; Y[P]=T;
279 if (N_UN_SYM==0) {
REPORT return; }
280 P_UN_SYM=PTS/
square(P_SYM); JL=(P_UN_SYM-3)*P_SYM; MS=P_UN_SYM*P_SYM;
282 for (J = P_SYM; J<=JL; J+=P_SYM)
285 do K = P_SYM *
BitReverse(K / P_SYM, P_UN_SYM, N_UN_SYM, UN_SYM);
291 for (L=0; L<P_SYM; L++)
for (M=L; M<PTS; M+=MS)
294 T=X[JJ]; X[JJ]=X[KK]; X[KK]=T; T=Y[JJ]; Y[JJ]=Y[KK]; Y[KK]=T;
311 for (
int i = 0; i < N_FACTOR; i++)
313 int P = FACTOR[i]; M /= P;
320 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;
323 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);
327 R_8_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,
328 X+3*M,Y+3*M,X+4*M,Y+4*M,X+5*M,Y+5*M,
329 X+6*M,Y+6*M,X+7*M,Y+7*M);
333 R_16_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,
334 X+3*M,Y+3*M,X+4*M,Y+4*M,X+5*M,Y+5*M,
335 X+6*M,Y+6*M,X+7*M,Y+7*M,X+8*M,Y+8*M,
336 X+9*M,Y+9*M,X+10*M,Y+10*M,X+11*M,Y+11*M,
337 X+12*M,Y+12*M,X+13*M,Y+13*M,X+14*M,Y+14*M,
352 Real ANGLE,IS,IU,RS,RU,T,TWOPI,XT,YT;
353 int J,JJ,K0,
K,M_OVER_2,MP,PM,PP,U,V;
355 Real AA [9][9], BB [9][9];
356 Real A [18], B [18], C [18], S [18];
357 Real IA [9], IB [9], RA [9], RB [9];
360 M_OVER_2=M/2+1; MP=M*P; PP=P/2; PM=P-1;
366 A[U]=cos(ANGLE); B[U]=sin(ANGLE);
367 A[JJ]=A[U]; B[JJ]= -B[U];
370 for (U=1; U<=PP; U++)
372 for (V=1; V<=PP; V++)
373 { JJ=U*V-U*V/P*P; AA[V-1][U-1]=A[JJ-1]; BB[V-1][U-1]=B[JJ-1]; }
376 for (J=0; J<M_OVER_2; J++)
378 NO_FOLD = (J==0 || 2*J==M);
380 ANGLE=TWOPI*
Real(J)/
Real(MP); ZERO=ANGLE==0.0;
381 C[0]=cos(ANGLE); S[0]=sin(ANGLE);
384 C[U]=C[U-1]*C[0]-S[U-1]*S[0];
385 S[U]=S[U-1]*C[0]+C[U-1]*S[0];
390 if (NO_FOLD) {
REPORT goto L1500; }
392 NO_FOLD=
true; K0=M-J;
394 { T=C[U]*A[U]+S[U]*B[U]; S[U]= -S[U]*A[U]+C[U]*B[U]; C[U]=T; }
397 for (K=K0; K<N; K+=MP)
400 for (U=1; U<=PP; U++)
402 RA[U-1]=XT; IA[U-1]=YT;
403 RB[U-1]=0.0; IB[U-1]=0.0;
405 for (U=1; U<=PP; U++)
408 RS=X[K+M*U]+X[K+M*JJ]; IS=Y[K+M*U]+Y[K+M*JJ];
409 RU=X[K+M*U]-X[K+M*JJ]; IU=Y[K+M*U]-Y[K+M*JJ];
413 RA[V]=RA[V]+RS*AA[V][U-1]; IA[V]=IA[V]+IS*AA[V][U-1];
414 RB[V]=RB[V]+RU*BB[V][U-1]; IB[V]=IB[V]+IU*BB[V][U-1];
418 for (U=1; U<=PP; U++)
423 XT=RA[U-1]+IB[U-1]; YT=IA[U-1]-RB[U-1];
424 X[K+M*U]=XT*C[U-1]+YT*S[U-1]; Y[K+M*U]=YT*C[U-1]-XT*S[U-1];
426 XT=RA[U-1]-IB[U-1]; YT=IA[U-1]+RB[U-1];
427 X[K+M*JJ]=XT*C[JJ-1]+YT*S[JJ-1];
428 Y[K+M*JJ]=YT*C[JJ-1]-XT*S[JJ-1];
433 X[K+M*U]=RA[U-1]+IB[U-1]; Y[K+M*U]=IA[U-1]-RB[U-1];
435 X[K+M*JJ]=RA[U-1]-IB[U-1]; Y[K+M*JJ]=IA[U-1]+RB[U-1];
450 int J,
K,K0,M2,M_OVER_2;
451 Real ANGLE,C,IS,IU,RS,RU,S,TWOPI;
453 M2=M*2; M_OVER_2=M/2+1;
456 for (J=0; J<M_OVER_2; J++)
458 NO_FOLD = (J==0 || 2*J==M);
460 ANGLE=TWOPI*
Real(J)/
Real(M2); ZERO=ANGLE==0.0;
461 C=cos(ANGLE); S=sin(ANGLE);
465 if (NO_FOLD) {
REPORT goto L600; }
467 NO_FOLD=
true; K0=M-J; C= -C;
470 for (K=K0; K<N; K+=M2)
472 RS=X0[
K]+X1[
K]; IS=Y0[
K]+Y1[
K];
473 RU=X0[
K]-X1[
K]; IU=Y0[
K]-Y1[
K];
475 if (!ZERO) { X1[
K]=RU*C+IU*S; Y1[
K]=IU*C-RU*S; }
476 else { X1[
K]=RU; Y1[
K]=IU; }
491 int J,
K,K0,M3,M_OVER_2;
492 Real ANGLE,A,B,C1,C2,S1,S2,T,TWOPI;
493 Real I0,I1,I2,IA,IB,IS,R0,R1,R2,RA,RB,RS;
495 M3=M*3; M_OVER_2=M/2+1; TWOPI=8.0*atan(1.0);
496 A=cos(TWOPI/3.0); B=sin(TWOPI/3.0);
498 for (J=0; J<M_OVER_2; J++)
500 NO_FOLD = (J==0 || 2*J==M);
502 ANGLE=TWOPI*
Real(J)/
Real(M3); ZERO=ANGLE==0.0;
503 C1=cos(ANGLE); S1=sin(ANGLE);
504 C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
508 if (NO_FOLD) {
REPORT goto L600; }
510 NO_FOLD=
true; K0=M-J;
511 T=C1*A+S1*B; S1=C1*B-S1*A; C1=T;
512 T=C2*A-S2*B; S2= -C2*B-S2*A; C2=T;
515 for (K=K0; K<N; K+=M3)
517 R0 = X0[
K]; I0 = Y0[
K];
518 RS=X1[
K]+X2[
K]; IS=Y1[
K]+Y2[
K];
519 X0[
K]=R0+RS; Y0[
K]=I0+IS;
520 RA=R0+RS*A; IA=I0+IS*A;
521 RB=(X1[
K]-X2[
K])*B; IB=(Y1[
K]-Y2[
K])*B;
525 R1=RA+IB; I1=IA-RB; R2=RA-IB; I2=IA+RB;
526 X1[
K]=R1*C1+I1*S1; Y1[
K]=I1*C1-R1*S1;
527 X2[
K]=R2*C2+I2*S2; Y2[
K]=I2*C2-R2*S2;
529 else {
REPORT X1[
K]=RA+IB; Y1[
K]=IA-RB; X2[
K]=RA-IB; Y2[
K]=IA+RB; }
545 int J,
K,K0,M4,M_OVER_2;
546 Real ANGLE,C1,C2,C3,S1,S2,S3,T,TWOPI;
547 Real I1,I2,I3,IS0,IS1,IU0,IU1,R1,R2,R3,RS0,RS1,RU0,RU1;
549 M4=M*4; M_OVER_2=M/2+1;
552 for (J=0; J<M_OVER_2; J++)
554 NO_FOLD = (J==0 || 2*J==M);
556 ANGLE=TWOPI*
Real(J)/
Real(M4); ZERO=ANGLE==0.0;
557 C1=cos(ANGLE); S1=sin(ANGLE);
558 C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
559 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
563 if (NO_FOLD) {
REPORT goto L600; }
565 NO_FOLD=
true; K0=M-J;
568 T=C3; C3= -S3; S3= -T;
571 for (K=K0; K<N; K+=M4)
573 RS0=X0[
K]+X2[
K]; IS0=Y0[
K]+Y2[
K];
574 RU0=X0[
K]-X2[
K]; IU0=Y0[
K]-Y2[
K];
575 RS1=X1[
K]+X3[
K]; IS1=Y1[
K]+Y3[
K];
576 RU1=X1[
K]-X3[
K]; IU1=Y1[
K]-Y3[
K];
577 X0[
K]=RS0+RS1; Y0[
K]=IS0+IS1;
581 R1=RU0+IU1; I1=IU0-RU1;
582 R2=RS0-RS1; I2=IS0-IS1;
583 R3=RU0-IU1; I3=IU0+RU1;
584 X2[
K]=R1*C1+I1*S1; Y2[
K]=I1*C1-R1*S1;
585 X1[
K]=R2*C2+I2*S2; Y1[
K]=I2*C2-R2*S2;
586 X3[
K]=R3*C3+I3*S3; Y3[
K]=I3*C3-R3*S3;
591 X2[
K]=RU0+IU1; Y2[
K]=IU0-RU1;
592 X1[
K]=RS0-RS1; Y1[
K]=IS0-IS1;
593 X3[
K]=RU0-IU1; Y3[
K]=IU0+RU1;
611 int J,
K,K0,M5,M_OVER_2;
612 Real ANGLE,A1,A2,B1,B2,C1,C2,C3,C4,S1,S2,S3,S4,T,TWOPI;
613 Real R0,R1,R2,R3,R4,RA1,RA2,RB1,RB2,RS1,RS2,RU1,RU2;
614 Real I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2;
616 M5=M*5; M_OVER_2=M/2+1;
618 A1=cos(TWOPI/5.0); B1=sin(TWOPI/5.0);
619 A2=cos(2.0*TWOPI/5.0); B2=sin(2.0*TWOPI/5.0);
621 for (J=0; J<M_OVER_2; J++)
623 NO_FOLD = (J==0 || 2*J==M);
625 ANGLE=TWOPI*
Real(J)/
Real(M5); ZERO=ANGLE==0.0;
626 C1=cos(ANGLE); S1=sin(ANGLE);
627 C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
628 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
629 C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
633 if (NO_FOLD) {
REPORT goto L600; }
635 NO_FOLD=
true; K0=M-J;
636 T=C1*A1+S1*B1; S1=C1*B1-S1*A1; C1=T;
637 T=C2*A2+S2*B2; S2=C2*B2-S2*A2; C2=T;
638 T=C3*A2-S3*B2; S3= -C3*B2-S3*A2; C3=T;
639 T=C4*A1-S4*B1; S4= -C4*B1-S4*A1; C4=T;
642 for (K=K0; K<N; K+=M5)
645 RS1=X1[
K]+X4[
K]; IS1=Y1[
K]+Y4[
K];
646 RU1=X1[
K]-X4[
K]; IU1=Y1[
K]-Y4[
K];
647 RS2=X2[
K]+X3[
K]; IS2=Y2[
K]+Y3[
K];
648 RU2=X2[
K]-X3[
K]; IU2=Y2[
K]-Y3[
K];
649 X0[
K]=R0+RS1+RS2; Y0[
K]=I0+IS1+IS2;
650 RA1=R0+RS1*A1+RS2*A2; IA1=I0+IS1*A1+IS2*A2;
651 RA2=R0+RS1*A2+RS2*A1; IA2=I0+IS1*A2+IS2*A1;
652 RB1=RU1*B1+RU2*B2; IB1=IU1*B1+IU2*B2;
653 RB2=RU1*B2-RU2*B1; IB2=IU1*B2-IU2*B1;
657 R1=RA1+IB1; I1=IA1-RB1;
658 R2=RA2+IB2; I2=IA2-RB2;
659 R3=RA2-IB2; I3=IA2+RB2;
660 R4=RA1-IB1; I4=IA1+RB1;
661 X1[
K]=R1*C1+I1*S1; Y1[
K]=I1*C1-R1*S1;
662 X2[
K]=R2*C2+I2*S2; Y2[
K]=I2*C2-R2*S2;
663 X3[
K]=R3*C3+I3*S3; Y3[
K]=I3*C3-R3*S3;
664 X4[
K]=R4*C4+I4*S4; Y4[
K]=I4*C4-R4*S4;
669 X1[
K]=RA1+IB1; Y1[
K]=IA1-RB1;
670 X2[
K]=RA2+IB2; Y2[
K]=IA2-RB2;
671 X3[
K]=RA2-IB2; Y3[
K]=IA2+RB2;
672 X4[
K]=RA1-IB1; Y4[
K]=IA1+RB1;
691 int J,
K,K0,M8,M_OVER_2;
692 Real ANGLE,C1,C2,C3,C4,C5,C6,C7,E,S1,S2,S3,S4,S5,S6,S7,T,TWOPI;
693 Real R1,R2,R3,R4,R5,R6,R7,RS0,RS1,RS2,RS3,RU0,RU1,RU2,RU3;
694 Real I1,I2,I3,I4,I5,I6,I7,IS0,IS1,IS2,IS3,IU0,IU1,IU2,IU3;
695 Real RSS0,RSS1,RSU0,RSU1,RUS0,RUS1,RUU0,RUU1;
696 Real ISS0,ISS1,ISU0,ISU1,IUS0,IUS1,IUU0,IUU1;
698 M8=M*8; M_OVER_2=M/2+1;
699 TWOPI=8.0*atan(1.0); E=cos(TWOPI/8.0);
701 for (J=0;J<M_OVER_2;J++)
703 NO_FOLD= (J==0 || 2*J==M);
705 ANGLE=TWOPI*
Real(J)/
Real(M8); ZERO=ANGLE==0.0;
706 C1=cos(ANGLE); S1=sin(ANGLE);
707 C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
708 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
709 C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
710 C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
711 C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
712 C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
716 if (NO_FOLD) {
REPORT goto L600; }
718 NO_FOLD=
true; K0=M-J;
719 T=(C1+S1)*E; S1=(C1-S1)*E; C1=T;
721 T=(-C3+S3)*E; S3=(C3+S3)*E; C3=T;
723 T= -(C5+S5)*E; S5=(-C5+S5)*E; C5=T;
724 T= -S6; S6= -C6; C6=T;
725 T=(C7-S7)*E; S7= -(C7+S7)*E; C7=T;
728 for (K=K0; K<N; K+=M8)
730 RS0=X0[
K]+X4[
K]; IS0=Y0[
K]+Y4[
K];
731 RU0=X0[
K]-X4[
K]; IU0=Y0[
K]-Y4[
K];
732 RS1=X1[
K]+X5[
K]; IS1=Y1[
K]+Y5[
K];
733 RU1=X1[
K]-X5[
K]; IU1=Y1[
K]-Y5[
K];
734 RS2=X2[
K]+X6[
K]; IS2=Y2[
K]+Y6[
K];
735 RU2=X2[
K]-X6[
K]; IU2=Y2[
K]-Y6[
K];
736 RS3=X3[
K]+X7[
K]; IS3=Y3[
K]+Y7[
K];
737 RU3=X3[
K]-X7[
K]; IU3=Y3[
K]-Y7[
K];
738 RSS0=RS0+RS2; ISS0=IS0+IS2;
739 RSU0=RS0-RS2; ISU0=IS0-IS2;
740 RSS1=RS1+RS3; ISS1=IS1+IS3;
741 RSU1=RS1-RS3; ISU1=IS1-IS3;
742 RUS0=RU0-IU2; IUS0=IU0+RU2;
743 RUU0=RU0+IU2; IUU0=IU0-RU2;
744 RUS1=RU1-IU3; IUS1=IU1+RU3;
745 RUU1=RU1+IU3; IUU1=IU1-RU3;
746 T=(RUS1+IUS1)*E; IUS1=(IUS1-RUS1)*E; RUS1=T;
747 T=(RUU1+IUU1)*E; IUU1=(IUU1-RUU1)*E; RUU1=T;
748 X0[
K]=RSS0+RSS1; Y0[
K]=ISS0+ISS1;
752 R1=RUU0+RUU1; I1=IUU0+IUU1;
753 R2=RSU0+ISU1; I2=ISU0-RSU1;
754 R3=RUS0+IUS1; I3=IUS0-RUS1;
755 R4=RSS0-RSS1; I4=ISS0-ISS1;
756 R5=RUU0-RUU1; I5=IUU0-IUU1;
757 R6=RSU0-ISU1; I6=ISU0+RSU1;
758 R7=RUS0-IUS1; I7=IUS0+RUS1;
759 X4[
K]=R1*C1+I1*S1; Y4[
K]=I1*C1-R1*S1;
760 X2[
K]=R2*C2+I2*S2; Y2[
K]=I2*C2-R2*S2;
761 X6[
K]=R3*C3+I3*S3; Y6[
K]=I3*C3-R3*S3;
762 X1[
K]=R4*C4+I4*S4; Y1[
K]=I4*C4-R4*S4;
763 X5[
K]=R5*C5+I5*S5; Y5[
K]=I5*C5-R5*S5;
764 X3[
K]=R6*C6+I6*S6; Y3[
K]=I6*C6-R6*S6;
765 X7[
K]=R7*C7+I7*S7; Y7[
K]=I7*C7-R7*S7;
770 X4[
K]=RUU0+RUU1; Y4[
K]=IUU0+IUU1;
771 X2[
K]=RSU0+ISU1; Y2[
K]=ISU0-RSU1;
772 X6[
K]=RUS0+IUS1; Y6[
K]=IUS0-RUS1;
773 X1[
K]=RSS0-RSS1; Y1[
K]=ISS0-ISS1;
774 X5[
K]=RUU0-RUU1; Y5[
K]=IUU0-IUU1;
775 X3[
K]=RSU0-ISU1; Y3[
K]=ISU0+RSU1;
776 X7[
K]=RUS0-IUS1; Y7[
K]=IUS0+RUS1;
799 int J,
K,K0,M16,M_OVER_2;
800 Real ANGLE,EI1,ER1,E2,EI3,ER3,EI5,ER5,T,TWOPI;
801 Real RS0,RS1,RS2,RS3,RS4,RS5,RS6,RS7;
802 Real IS0,IS1,IS2,IS3,IS4,IS5,IS6,IS7;
803 Real RU0,RU1,RU2,RU3,RU4,RU5,RU6,RU7;
804 Real IU0,IU1,IU2,IU3,IU4,IU5,IU6,IU7;
805 Real RUS0,RUS1,RUS2,RUS3,RUU0,RUU1,RUU2,RUU3;
806 Real ISS0,ISS1,ISS2,ISS3,ISU0,ISU1,ISU2,ISU3;
807 Real RSS0,RSS1,RSS2,RSS3,RSU0,RSU1,RSU2,RSU3;
808 Real IUS0,IUS1,IUS2,IUS3,IUU0,IUU1,IUU2,IUU3;
809 Real RSSS0,RSSS1,RSSU0,RSSU1,RSUS0,RSUS1,RSUU0,RSUU1;
810 Real ISSS0,ISSS1,ISSU0,ISSU1,ISUS0,ISUS1,ISUU0,ISUU1;
811 Real RUSS0,RUSS1,RUSU0,RUSU1,RUUS0,RUUS1,RUUU0,RUUU1;
812 Real IUSS0,IUSS1,IUSU0,IUSU1,IUUS0,IUUS1,IUUU0,IUUU1;
813 Real R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15;
814 Real I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14,I15;
815 Real C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15;
816 Real S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,S15;
818 M16=M*16; M_OVER_2=M/2+1;
820 ER1=cos(TWOPI/16.0); EI1=sin(TWOPI/16.0);
822 ER3=cos(3.0*TWOPI/16.0); EI3=sin(3.0*TWOPI/16.0);
823 ER5=cos(5.0*TWOPI/16.0); EI5=sin(5.0*TWOPI/16.0);
825 for (J=0; J<M_OVER_2; J++)
827 NO_FOLD = (J==0 || 2*J==M);
831 C1=cos(ANGLE); S1=sin(ANGLE);
832 C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
833 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
834 C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
835 C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
836 C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
837 C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
838 C8=C4*C4-S4*S4; S8=C4*S4+S4*C4;
839 C9=C8*C1-S8*S1; S9=S8*C1+C8*S1;
840 C10=C8*C2-S8*S2; S10=S8*C2+C8*S2;
841 C11=C8*C3-S8*S3; S11=S8*C3+C8*S3;
842 C12=C8*C4-S8*S4; S12=S8*C4+C8*S4;
843 C13=C8*C5-S8*S5; S13=S8*C5+C8*S5;
844 C14=C8*C6-S8*S6; S14=S8*C6+C8*S6;
845 C15=C8*C7-S8*S7; S15=S8*C7+C8*S7;
849 if (NO_FOLD) {
REPORT goto L600; }
851 NO_FOLD=
true; K0=M-J;
852 T=C1*ER1+S1*EI1; S1= -S1*ER1+C1*EI1; C1=T;
853 T=(C2+S2)*E2; S2=(C2-S2)*E2; C2=T;
854 T=C3*ER3+S3*EI3; S3= -S3*ER3+C3*EI3; C3=T;
856 T=S5*ER1-C5*EI1; S5=C5*ER1+S5*EI1; C5=T;
857 T=(-C6+S6)*E2; S6=(C6+S6)*E2; C6=T;
858 T=S7*ER3-C7*EI3; S7=C7*ER3+S7*EI3; C7=T;
860 T= -(C9*ER1+S9*EI1); S9=S9*ER1-C9*EI1; C9=T;
861 T= -(C10+S10)*E2; S10=(-C10+S10)*E2; C10=T;
862 T= -(C11*ER3+S11*EI3); S11=S11*ER3-C11*EI3; C11=T;
863 T= -S12; S12= -C12; C12=T;
864 T= -S13*ER1+C13*EI1; S13= -(C13*ER1+S13*EI1); C13=T;
865 T=(C14-S14)*E2; S14= -(C14+S14)*E2; C14=T;
866 T= -S15*ER3+C15*EI3; S15= -(C15*ER3+S15*EI3); C15=T;
869 for (K=K0; K<N; K+=M16)
871 RS0=X0[
K]+X8[
K]; IS0=Y0[
K]+Y8[
K];
872 RU0=X0[
K]-X8[
K]; IU0=Y0[
K]-Y8[
K];
873 RS1=X1[
K]+X9[
K]; IS1=Y1[
K]+Y9[
K];
874 RU1=X1[
K]-X9[
K]; IU1=Y1[
K]-Y9[
K];
875 RS2=X2[
K]+X10[
K]; IS2=Y2[
K]+Y10[
K];
876 RU2=X2[
K]-X10[
K]; IU2=Y2[
K]-Y10[
K];
877 RS3=X3[
K]+X11[
K]; IS3=Y3[
K]+Y11[
K];
878 RU3=X3[
K]-X11[
K]; IU3=Y3[
K]-Y11[
K];
879 RS4=X4[
K]+X12[
K]; IS4=Y4[
K]+Y12[
K];
880 RU4=X4[
K]-X12[
K]; IU4=Y4[
K]-Y12[
K];
881 RS5=X5[
K]+X13[
K]; IS5=Y5[
K]+Y13[
K];
882 RU5=X5[
K]-X13[
K]; IU5=Y5[
K]-Y13[
K];
883 RS6=X6[
K]+X14[
K]; IS6=Y6[
K]+Y14[
K];
884 RU6=X6[
K]-X14[
K]; IU6=Y6[
K]-Y14[
K];
885 RS7=X7[
K]+X15[
K]; IS7=Y7[
K]+Y15[
K];
886 RU7=X7[
K]-X15[
K]; IU7=Y7[
K]-Y15[
K];
887 RSS0=RS0+RS4; ISS0=IS0+IS4;
888 RSS1=RS1+RS5; ISS1=IS1+IS5;
889 RSS2=RS2+RS6; ISS2=IS2+IS6;
890 RSS3=RS3+RS7; ISS3=IS3+IS7;
891 RSU0=RS0-RS4; ISU0=IS0-IS4;
892 RSU1=RS1-RS5; ISU1=IS1-IS5;
893 RSU2=RS2-RS6; ISU2=IS2-IS6;
894 RSU3=RS3-RS7; ISU3=IS3-IS7;
895 RUS0=RU0-IU4; IUS0=IU0+RU4;
896 RUS1=RU1-IU5; IUS1=IU1+RU5;
897 RUS2=RU2-IU6; IUS2=IU2+RU6;
898 RUS3=RU3-IU7; IUS3=IU3+RU7;
899 RUU0=RU0+IU4; IUU0=IU0-RU4;
900 RUU1=RU1+IU5; IUU1=IU1-RU5;
901 RUU2=RU2+IU6; IUU2=IU2-RU6;
902 RUU3=RU3+IU7; IUU3=IU3-RU7;
903 T=(RSU1+ISU1)*E2; ISU1=(ISU1-RSU1)*E2; RSU1=T;
904 T=(RSU3+ISU3)*E2; ISU3=(ISU3-RSU3)*E2; RSU3=T;
905 T=RUS1*ER3+IUS1*EI3; IUS1=IUS1*ER3-RUS1*EI3; RUS1=T;
906 T=(RUS2+IUS2)*E2; IUS2=(IUS2-RUS2)*E2; RUS2=T;
907 T=RUS3*ER5+IUS3*EI5; IUS3=IUS3*ER5-RUS3*EI5; RUS3=T;
908 T=RUU1*ER1+IUU1*EI1; IUU1=IUU1*ER1-RUU1*EI1; RUU1=T;
909 T=(RUU2+IUU2)*E2; IUU2=(IUU2-RUU2)*E2; RUU2=T;
910 T=RUU3*ER3+IUU3*EI3; IUU3=IUU3*ER3-RUU3*EI3; RUU3=T;
911 RSSS0=RSS0+RSS2; ISSS0=ISS0+ISS2;
912 RSSS1=RSS1+RSS3; ISSS1=ISS1+ISS3;
913 RSSU0=RSS0-RSS2; ISSU0=ISS0-ISS2;
914 RSSU1=RSS1-RSS3; ISSU1=ISS1-ISS3;
915 RSUS0=RSU0-ISU2; ISUS0=ISU0+RSU2;
916 RSUS1=RSU1-ISU3; ISUS1=ISU1+RSU3;
917 RSUU0=RSU0+ISU2; ISUU0=ISU0-RSU2;
918 RSUU1=RSU1+ISU3; ISUU1=ISU1-RSU3;
919 RUSS0=RUS0-IUS2; IUSS0=IUS0+RUS2;
920 RUSS1=RUS1-IUS3; IUSS1=IUS1+RUS3;
921 RUSU0=RUS0+IUS2; IUSU0=IUS0-RUS2;
922 RUSU1=RUS1+IUS3; IUSU1=IUS1-RUS3;
923 RUUS0=RUU0+RUU2; IUUS0=IUU0+IUU2;
924 RUUS1=RUU1+RUU3; IUUS1=IUU1+IUU3;
925 RUUU0=RUU0-RUU2; IUUU0=IUU0-IUU2;
926 RUUU1=RUU1-RUU3; IUUU1=IUU1-IUU3;
927 X0[
K]=RSSS0+RSSS1; Y0[
K]=ISSS0+ISSS1;
931 R1=RUUS0+RUUS1; I1=IUUS0+IUUS1;
932 R2=RSUU0+RSUU1; I2=ISUU0+ISUU1;
933 R3=RUSU0+RUSU1; I3=IUSU0+IUSU1;
934 R4=RSSU0+ISSU1; I4=ISSU0-RSSU1;
935 R5=RUUU0+IUUU1; I5=IUUU0-RUUU1;
936 R6=RSUS0+ISUS1; I6=ISUS0-RSUS1;
937 R7=RUSS0+IUSS1; I7=IUSS0-RUSS1;
938 R8=RSSS0-RSSS1; I8=ISSS0-ISSS1;
939 R9=RUUS0-RUUS1; I9=IUUS0-IUUS1;
940 R10=RSUU0-RSUU1; I10=ISUU0-ISUU1;
941 R11=RUSU0-RUSU1; I11=IUSU0-IUSU1;
942 R12=RSSU0-ISSU1; I12=ISSU0+RSSU1;
943 R13=RUUU0-IUUU1; I13=IUUU0+RUUU1;
944 R14=RSUS0-ISUS1; I14=ISUS0+RSUS1;
945 R15=RUSS0-IUSS1; I15=IUSS0+RUSS1;
946 X8[
K]=R1*C1+I1*S1; Y8[
K]=I1*C1-R1*S1;
947 X4[
K]=R2*C2+I2*S2; Y4[
K]=I2*C2-R2*S2;
948 X12[
K]=R3*C3+I3*S3; Y12[
K]=I3*C3-R3*S3;
949 X2[
K]=R4*C4+I4*S4; Y2[
K]=I4*C4-R4*S4;
950 X10[
K]=R5*C5+I5*S5; Y10[
K]=I5*C5-R5*S5;
951 X6[
K]=R6*C6+I6*S6; Y6[
K]=I6*C6-R6*S6;
952 X14[
K]=R7*C7+I7*S7; Y14[
K]=I7*C7-R7*S7;
953 X1[
K]=R8*C8+I8*S8; Y1[
K]=I8*C8-R8*S8;
954 X9[
K]=R9*C9+I9*S9; Y9[
K]=I9*C9-R9*S9;
955 X5[
K]=R10*C10+I10*S10; Y5[
K]=I10*C10-R10*S10;
956 X13[
K]=R11*C11+I11*S11; Y13[
K]=I11*C11-R11*S11;
957 X3[
K]=R12*C12+I12*S12; Y3[
K]=I12*C12-R12*S12;
958 X11[
K]=R13*C13+I13*S13; Y11[
K]=I13*C13-R13*S13;
959 X7[
K]=R14*C14+I14*S14; Y7[
K]=I14*C14-R14*S14;
960 X15[
K]=R15*C15+I15*S15; Y15[
K]=I15*C15-R15*S15;
965 X8[
K]=RUUS0+RUUS1; Y8[
K]=IUUS0+IUUS1;
966 X4[
K]=RSUU0+RSUU1; Y4[
K]=ISUU0+ISUU1;
967 X12[
K]=RUSU0+RUSU1; Y12[
K]=IUSU0+IUSU1;
968 X2[
K]=RSSU0+ISSU1; Y2[
K]=ISSU0-RSSU1;
969 X10[
K]=RUUU0+IUUU1; Y10[
K]=IUUU0-RUUU1;
970 X6[
K]=RSUS0+ISUS1; Y6[
K]=ISUS0-RSUS1;
971 X14[
K]=RUSS0+IUSS1; Y14[
K]=IUSS0-RUSS1;
972 X1[
K]=RSSS0-RSSS1; Y1[
K]=ISSS0-ISSS1;
973 X9[
K]=RUUS0-RUUS1; Y9[
K]=IUUS0-IUUS1;
974 X5[
K]=RSUU0-RSUU1; Y5[
K]=ISUU0-ISUU1;
975 X13[
K]=RUSU0-RUSU1; Y13[
K]=IUSU0-IUSU1;
976 X3[
K]=RSSU0-ISSU1; Y3[
K]=ISSU0+RSSU1;
977 X11[
K]=RUUU0-IUUU1; Y11[
K]=IUUU0+RUUU1;
978 X7[
K]=RSUS0-ISUS1; Y7[
K]=ISUS0+RSUS1;
979 X15[
K]=RUSS0-IUSS1; Y15[
K]=IUSS0+RUSS1;
995 const int NP = 16, NQ = 10, PMAX=19;
997 if (PTS<=1) {
REPORT return true; }
999 int N = PTS, F = 2, P = 0, Q = 0;
1004 for (
int J = F; J <= PMAX; J++)
1005 if (N % J == 0) { fail =
false; F=J;
break; }
1006 if (fail || P >= NP || Q >= NQ) {
REPORT return false; }
1008 if (N % F != 0) Q++;
else { N /= F; P++; }
1021 : Radix(rx), Value(vx), n(nx), reverse(0),
1022 product(1), counter(0), finish(false)
1031 for (
int k = 0; k <
n; k++)
1049 const int*
d = f.
Data() +
n;
int sum = 0;
int q = 1;
1053 int c = x / prod; x-= c * prod;
1054 sum += q * c; q *= *
d;
1060 #ifdef use_namespace const SimpleIntArray & Radix
static void R_3_FTK(int N, int M, Real *X0, Real *Y0, Real *X1, Real *Y1, Real *X2, Real *Y2)
int * Data()
pointer to the data
static void R_P_FTK(int N, int M, int P, Real *X, Real *Y)
static void GR_1D_FS(int PTS, int N_SYM, int N_UN_SYM, const SimpleIntArray &SYM, int P_SYM, const SimpleIntArray &UN_SYM, Real *X, Real *Y)
MultiRadixCounter(int nx, const SimpleIntArray &rx, SimpleIntArray &vx)
static void R_2_FTK(int N, int M, Real *X0, Real *Y0, Real *X1, Real *Y1)
static void R_16_FTK(int N, int M, Real *X0, Real *Y0, Real *X1, Real *Y1, Real *X2, Real *Y2, Real *X3, Real *Y3, Real *X4, Real *Y4, Real *X5, Real *Y5, Real *X6, Real *Y6, Real *X7, Real *Y7, Real *X8, Real *Y8, Real *X9, Real *Y9, Real *X10, Real *Y10, Real *X11, Real *Y11, Real *X12, Real *Y12, Real *X13, Real *Y13, Real *X14, Real *Y14, Real *X15, Real *Y15)
static void R_5_FTK(int N, int M, Real *X0, Real *Y0, Real *X1, Real *Y1, Real *X2, Real *Y2, Real *X3, Real *Y3, Real *X4, Real *Y4)
static bool ar_1d_ft(int PTS, Real *X, Real *Y)
static int BitReverse(int x, int prod, int n, const SimpleIntArray &f)
Real sum(const BaseMatrix &B)
static void R_4_FTK(int N, int M, Real *X0, Real *Y0, Real *X1, Real *Y1, Real *X2, Real *Y2, Real *X3, Real *Y3)
static bool CanFactor(int PTS)
static void R_8_FTK(int N, int M, Real *X0, Real *Y0, Real *X1, Real *Y1, Real *X2, Real *Y2, Real *X3, Real *Y3, Real *X4, Real *Y4, Real *X5, Real *Y5, Real *X6, Real *Y6, Real *X7, Real *Y7)
static void GR_1D_FT(int N, int N_FACTOR, const SimpleIntArray &FACTOR, Real *X, Real *Y)