Skip to content

Instantly share code, notes, and snippets.

@cadebrown
Created August 6, 2020 19:12
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save cadebrown/52d316379ca6335ad8614991215dc335 to your computer and use it in GitHub Desktop.
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