| 1119 | } |
| 1120 | |
| 1121 | NOINLINE static void radf3(size_t ido, size_t l1, const double * restrict cc, |
| 1122 | double * restrict ch, const double * restrict wa) |
| 1123 | { |
| 1124 | const size_t cdim=3; |
| 1125 | static const double taur=-0.5, taui=0.86602540378443864676; |
| 1126 | |
| 1127 | for (size_t k=0; k<l1; k++) |
| 1128 | { |
| 1129 | double cr2=CC(0,k,1)+CC(0,k,2); |
| 1130 | CH(0,0,k) = CC(0,k,0)+cr2; |
| 1131 | CH(0,2,k) = taui*(CC(0,k,2)-CC(0,k,1)); |
| 1132 | CH(ido-1,1,k) = CC(0,k,0)+taur*cr2; |
| 1133 | } |
| 1134 | if (ido==1) return; |
| 1135 | for (size_t k=0; k<l1; k++) |
| 1136 | for (size_t i=2; i<ido; i+=2) |
| 1137 | { |
| 1138 | size_t ic=ido-i; |
| 1139 | double di2, di3, dr2, dr3; |
| 1140 | MULPM (dr2,di2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1)) // d2=conj(WA0)*CC1 |
| 1141 | MULPM (dr3,di3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2)) // d3=conj(WA1)*CC2 |
| 1142 | double cr2=dr2+dr3; // c add |
| 1143 | double ci2=di2+di3; |
| 1144 | CH(i-1,0,k) = CC(i-1,k,0)+cr2; // c add |
| 1145 | CH(i ,0,k) = CC(i ,k,0)+ci2; |
| 1146 | double tr2 = CC(i-1,k,0)+taur*cr2; // c add |
| 1147 | double ti2 = CC(i ,k,0)+taur*ci2; |
| 1148 | double tr3 = taui*(di2-di3); // t3 = taui*i*(d3-d2)? |
| 1149 | double ti3 = taui*(dr3-dr2); |
| 1150 | PM(CH(i-1,2,k),CH(ic-1,1,k),tr2,tr3) // PM(i) = t2+t3 |
| 1151 | PM(CH(i ,2,k),CH(ic ,1,k),ti3,ti2) // PM(ic) = conj(t2-t3) |
| 1152 | } |
| 1153 | } |
| 1154 | |
| 1155 | NOINLINE static void radf4(size_t ido, size_t l1, const double * restrict cc, |
| 1156 | double * restrict ch, const double * restrict wa) |