Skip to content

GitLab

Commit 933310b9 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

switch to new libsharp

parent 09727e86
Expand all Hide whitespace changes
Inline Side-by-side
c_utils/trig_utils.c deleted 100644 → 0
View file @ 09727e86
/*
* This file is part of libc_utils.
*
* libc_utils is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* libc_utils is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libc_utils; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libc_utils is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*! \file trig_utils.c
*
* Copyright (C) 2016-2017 Max-Planck-Society
* \author Martin Reinecke
*
* Many inspirations for this code come from Tasche and Zeuner: "Improved
* Roundoff Error Analysis for Precomputed Twiddle Factors", Journal of
* Computational Analysis and Applications, 4, 2002.
*/
#include <math.h>
#include "c_utils.h"
#include "trig_utils.h"
/* Code for accurate calculation of sin/cos(2*pi*m/n). Adapted from FFTW. */
void fracsincos(int m, int n, double *s, double *c)
{
static const double twopi=6.28318530717958647692;
UTIL_ASSERT (n>0,"denominator must be positive");
int quarter_n = n;
unsigned octant = 0;
m%=n;
if (m<0) m+=n;
n<<=2;
m<<=2;
if (m > n - m) { m = n - m; octant |= 4; }
if (m - quarter_n > 0) { m = m - quarter_n; octant |= 2; }
if (m > quarter_n - m) { m = quarter_n - m; octant |= 1; }
double theta = (twopi*m)/n;
*c = cos(theta); *s = sin(theta);
if (octant & 1) { double t = *c; *c = *s; *s = t; }
if (octant & 2) { double t = *c; *c = -(*s); *s = t; }
if (octant & 4) { *s = -(*s); }
}
void sincos_multi (size_t n, double alpha, double *s, double *c,
int stride)
{
if (n==0) return;
s[0] = 0.; c[0]=1.;
if (n==1) return;
size_t l1=(size_t)sqrt(n);
double scur=0.,ccur=1.;
for (size_t i=1,m=0,a=1; i<n; ++i,++a)
{
if (a==l1)
{
a=0;
++m;
scur=sin(i*alpha); ccur=cos(i*alpha);
}
if (m==0)
{
s[i*stride]=sin(i*alpha);
c[i*stride]=cos(i*alpha);
}
else
{
c[i*stride]=ccur*c[a*stride] - scur*s[a*stride];
s[i*stride]=ccur*s[a*stride] + scur*c[a*stride];
}
}
}
static void fracsincos_multi_priv (size_t n, int num, int den, double *s,
double *c, int stride, int exact)
{
if (n==0) return;
s[0] = 0.; c[0]=1.;
if (n==1) return;
if (exact)
{
for (size_t i=1; i<n; ++i)
fracsincos(i*num,den,&s[i*stride],&c[i*stride]);
}
else
{
size_t l1=(size_t)sqrt(n);
double scur=0.,ccur=1.;
for (size_t i=1,m=0,a=1; i<n; ++i,++a)
{
if (a==l1)
{
a=0;
++m;
fracsincos(i*num,den,&scur,&ccur);
s[i*stride]=scur; c[i*stride]=ccur;
}
else
{
if (m==0)
fracsincos(i*num,den,&s[i*stride],&c[i*stride]);
else
{
c[i*stride]=ccur*c[a*stride] - scur*s[a*stride];
s[i*stride]=ccur*s[a*stride] + scur*c[a*stride];
}
}
}
}
}
void fracsincos_multi (size_t n, int num, int den, double *s, double *c,
int stride)
{ fracsincos_multi_priv (n,num,den,s,c,stride,0); }
static void sincos_2pibyn_priv (size_t n, size_t nang, double *s, double *c,
int stride, int exact)
{
// nmax: number of sin/cos pairs that must be genuinely computed; the rest
// can be obtained via symmetries
size_t nmax = ((n&3)==0) ? n/8+1 : ( ((n&1)==0) ? n/4+1 : n/2+1 );
size_t ngoal = (nang<nmax) ? nang : nmax;
fracsincos_multi_priv (ngoal, 1, n, s, c, stride, exact);
size_t ndone=ngoal;
if (ndone==nang) return;
if ((n&3)==0)
{
ngoal=n/4+1;
if (nang<ngoal) ngoal=nang;
for (size_t i=ndone; i<ngoal; ++i)
{
s[i*stride]=c[(n/4-i)*stride];
c[i*stride]=s[(n/4-i)*stride];
}
ndone=ngoal;
if (ngoal==nang) return;
}
if ((n&1)==0)
{
ngoal=n/2+1;
if (nang<ngoal) ngoal=nang;
for (size_t i=ndone; i<ngoal; ++i)
{
c[i*stride]=-c[(n/2-i)*stride];
s[i*stride]= s[(n/2-i)*stride];
}
ndone=ngoal;
if (ngoal==nang) return;
}
ngoal=n;
if (nang<ngoal) ngoal=nang;
for (size_t i=ndone; i<ngoal; ++i)
{
c[i*stride]= c[(n-i)*stride];
s[i*stride]=-s[(n-i)*stride];
}
ndone=ngoal;
if (ngoal==nang) return;
for (size_t i=ndone; i<nang; ++i)
{
c[i*stride]= c[(i-n)*stride];
s[i*stride]= s[(i-n)*stride];
}
}
void sincos_2pibyn (size_t n, size_t nang, double *s, double *c, int stride)
{ sincos_2pibyn_priv (n,nang,s,c,stride,0); }
void triggen_init (struct triggen *tg, size_t n)
{
tg->n=n;
tg->ilg=1;
while ((((size_t)1)<<(2*tg->ilg))<n) ++tg->ilg;
size_t s=((size_t)1)<<tg->ilg;
tg->mask=s-1;
size_t s1=n/s+1;
tg->t1=RALLOC(double,2*(s1+s));
tg->t2=tg->t1+2*s1;
fracsincos_multi_priv(s1,s,n,&(tg->t1[1]),&(tg->t1[0]),2,1);
sincos_2pibyn_priv(n,s,&(tg->t2[1]),&(tg->t2[0]),2,1);
}
void triggen_get (const struct triggen *tg,size_t i, double *s, double *c)
{
if (i>=tg->n) i%=tg->n;
size_t i1=i>>tg->ilg,
i2=i&tg->mask;
*c = tg->t1[2*i1]*tg->t2[2*i2 ] - tg->t1[2*i1+1]*tg->t2[2*i2+1];
*s = tg->t1[2*i1]*tg->t2[2*i2+1] + tg->t1[2*i1+1]*tg->t2[2*i2 ];
}
void triggen_destroy (struct triggen *tg)
{ DEALLOC(tg->t1); }
int trigtest(int argc, const char **argv);
int trigtest(int argc, const char **argv)
{
UTIL_ASSERT((argc==1)||(argv[0]==NULL),"problem with args");
#define LENGTH 12345
static const double twopi=6.28318530717958647692;
double *sc=RALLOC(double,2*LENGTH+34);
for (int i=1; i<LENGTH; ++i)
{
sc[0]=sc[1]=sc[2*i+30+2]=sc[2*i+30+3]=10;
sincos_2pibyn(i,i+15,&sc[2],&sc[3],2);
UTIL_ASSERT(fabs(sc[0]-10.)<1e-16,"bad memory access");
UTIL_ASSERT(fabs(sc[1]-10.)<1e-16,"bad memory access");
UTIL_ASSERT(fabs(sc[2*i+30+2]-10.)<1e-16,"bad memory access");
UTIL_ASSERT(fabs(sc[2*i+30+3]-10.)<1e-16,"bad memory access");
struct triggen tg;
triggen_init(&tg,i);
for (int j=0; j<i; ++j)
{
double c, s, c2, s2;
fracsincos(j,i,&s,&c);
triggen_get(&tg,j,&s2,&c2);
UTIL_ASSERT(fabs(s2-s)<4e-16,"bad sin");
UTIL_ASSERT(fabs(c2-c)<4e-16,"bad cos");
UTIL_ASSERT(fabs(sc[2*j+2]-s)<4e-16,"bad sin");
UTIL_ASSERT(fabs(sc[2*j+3]-c)<4e-16,"bad cos");
}
triggen_destroy(&tg);
sc[0]=sc[1]=sc[2*i+2]=sc[2*i+3]=10;
sincos_multi(i,twopi*1.1/i,&sc[2],&sc[3],2);
UTIL_ASSERT(fabs(sc[0]-10.)<1e-16,"bad memory access");
UTIL_ASSERT(fabs(sc[1]-10.)<1e-16,"bad memory access");
UTIL_ASSERT(fabs(sc[2*i+2]-10.)<1e-16,"bad memory access");
UTIL_ASSERT(fabs(sc[2*i+3]-10.)<1e-16,"bad memory access");
for (int j=0; j<i; ++j)
{
double ang=twopi*1.1/i*j;
UTIL_ASSERT(fabs(sc[2*j+2]-sin(ang))<1e-15,"bad sin");
UTIL_ASSERT(fabs(sc[2*j+3]-cos(ang))<1e-15,"bad cos");
}
}
DEALLOC(sc);
return 0;
}
c_utils/trig_utils.h deleted 100644 → 0
View file @ 09727e86
/*
* This file is part of libc_utils.
*
* libc_utils is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* libc_utils is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libc_utils; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libc_utils is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*! \file trig_utils.h
*
* Copyright (C) 2016-2017 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef PLANCK_TRIGHELPER_H
#define PLANCK_TRIGHELPER_H
#include <stdlib.h>
#ifdef __cplusplus
extern "C" {
#endif
/*! Computes sin/cos(2*pi*m/n). Adapted from FFTW. */
void fracsincos(int m, int n, double *s, double *c);
/*! Computes sine and cosine of \a i*alpha for \a i=[0;n[. Stores the sines
in \a s[i*stride] and the cosines in c[i*stride]. */
void sincos_multi (size_t n, double alpha, double *s, double *c, int stride);
/*! Computes sine and cosine of \a i*2*pi*num/den for \a i=[0;n[.
Stores the sines in \a s[i*stride] and the cosines in c[i*stride]. */
void fracsincos_multi (size_t n, int num, int den, double *s, double *c,
int stride);
/*! Computes sine and cosine of \a i*2pi/n for \a i=[0;nang[. Stores the sines
in \a s[i*stride] and the cosines in c[i*stride]. */
void sincos_2pibyn (size_t n, size_t nang, double *s, double *c, int stride);
typedef struct triggen
{
size_t n, ilg, mask;
double *t1, *t2;
} triggen;
void triggen_init (struct triggen *tg, size_t n);
void triggen_get (const struct triggen *tg,size_t i, double *s, double *c);
void triggen_destroy (struct triggen *tg);
#ifdef __cplusplus
}
#endif
#endif
libfftpack/bluestein.c deleted 100644 → 0
View file @ 09727e86
/*
* This file is part of libfftpack.
*
* libfftpack is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* libfftpack is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libfftpack; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*
* Copyright (C) 2005-2016 Max-Planck-Society
* \author Martin Reinecke
*/
#include <math.h>
#include <stdlib.h>
#include "trig_utils.h"
#include "fftpack.h"
#include "bluestein.h"
/* returns the sum of all prime factors of n */
size_t prime_factor_sum (size_t n)
{
size_t result=0,tmp;
while (((tmp=(n>>1))<<1)==n)
{ result+=2; n=tmp; }
size_t limit=(size_t)sqrt(n+0.01);
for (size_t x=3; x<=limit; x+=2)
while ((tmp=(n/x))*x==n)
{
result+=x;
n=tmp;
limit=(size_t)sqrt(n+0.01);
}
if (n>1) result+=n;
return result;
}
/* returns the smallest composite of 2, 3 and 5 which is >= n */
static size_t good_size(size_t n)
{
if (n<=6) return n;
size_t bestfac=2*n;
for (size_t f2=1; f2<bestfac; f2*=2)
for (size_t f23=f2; f23<bestfac; f23*=3)
for (size_t f235=f23; f235<bestfac; f235*=5)
if (f235>=n) bestfac=f235;
return bestfac;
}
void bluestein_i (size_t n, double **tstorage, size_t *worksize)
{
size_t n2=good_size(n*2-1);
*worksize=2+2*n+8*n2+16;
*tstorage = RALLOC(double,2+2*n+8*n2+16);
((size_t *)(*tstorage))[0]=n2;
double *bk = *tstorage+2;
double *bkf = *tstorage+2+2*n;
double *work= *tstorage+2+2*(n+n2);
/* initialize b_k */
double *tmp=RALLOC(double,4*n);
sincos_2pibyn(2*n,2*n,&tmp[1],&tmp[0],2);
bk[0] = 1;
bk[1] = 0;
size_t coeff=0;
for (size_t m=1; m<n; ++m)
{
coeff+=2*m-1;
if (coeff>=2*n) coeff-=2*n;
bk[2*m ] = tmp[2*coeff ];
bk[2*m+1] = tmp[2*coeff+1];
}
DEALLOC(tmp);
/* initialize the zero-padded, Fourier transformed b_k. Add normalisation. */
double xn2 = 1./n2;
bkf[0] = bk[0]*xn2;
bkf[1] = bk[1]*xn2;
for (size_t m=2; m<2*n; m+=2)
{
bkf[m] = bkf[2*n2-m] = bk[m] *xn2;
bkf[m+1] = bkf[2*n2-m+1] = bk[m+1] *xn2;
}
for (size_t m=2*n;m<=(2*n2-2*n+1);++m)
bkf[m]=0.;
cffti (n2,work);
cfftf (n2,bkf,work);
}
void bluestein (size_t n, double *data, double *tstorage, int isign)
{
size_t n2=*((size_t *)tstorage);
double *bk = tstorage+2;
double *bkf = tstorage+2+2*n;
double *work= tstorage+2+2*(n+n2);
double *akf = tstorage+2+2*n+6*n2+16;
/* initialize a_k and FFT it */
if (isign>0)
for (size_t m=0; m<2*n; m+=2)
{
akf[m] = data[m]*bk[m] - data[m+1]*bk[m+1];
akf[m+1] = data[m]*bk[m+1] + data[m+1]*bk[m];
}
else
for (size_t m=0; m<2*n; m+=2)
{
akf[m] = data[m]*bk[m] + data[m+1]*bk[m+1];
akf[m+1] =-data[m]*bk[m+1] + data[m+1]*bk[m];
}
for (size_t m=2*n; m<2*n2; ++m)
akf[m]=0;
cfftf (n2,akf,work);
/* do the convolution */
if (isign>0)
for (size_t m=0; m<2*n2; m+=2)
{
double im = -akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
akf[m ] = akf[m]*bkf[m] + akf[m+1]*bkf[m+1];
akf[m+1] = im;
}
else
for (size_t m=0; m<2*n2; m+=2)
{
double im = akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
akf[m ] = akf[m]*bkf[m] - akf[m+1]*bkf[m+1];
akf[m+1] = im;
}
/* inverse FFT */
cfftb (n2,akf,work);
/* multiply by b_k* */
if (isign>0)
for (size_t m=0; m<2*n; m+=2)
{
data[m] = bk[m] *akf[m] - bk[m+1]*akf[m+1];
data[m+1] = bk[m+1]*akf[m] + bk[m] *akf[m+1];
}
else
for (size_t m=0; m<2*n; m+=2)
{
data[m] = bk[m] *akf[m] + bk[m+1]*akf[m+1];
data[m+1] =-bk[m+1]*akf[m] + bk[m] *akf[m+1];
}
}
libfftpack/bluestein.h deleted 100644 → 0
View file @ 09727e86
/*
* This file is part of libfftpack.
*
* libfftpack is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* libfftpack is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libfftpack; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*! \file bluestein.h
* Interface for the Bluestein algorithm
*
* Copyright (C) 2005 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef PLANCK_BLUESTEIN_H
#define PLANCK_BLUESTEIN_H
#include "c_utils.h"
#ifdef __cplusplus
extern "C" {
#endif
size_t prime_factor_sum (size_t n);
void bluestein_i (size_t n, double **tstorage, size_t *worksize);
void bluestein (size_t n, double *data, double *tstorage, int isign);
#ifdef __cplusplus
}
#endif
#endif
libfftpack/fftpack.c deleted 100644 → 0
View file @ 09727e86
This diff is collapsed.
libfftpack/fftpack.h deleted 100644 → 0
View file @ 09727e86
/*
* This file is part of libfftpack.
*
* libfftpack is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* libfftpack is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libfftpack; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*
fftpack.h : function declarations for fftpack.c
Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber
(Version 4, 1985).
Pekka Janhunen 23.2.1995
(reformatted by joerg arndt)
reformatted and slightly enhanced by Martin Reinecke (2004)
*/
#ifndef PLANCK_FFTPACK_H
#define PLANCK_FFTPACK_H
#include "c_utils.h"
#ifdef __cplusplus
extern "C" {
#endif
/*! forward complex transform */
void cfftf(size_t N, double complex_data[], double wrk[]);
/*! backward complex transform */
void cfftb(size_t N, double complex_data[], double wrk[]);
/*! initializer for complex transforms */
void cffti(size_t N, double wrk[]);
/*! forward real transform */