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