| 777 | #define CH2(a,b) ch[(a)+idl1*(b)] |
| 778 | |
| 779 | NOINLINE static int passg (size_t ido, size_t ip, size_t l1, |
| 780 | cmplx * restrict cc, cmplx * restrict ch, const cmplx * restrict wa, |
| 781 | const cmplx * restrict csarr, const int sign) |
| 782 | { |
| 783 | const size_t cdim=ip; |
| 784 | size_t ipph = (ip+1)/2; |
| 785 | size_t idl1 = ido*l1; |
| 786 | |
| 787 | cmplx * restrict wal=RALLOC(cmplx,ip); |
| 788 | if (!wal) return -1; |
| 789 | wal[0]=(cmplx){1.,0.}; |
| 790 | for (size_t i=1; i<ip; ++i) |
| 791 | wal[i]=(cmplx){csarr[i].r,sign*csarr[i].i}; |
| 792 | |
| 793 | for (size_t k=0; k<l1; ++k) |
| 794 | for (size_t i=0; i<ido; ++i) |
| 795 | CH(i,k,0) = CC(i,0,k); |
| 796 | for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc) |
| 797 | for (size_t k=0; k<l1; ++k) |
| 798 | for (size_t i=0; i<ido; ++i) |
| 799 | PMC(CH(i,k,j),CH(i,k,jc),CC(i,j,k),CC(i,jc,k)) |
| 800 | for (size_t k=0; k<l1; ++k) |
| 801 | for (size_t i=0; i<ido; ++i) |
| 802 | { |
| 803 | cmplx tmp = CH(i,k,0); |
| 804 | for (size_t j=1; j<ipph; ++j) |
| 805 | ADDC(tmp,tmp,CH(i,k,j)) |
| 806 | CX(i,k,0) = tmp; |
| 807 | } |
| 808 | for (size_t l=1, lc=ip-1; l<ipph; ++l, --lc) |
| 809 | { |
| 810 | // j=0 |
| 811 | for (size_t ik=0; ik<idl1; ++ik) |
| 812 | { |
| 813 | CX2(ik,l).r = CH2(ik,0).r+wal[l].r*CH2(ik,1).r+wal[2*l].r*CH2(ik,2).r; |
| 814 | CX2(ik,l).i = CH2(ik,0).i+wal[l].r*CH2(ik,1).i+wal[2*l].r*CH2(ik,2).i; |
| 815 | CX2(ik,lc).r=-wal[l].i*CH2(ik,ip-1).i-wal[2*l].i*CH2(ik,ip-2).i; |
| 816 | CX2(ik,lc).i=wal[l].i*CH2(ik,ip-1).r+wal[2*l].i*CH2(ik,ip-2).r; |
| 817 | } |
| 818 | |
| 819 | size_t iwal=2*l; |
| 820 | size_t j=3, jc=ip-3; |
| 821 | for (; j<ipph-1; j+=2, jc-=2) |
| 822 | { |
| 823 | iwal+=l; if (iwal>ip) iwal-=ip; |
| 824 | cmplx xwal=wal[iwal]; |
| 825 | iwal+=l; if (iwal>ip) iwal-=ip; |
| 826 | cmplx xwal2=wal[iwal]; |
| 827 | for (size_t ik=0; ik<idl1; ++ik) |
| 828 | { |
| 829 | CX2(ik,l).r += CH2(ik,j).r*xwal.r+CH2(ik,j+1).r*xwal2.r; |
| 830 | CX2(ik,l).i += CH2(ik,j).i*xwal.r+CH2(ik,j+1).i*xwal2.r; |
| 831 | CX2(ik,lc).r -= CH2(ik,jc).i*xwal.i+CH2(ik,jc-1).i*xwal2.i; |
| 832 | CX2(ik,lc).i += CH2(ik,jc).r*xwal.i+CH2(ik,jc-1).r*xwal2.i; |
| 833 | } |
| 834 | } |
| 835 | for (; j<ipph; ++j, --jc) |
| 836 | { |