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 "rdft.h" cannam@127: cannam@127: /* cannam@127: * Compute DHTs of prime sizes using Rader's trick: turn them cannam@127: * into convolutions of size n - 1, which we then perform via a pair cannam@127: * of FFTs. (We can then do prime real FFTs via rdft-dht.c.) cannam@127: * cannam@127: * Optionally (determined by the "pad" field of the solver), we can cannam@127: * perform the (cyclic) convolution by zero-padding to a size cannam@127: * >= 2*(n-1) - 1. This is advantageous if n-1 has large prime factors. cannam@127: * cannam@127: */ cannam@127: cannam@127: typedef struct { cannam@127: solver super; cannam@127: int pad; cannam@127: } S; cannam@127: cannam@127: typedef struct { cannam@127: plan_rdft super; cannam@127: cannam@127: plan *cld1, *cld2; cannam@127: R *omega; cannam@127: INT n, npad, 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: /***************************************************************************/ cannam@127: cannam@127: /* If R2HC_ONLY_CONV is 1, we use a trick to perform the convolution cannam@127: purely in terms of R2HC transforms, as opposed to R2HC followed by H2RC. cannam@127: This requires a few more operations, but allows us to share the same cannam@127: plan/codelets for both Rader children. */ cannam@127: #define R2HC_ONLY_CONV 1 cannam@127: cannam@127: static void apply(const plan *ego_, R *I, R *O) cannam@127: { cannam@127: const P *ego = (const P *) ego_; cannam@127: INT n = ego->n; /* prime */ cannam@127: INT npad = ego->npad; /* == n - 1 for unpadded Rader; always even */ cannam@127: INT is = ego->is, os; cannam@127: INT k, gpower, g; cannam@127: R *buf, *omega; cannam@127: R r0; cannam@127: cannam@127: buf = (R *) MALLOC(sizeof(R) * npad, BUFFERS); cannam@127: cannam@127: /* First, permute the input, storing in buf: */ cannam@127: g = ego->g; cannam@127: for (gpower = 1, k = 0; k < n - 1; ++k, gpower = MULMOD(gpower, g, n)) { cannam@127: buf[k] = I[gpower * is]; cannam@127: } cannam@127: /* gpower == g^(n-1) mod n == 1 */; cannam@127: cannam@127: A(n - 1 <= npad); cannam@127: for (k = n - 1; k < npad; ++k) /* optionally, zero-pad convolution */ cannam@127: buf[k] = 0; cannam@127: cannam@127: os = ego->os; cannam@127: cannam@127: /* compute RDFT of buf, storing in buf (i.e., in-place): */ cannam@127: { cannam@127: plan_rdft *cld = (plan_rdft *) ego->cld1; cannam@127: cld->apply((plan *) cld, buf, buf); cannam@127: } cannam@127: cannam@127: /* set output DC component: */ cannam@127: O[0] = (r0 = I[0]) + buf[0]; cannam@127: cannam@127: /* now, multiply by omega: */ cannam@127: omega = ego->omega; cannam@127: buf[0] *= omega[0]; cannam@127: for (k = 1; k < npad/2; ++k) { cannam@127: E rB, iB, rW, iW, a, b; cannam@127: rW = omega[k]; cannam@127: iW = omega[npad - k]; cannam@127: rB = buf[k]; cannam@127: iB = buf[npad - k]; cannam@127: a = rW * rB - iW * iB; cannam@127: b = rW * iB + iW * rB; cannam@127: #if R2HC_ONLY_CONV cannam@127: buf[k] = a + b; cannam@127: buf[npad - k] = a - b; cannam@127: #else cannam@127: buf[k] = a; cannam@127: buf[npad - k] = b; cannam@127: #endif cannam@127: } cannam@127: /* Nyquist component: */ cannam@127: A(k + k == npad); /* since npad is even */ cannam@127: buf[k] *= omega[k]; cannam@127: cannam@127: /* this will add input[0] to all of the outputs after the ifft */ cannam@127: buf[0] += r0; cannam@127: cannam@127: /* inverse FFT: */ cannam@127: { cannam@127: plan_rdft *cld = (plan_rdft *) ego->cld2; cannam@127: cld->apply((plan *) cld, buf, buf); cannam@127: } cannam@127: cannam@127: /* do inverse permutation to unshuffle the output: */ cannam@127: A(gpower == 1); cannam@127: #if R2HC_ONLY_CONV cannam@127: O[os] = buf[0]; cannam@127: gpower = g = ego->ginv; cannam@127: A(npad == n - 1 || npad/2 >= n - 1); cannam@127: if (npad == n - 1) { cannam@127: for (k = 1; k < npad/2; ++k, gpower = MULMOD(gpower, g, n)) { cannam@127: O[gpower * os] = buf[k] + buf[npad - k]; cannam@127: } cannam@127: O[gpower * os] = buf[k]; cannam@127: ++k, gpower = MULMOD(gpower, g, n); cannam@127: for (; k < npad; ++k, gpower = MULMOD(gpower, g, n)) { cannam@127: O[gpower * os] = buf[npad - k] - buf[k]; cannam@127: } cannam@127: } cannam@127: else { cannam@127: for (k = 1; k < n - 1; ++k, gpower = MULMOD(gpower, g, n)) { cannam@127: O[gpower * os] = buf[k] + buf[npad - k]; cannam@127: } cannam@127: } cannam@127: #else cannam@127: g = ego->ginv; cannam@127: for (k = 0; k < n - 1; ++k, gpower = MULMOD(gpower, g, n)) { cannam@127: O[gpower * os] = buf[k]; cannam@127: } cannam@127: #endif cannam@127: A(gpower == 1); cannam@127: cannam@127: X(ifree)(buf); cannam@127: } cannam@127: cannam@127: static R *mkomega(enum wakefulness wakefulness, cannam@127: plan *p_, INT n, INT npad, INT ginv) cannam@127: { cannam@127: plan_rdft *p = (plan_rdft *) 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, npad + 1, ginv, omegas))) cannam@127: return omega; cannam@127: cannam@127: omega = (R *)MALLOC(sizeof(R) * npad, TWIDDLES); cannam@127: cannam@127: scale = npad; /* 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[i] = (w[0] + w[1]) / scale; cannam@127: } cannam@127: X(triggen_destroy)(t); cannam@127: A(gpower == 1); cannam@127: cannam@127: A(npad == n - 1 || npad >= 2*(n - 1) - 1); cannam@127: cannam@127: for (; i < npad; ++i) cannam@127: omega[i] = K(0.0); cannam@127: if (npad > n - 1) cannam@127: for (i = 1; i < n-1; ++i) cannam@127: omega[npad - i] = omega[n - 1 - i]; cannam@127: cannam@127: p->apply(p_, omega, omega); cannam@127: cannam@127: X(rader_tl_insert)(n, npad + 1, 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: 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: A(!ego->omega); cannam@127: ego->omega = mkomega(wakefulness, cannam@127: ego->cld_omega,ego->n,ego->npad,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: cannam@127: p->print(p, "(dht-rader-%D/%D%ois=%oos=%(%p%)", cannam@127: ego->n, ego->npad, 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_, const planner *plnr) cannam@127: { cannam@127: const problem_rdft *p = (const problem_rdft *) p_; cannam@127: UNUSED(ego); cannam@127: return (1 cannam@127: && p->sz->rnk == 1 cannam@127: && p->vecsz->rnk == 0 cannam@127: && p->kind[0] == DHT cannam@127: && X(is_prime)(p->sz->dims[0].n) cannam@127: && p->sz->dims[0].n > 2 cannam@127: && CIMPLIES(NO_SLOWP(plnr), p->sz->dims[0].n > RADER_MAX_SLOW) cannam@127: /* proclaim the solver SLOW if p-1 is not easily cannam@127: factorizable. Unlike in the complex case where cannam@127: Bluestein can solve the problem, in the DHT case we cannam@127: may have no other choice */ 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 choose_transform_size(INT minsz) cannam@127: { cannam@127: static const INT primes[] = { 2, 3, 5, 0 }; cannam@127: while (!X(factors_into)(minsz, primes) || minsz % 2) cannam@127: ++minsz; cannam@127: return minsz; cannam@127: } cannam@127: cannam@127: static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr) cannam@127: { cannam@127: const S *ego = (const S *) ego_; cannam@127: const problem_rdft *p = (const problem_rdft *) p_; cannam@127: P *pln; cannam@127: INT n, npad; cannam@127: INT is, os; 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: problem *cldp; cannam@127: cannam@127: static const plan_adt padt = { cannam@127: X(rdft_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: if (ego->pad) cannam@127: npad = choose_transform_size(2 * (n - 1) - 1); cannam@127: else cannam@127: npad = n - 1; cannam@127: cannam@127: /* initial allocation for the purpose of planning */ cannam@127: buf = (R *) MALLOC(sizeof(R) * npad, BUFFERS); cannam@127: cannam@127: cld1 = X(mkplan_f_d)(plnr, cannam@127: X(mkproblem_rdft_1_d)(X(mktensor_1d)(npad, 1, 1), cannam@127: X(mktensor_1d)(1, 0, 0), cannam@127: buf, buf, cannam@127: R2HC), cannam@127: NO_SLOW, 0, 0); cannam@127: if (!cld1) goto nada; cannam@127: cannam@127: cldp = cannam@127: X(mkproblem_rdft_1_d)( cannam@127: X(mktensor_1d)(npad, 1, 1), cannam@127: X(mktensor_1d)(1, 0, 0), cannam@127: buf, buf, cannam@127: #if R2HC_ONLY_CONV cannam@127: R2HC cannam@127: #else cannam@127: HC2R cannam@127: #endif cannam@127: ); cannam@127: if (!(cld2 = X(mkplan_f_d)(plnr, cldp, NO_SLOW, 0, 0))) cannam@127: goto nada; cannam@127: cannam@127: /* plan for omega */ cannam@127: cld_omega = X(mkplan_f_d)(plnr, cannam@127: X(mkproblem_rdft_1_d)( cannam@127: X(mktensor_1d)(npad, 1, 1), cannam@127: X(mktensor_1d)(1, 0, 0), cannam@127: buf, buf, R2HC), 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 = MKPLAN_RDFT(P, &padt, apply); 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->npad = npad; 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 += (npad/2-1)*6 + npad + n + (n-1) * ego->pad; cannam@127: pln->super.super.ops.add += (npad/2-1)*2 + 2 + (n-1) * ego->pad; cannam@127: pln->super.super.ops.mul += (npad/2-1)*4 + 2 + ego->pad; cannam@127: #if R2HC_ONLY_CONV cannam@127: pln->super.super.ops.other += n-2 - ego->pad; cannam@127: pln->super.super.ops.add += (npad/2-1)*2 + (n-2) - ego->pad; cannam@127: #endif cannam@127: cannam@127: return &(pln->super.super); 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: /* constructors */ cannam@127: cannam@127: static solver *mksolver(int pad) cannam@127: { cannam@127: static const solver_adt sadt = { PROBLEM_RDFT, mkplan, 0 }; cannam@127: S *slv = MKSOLVER(S, &sadt); cannam@127: slv->pad = pad; cannam@127: return &(slv->super); cannam@127: } cannam@127: cannam@127: void X(dht_rader_register)(planner *p) cannam@127: { cannam@127: REGISTER_SOLVER(p, mksolver(0)); cannam@127: REGISTER_SOLVER(p, mksolver(1)); cannam@127: }