Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Testing out my generated code for the Gamma and Zeta functions. Code generated from https://github.com/ChemicalDevelopment/kscript/blob/master/modules/m/tools/table_zeta.py, check my blog for more info: http://cade.site/Implementing-Zeta-Gamma/
$ gcc -std=c99 -Ofast -fno-math-errno t.c -lm -o test_gz && ./test_gz
# -- ACCURACY
|my_gamma(x)-tgamma(x)| <= 0.000000 , at x=0.000480, accurate to 14.76 digits
# -- SPEED
my_gamma(x), x in [0, 1) : 0.039 us/iter
my_gamma(x), x in [0, 4) : 0.038 us/iter
my_gamma(x), x in [0, 16) : 0.037 us/iter
my_gamma(x), x in [0, 256) : 0.052 us/iter
my_cgamma(x), x in [0i, 1i) : 0.137 us/iter
my_cgamma(x), x in [0i, 4i) : 0.135 us/iter
my_cgamma(x), x in [0i, 16i) : 0.139 us/iter
my_cgamma(x), x in [0i, 256i) : 0.155 us/iter
my_zeta(x), x in [0, 1) : 0.395 us/iter
my_zeta(x), x in [0, 4) : 0.396 us/iter
my_zeta(x), x in [0, 16) : 0.397 us/iter
my_zeta(x), x in [0, 256) : 0.416 us/iter
my_czeta(x), x in [0i, 1i) : 1.356 us/iter
my_czeta(x), x in [0i, 4i) : 1.494 us/iter
my_czeta(x), x in [0i, 16i) : 2.005 us/iter
my_czeta(x), x in [0i, 256i) : 11.456 us/iter
tgamma(x), x in [0, 1) : 0.029 us/iter
tgamma(x), x in [0, 4) : 0.036 us/iter
tgamma(x), x in [0, 16) : 0.058 us/iter
tgamma(x), x in [0, 256) : 0.050 us/iter
/* test_gz.c - implementations of the Gamma and Zeta functions
*
* NOTE: DO NOT EDIT THIS FILE, it was generated via the `gen_gz.py` script, available here: https://github.com/ChemicalDevelopment/kscript/blob/master/modules/m/tools/gen_gz.py
*
* Generation Arguments:
* $ ./tools/gen_gz.py --lgamma
*
*
* To embed this in your own, compile it like so:
* $ cc -std=c99 -Ofast -fno-math-errno my_gz.c -o my_gz.o
* Then,
* $ cc other_objs.o... my_gz.o -lm -o total.o
*
* The flags `-Ofast` and `-fno-math-errno` are not required; they are just mainly included to generate more efficient code
*
*
* You can embed all this data & use these methods in any non-commercial project for free; for commercial projects,
* contact the @author below.
*
* @author: Cade Brown <brown.cade@gmail.com>
*/
// std includes
#include <complex.h>
#include <math.h>
/* Constants */
// PI, 3.1415... constant
#define MY_PI 3.14159265358979323846264338327950288419716940l
// log(PI), a constant that should be precomputed
#define MY_LOG_PI 1.14472988584940017414342735135305871164729481l
// sqrt(2 * PI), a constant that should be precomputed
#define MY_SQRT_2PI 2.50662827463100050241576528481104525300698674l
// log(sqrt(2 * PI)), a constant that should be precomputed
#define MY_LOG_SQRT_2PI 0.91893853320467274178032973640561763986139747l
/* Forward Declarations */
// gamma(x) - compute the Gamma function evaluated at 'x'
double my_gamma(double x);
// cgamma(x) - compute the Gamma function evaluated at 'x'
double complex my_cgamma(double complex x);
// zeta(x) - compute the Riemann Zeta function evaluated at 'x'
double my_zeta(double x);
// czeta(x) - compute the Riemann Zeta function evaluated at 'x'
double complex my_czeta(double complex x);
// lgamma(x) - compute the Gamma function evaluated at 'x'
double my_lgamma(double x);
// clgamma(x) - compute the Gamma function evaluated at 'x'
double complex my_lcgamma(double complex x);
/* -- GAMMA -- */
// NOTE: correct to 14.188319365090821 digits
double my_gamma(double x) {
if (x < 0) {
int ix = (int)x;
if (x == ix) return INFINITY;
// use reflection formula, since it won't converge otherwise
// Gamma(x) = pi / (sin(pi * x) * Gamma(1 - x))
return MY_PI / (sin(MY_PI * x) * my_gamma(1 - x));
} else {
// shift off by 1 to make indexing cleaner
x -= 1.0;
// constant
static const double g = 4.7421875;
// length (n) used
static const int a_n = 15;
// array of data
static const long double a[] = {
0.99999999999999709182046422698042381484509186l,
57.15623566586292351657939348597744427325444289l,
-59.59796035547549124814226613160283600191010086l,
14.13609797474174717386341954087672601565916414l,
-0.49191381609762019978284002853030784416348279l,
0.00003399464998481188869891934155234155285737l,
0.00004652362892704857566523022496582288935897l,
-0.00009837447530487956467653837063470730301913l,
0.00015808870322491248883607241344366395885874l,
-0.00021026444172410488319269928300571190544166l,
0.00021743961811521264319614464960383469237637l,
-0.00016431810653676389021706956228018413242085l,
0.00008441822398385274329281181534545498133891l,
-0.00002619083840158140866966503624795490184143l,
0.00000368991826595316227036759674456787362014l,
};
// keep track of sum
long double sum = a[0];
// loop var
int i;
for (i = 1; i < a_n; ++i) {
sum += a[i] / (x + i);
}
// temporary variable
double tmp = x + g + 0.5;
return MY_SQRT_2PI * pow(tmp, x + 0.5) * exp(-tmp) * sum;
}
}
double complex my_cgamma(double complex x) {
double x_re = creal(x), x_im = cimag(x);
// short circuit for real only
if (x_im == 0.0) return my_gamma(x_re);
if (x_re < 0) {
// use reflection formula, since it won't converge otherwise
// Gamma(x) = pi / (sin(pi * x) * Gamma(1 - x))
return MY_PI / (csin(MY_PI * x) * my_cgamma(1 - x));
} else {
// shift off by 1 to make indexing cleaner
x -= 1.0;
// constant
static const double g = 4.7421875;
// length (n) used
static const int a_n = 15;
// array of data
static const long double a[] = {
0.99999999999999709182046422698042381484509186l,
57.15623566586292351657939348597744427325444289l,
-59.59796035547549124814226613160283600191010086l,
14.13609797474174717386341954087672601565916414l,
-0.49191381609762019978284002853030784416348279l,
0.00003399464998481188869891934155234155285737l,
0.00004652362892704857566523022496582288935897l,
-0.00009837447530487956467653837063470730301913l,
0.00015808870322491248883607241344366395885874l,
-0.00021026444172410488319269928300571190544166l,
0.00021743961811521264319614464960383469237637l,
-0.00016431810653676389021706956228018413242085l,
0.00008441822398385274329281181534545498133891l,
-0.00002619083840158140866966503624795490184143l,
0.00000368991826595316227036759674456787362014l,
};
// keep track of sum
long double complex sum = a[0];
// loop var
int i;
for (i = 1; i < a_n; ++i) {
sum += a[i] / (x + i);
}
// temporary variable
double complex tmp = x + g + 0.5;
return MY_SQRT_2PI * pow(tmp, x + 0.5) * exp(-tmp) * sum;
}
}
double my_lgamma(double x) {
if (x < 0) {
int ix = (int)x;
if (x == ix) return INFINITY;
// use reflection formula, since it won't converge otherwise
// log(Gamma(x)) = log(pi / (sin(pi * x) * Gamma(1 - x)))
// log(Gamma(x)) = log(pi) - log(sin(pi * x)) - log(Gamma(1 - x))
return MY_LOG_PI - log(sin(MY_PI * x)) - my_lgamma(1 - x);
} else {
// shift off by 1 to make indexing cleaner
x -= 1.0;
// constant
static const double g = 4.7421875;
// length (n) used
static const int a_n = 15;
// array of data
static const long double a[] = {
0.99999999999999709182046422698042381484509186l,
57.15623566586292351657939348597744427325444289l,
-59.59796035547549124814226613160283600191010086l,
14.13609797474174717386341954087672601565916414l,
-0.49191381609762019978284002853030784416348279l,
0.00003399464998481188869891934155234155285737l,
0.00004652362892704857566523022496582288935897l,
-0.00009837447530487956467653837063470730301913l,
0.00015808870322491248883607241344366395885874l,
-0.00021026444172410488319269928300571190544166l,
0.00021743961811521264319614464960383469237637l,
-0.00016431810653676389021706956228018413242085l,
0.00008441822398385274329281181534545498133891l,
-0.00002619083840158140866966503624795490184143l,
0.00000368991826595316227036759674456787362014l,
};
// keep track of sum
long double sum = a[0];
// loop var
int i;
for (i = 1; i < a_n; ++i) {
sum += a[i] / (x + i);
}
// temporary variable
double tmp = x + g + 0.5;
// log(sqrt(2pi) * pow(tmp, (x + 0.5)) * exp(-tmp) * sum)
// = log(sqrt(2pi)) + log(tmp) * (x + 0.5) - tmp + log(sum)
return MY_LOG_SQRT_2PI + log(sum) + log(tmp) * (x + 0.5) - tmp;
}
}
double complex my_clgamma(double complex x) {
double x_re = creal(x), x_im = cimag(x);
// short circuit for real only
if (x_im == 0.0) return my_gamma(x_re);
if (x_re < 0) {
// use reflection formula, since it won't converge otherwise
// log(Gamma(x)) = log(pi / (sin(pi * x) * Gamma(1 - x)))
// log(Gamma(x)) = log(pi) - log(sin(pi * x)) - log(Gamma(1 - x))
return MY_LOG_PI - clog(csin(MY_PI * x)) - my_clgamma(1 - x);
} else {
// shift off by 1 to make indexing cleaner
x -= 1.0;
// constant
static const double g = 4.7421875;
// length (n) used
static const int a_n = 15;
// array of data
static const long double a[] = {
0.99999999999999709182046422698042381484509186l,
57.15623566586292351657939348597744427325444289l,
-59.59796035547549124814226613160283600191010086l,
14.13609797474174717386341954087672601565916414l,
-0.49191381609762019978284002853030784416348279l,
0.00003399464998481188869891934155234155285737l,
0.00004652362892704857566523022496582288935897l,
-0.00009837447530487956467653837063470730301913l,
0.00015808870322491248883607241344366395885874l,
-0.00021026444172410488319269928300571190544166l,
0.00021743961811521264319614464960383469237637l,
-0.00016431810653676389021706956228018413242085l,
0.00008441822398385274329281181534545498133891l,
-0.00002619083840158140866966503624795490184143l,
0.00000368991826595316227036759674456787362014l,
};
// keep track of sum
long double complex sum = a[0];
// loop var
int i;
for (i = 1; i < a_n; ++i) {
sum += a[i] / (x + i);
}
// temporary variable
double complex tmp = x + g + 0.5;
// log(sqrt(2pi) * pow(tmp, (x + 0.5)) * exp(-tmp) * sum)
// = log(sqrt(2pi)) + log(tmp) * (x + 0.5) - tmp + log(sum)
return MY_LOG_SQRT_2PI + clog(sum) + clog(tmp) * (x + 0.5) - tmp;
}
}
/* -- ZETA -- */
double my_zeta(double x) {
if (x < 0) {
// check for negative even integers (which are exactly 0)
int ix = (int)x;
if (ix == x && ix % 2 == 0) return 0.0;
// use reflection formula with the functional equation
// Zeta(x) = 2 * (2*PI)^(x-1) * sin(x * PI/2) * Gamma(1-x) * Zeta(1 - x)
return 2 * pow(2 *MY_PI, x - 1) * sin(x * MY_PI / 2) * my_gamma(1 - x) * my_zeta(1 - x);
} else {
// for x >= 0, summation with the coefficients will work fine
// the 'n' used in computation
static const int n = 24;
// d_k for k == 0 through n-1 (0-indexing)
// NOTE: we bake in the alternating sign here, instead of at runtime (cheaper this way)
static const long double d[] = {
0.99999999999999999915316831449165727118563471l,
-0.99999999999999902360306660888083367703682290l,
0.99999999999981204316690636680631146518122476l,
-0.99999999998555166856908523792275077432760672l,
0.99999999940800649735732951813854279475607626l,
-0.99999998503354890275363160350936641430909587l,
0.99999974502366603527976642288281533472466576l,
-0.99999689655472650921631153193034098361274796l,
0.99997187750254100529229940639777459967973657l,
-0.99980442972843671759172152048046755950062770l,
0.99893193869494595536239464122713087646211307l,
-0.99533621807207493526577477521337969787914367l,
0.98348076239521758639648463726948994240630979l,
-0.95196348945735681894089485516634917708776065l,
0.88409295990333918743600045666910953960284266l,
-0.76551456344114746342744932366243798905448018l,
0.59768788135151320888308864751186656065740265l,
-0.40622785187314257981287325939000254252169923l,
0.23178649168173822888223257243452643710916944l,
-0.10672469148761619066168789786047244403958905l,
0.03778036573957455420677224392862216580892294l,
-0.00959406764049133001149894522477629131740788l,
0.00154935253821599118198119612325579859784227l,
-0.00011918096447815316784470739409659989214171l,
};
// use `long double` will prevent some rounding errors
long double sum = 0.0;
int i;
// compute elementwise sum
for (i = 0; i < n; ++i) {
sum += d[i] * pow(i + 1, -x);
}
// divide by the normalization factor in Proposition 1 (without d0, since everything has been normalized by that)
sum /= 1 - pow(2, 1 - x);
return (double)sum;
}
}
double complex my_czeta(double complex x) {
// get real and imaginary components
double x_re = creal(x), x_im = cimag(x);
// use the real-only version of the function if it is a real argument
if (x_im == 0.0) return my_zeta(x_re);
if (x_re < 0) {
// use reflection formula with the functional equation
// Zeta(x) = 2 * (2*PI)^(x-1) * sin(x * PI/2) * Gamma(1-x) * Zeta(1 - x)
return 2 * cpow(2 * MY_PI, x - 1) * csin(x * MY_PI / 2) * my_cgamma(1 - x) * my_czeta(1 - x);
} else {
// when the real component is >= 0, summation with the coefficients will work fine
/* Generated Tables */
// the 'n' used in computation
// NOTE: this table is useful for abs(imag(x)) <= 1
static const int n_0 = 28;
// d_k for k == 0 through n-1 (0-indexing)
// NOTE: we bake in the alternating sign here, instead of at runtime (cheaper this way)
static const long double d_0[] = {
0.99999999999999999999926617649791646219655661l,
-0.99999999999999999884863092523092918639732831l,
0.99999999999999969853283034454725020279874115l,
-0.99999999999996846568956995344463590854567626l,
0.99999999999823950472336973169277319096529842l,
-0.99999999993922430374373549589585909755506829l,
0.99999999858187468121214807256683494911977531l,
-0.99999997626764352398956801432265773967803573l,
0.99999970291831184801296230083148692401672589l,
-0.99999713021871960352726146797340865896910381l,
0.99997809224173699433307530482362949761670048l,
-0.99986534785856725676750529967688537324402619l,
0.99932368462638221324596244886317990614922143l,
-0.99719036481962265722265522104304760497583652l,
0.99024860989286537174998884480293456147513962l,
-0.97148193450466636550664250351931860939049695l,
0.92918124280304037965974490360987781416745159l,
-0.84955641136468558277146706848622455257113091l,
0.72443167624727090194703047043476942720548413l,
-0.56068380098692736516398399787809984635712702l,
0.38308033628147783680698743918201976251390894l,
-0.22466052107382947701816820912558595950393043l,
0.10978103562726840427490183299586757804003482l,
-0.04318423246984169543822567292066851632183446l,
0.01307397572313281112377102607815830224711621l,
-0.00284877832996473285861744804837524179398741l,
0.00039658167006469598959871666566708702468502l,
-0.00002643877800431306597324777771113913497900l,
};
// the 'n' used in computation
// NOTE: this table is useful for abs(imag(x)) <= 2
static const int n_1 = 28;
// d_k for k == 0 through n-1 (0-indexing)
// NOTE: we bake in the alternating sign here, instead of at runtime (cheaper this way)
static const long double d_1[] = {
0.99999999999999999999926617649791646219655661l,
-0.99999999999999999884863092523092918639732831l,
0.99999999999999969853283034454725020279874115l,
-0.99999999999996846568956995344463590854567626l,
0.99999999999823950472336973169277319096529842l,
-0.99999999993922430374373549589585909755506829l,
0.99999999858187468121214807256683494911977531l,
-0.99999997626764352398956801432265773967803573l,
0.99999970291831184801296230083148692401672589l,
-0.99999713021871960352726146797340865896910381l,
0.99997809224173699433307530482362949761670048l,
-0.99986534785856725676750529967688537324402619l,
0.99932368462638221324596244886317990614922143l,
-0.99719036481962265722265522104304760497583652l,
0.99024860989286537174998884480293456147513962l,
-0.97148193450466636550664250351931860939049695l,
0.92918124280304037965974490360987781416745159l,
-0.84955641136468558277146706848622455257113091l,
0.72443167624727090194703047043476942720548413l,
-0.56068380098692736516398399787809984635712702l,
0.38308033628147783680698743918201976251390894l,
-0.22466052107382947701816820912558595950393043l,
0.10978103562726840427490183299586757804003482l,
-0.04318423246984169543822567292066851632183446l,
0.01307397572313281112377102607815830224711621l,
-0.00284877832996473285861744804837524179398741l,
0.00039658167006469598959871666566708702468502l,
-0.00002643877800431306597324777771113913497900l,
};
// the 'n' used in computation
// NOTE: this table is useful for abs(imag(x)) <= 4
static const int n_2 = 28;
// d_k for k == 0 through n-1 (0-indexing)
// NOTE: we bake in the alternating sign here, instead of at runtime (cheaper this way)
static const long double d_2[] = {
0.99999999999999999999926617649791646219655661l,
-0.99999999999999999884863092523092918639732831l,
0.99999999999999969853283034454725020279874115l,
-0.99999999999996846568956995344463590854567626l,
0.99999999999823950472336973169277319096529842l,
-0.99999999993922430374373549589585909755506829l,
0.99999999858187468121214807256683494911977531l,
-0.99999997626764352398956801432265773967803573l,
0.99999970291831184801296230083148692401672589l,
-0.99999713021871960352726146797340865896910381l,
0.99997809224173699433307530482362949761670048l,
-0.99986534785856725676750529967688537324402619l,
0.99932368462638221324596244886317990614922143l,
-0.99719036481962265722265522104304760497583652l,
0.99024860989286537174998884480293456147513962l,
-0.97148193450466636550664250351931860939049695l,
0.92918124280304037965974490360987781416745159l,
-0.84955641136468558277146706848622455257113091l,
0.72443167624727090194703047043476942720548413l,
-0.56068380098692736516398399787809984635712702l,
0.38308033628147783680698743918201976251390894l,
-0.22466052107382947701816820912558595950393043l,
0.10978103562726840427490183299586757804003482l,
-0.04318423246984169543822567292066851632183446l,
0.01307397572313281112377102607815830224711621l,
-0.00284877832996473285861744804837524179398741l,
0.00039658167006469598959871666566708702468502l,
-0.00002643877800431306597324777771113913497900l,
};
// the 'n' used in computation
// NOTE: this table is useful for abs(imag(x)) <= 8
static const int n_3 = 32;
// d_k for k == 0 through n-1 (0-indexing)
// NOTE: we bake in the alternating sign here, instead of at runtime (cheaper this way)
static const long double d_3[] = {
0.99999999999999999999999936410394010364069854l,
-0.99999999999999999999869704897327235979130703l,
0.99999999999999999955460758941557170714880244l,
-0.99999999999999993915844108956828294776817835l,
0.99999999999999556043636985063984789267293190l,
-0.99999999999979939368757834664595742440589070l,
0.99999999999386089120143554392181870323091616l,
-0.99999999986491055150233468476909218628861191l,
0.99999999776946753139194572353728628597616788l,
-0.99999997147371198490863326886364361538863497l,
0.99999971045373850771112121878738005366186089l,
-0.99999762229395069013102481817727155984766830l,
0.99998395846577388379169837070677902423740805l,
-0.99990996358087794792334591671272713908646026l,
0.99957522481587252375698957721582575387979167l,
-0.99830090896564497796506751237244944605854296l,
0.99419535104495219938544763410689432529572957l,
-0.98295446518722651685729053425681827272588758l,
0.95672573151919992429159063460664081672958959l,
-0.90449212250748267024894503075137988017935178l,
0.81569498718756333837644750419743628804394751l,
-0.68698555062628653677728732633527414752712459l,
0.52834368695773606038762478152842313712266843l,
-0.36280435095577034589406386520823077843975766l,
0.21751716776255575604599178439529599555316044l,
-0.11124997091266165604283049100069226864182076l,
0.04729731398489733341649360402565925651870683l,
-0.01619245778522998991689159610488848332675205l,
0.00427566222821457909561550216121656372723691l,
-0.00081524972526999519107071989082181393427425l,
0.00009970680093230157013095135355374702793282l,
-0.00000586510593719421000770302079727923693722l,
};
// the 'n' used in computation
// NOTE: this table is useful for abs(imag(x)) <= 16
static const int n_4 = 40;
// d_k for k == 0 through n-1 (0-indexing)
// NOTE: we bake in the alternating sign here, instead of at runtime (cheaper this way)
static const long double d_4[] = {
0.99999999999999999999999999999952249802893360l,
-0.99999999999999999999999999847151619061645828l,
0.99999999999999999999999918404415433976711087l,
-0.99999999999999999999982587390155247902598209l,
0.99999999999999999998013041412444101908770422l,
-0.99999999999999999859357004581579423943294215l,
0.99999999999999993241682519472128884681929789l,
-0.99999999999999765768212306039587271390194357l,
0.99999999999993885579007288808386567798833457l,
-0.99999999999875820603596746754474401572293016l,
0.99999999997988023786242921702973554223714809l,
-0.99999999973471182002427011943222289956465373l,
0.99999999710714595058486935605018522918335542l,
-0.99999997356415576040783851614712770256692264l,
0.99999979531008717763889072830969214389964587l,
-0.99999864464934101576485535302803915829570753l,
0.99999226497786531989066728443016474214483974l,
-0.99996169714020637067530627125532261983907745l,
0.99983447671109245822661329261336064414747632l,
-0.99937264664694067363334944980498658729830693l,
0.99790544805852000411782662611376746823133034l,
-0.99381569589567493229755394683510790288784157l,
0.98379450613513066540601898851805412649650228l,
-0.96218359256456563767470881753866859132213581l,
0.92114584711407779246152939711504818676230162l,
-0.85253743676093566430104657831702811448594208l,
0.75164271565337371112386596243758683172658982l,
-0.62134680747253898752398491551570878103526699l,
0.47396013730954282615321038062096101200651090l,
-0.32844589308328164323905185615137337797449038l,
0.20364893108584408636690912160626940709279141l,
-0.11125562236219914738171619280502902674674247l,
0.05268486415345994498931710401138557134880073l,
-0.02122868072387133838976290807186237292529077l,
0.00711620510269506273359768056958286950613923l,
-0.00192702152025012162275928635756456762779283l,
0.00040437375544820228272454627187844618461607l,
-0.00006162298129064324579706683015421003598539l,
0.00000606127684825999138987542591680754452315l,
-0.00000028863223086952339951787742460988307253l,
};
// the 'n' used in computation
// NOTE: this table is useful for abs(imag(x)) <= 32
static const int n_5 = 56;
// d_k for k == 0 through n-1 (0-indexing)
// NOTE: we bake in the alternating sign here, instead of at runtime (cheaper this way)
static const long double d_5[] = {
0.99999999999999999999999999999999999999999973l,
-0.99999999999999999999999999999999999999831100l,
0.99999999999999999999999999999999999823359194l,
-0.99999999999999999999999999999999926128690977l,
0.99999999999999999999999999999983465898154359l,
-0.99999999999999999999999999997700980597076557l,
0.99999999999999999999999999782524866304560153l,
-0.99999999999999999999999985122174222199046902l,
0.99999999999999999999999230847133582970390137l,
-0.99999999999999999999968941527854580219624376l,
0.99999999999999999998994900786829807550262264l,
-0.99999999999999999973391544165607211669886744l,
0.99999999999999999414013861462809192979073756l,
-0.99999999999999989114579851328829673453458629l,
0.99999999999999827429814200177944835186262424l,
-0.99999999999997641897809536345294607574506819l,
0.99999999999971988334641889503178427825230352l,
-0.99999999999708593461476745455782357672231323l,
0.99999999997328005512584157827412161717982961l,
-0.99999999978283301921443456800450594083996066l,
0.99999999842772910984480776416300978226781619l,
-0.99999998981550077991379978434806555496610359l,
0.99999994074582773728363803889082519010750856l,
-0.99999968928153226084849031144618471164375202l,
0.99999852692795370222006810320473484044691989l,
-0.99999366876442535921800956753434844001607868l,
0.99997526929441077508804126728044157503592443l,
-0.99991200906207551728269952217266954449619072l,
0.99971425921891841918210526181629645682198416l,
-0.99915151556361836505918003452266124543694627l,
0.99769219727106059758786546204933603625201751l,
-0.99424107542848083231461386601511586942229704l,
0.98679445835743818403155610671508723563792995l,
-0.97213035089446189202799621147810777218563784l,
0.94577615600709361307164229966941860476962207l,
-0.90256182401786860906867687906386618118062727l,
0.83794321022177394345860886750110351961680171l,
-0.74990304037000483222386237266549944906705049l,
0.64073322975381113429277671906935040158535898l,
-0.51771270390459386030050437415780582060698932l,
0.39196703982453949479574498363113208751833932l,
-0.27564946076584421902621145526680155212287626l,
0.17855130957226898882272821730346194856556201l,
-0.10565464011257803624315010376819413128477670l,
0.05663791409657894743964068259792784035459347l,
-0.02726459513568061707049646017379822781216158l,
0.01167273781816602316695167703323062172776510l,
-0.00439582431376662450229110562925946720656427l,
0.00143718448105686899258042594001330043412867l,
-0.00040138042432869033317403608036557867748385l,
0.00009377800748214030704729000083383100429841l,
-0.00001781794647868240892397006040199581884086l,
0.00000264295296230825825332922092065458500225l,
-0.00000028694318725933351308956497961616163539l,
0.00000002027126048847910323897048577069107777l,
-0.00000000069900898236134838755070640588589923l,
};
// the 'n' used in computation
// NOTE: this table is useful for abs(imag(x)) <= 64
static const int n_6 = 84;
// d_k for k == 0 through n-1 (0-indexing)
// NOTE: we bake in the alternating sign here, instead of at runtime (cheaper this way)
static const long double d_6[] = {
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-0.99999999999999999999999999999999999999999998l,
0.99999999999999999999999999999999999999999812l,
-0.99999999999999999999999999999999999999982777l,
0.99999999999999999999999999999999999998732037l,
-0.99999999999999999999999999999999999923406103l,
0.99999999999999999999999999999999996138004976l,
-0.99999999999999999999999999999999835124512139l,
0.99999999999999999999999999999993967929140860l,
-0.99999999999999999999999999999808915238350192l,
0.99999999999999999999999999994711758259353996l,
-0.99999999999999999999999999871144316344294696l,
0.99999999999999999999999997216601841934116025l,
-0.99999999999999999999999946376260625065487830l,
0.99999999999999999999999073617069735487370473l,
-0.99999999999999999999985579807140441568789463l,
0.99999999999999999999796866165106024504490551l,
-0.99999999999999999997400294077279572433819155l,
0.99999999999999999969665496174366003878371995l,
-0.99999999999999999676242654687215107161151427l,
0.99999999999999996830085349100688083347508175l,
-0.99999999999999971451352142053808402102163736l,
0.99999999999999762917267725189382761531651705l,
-0.99999999999998180429938671863699013608782424l,
0.99999999999987067256209783254518795144790809l,
-0.99999999999914711151000858581738864034973377l,
0.99999999999477200772778243620832435350114992l,
-0.99999999997016536342532263719492036443759955l,
0.99999999984125592834825712296052036203004463l,
-0.99999999921138788118081485879243753453143260l,
0.99999999633755330916723438745903302189887646l,
-0.99999998408039882208969013356757801059674882l,
0.99999993516360052874934334628287299356152085l,
-0.99999975233238373073001189757040908706688255l,
0.99999911184454513766861976553022698425180473l,
-0.99999700735695262713386790018169729261731836l,
0.99999051761466361730926016854179649368710366l,
-0.99997172488350441621934273169454340042379522l,
0.99992059962167267114348061379295152131091797l,
-0.99978988234797552518172578550024219743679487l,
0.99947567374598300662769399326869866737912009l,
-0.99876545160283816944768161458014911877014574l,
0.99725560654677543882859389549971067047729545l,
-0.99423655175762749188901685318107254454099225l,
0.98855828911379569673599417358783598848705632l,
-0.97851357367784242214186824533049343388339408l,
0.96180363817852508147090819664249452396580063l,
-0.93566827814509145466047939811952192523754664l,
0.89724776098967691223614879676020568197036261l,
-0.84418322687010520059754039321472803488908090l,
0.77535982370537110483974423018647179441783041l,
-0.69158781923222791112120333309882563306173603l,
0.59595962701955740482840417411451710644025138l,
-0.49366807678076861466962842589778030206924709l,
0.39123326106545491023892161921015170945682682l,
-0.29530730980895195331912139595789743061425144l,
0.21140669416760987062110572731214132782380065l,
-0.14296509989713602436687681869241596854432688l,
0.09097710739887008849413945626399761615856519l,
-0.05427240667260481236216507515770821292377577l,
0.03023562486458792669720567823152146747901778l,
-0.01566951311918436002783290936142737866724657l,
0.00752328025416236533496147195630068455399675l,
-0.00333167159495714136520301733464234961621540l,
0.00135433513598982604956156751825881854073581l,
-0.00050264451456668514045698939337161586393349l,
0.00016928076123186988563070250105465522023562l,
-0.00005136730894225738548241547934141000058471l,
0.00001392724052818112371442656842051011194027l,
-0.00000334016077977119197511158331894021054282l,
0.00000069984862920050317780246103236883917654l,
-0.00000012611133139378360212031579396609656575l,
0.00000001914657793999154888726737648389790200l,
-0.00000000238132108969954280123608376557215804l,
0.00000000023294304992155774461414924270494400l,
-0.00000000001680379770250120030865287785548577l,
0.00000000000079468056222991653920492026151871l,
-0.00000000000001848094330767247765592837817485l,
};
// the 'n' used in computation
// NOTE: this table is useful for abs(imag(x)) <= 128
static const int n_7 = 144;
// d_k for k == 0 through n-1 (0-indexing)
// NOTE: we bake in the alternating sign here, instead of at runtime (cheaper this way)
static const long double d_7[] = {
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-0.99999999999999999999999999999999999999999993l,
0.99999999999999999999999999999999999999999855l,
-0.99999999999999999999999999999999999999997330l,
0.99999999999999999999999999999999999999953770l,
-0.99999999999999999999999999999999999999247435l,
0.99999999999999999999999999999999999988463944l,
-0.99999999999999999999999999999999999833239166l,
0.99999999999999999999999999999999997723598271l,
-0.99999999999999999999999999999999970617909193l,
0.99999999999999999999999999999999640974987909l,
-0.99999999999999999999999999999995842090951962l,
0.99999999999999999999999999999954311548608313l,
-0.99999999999999999999999999999523166569487266l,
0.99999999999999999999999999995268707206551424l,
-0.99999999999999999999999999955326716683108809l,
0.99999999999999999999999999598253336518687039l,
-0.99999999999999999999999996556063293662803555l,
0.99999999999999999999999971835600607701859118l,
-0.99999999999999999999999780105085889573465320l,
0.99999999999999999999998359749939482751606621l,
-0.99999999999999999999988302830331688224945688l,
0.99999999999999999999920198932987642765275771l,
-0.99999999999999999999478851412086356968773369l,
0.99999999999999999996740174704852936936271271l,
-0.99999999999999999980458867307722146284332191l,
0.99999999999999999887679512317186936105280376l,
-0.99999999999999999380641273586114214741730635l,
0.99999999999999996721997291370141597505030689l,
-0.99999999999999983340579571811724588484039392l,
0.99999999999999918663727260612709044882581455l,
-0.99999999999999618351427787721141410373820254l,
0.99999999999998278248155946323285837028272961l,
-0.99999999999992529162576785715647400214816487l,
0.99999999999968809940651149419320217623471823l,
-0.99999999999874668532698236367230934780305406l,
0.99999999999515110636165609163195276898389903l,
-0.99999999998193248141450216142859943915385544l,
0.99999999993514315690112309737328755467701048l,
-0.99999999977564467037630272710793579699353943l,
0.99999999925190606354505803749046522577186421l,
-0.99999999759494380647192030803641189126653029l,
0.99999999254326885993678109595766925625075411l,
-0.99999997769893276763361906926922231061724420l,
0.99999993564953284593573023948637581468499797l,
-0.99999982080857663920467278221814538463111388l,
0.99999951837537879328048155935648103667986911l,
-0.99999875029106680363174194573955570854972367l,
0.99999686889745098382547437604132356072125987l,
-0.99999242381492917590780406543606718242848369l,
0.99998229308047106179355388426889875422820212l,
-0.99996001960310099447634038879883869952699261l,
0.99991277636776945112432762532613164051312987l,
-0.99981610280610963579779928041491275611703716l,
0.99962525039545901378091333915662249489195936l,
-0.99926174730075654121690363829850752467837573l,
0.99859381716169423918884229350128307716637678l,
-0.99740981768589979243389732777209908543027073l,
0.99538517858229128848294143637519445956152939l,
-0.99204561008942241482956507110780603839899686l,
0.98673256598810464749647778111381688666741965l,
-0.97858058671978441095699736142296060334762005l,
0.96651918237351112204048172501887233906297952l,
-0.94931284071266069465501472240905296338136270l,
0.92564996078835768910169356408104217786950670l,
-0.89428396010896968953783705791525790077776774l,
0.85421772603869123590647541484070514455816097l,
-0.80490807201490962318992329388361441739082269l,
0.74645381234694789349443549963128928221296390l,
-0.67972430215633619529721762283111665186892205l,
0.60638890578102575109756159440801235912398762l,
-0.52882338531549864120863916485769979453016308l,
0.44989430600473576062525291977776713076159922l,
-0.37265109052391150467568657507818254330653867l,
0.29997903151477853923521292469052989399797325l,
-0.23427779429713451853168149503953747830957326l,
0.17722470347887992807724228254371459100294244l,
-0.12966168625511452863736052933320891512974752l,
0.09161536920454525520082462562382336716813157l,
-0.06243038731787003864341087527912749847428471l,
0.04097476136465218537517519580197665436602286l,
-0.02586809749183906069012498766006032892319955l,
0.01568851897728729545752589240747978130938163l,
-0.00912874653353444898472732498886504415678407l,
0.00508966259664760086286426087076444821469081l,
-0.00271544504042932066820163575027091393773605l,
0.00138444428921604237725440651605484169156444l,
-0.00067356167531100103921869037760848932355671l,
0.00031224975463384082939988824074361769446545l,
-0.00013771110245834068053021620895100658661382l,
0.00005768379157591993584276198237942931784642l,
-0.00002290768364059413505817229206993968615271l,
0.00000860825856700506700298277831183813416099l,
-0.00000305456780671681571650731559324445357327l,
0.00000102116411774357391402149089276418930166l,
-0.00000032082395562973770124565022483974165535l,
0.00000009446345676400590566623318482444715586l,
-0.00000002598678162240081858422722083488372693l,
0.00000000665638686888521587670262058366321731l,
-0.00000000158137839500352173239308236429775066l,
0.00000000034691687432959612972319468931696147l,
-0.00000000006992033177131011224247232839003291l,
0.00000000001287098048395716187413597395324471l,
-0.00000000000214903686280732662828947056867143l,
0.00000000000032278820644623612607613858740949l,
-0.00000000000004318081199855161024779972424506l,
0.00000000000000508149750657549676587228310083l,
-0.00000000000000051785137614975979629834293008l,
0.00000000000000004477143208928946063816162226l,
-0.00000000000000000319278806932655863379693521l,
0.00000000000000000018030566849463306614839286l,
-0.00000000000000000000755990198974588789197610l,
0.00000000000000000000020920123019387786812519l,
-0.00000000000000000000000286577027662846394692l,
};
// the 'n' used in computation
// NOTE: this table is useful for abs(imag(x)) <= 256
static const int n_8 = 256;
// d_k for k == 0 through n-1 (0-indexing)
// NOTE: we bake in the alternating sign here, instead of at runtime (cheaper this way)
static const long double d_8[] = {
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
1.00000000000000000000000000000000000000000000l,
-1.00000000000000000000000000000000000000000000l,
0.99999999999999999999999999999999999999999999l,
-0.99999999999999999999999999999999999999999995l,
0.99999999999999999999999999999999999999999959l,
-0.99999999999999999999999999999999999999999681l,
0.99999999999999999999999999999999999999997591l,
-0.99999999999999999999999999999999999999982253l,
0.99999999999999999999999999999999999999872550l,
-0.99999999999999999999999999999999999999107453l,
0.99999999999999999999999999999999999993903649l,
-0.99999999999999999999999999999999999959379793l,
0.99999999999999999999999999999999999735923359l,
-0.99999999999999999999999999999999998324619563l,
0.99999999999999999999999999999999989625427936l,
-0.99999999999999999999999999999999937284370046l,
0.99999999999999999999999999999999629827901578l,
-0.99999999999999999999999999999997866339700344l,
0.99999999999999999999999999999987988147248504l,
-0.99999999999999999999999999999933942011983700l,
0.99999999999999999999999999999645074737292506l,
-0.99999999999999999999999999998136589718275228l,
0.99999999999999999999999999990439027785948768l,
-0.99999999999999999999999999952051186461099927l,
0.99999999999999999999999999764931485266533608l,
-0.99999999999999999999999998873298843007312764l,
0.99999999999999999999999994719488789443619588l,
-0.99999999999999999999999975798078802923239191l,
0.99999999999999999999999891511797953877908331l,
-0.99999999999999999999999524309343634886729832l,
0.99999999999999999999997959552888911543074821l,
-0.99999999999999999999991436872531098412143347l,
0.99999999999999999999964836495806411382244506l,
-0.99999999999999999999858697608917096284398666l,
0.99999999999999999999444290136752074483875177l,
-0.99999999999999999997860904263133767351248838l,
0.99999999999999999991939966117480080376472220l,
-0.99999999999999999970269118365575808338591797l,
0.99999999999999999892630359421776507366059944l,
-0.99999999999999999620339589053811001757103541l,
0.99999999999999998685425013200231806724530099l,
-0.99999999999999995542548418135656171297955686l,
0.99999999999999985197571238599280041289640651l,
-0.99999999999999951854338767881553025629105251l,
0.99999999999999846612191430501392835998233035l,
-0.99999999999999521298629028192174011992188142l,
0.99999999999998536447699027524798127331240836l,
-0.99999999999995616205614667569267135695196919l,
0.99999999999987134754072358093431802398828203l,
-0.99999999999963005370479496599662841549353848l,
0.99999999999895759205523587234466765901622430l,
-0.99999999999712167349183186532051467756429979l,
0.99999999999221121480195587873184858647311346l,
-0.99999999997934382532299129429628572737254096l,
0.99999999994630882078428424772417185571489034l,
-0.99999999986321122374138274932828743964174335l,
0.99999999965840356519564496932932425194778470l,
-0.99999999916379703266007938390939351812093873l,
0.99999999799338019184435160458619337148460044l,
-0.99999999527945232532350300051307447312833302l,
0.99999998911290791951006832082570855940437995l,
-0.99999997538264401174085432073859452493090682l,
0.99999994542443607476879897705258453304620478l,
-0.99999988136870451164515839014728694867991669l,
0.99999974715109552957123445455713015454877470l,
-0.99999947155486283530785926327776661822600073l,
0.99999891699161439000747987586005694987737534l,
-0.99999782342943075656585316436888083755110489l,
0.99999571018451536057062937206252111279485384l,
-0.99999170826560581041549742223978906665241384l,
0.99998428156747311865418233297131355892715879l,
-0.99997077552846459505764866243279635143554743l,
0.99994670645171379876545466358444767340351933l,
-0.99990467413979320993824710405472884198899109l,
0.99983274721671069848609882292528047811528009l,
-0.99971213960876623370807195401608967647929817l,
0.99951397672887773335126325041592554923435574l,
-0.99919495348173624429400702668059793101394012l,
0.99869172931192117211532954680341872758671321l,
-0.99791400595340509376825174891347740802750445l,
0.99673640740401800346058493075768390673414914l,
-0.99498953804290360790159126888029124027568696l,
0.99245092916714674015536519211530429189473314l,
-0.98883697146394160541412265343246020290366541l,
0.98379732040408312919850212056461031780422977l,
-0.97691357627959606992356486576534619250252450l,
0.96770418342051418220062375590650125557170367l,
-0.95563736108075812388171593046962445235632303l,
0.94015338551204459772028553499752589955782574l,
-0.92069664736879928399045191977858598838066375l,
0.89675664089686597794176010798214473554402755l,
-0.86791552126282237505786740402953835464488222l,
0.83389830483268588269684956608722640642944389l,
-0.79462046442903094301123388157323805459722594l,
0.75022689131545940716027139898880891218546212l,
-0.70111622129401185417586771152804896398519139l,
0.64794550790823383724478694124888404712851595l,
-0.59161215925369938881041152114371556472762155l,
0.53321273142381293414283477411732688264721432l,
-0.47398120939105807695164323773231853048191614l,
0.41521230754673450826453329903408801123698993l,
-0.35817757001912868422720693229744462206814350l,
0.30404321883007573306845210245473609545121359l,
-0.25379854745618874873652935992439093693002881l,
0.20820219625809654323507055949058410104744421l,
-0.16775113158321094049919192591050686911651368l,
0.13267403388343141656082023250221691880636464l,
-0.10294762983723705599362114846693166261324941l,
0.07833180601728746953865543008026231353957269l,
-0.05841751625999305206772962484383574363307839l,
0.04268074273686997494583696625091869650967784l,
-0.03053607493645318287965282501745883310554243l,
0.02138462662748411611772121247927510002493883l,
-0.01465269312751082595470250600467682297138032l,
0.00981939748270250647887371447529968428928710l,
-0.00643326250931798227524739018920914035267124l,
0.00411894900890503819161816709256313682898036l,
-0.00257621004727386865096465513527334872508992l,
0.00157342972221360844953987236303498645756113l,
-0.00093802821650676102820308366632207175765380l,
0.00054565377614753419147605657023875821854930l,
-0.00030958481892648211762255455381030816809299l,
0.00017124957003896143061012155401132699601970l,
-0.00009231870586364996699661097245746857545812l,
0.00004848189893700295825597725212314865185659l,
-0.00002479214820287782457152768574252869805572l,
0.00001233966036911746989360582956835991001400l,
-0.00000597517930630814332110286166797582117563l,
0.00000281355335727534795249160212796231369609l,
-0.00000128767581804991985821030484436245449533l,
0.00000057251643344707469559552580839182598594l,
-0.00000024715742659279604133286842294676258322l,
0.00000010354484826474689119793416103972601769l,
-0.00000004207333970084619788924218228882014024l,
0.00000001657114192913771871754520708776438275l,
-0.00000000632260277076656288136553910219733469l,
0.00000000233536662493775951302205422758475775l,
-0.00000000083451299529880553518221953982111950l,
0.00000000028828183893300275695618734604505940l,
-0.00000000009620055317799518659098921196996134l,
0.00000000003098596253516576013242490130386893l,
-0.00000000000962520217587188521320109565172334l,
0.00000000000288084417412695310933063592130896l,
-0.00000000000082999856349983789793840109905378l,
0.00000000000022995227697704852295910997736457l,
-0.00000000000006119637062352605673393372468540l,
0.00000000000001562544368514569844559665564208l,
-0.00000000000000382306435786543671351933421460l,
0.00000000000000089510836107212370070219820829l,
-0.00000000000000020025805493394670490053739985l,
0.00000000000000004274329964176853919179344471l,
-0.00000000000000000868888432067405994898926911l,
0.00000000000000000167905600909534130399542657l,
-0.00000000000000000030780982302595715495884130l,
0.00000000000000000005341205008412009783702636l,
-0.00000000000000000000875100181965008165378698l,
0.00000000000000000000135002810725219325770731l,
-0.00000000000000000000019550659026330658598919l,
0.00000000000000000000002648608634087302053150l,
-0.00000000000000000000000334365797796323239178l,
0.00000000000000000000000039160549185724073184l,
-0.00000000000000000000000004233376077869550503l,
0.00000000000000000000000000419921808662759366l,
-0.00000000000000000000000000037955686794043589l,
0.00000000000000000000000000003100400624810521l,
-0.00000000000000000000000000000226587704336024l,
0.00000000000000000000000000000014633466946189l,
-0.00000000000000000000000000000000822092916375l,
0.00000000000000000000000000000000039356127671l,
-0.00000000000000000000000000000000001560978678l,
0.00000000000000000000000000000000000049244280l,
-0.00000000000000000000000000000000000001158430l,
0.00000000000000000000000000000000000000018063l,
-0.00000000000000000000000000000000000000000140l,
};
// absolute value of the imaginary portion (used for discriminating amongst tables)
double x_im_abs = fabs(x_im);
// the sum of all elements in the series
long double complex sum = 0.0;
// tmp vars
int i;
if (0) {
// just used for easier code generation
} else if (x_im_abs <= 1) {
for (i = 0; i < n_0; ++i) {
sum += d_0[i] * cpow(i + 1, -x);
}
} else if (x_im_abs <= 2) {
for (i = 0; i < n_1; ++i) {
sum += d_1[i] * cpow(i + 1, -x);
}
} else if (x_im_abs <= 4) {
for (i = 0; i < n_2; ++i) {
sum += d_2[i] * cpow(i + 1, -x);
}
} else if (x_im_abs <= 8) {
for (i = 0; i < n_3; ++i) {
sum += d_3[i] * cpow(i + 1, -x);
}
} else if (x_im_abs <= 16) {
for (i = 0; i < n_4; ++i) {
sum += d_4[i] * cpow(i + 1, -x);
}
} else if (x_im_abs <= 32) {
for (i = 0; i < n_5; ++i) {
sum += d_5[i] * cpow(i + 1, -x);
}
} else if (x_im_abs <= 64) {
for (i = 0; i < n_6; ++i) {
sum += d_6[i] * cpow(i + 1, -x);
}
} else if (x_im_abs <= 128) {
for (i = 0; i < n_7; ++i) {
sum += d_7[i] * cpow(i + 1, -x);
}
} else {
for (i = 0; i < n_8; ++i) {
sum += d_8[i] * cpow(i + 1, -x);
}
}
// transform by the normalization factor
sum /= 1 - cpow(2, 1 - x);
return sum;
}
}
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
// return a random number in [0, 1)
static double my_rand() {
return (double)rand() / (double)(RAND_MAX);
}
// start time
static struct timeval my_stime;
// return time (in us) since it began
static double my_time() {
struct timeval curtime;
gettimeofday(&curtime, NULL);
return (curtime.tv_sec - my_stime.tv_sec) * 1000000.0 + (curtime.tv_usec - my_stime.tv_usec);
}
int main(int argc, char** argv) {
// number of trials
int N = 0x100000;
// temp vars
int i;
// initialize the timer
gettimeofday(&my_stime, NULL);
// calculate some random numbers
// in the range [a, b)
double* X_0_1 = malloc(sizeof(*X_0_1) * N);
double* X_0_4 = malloc(sizeof(*X_0_4) * N);
double* X_0_16 = malloc(sizeof(*X_0_16) * N);
double* X_0_256 = malloc(sizeof(*X_0_256) * N);
// random double complex
double complex* Y_0_1 = malloc(sizeof(*Y_0_1) * N);
double complex* Y_0_4 = malloc(sizeof(*Y_0_4) * N);
double complex* Y_0_16 = malloc(sizeof(*Y_0_16) * N);
double complex* Y_0_256 = malloc(sizeof(*Y_0_256) * N);
// result storage
double* Z = malloc(sizeof(*Z) * N);
double complex* W = malloc(sizeof(*W) * N);
// to seed the random number generator
//srand(time(NULL));
for (i = 0; i < N; ++i) {
X_0_1[i] = my_rand();
X_0_4[i] = my_rand() * 4;
X_0_16[i] = my_rand() * 16;
X_0_256[i] = my_rand() * 256;
Y_0_1[i] = my_rand() + my_rand() * I;
Y_0_4[i] = (my_rand() + my_rand() * I) * 4;
Y_0_16[i] = (my_rand() + my_rand() * I) * 16;
Y_0_256[i] = (my_rand() + my_rand() * I) * 256;
}
// start & elapsed time
double st, et;
// macro for a specific test
#define MY_TEST(msg, ...) { \
st = my_time(); \
for (i = 0; i < N; ++i) { __VA_ARGS__; } \
et = my_time() - st; \
printf(msg " %8.3f us/iter\n", et / N); \
}
printf("# -- ACCURACY\n");
double amaxerr = 0.0;
int amaxerri = -1;
for (i = 0; i < N; ++i) {
double x = X_0_16[i];
if (x > 0.0) {
double my_val = my_gamma(x);
double c_val = tgamma(x);
double this_err = fabs(my_val - c_val) / c_val;
if (this_err > amaxerr) {
amaxerr = this_err;
amaxerri = i;
}
}
}
printf("|my_gamma(x)-tgamma(x)| <= %f , at x=%f, accurate to %.2f digits\n", amaxerr, X_0_16[amaxerri], -log10(amaxerr));
printf("\n# -- SPEED\n");
// run tests on our functions
MY_TEST("my_gamma(x), x in [0, 1) :", Z[i] = my_gamma(X_0_1[i]))
MY_TEST("my_gamma(x), x in [0, 4) :", Z[i] = my_gamma(X_0_4[i]))
MY_TEST("my_gamma(x), x in [0, 16) :", Z[i] = my_gamma(X_0_16[i]))
MY_TEST("my_gamma(x), x in [0, 256) :", Z[i] = my_gamma(X_0_256[i]))
MY_TEST("my_cgamma(x), x in [0i, 1i) :", W[i] = my_cgamma(Y_0_1[i]))
MY_TEST("my_cgamma(x), x in [0i, 4i) :", W[i] = my_cgamma(Y_0_4[i]))
MY_TEST("my_cgamma(x), x in [0i, 16i) :", W[i] = my_cgamma(Y_0_16[i]))
MY_TEST("my_cgamma(x), x in [0i, 256i) :", W[i] = my_cgamma(Y_0_256[i]))
MY_TEST("my_zeta(x), x in [0, 1) :", Z[i] = my_zeta(X_0_1[i]))
MY_TEST("my_zeta(x), x in [0, 4) :", Z[i] = my_zeta(X_0_4[i]))
MY_TEST("my_zeta(x), x in [0, 16) :", Z[i] = my_zeta(X_0_16[i]))
MY_TEST("my_zeta(x), x in [0, 256) :", Z[i] = my_zeta(X_0_256[i]))
MY_TEST("my_czeta(x), x in [0i, 1i) :", W[i] = my_czeta(Y_0_1[i]))
MY_TEST("my_czeta(x), x in [0i, 4i) :", W[i] = my_czeta(Y_0_4[i]))
MY_TEST("my_czeta(x), x in [0i, 16i) :", W[i] = my_czeta(Y_0_16[i]))
MY_TEST("my_czeta(x), x in [0i, 256i) :", W[i] = my_czeta(Y_0_256[i]))
// compare to std C functions
MY_TEST("tgamma(x), x in [0, 1) :", Z[i] = tgamma(X_0_1[i]))
MY_TEST("tgamma(x), x in [0, 4) :", Z[i] = tgamma(X_0_4[i]))
MY_TEST("tgamma(x), x in [0, 16) :", Z[i] = tgamma(X_0_16[i]))
MY_TEST("tgamma(x), x in [0, 256) :", Z[i] = tgamma(X_0_256[i]))
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment