Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/kernel/trig.c @ 84:08ae793730bd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add null config files
author Chris Cannam
date Mon, 02 Mar 2020 14:03:47 +0000
parents 2cd0e3b3e1fd
children
rev   line source
Chris@42 1 /*
Chris@42 2 * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@42 3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@42 4 *
Chris@42 5 * This program is free software; you can redistribute it and/or modify
Chris@42 6 * it under the terms of the GNU General Public License as published by
Chris@42 7 * the Free Software Foundation; either version 2 of the License, or
Chris@42 8 * (at your option) any later version.
Chris@42 9 *
Chris@42 10 * This program is distributed in the hope that it will be useful,
Chris@42 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@42 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Chris@42 13 * GNU General Public License for more details.
Chris@42 14 *
Chris@42 15 * You should have received a copy of the GNU General Public License
Chris@42 16 * along with this program; if not, write to the Free Software
Chris@42 17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Chris@42 18 *
Chris@42 19 */
Chris@42 20
Chris@42 21
Chris@42 22 /* trigonometric functions */
Chris@42 23 #include "ifftw.h"
Chris@42 24 #include <math.h>
Chris@42 25
Chris@42 26 #if defined(TRIGREAL_IS_LONG_DOUBLE)
Chris@42 27 # define COS cosl
Chris@42 28 # define SIN sinl
Chris@42 29 # define KTRIG(x) (x##L)
Chris@42 30 # if defined(HAVE_DECL_SINL) && !HAVE_DECL_SINL
Chris@42 31 extern long double sinl(long double x);
Chris@42 32 # endif
Chris@42 33 # if defined(HAVE_DECL_COSL) && !HAVE_DECL_COSL
Chris@42 34 extern long double cosl(long double x);
Chris@42 35 # endif
Chris@42 36 #elif defined(TRIGREAL_IS_QUAD)
Chris@42 37 # define COS cosq
Chris@42 38 # define SIN sinq
Chris@42 39 # define KTRIG(x) (x##Q)
Chris@42 40 extern __float128 sinq(__float128 x);
Chris@42 41 extern __float128 cosq(__float128 x);
Chris@42 42 #else
Chris@42 43 # define COS cos
Chris@42 44 # define SIN sin
Chris@42 45 # define KTRIG(x) (x)
Chris@42 46 #endif
Chris@42 47
Chris@42 48 static const trigreal K2PI =
Chris@42 49 KTRIG(6.2831853071795864769252867665590057683943388);
Chris@42 50 #define by2pi(m, n) ((K2PI * (m)) / (n))
Chris@42 51
Chris@42 52 /*
Chris@42 53 * Improve accuracy by reducing x to range [0..1/8]
Chris@42 54 * before multiplication by 2 * PI.
Chris@42 55 */
Chris@42 56
Chris@42 57 static void real_cexp(INT m, INT n, trigreal *out)
Chris@42 58 {
Chris@42 59 trigreal theta, c, s, t;
Chris@42 60 unsigned octant = 0;
Chris@42 61 INT quarter_n = n;
Chris@42 62
Chris@42 63 n += n; n += n;
Chris@42 64 m += m; m += m;
Chris@42 65
Chris@42 66 if (m < 0) m += n;
Chris@42 67 if (m > n - m) { m = n - m; octant |= 4; }
Chris@42 68 if (m - quarter_n > 0) { m = m - quarter_n; octant |= 2; }
Chris@42 69 if (m > quarter_n - m) { m = quarter_n - m; octant |= 1; }
Chris@42 70
Chris@42 71 theta = by2pi(m, n);
Chris@42 72 c = COS(theta); s = SIN(theta);
Chris@42 73
Chris@42 74 if (octant & 1) { t = c; c = s; s = t; }
Chris@42 75 if (octant & 2) { t = c; c = -s; s = t; }
Chris@42 76 if (octant & 4) { s = -s; }
Chris@42 77
Chris@42 78 out[0] = c;
Chris@42 79 out[1] = s;
Chris@42 80 }
Chris@42 81
Chris@42 82 static INT choose_twshft(INT n)
Chris@42 83 {
Chris@42 84 INT log2r = 0;
Chris@42 85 while (n > 0) {
Chris@42 86 ++log2r;
Chris@42 87 n /= 4;
Chris@42 88 }
Chris@42 89 return log2r;
Chris@42 90 }
Chris@42 91
Chris@42 92 static void cexpl_sqrtn_table(triggen *p, INT m, trigreal *res)
Chris@42 93 {
Chris@42 94 m += p->n * (m < 0);
Chris@42 95
Chris@42 96 {
Chris@42 97 INT m0 = m & p->twmsk;
Chris@42 98 INT m1 = m >> p->twshft;
Chris@42 99 trigreal wr0 = p->W0[2 * m0];
Chris@42 100 trigreal wi0 = p->W0[2 * m0 + 1];
Chris@42 101 trigreal wr1 = p->W1[2 * m1];
Chris@42 102 trigreal wi1 = p->W1[2 * m1 + 1];
Chris@42 103
Chris@42 104 res[0] = wr1 * wr0 - wi1 * wi0;
Chris@42 105 res[1] = wi1 * wr0 + wr1 * wi0;
Chris@42 106 }
Chris@42 107 }
Chris@42 108
Chris@42 109 /* multiply (xr, xi) by exp(FFT_SIGN * 2*pi*i*m/n) */
Chris@42 110 static void rotate_sqrtn_table(triggen *p, INT m, R xr, R xi, R *res)
Chris@42 111 {
Chris@42 112 m += p->n * (m < 0);
Chris@42 113
Chris@42 114 {
Chris@42 115 INT m0 = m & p->twmsk;
Chris@42 116 INT m1 = m >> p->twshft;
Chris@42 117 trigreal wr0 = p->W0[2 * m0];
Chris@42 118 trigreal wi0 = p->W0[2 * m0 + 1];
Chris@42 119 trigreal wr1 = p->W1[2 * m1];
Chris@42 120 trigreal wi1 = p->W1[2 * m1 + 1];
Chris@42 121 trigreal wr = wr1 * wr0 - wi1 * wi0;
Chris@42 122 trigreal wi = wi1 * wr0 + wr1 * wi0;
Chris@42 123
Chris@42 124 #if FFT_SIGN == -1
Chris@42 125 res[0] = xr * wr + xi * wi;
Chris@42 126 res[1] = xi * wr - xr * wi;
Chris@42 127 #else
Chris@42 128 res[0] = xr * wr - xi * wi;
Chris@42 129 res[1] = xi * wr + xr * wi;
Chris@42 130 #endif
Chris@42 131 }
Chris@42 132 }
Chris@42 133
Chris@42 134 static void cexpl_sincos(triggen *p, INT m, trigreal *res)
Chris@42 135 {
Chris@42 136 real_cexp(m, p->n, res);
Chris@42 137 }
Chris@42 138
Chris@42 139 static void cexp_zero(triggen *p, INT m, R *res)
Chris@42 140 {
Chris@42 141 UNUSED(p); UNUSED(m);
Chris@42 142 res[0] = 0;
Chris@42 143 res[1] = 0;
Chris@42 144 }
Chris@42 145
Chris@42 146 static void cexpl_zero(triggen *p, INT m, trigreal *res)
Chris@42 147 {
Chris@42 148 UNUSED(p); UNUSED(m);
Chris@42 149 res[0] = 0;
Chris@42 150 res[1] = 0;
Chris@42 151 }
Chris@42 152
Chris@42 153 static void cexp_generic(triggen *p, INT m, R *res)
Chris@42 154 {
Chris@42 155 trigreal resl[2];
Chris@42 156 p->cexpl(p, m, resl);
Chris@42 157 res[0] = (R)resl[0];
Chris@42 158 res[1] = (R)resl[1];
Chris@42 159 }
Chris@42 160
Chris@42 161 static void rotate_generic(triggen *p, INT m, R xr, R xi, R *res)
Chris@42 162 {
Chris@42 163 trigreal w[2];
Chris@42 164 p->cexpl(p, m, w);
Chris@42 165 res[0] = xr * w[0] - xi * (FFT_SIGN * w[1]);
Chris@42 166 res[1] = xi * w[0] + xr * (FFT_SIGN * w[1]);
Chris@42 167 }
Chris@42 168
Chris@42 169 triggen *X(mktriggen)(enum wakefulness wakefulness, INT n)
Chris@42 170 {
Chris@42 171 INT i, n0, n1;
Chris@42 172 triggen *p = (triggen *)MALLOC(sizeof(*p), TWIDDLES);
Chris@42 173
Chris@42 174 p->n = n;
Chris@42 175 p->W0 = p->W1 = 0;
Chris@42 176 p->cexp = 0;
Chris@42 177 p->rotate = 0;
Chris@42 178
Chris@42 179 switch (wakefulness) {
Chris@42 180 case SLEEPY:
Chris@42 181 A(0 /* can't happen */);
Chris@42 182 break;
Chris@42 183
Chris@42 184 case AWAKE_SQRTN_TABLE: {
Chris@42 185 INT twshft = choose_twshft(n);
Chris@42 186
Chris@42 187 p->twshft = twshft;
Chris@42 188 p->twradix = ((INT)1) << twshft;
Chris@42 189 p->twmsk = p->twradix - 1;
Chris@42 190
Chris@42 191 n0 = p->twradix;
Chris@42 192 n1 = (n + n0 - 1) / n0;
Chris@42 193
Chris@42 194 p->W0 = (trigreal *)MALLOC(n0 * 2 * sizeof(trigreal), TWIDDLES);
Chris@42 195 p->W1 = (trigreal *)MALLOC(n1 * 2 * sizeof(trigreal), TWIDDLES);
Chris@42 196
Chris@42 197 for (i = 0; i < n0; ++i)
Chris@42 198 real_cexp(i, n, p->W0 + 2 * i);
Chris@42 199
Chris@42 200 for (i = 0; i < n1; ++i)
Chris@42 201 real_cexp(i * p->twradix, n, p->W1 + 2 * i);
Chris@42 202
Chris@42 203 p->cexpl = cexpl_sqrtn_table;
Chris@42 204 p->rotate = rotate_sqrtn_table;
Chris@42 205 break;
Chris@42 206 }
Chris@42 207
Chris@42 208 case AWAKE_SINCOS:
Chris@42 209 p->cexpl = cexpl_sincos;
Chris@42 210 break;
Chris@42 211
Chris@42 212 case AWAKE_ZERO:
Chris@42 213 p->cexp = cexp_zero;
Chris@42 214 p->cexpl = cexpl_zero;
Chris@42 215 break;
Chris@42 216 }
Chris@42 217
Chris@42 218 if (!p->cexp) {
Chris@42 219 if (sizeof(trigreal) == sizeof(R))
Chris@42 220 p->cexp = (void (*)(triggen *, INT, R *))p->cexpl;
Chris@42 221 else
Chris@42 222 p->cexp = cexp_generic;
Chris@42 223 }
Chris@42 224 if (!p->rotate)
Chris@42 225 p->rotate = rotate_generic;
Chris@42 226 return p;
Chris@42 227 }
Chris@42 228
Chris@42 229 void X(triggen_destroy)(triggen *p)
Chris@42 230 {
Chris@42 231 X(ifree0)(p->W0);
Chris@42 232 X(ifree0)(p->W1);
Chris@42 233 X(ifree)(p);
Chris@42 234 }