Created
August 6, 2020 19:12
Star
You must be signed in to star a gist
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/
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
$ 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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* 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