| 1419 | } |
| 1420 | |
| 1421 | NOINLINE static void radb3(size_t ido, size_t l1, const double * restrict cc, |
| 1422 | double * restrict ch, const double * restrict wa) |
| 1423 | { |
| 1424 | const size_t cdim=3; |
| 1425 | static const double taur=-0.5, taui=0.86602540378443864676; |
| 1426 | |
| 1427 | for (size_t k=0; k<l1; k++) |
| 1428 | { |
| 1429 | double tr2=2.*CC(ido-1,1,k); |
| 1430 | double cr2=CC(0,0,k)+taur*tr2; |
| 1431 | CH(0,k,0)=CC(0,0,k)+tr2; |
| 1432 | double ci3=2.*taui*CC(0,2,k); |
| 1433 | PM (CH(0,k,2),CH(0,k,1),cr2,ci3); |
| 1434 | } |
| 1435 | if (ido==1) return; |
| 1436 | for (size_t k=0; k<l1; k++) |
| 1437 | for (size_t i=2; i<ido; i+=2) |
| 1438 | { |
| 1439 | size_t ic=ido-i; |
| 1440 | double tr2=CC(i-1,2,k)+CC(ic-1,1,k); // t2=CC(I) + conj(CC(ic)) |
| 1441 | double ti2=CC(i ,2,k)-CC(ic ,1,k); |
| 1442 | double cr2=CC(i-1,0,k)+taur*tr2; // c2=CC +taur*t2 |
| 1443 | double ci2=CC(i ,0,k)+taur*ti2; |
| 1444 | CH(i-1,k,0)=CC(i-1,0,k)+tr2; // CH=CC+t2 |
| 1445 | CH(i ,k,0)=CC(i ,0,k)+ti2; |
| 1446 | double cr3=taui*(CC(i-1,2,k)-CC(ic-1,1,k));// c3=taui*(CC(i)-conj(CC(ic))) |
| 1447 | double ci3=taui*(CC(i ,2,k)+CC(ic ,1,k)); |
| 1448 | double di2, di3, dr2, dr3; |
| 1449 | PM(dr3,dr2,cr2,ci3) // d2= (cr2-ci3, ci2+cr3) = c2+i*c3 |
| 1450 | PM(di2,di3,ci2,cr3) // d3= (cr2+ci3, ci2-cr3) = c2-i*c3 |
| 1451 | MULPM(CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),di2,dr2) // ch = WA*d2 |
| 1452 | MULPM(CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),di3,dr3) |
| 1453 | } |
| 1454 | } |
| 1455 | |
| 1456 | NOINLINE static void radb4(size_t ido, size_t l1, const double * restrict cc, |
| 1457 | double * restrict ch, const double * restrict wa) |