LevelS C support library  3.50
trig_utils.c
Go to the documentation of this file.
1 /*
2  * This file is part of libc_utils.
3  *
4  * libc_utils is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * libc_utils is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with libc_utils; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18 
19 /*
20  * libc_utils is being developed at the Max-Planck-Institut fuer Astrophysik
21  * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
22  * (DLR).
23  */
24 
25 /*! \file trig_utils.c
26  *
27  * Copyright (C) 2016-2017 Max-Planck-Society
28  * \author Martin Reinecke
29  *
30  * Many inspirations for this code come from Tasche and Zeuner: "Improved
31  * Roundoff Error Analysis for Precomputed Twiddle Factors", Journal of
32  * Computational Analysis and Applications, 4, 2002.
33  */
34 
35 #include <math.h>
36 #include "c_utils.h"
37 #include "trig_utils.h"
38 
39 /* Code for accurate calculation of sin/cos(2*pi*m/n). Adapted from FFTW. */
40 void fracsincos(int m, int n, double *s, double *c)
41  {
42  static const double twopi=6.28318530717958647692;
43  UTIL_ASSERT (n>0,"denominator must be positive");
44  int quarter_n = n;
45  unsigned octant = 0;
46  m%=n;
47  if (m<0) m+=n;
48 
49  n<<=2;
50  m<<=2;
51 
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; }
55 
56  double theta = (twopi*m)/n;
57  *c = cos(theta); *s = sin(theta);
58 
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); }
62  }
63 
64 void sincos_multi (size_t n, double alpha, double *s, double *c,
65  int stride)
66  {
67  if (n==0) return;
68  s[0] = 0.; c[0]=1.;
69  if (n==1) return;
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)
73  {
74  if (a==l1)
75  {
76  a=0;
77  ++m;
78  scur=sin(i*alpha); ccur=cos(i*alpha);
79  }
80  if (m==0)
81  {
82  s[i*stride]=sin(i*alpha);
83  c[i*stride]=cos(i*alpha);
84  }
85  else
86  {
87  c[i*stride]=ccur*c[a*stride] - scur*s[a*stride];
88  s[i*stride]=ccur*s[a*stride] + scur*c[a*stride];
89  }
90  }
91  }
92 
93 static void fracsincos_multi_priv (size_t n, int num, int den, double *s,
94  double *c, int stride, int exact)
95  {
96  if (n==0) return;
97  s[0] = 0.; c[0]=1.;
98  if (n==1) return;
99  if (exact)
100  {
101  for (size_t i=1; i<n; ++i)
102  fracsincos(i*num,den,&s[i*stride],&c[i*stride]);
103  }
104  else
105  {
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)
109  {
110  if (a==l1)
111  {
112  a=0;
113  ++m;
114  fracsincos(i*num,den,&scur,&ccur);
115  s[i*stride]=scur; c[i*stride]=ccur;
116  }
117  else
118  {
119  if (m==0)
120  fracsincos(i*num,den,&s[i*stride],&c[i*stride]);
121  else
122  {
123  c[i*stride]=ccur*c[a*stride] - scur*s[a*stride];
124  s[i*stride]=ccur*s[a*stride] + scur*c[a*stride];
125  }
126  }
127  }
128  }
129  }
130 
131 void fracsincos_multi (size_t n, int num, int den, double *s, double *c,
132  int stride)
133  { fracsincos_multi_priv (n,num,den,s,c,stride,0); }
134 
135 static void sincos_2pibyn_priv (size_t n, size_t nang, double *s, double *c,
136  int stride, int exact)
137  {
138  // nmax: number of sin/cos pairs that must be genuinely computed; the rest
139  // can be obtained via symmetries
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);
143  size_t ndone=ngoal;
144  if (ndone==nang) return;
145  if ((n&3)==0)
146  {
147  ngoal=n/4+1;
148  if (nang<ngoal) ngoal=nang;
149  for (size_t i=ndone; i<ngoal; ++i)
150  {
151  s[i*stride]=c[(n/4-i)*stride];
152  c[i*stride]=s[(n/4-i)*stride];
153  }
154  ndone=ngoal;
155  if (ngoal==nang) return;
156  }
157  if ((n&1)==0)
158  {
159  ngoal=n/2+1;
160  if (nang<ngoal) ngoal=nang;
161  for (size_t i=ndone; i<ngoal; ++i)
162  {
163  c[i*stride]=-c[(n/2-i)*stride];
164  s[i*stride]= s[(n/2-i)*stride];
165  }
166  ndone=ngoal;
167  if (ngoal==nang) return;
168  }
169  ngoal=n;
170  if (nang<ngoal) ngoal=nang;
171  for (size_t i=ndone; i<ngoal; ++i)
172  {
173  c[i*stride]= c[(n-i)*stride];
174  s[i*stride]=-s[(n-i)*stride];
175  }
176  ndone=ngoal;
177  if (ngoal==nang) return;
178  for (size_t i=ndone; i<nang; ++i)
179  {
180  c[i*stride]= c[(i-n)*stride];
181  s[i*stride]= s[(i-n)*stride];
182  }
183  }
184 
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); }
187 
188 void triggen_init (struct triggen *tg, size_t n)
189  {
190  tg->n=n;
191  tg->ilg=1;
192  while ((((size_t)1)<<(2*tg->ilg))<n) ++tg->ilg;
193  size_t s=((size_t)1)<<tg->ilg;
194  tg->mask=s-1;
195  size_t s1=n/s+1;
196  tg->t1=RALLOC(double,2*(s1+s));
197  tg->t2=tg->t1+2*s1;
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);
200  }
201 void triggen_get (const struct triggen *tg,size_t i, double *s, double *c)
202  {
203  if (i>=tg->n) i%=tg->n;
204  size_t i1=i>>tg->ilg,
205  i2=i&tg->mask;
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 ];
208  }
209 void triggen_destroy (struct triggen *tg)
210  { DEALLOC(tg->t1); }
211 
212 int trigtest(int argc, const char **argv);
213 int trigtest(int argc, const char **argv)
214  {
215  UTIL_ASSERT((argc==1)||(argv[0]==NULL),"problem with args");
216 #define LENGTH 12345
217  static const double twopi=6.28318530717958647692;
218  double *sc=RALLOC(double,2*LENGTH+34);
219  for (int i=1; i<LENGTH; ++i)
220  {
221  sc[0]=sc[1]=sc[2*i+30+2]=sc[2*i+30+3]=10;
222  sincos_2pibyn(i,i+15,&sc[2],&sc[3],2);
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");
227  struct triggen tg;
228  triggen_init(&tg,i);
229  for (int j=0; j<i; ++j)
230  {
231  double c, s, c2, s2;
232  fracsincos(j,i,&s,&c);
233  triggen_get(&tg,j,&s2,&c2);
234  UTIL_ASSERT(fabs(s2-s)<4e-16,"bad sin");
235  UTIL_ASSERT(fabs(c2-c)<4e-16,"bad cos");
236  UTIL_ASSERT(fabs(sc[2*j+2]-s)<4e-16,"bad sin");
237  UTIL_ASSERT(fabs(sc[2*j+3]-c)<4e-16,"bad cos");
238  }
239  triggen_destroy(&tg);
240  sc[0]=sc[1]=sc[2*i+2]=sc[2*i+3]=10;
241  sincos_multi(i,twopi*1.1/i,&sc[2],&sc[3],2);
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)
247  {
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");
251  }
252  }
253  DEALLOC(sc);
254  return 0;
255  }
#define RALLOC(type, num)
Definition: c_utils.h:83
void fracsincos(int m, int n, double *s, double *c)
Definition: trig_utils.c:40
void sincos_2pibyn(size_t n, size_t nang, double *s, double *c, int stride)
Definition: trig_utils.c:185
void sincos_multi(size_t n, double alpha, double *s, double *c, int stride)
Definition: trig_utils.c:64
void fracsincos_multi(size_t n, int num, int den, double *s, double *c, int stride)
Definition: trig_utils.c:131
#define DEALLOC(ptr)
Definition: c_utils.h:88
#define UTIL_ASSERT(cond, msg)
Definition: c_utils.h:59

Generated on Mon Dec 10 2018 10:24:19 for LevelS C support library