cannam@95: /* cannam@95: * Copyright (c) 2003, 2007-11 Matteo Frigo cannam@95: * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology cannam@95: * cannam@95: * This program is free software; you can redistribute it and/or modify cannam@95: * it under the terms of the GNU General Public License as published by cannam@95: * the Free Software Foundation; either version 2 of the License, or cannam@95: * (at your option) any later version. cannam@95: * cannam@95: * This program is distributed in the hope that it will be useful, cannam@95: * but WITHOUT ANY WARRANTY; without even the implied warranty of cannam@95: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the cannam@95: * GNU General Public License for more details. cannam@95: * cannam@95: * You should have received a copy of the GNU General Public License cannam@95: * along with this program; if not, write to the Free Software cannam@95: * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA cannam@95: * cannam@95: */ cannam@95: cannam@95: #include "dft.h" cannam@95: cannam@95: /* cannam@95: * Compute transforms of prime sizes using Rader's trick: turn them cannam@95: * into convolutions of size n - 1, which you then perform via a pair cannam@95: * of FFTs. cannam@95: */ cannam@95: cannam@95: typedef struct { cannam@95: solver super; cannam@95: } S; cannam@95: cannam@95: typedef struct { cannam@95: plan_dft super; cannam@95: cannam@95: plan *cld1, *cld2; cannam@95: R *omega; cannam@95: INT n, g, ginv; cannam@95: INT is, os; cannam@95: plan *cld_omega; cannam@95: } P; cannam@95: cannam@95: static rader_tl *omegas = 0; cannam@95: cannam@95: static R *mkomega(enum wakefulness wakefulness, plan *p_, INT n, INT ginv) cannam@95: { cannam@95: plan_dft *p = (plan_dft *) p_; cannam@95: R *omega; cannam@95: INT i, gpower; cannam@95: trigreal scale; cannam@95: triggen *t; cannam@95: cannam@95: if ((omega = X(rader_tl_find)(n, n, ginv, omegas))) cannam@95: return omega; cannam@95: cannam@95: omega = (R *)MALLOC(sizeof(R) * (n - 1) * 2, TWIDDLES); cannam@95: cannam@95: scale = n - 1.0; /* normalization for convolution */ cannam@95: cannam@95: t = X(mktriggen)(wakefulness, n); cannam@95: for (i = 0, gpower = 1; i < n-1; ++i, gpower = MULMOD(gpower, ginv, n)) { cannam@95: trigreal w[2]; cannam@95: t->cexpl(t, gpower, w); cannam@95: omega[2*i] = w[0] / scale; cannam@95: omega[2*i+1] = FFT_SIGN * w[1] / scale; cannam@95: } cannam@95: X(triggen_destroy)(t); cannam@95: A(gpower == 1); cannam@95: cannam@95: p->apply(p_, omega, omega + 1, omega, omega + 1); cannam@95: cannam@95: X(rader_tl_insert)(n, n, ginv, omega, &omegas); cannam@95: return omega; cannam@95: } cannam@95: cannam@95: static void free_omega(R *omega) cannam@95: { cannam@95: X(rader_tl_delete)(omega, &omegas); cannam@95: } cannam@95: cannam@95: cannam@95: /***************************************************************************/ cannam@95: cannam@95: /* Below, we extensively use the identity that fft(x*)* = ifft(x) in cannam@95: order to share data between forward and backward transforms and to cannam@95: obviate the necessity of having separate forward and backward cannam@95: plans. (Although we often compute separate plans these days anyway cannam@95: due to the differing strides, etcetera.) cannam@95: cannam@95: Of course, since the new FFTW gives us separate pointers to cannam@95: the real and imaginary parts, we could have instead used the cannam@95: fft(r,i) = ifft(i,r) form of this identity, but it was easier to cannam@95: reuse the code from our old version. */ cannam@95: cannam@95: static void apply(const plan *ego_, R *ri, R *ii, R *ro, R *io) cannam@95: { cannam@95: const P *ego = (const P *) ego_; cannam@95: INT is, os; cannam@95: INT k, gpower, g, r; cannam@95: R *buf; cannam@95: R r0 = ri[0], i0 = ii[0]; cannam@95: cannam@95: r = ego->n; is = ego->is; os = ego->os; g = ego->g; cannam@95: buf = (R *) MALLOC(sizeof(R) * (r - 1) * 2, BUFFERS); cannam@95: cannam@95: /* First, permute the input, storing in buf: */ cannam@95: for (gpower = 1, k = 0; k < r - 1; ++k, gpower = MULMOD(gpower, g, r)) { cannam@95: R rA, iA; cannam@95: rA = ri[gpower * is]; cannam@95: iA = ii[gpower * is]; cannam@95: buf[2*k] = rA; buf[2*k + 1] = iA; cannam@95: } cannam@95: /* gpower == g^(r-1) mod r == 1 */; cannam@95: cannam@95: cannam@95: /* compute DFT of buf, storing in output (except DC): */ cannam@95: { cannam@95: plan_dft *cld = (plan_dft *) ego->cld1; cannam@95: cld->apply(ego->cld1, buf, buf+1, ro+os, io+os); cannam@95: } cannam@95: cannam@95: /* set output DC component: */ cannam@95: { cannam@95: ro[0] = r0 + ro[os]; cannam@95: io[0] = i0 + io[os]; cannam@95: } cannam@95: cannam@95: /* now, multiply by omega: */ cannam@95: { cannam@95: const R *omega = ego->omega; cannam@95: for (k = 0; k < r - 1; ++k) { cannam@95: E rB, iB, rW, iW; cannam@95: rW = omega[2*k]; cannam@95: iW = omega[2*k+1]; cannam@95: rB = ro[(k+1)*os]; cannam@95: iB = io[(k+1)*os]; cannam@95: ro[(k+1)*os] = rW * rB - iW * iB; cannam@95: io[(k+1)*os] = -(rW * iB + iW * rB); cannam@95: } cannam@95: } cannam@95: cannam@95: /* this will add input[0] to all of the outputs after the ifft */ cannam@95: ro[os] += r0; cannam@95: io[os] -= i0; cannam@95: cannam@95: /* inverse FFT: */ cannam@95: { cannam@95: plan_dft *cld = (plan_dft *) ego->cld2; cannam@95: cld->apply(ego->cld2, ro+os, io+os, buf, buf+1); cannam@95: } cannam@95: cannam@95: /* finally, do inverse permutation to unshuffle the output: */ cannam@95: { cannam@95: INT ginv = ego->ginv; cannam@95: gpower = 1; cannam@95: for (k = 0; k < r - 1; ++k, gpower = MULMOD(gpower, ginv, r)) { cannam@95: ro[gpower * os] = buf[2*k]; cannam@95: io[gpower * os] = -buf[2*k+1]; cannam@95: } cannam@95: A(gpower == 1); cannam@95: } cannam@95: cannam@95: cannam@95: X(ifree)(buf); cannam@95: } cannam@95: cannam@95: /***************************************************************************/ cannam@95: cannam@95: static void awake(plan *ego_, enum wakefulness wakefulness) cannam@95: { cannam@95: P *ego = (P *) ego_; cannam@95: cannam@95: X(plan_awake)(ego->cld1, wakefulness); cannam@95: X(plan_awake)(ego->cld2, wakefulness); cannam@95: X(plan_awake)(ego->cld_omega, wakefulness); cannam@95: cannam@95: switch (wakefulness) { cannam@95: case SLEEPY: cannam@95: free_omega(ego->omega); cannam@95: ego->omega = 0; cannam@95: break; cannam@95: default: cannam@95: ego->g = X(find_generator)(ego->n); cannam@95: ego->ginv = X(power_mod)(ego->g, ego->n - 2, ego->n); cannam@95: A(MULMOD(ego->g, ego->ginv, ego->n) == 1); cannam@95: cannam@95: ego->omega = mkomega(wakefulness, cannam@95: ego->cld_omega, ego->n, ego->ginv); cannam@95: break; cannam@95: } cannam@95: } cannam@95: cannam@95: static void destroy(plan *ego_) cannam@95: { cannam@95: P *ego = (P *) ego_; cannam@95: X(plan_destroy_internal)(ego->cld_omega); cannam@95: X(plan_destroy_internal)(ego->cld2); cannam@95: X(plan_destroy_internal)(ego->cld1); cannam@95: } cannam@95: cannam@95: static void print(const plan *ego_, printer *p) cannam@95: { cannam@95: const P *ego = (const P *)ego_; cannam@95: p->print(p, "(dft-rader-%D%ois=%oos=%(%p%)", cannam@95: ego->n, ego->is, ego->os, ego->cld1); cannam@95: if (ego->cld2 != ego->cld1) cannam@95: p->print(p, "%(%p%)", ego->cld2); cannam@95: if (ego->cld_omega != ego->cld1 && ego->cld_omega != ego->cld2) cannam@95: p->print(p, "%(%p%)", ego->cld_omega); cannam@95: p->putchr(p, ')'); cannam@95: } cannam@95: cannam@95: static int applicable(const solver *ego_, const problem *p_, cannam@95: const planner *plnr) cannam@95: { cannam@95: const problem_dft *p = (const problem_dft *) p_; cannam@95: UNUSED(ego_); cannam@95: return (1 cannam@95: && p->sz->rnk == 1 cannam@95: && p->vecsz->rnk == 0 cannam@95: && CIMPLIES(NO_SLOWP(plnr), p->sz->dims[0].n > RADER_MAX_SLOW) cannam@95: && X(is_prime)(p->sz->dims[0].n) cannam@95: cannam@95: /* proclaim the solver SLOW if p-1 is not easily factorizable. cannam@95: Bluestein should take care of this case. */ cannam@95: && CIMPLIES(NO_SLOWP(plnr), X(factors_into_small_primes)(p->sz->dims[0].n - 1)) cannam@95: ); cannam@95: } cannam@95: cannam@95: static int mkP(P *pln, INT n, INT is, INT os, R *ro, R *io, cannam@95: planner *plnr) cannam@95: { cannam@95: plan *cld1 = (plan *) 0; cannam@95: plan *cld2 = (plan *) 0; cannam@95: plan *cld_omega = (plan *) 0; cannam@95: R *buf = (R *) 0; cannam@95: cannam@95: /* initial allocation for the purpose of planning */ cannam@95: buf = (R *) MALLOC(sizeof(R) * (n - 1) * 2, BUFFERS); cannam@95: cannam@95: cld1 = X(mkplan_f_d)(plnr, cannam@95: X(mkproblem_dft_d)(X(mktensor_1d)(n - 1, 2, os), cannam@95: X(mktensor_1d)(1, 0, 0), cannam@95: buf, buf + 1, ro + os, io + os), cannam@95: NO_SLOW, 0, 0); cannam@95: if (!cld1) goto nada; cannam@95: cannam@95: cld2 = X(mkplan_f_d)(plnr, cannam@95: X(mkproblem_dft_d)(X(mktensor_1d)(n - 1, os, 2), cannam@95: X(mktensor_1d)(1, 0, 0), cannam@95: ro + os, io + os, buf, buf + 1), cannam@95: NO_SLOW, 0, 0); cannam@95: cannam@95: if (!cld2) goto nada; cannam@95: cannam@95: /* plan for omega array */ cannam@95: cld_omega = X(mkplan_f_d)(plnr, cannam@95: X(mkproblem_dft_d)(X(mktensor_1d)(n - 1, 2, 2), cannam@95: X(mktensor_1d)(1, 0, 0), cannam@95: buf, buf + 1, buf, buf + 1), cannam@95: NO_SLOW, ESTIMATE, 0); cannam@95: if (!cld_omega) goto nada; cannam@95: cannam@95: /* deallocate buffers; let awake() or apply() allocate them for real */ cannam@95: X(ifree)(buf); cannam@95: buf = 0; cannam@95: cannam@95: pln->cld1 = cld1; cannam@95: pln->cld2 = cld2; cannam@95: pln->cld_omega = cld_omega; cannam@95: pln->omega = 0; cannam@95: pln->n = n; cannam@95: pln->is = is; cannam@95: pln->os = os; cannam@95: cannam@95: X(ops_add)(&cld1->ops, &cld2->ops, &pln->super.super.ops); cannam@95: pln->super.super.ops.other += (n - 1) * (4 * 2 + 6) + 6; cannam@95: pln->super.super.ops.add += (n - 1) * 2 + 4; cannam@95: pln->super.super.ops.mul += (n - 1) * 4; cannam@95: cannam@95: return 1; cannam@95: cannam@95: nada: cannam@95: X(ifree0)(buf); cannam@95: X(plan_destroy_internal)(cld_omega); cannam@95: X(plan_destroy_internal)(cld2); cannam@95: X(plan_destroy_internal)(cld1); cannam@95: return 0; cannam@95: } cannam@95: cannam@95: static plan *mkplan(const solver *ego, const problem *p_, planner *plnr) cannam@95: { cannam@95: const problem_dft *p = (const problem_dft *) p_; cannam@95: P *pln; cannam@95: INT n; cannam@95: INT is, os; cannam@95: cannam@95: static const plan_adt padt = { cannam@95: X(dft_solve), awake, print, destroy cannam@95: }; cannam@95: cannam@95: if (!applicable(ego, p_, plnr)) cannam@95: return (plan *) 0; cannam@95: cannam@95: n = p->sz->dims[0].n; cannam@95: is = p->sz->dims[0].is; cannam@95: os = p->sz->dims[0].os; cannam@95: cannam@95: pln = MKPLAN_DFT(P, &padt, apply); cannam@95: if (!mkP(pln, n, is, os, p->ro, p->io, plnr)) { cannam@95: X(ifree)(pln); cannam@95: return (plan *) 0; cannam@95: } cannam@95: return &(pln->super.super); cannam@95: } cannam@95: cannam@95: static solver *mksolver(void) cannam@95: { cannam@95: static const solver_adt sadt = { PROBLEM_DFT, mkplan, 0 }; cannam@95: S *slv = MKSOLVER(S, &sadt); cannam@95: return &(slv->super); cannam@95: } cannam@95: cannam@95: void X(dft_rader_register)(planner *p) cannam@95: { cannam@95: REGISTER_SOLVER(p, mksolver()); cannam@95: }