42 static const double twopi=6.28318530717958647692;
52 if (m > n - m) { m = n - m; octant |= 4; }
53 if (m - quarter_n > 0) { m = m - quarter_n; octant |= 2; }
54 if (m > quarter_n - m) { m = quarter_n - m; octant |= 1; }
56 double theta = (twopi*m)/n;
57 *c = cos(theta); *s = sin(theta);
59 if (octant & 1) {
double t = *c; *c = *s; *s = t; }
60 if (octant & 2) {
double t = *c; *c = -(*s); *s = t; }
61 if (octant & 4) { *s = -(*s); }
70 size_t l1=(size_t)sqrt(n);
71 double scur=0.,ccur=1.;
72 for (
size_t i=1,m=0,a=1; i<n; ++i,++a)
78 scur=sin(i*alpha); ccur=cos(i*alpha);
82 s[i*stride]=sin(i*alpha);
83 c[i*stride]=cos(i*alpha);
87 c[i*stride]=ccur*c[a*stride] - scur*s[a*stride];
88 s[i*stride]=ccur*s[a*stride] + scur*c[a*stride];
93 static void fracsincos_multi_priv (
size_t n,
int num,
int den,
double *s,
94 double *c,
int stride,
int exact)
101 for (
size_t i=1; i<n; ++i)
102 fracsincos(i*num,den,&s[i*stride],&c[i*stride]);
106 size_t l1=(size_t)sqrt(n);
107 double scur=0.,ccur=1.;
108 for (
size_t i=1,m=0,a=1; i<n; ++i,++a)
115 s[i*stride]=scur; c[i*stride]=ccur;
120 fracsincos(i*num,den,&s[i*stride],&c[i*stride]);
123 c[i*stride]=ccur*c[a*stride] - scur*s[a*stride];
124 s[i*stride]=ccur*s[a*stride] + scur*c[a*stride];
133 { fracsincos_multi_priv (n,num,den,s,c,stride,0); }
135 static void sincos_2pibyn_priv (
size_t n,
size_t nang,
double *s,
double *c,
136 int stride,
int exact)
140 size_t nmax = ((n&3)==0) ? n/8+1 : ( ((n&1)==0) ? n/4+1 : n/2+1 );
141 size_t ngoal = (nang<nmax) ? nang : nmax;
142 fracsincos_multi_priv (ngoal, 1, n, s, c, stride, exact);
144 if (ndone==nang)
return;
148 if (nang<ngoal) ngoal=nang;
149 for (
size_t i=ndone; i<ngoal; ++i)
151 s[i*stride]=c[(n/4-i)*stride];
152 c[i*stride]=s[(n/4-i)*stride];
155 if (ngoal==nang)
return;
160 if (nang<ngoal) ngoal=nang;
161 for (
size_t i=ndone; i<ngoal; ++i)
163 c[i*stride]=-c[(n/2-i)*stride];
164 s[i*stride]= s[(n/2-i)*stride];
167 if (ngoal==nang)
return;
170 if (nang<ngoal) ngoal=nang;
171 for (
size_t i=ndone; i<ngoal; ++i)
173 c[i*stride]= c[(n-i)*stride];
174 s[i*stride]=-s[(n-i)*stride];
177 if (ngoal==nang)
return;
178 for (
size_t i=ndone; i<nang; ++i)
180 c[i*stride]= c[(i-n)*stride];
181 s[i*stride]= s[(i-n)*stride];
185 void sincos_2pibyn (
size_t n,
size_t nang,
double *s,
double *c,
int stride)
186 { sincos_2pibyn_priv (n,nang,s,c,stride,0); }
188 void triggen_init (
struct triggen *tg,
size_t n)
192 while ((((
size_t)1)<<(2*tg->ilg))<n) ++tg->ilg;
193 size_t s=((size_t)1)<<tg->ilg;
196 tg->t1=
RALLOC(
double,2*(s1+s));
198 fracsincos_multi_priv(s1,s,n,&(tg->t1[1]),&(tg->t1[0]),2,1);
199 sincos_2pibyn_priv(n,s,&(tg->t2[1]),&(tg->t2[0]),2,1);
201 void triggen_get (
const struct triggen *tg,
size_t i,
double *s,
double *c)
203 if (i>=tg->n) i%=tg->n;
204 size_t i1=i>>tg->ilg,
206 *c = tg->t1[2*i1]*tg->t2[2*i2 ] - tg->t1[2*i1+1]*tg->t2[2*i2+1];
207 *s = tg->t1[2*i1]*tg->t2[2*i2+1] + tg->t1[2*i1+1]*tg->t2[2*i2 ];
209 void triggen_destroy (
struct triggen *tg)
212 int trigtest(
int argc,
const char **argv);
213 int trigtest(
int argc,
const char **argv)
215 UTIL_ASSERT((argc==1)||(argv[0]==NULL),
"problem with args");
217 static const double twopi=6.28318530717958647692;
218 double *sc=
RALLOC(
double,2*LENGTH+34);
219 for (
int i=1; i<LENGTH; ++i)
221 sc[0]=sc[1]=sc[2*i+30+2]=sc[2*i+30+3]=10;
223 UTIL_ASSERT(fabs(sc[0]-10.)<1e-16,
"bad memory access");
224 UTIL_ASSERT(fabs(sc[1]-10.)<1e-16,
"bad memory access");
225 UTIL_ASSERT(fabs(sc[2*i+30+2]-10.)<1e-16,
"bad memory access");
226 UTIL_ASSERT(fabs(sc[2*i+30+3]-10.)<1e-16,
"bad memory access");
229 for (
int j=0; j<i; ++j)
233 triggen_get(&tg,j,&s2,&c2);
239 triggen_destroy(&tg);
240 sc[0]=sc[1]=sc[2*i+2]=sc[2*i+3]=10;
242 UTIL_ASSERT(fabs(sc[0]-10.)<1e-16,
"bad memory access");
243 UTIL_ASSERT(fabs(sc[1]-10.)<1e-16,
"bad memory access");
244 UTIL_ASSERT(fabs(sc[2*i+2]-10.)<1e-16,
"bad memory access");
245 UTIL_ASSERT(fabs(sc[2*i+3]-10.)<1e-16,
"bad memory access");
246 for (
int j=0; j<i; ++j)
248 double ang=twopi*1.1/i*j;
249 UTIL_ASSERT(fabs(sc[2*j+2]-sin(ang))<1e-15,
"bad sin");
250 UTIL_ASSERT(fabs(sc[2*j+3]-cos(ang))<1e-15,
"bad cos");
#define RALLOC(type, num)
void fracsincos(int m, int n, double *s, double *c)
void sincos_2pibyn(size_t n, size_t nang, double *s, double *c, int stride)
void sincos_multi(size_t n, double alpha, double *s, double *c, int stride)
void fracsincos_multi(size_t n, int num, int den, double *s, double *c, int stride)
#define UTIL_ASSERT(cond, msg)