added gamma, zeta, and bernoulli numbers from my old code base
This commit is contained in:
parent
db964afe9f
commit
46ac5380eb
651
miscfunc.cc
651
miscfunc.cc
@ -18,8 +18,9 @@
|
||||
|
||||
#include "miscfunc.h"
|
||||
#include "laerror.h"
|
||||
#include <stdio.h>
|
||||
#include <iostream>
|
||||
#include <stdint.h>
|
||||
#include <math.h>
|
||||
|
||||
|
||||
namespace LA {
|
||||
@ -55,6 +56,92 @@ return a[n];
|
||||
}
|
||||
|
||||
|
||||
#define EPS 1e-15
|
||||
|
||||
static double stirfact(double n) /* n should not be less than 9 or 7 at worst */
|
||||
{
|
||||
double logz,cor,nn,z;
|
||||
|
||||
cor=0.0;
|
||||
z=n+1.0;
|
||||
logz=log(z);
|
||||
if(z>=2.0) /*does not have sense for smaller anyway and could make troubles*/
|
||||
{
|
||||
int i,ii,m;
|
||||
m=(int)ceil(0.5-log(EPS/10000.)/logz/2.0);
|
||||
if(m<2)m=2; if(m>16)m=16;
|
||||
nn=z*z;
|
||||
for(i=m;i>0;i--) {ii=2*i;cor = cor/nn + bernoulli_number(ii)/((double)ii*(ii-1));}
|
||||
cor /=z;
|
||||
}
|
||||
|
||||
/*0.5*ln(2*M_PI)=0.91893853320467275836*/
|
||||
return((z-0.5)*logz-z+0.918938533204672741780329736406+cor);
|
||||
/*cor=log((1.0/(24.0*n)+1.0)/(12.0*n)+1.0)*/
|
||||
}
|
||||
|
||||
|
||||
|
||||
#define LIMITFACT 180
|
||||
double factorialln(int n)
|
||||
{
|
||||
static int top=MAXFACT;
|
||||
static double a[LIMITFACT+1]={
|
||||
0.0,
|
||||
0.0,
|
||||
0.69314718055994530843,
|
||||
1.7917594692280549573,
|
||||
3.1780538303479457518,
|
||||
4.7874917427820458116,
|
||||
6.579251212010101213,
|
||||
8.5251613610654146669,
|
||||
10.604602902745250859,
|
||||
12.801827480081469091,
|
||||
15.10441257307551588,
|
||||
17.502307845873886549,
|
||||
19.987214495661884683,
|
||||
22.552163853123424531,
|
||||
25.191221182738679829,
|
||||
27.899271383840890337,
|
||||
30.671860106080671926,
|
||||
33.505073450136890756,
|
||||
36.395445208033052609,
|
||||
39.339884187199494647,
|
||||
42.335616460753485057,
|
||||
45.380138898476907627,
|
||||
48.471181351835227247,
|
||||
51.606675567764376922,
|
||||
54.784729398112318677,
|
||||
58.003605222980517908,
|
||||
61.261701761002001376,
|
||||
64.557538627006337606,
|
||||
67.889743137181540078,
|
||||
71.257038967168014665,
|
||||
74.658236348830158136};
|
||||
|
||||
|
||||
if(n < 0) laerror("negative argument of factln");
|
||||
if(n < 2) return(0.0);
|
||||
if(n <= LIMITFACT) /*supposed LIMITFACT > MAXFACT*/
|
||||
if(a[n]) return(a[n]); else
|
||||
{
|
||||
if(n<=MAXFACT) return(a[n]= log(factorial(n)));
|
||||
|
||||
{
|
||||
int j;
|
||||
if(!a[top]) a[top]= log(factorial(top));
|
||||
while (top<n) {
|
||||
j=top++;
|
||||
a[top]=a[j]+log((double)top);
|
||||
}
|
||||
return a[n];
|
||||
}
|
||||
}
|
||||
else return(stirfact((double)n));
|
||||
}
|
||||
|
||||
|
||||
|
||||
#define MAXBINOM 100
|
||||
#define ibidxmaxeven(n) ((n-2)*(n-2)/4)
|
||||
#define ibidx(n,k) (k-2+(n&1?(n-3)*(n-3)/4:(n/2-1)*(n/2-2)))
|
||||
@ -110,4 +197,566 @@ while(i);
|
||||
return y ;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
static void splintx(double ary[][3], int n,double x,double *y)
|
||||
{
|
||||
int klo,khi,k;
|
||||
double h,b,a;
|
||||
|
||||
klo=1;
|
||||
khi=n;
|
||||
while (khi-klo > 1) {
|
||||
k=(khi+klo) >> 1;
|
||||
if (ary[k-1][0] > x) khi=k;
|
||||
else klo=k;
|
||||
}
|
||||
h=ary[khi-1][0]-ary[klo-1][0];
|
||||
if (h == 0.0) laerror("error in interpolation for zeta");
|
||||
a=(ary[khi-1][0]-x)/h;
|
||||
b=(x-ary[klo-1][0])/h;
|
||||
*y=a*ary[klo-1][1]+b*ary[khi-1][1]+((a*a*a-a)*ary[klo-1][2]+(b*b*b-b)*ary[khi-1][2])*(h*h)/6.0;
|
||||
}
|
||||
|
||||
|
||||
double zeta(double x)
|
||||
{
|
||||
static int nzet1=87;
|
||||
static double zet1[87][3]= {{-1.08,-0.070868879815766575,0.},
|
||||
{-1.06,-0.073845053393926163,-0.28794916575541835},
|
||||
{-1.04,-0.076912648952586551,-0.21953304449032346},
|
||||
{-1.02,-0.080074446235881738,-0.24694452580526705},
|
||||
{-1.00,-0.083333333333333329,-0.24903606463466105},
|
||||
{-0.98,-0.086692312167914967,-0.25828727260680795},
|
||||
{-0.96,-0.090154504320474158,-0.26601461460137482},
|
||||
{-0.94,-0.093723157214775366,-0.27456539511795663},
|
||||
{-0.92,-0.097401650689449851,-0.28333251052626718},
|
||||
{-0.90,-0.10119350398535187,-0.29250188118973602},
|
||||
{-0.88,-0.105102383179246,-0.302048434596348},
|
||||
{-0.86,-0.10913210909741072,-0.31200524448380662},
|
||||
{-0.84,-0.11328666574566436,-0.32239153880210036},
|
||||
{-0.82,-0.11757020929552797,-0.33323212445747097},
|
||||
{-0.80,-0.12198707766977114,-0.34455232906170852},
|
||||
{-0.78,-0.12654180077447311,-0.35637951617729297},
|
||||
{-0.76,-0.13123911142901706,-0.36874285385897365},
|
||||
{-0.74,-0.13608395705016335,-0.38167356742177061},
|
||||
{-0.72,-0.14108151215157036,-0.39520508036475621},
|
||||
{-0.70,-0.14623719172590804,-0.40937320507934721},
|
||||
{-0.68,-0.15155666558310329,-0.4242163421817764},
|
||||
{-0.66,-0.15704587372534554,-0.43977570189797999},
|
||||
{-0.64,-0.16271104284734619,-0.45609554660252544},
|
||||
{-0.62,-0.16855870405908679,-0.47322345779107611},
|
||||
{-0.60,-0.17459571193801338,-0.49121063002286869},
|
||||
{-0.58,-0.18082926502846644,-0.5101121950146944},
|
||||
{-0.56,-0.187266927918217,-0.52998757938134999},
|
||||
{-0.54,-0.19391665503547184,-0.5509009000235372},
|
||||
{-0.52,-0.2007868163248025,-0.5729214016617592},
|
||||
{-0.50,-0.20788622497735457,-0.59612394165060356},
|
||||
{-0.48,-0.21522416740965075,-0.62058952789755717},
|
||||
{-0.46,-0.22281043570659706,-0.64640591651133084},
|
||||
{-0.44,-0.23065536276825693,-0.67366827676028285},
|
||||
{-0.42,-0.23876986042695231,-0.70247993198005143},
|
||||
{-0.40,-0.24716546083171484,-0.73295318632721385},
|
||||
{-0.38,-0.25585436143154672,-0.76521024875093424},
|
||||
{-0.36,-0.26484947392794789,-0.7993842672084216},
|
||||
{-0.34,-0.27416447761139873,-0.83562048816081747},
|
||||
{-0.32,-0.28381387754675275,-0.87407755869560466},
|
||||
{-0.30,-0.29381306812972124,-0.91492899127393734},
|
||||
{-0.28,-0.30417840260190093,-0.95836481437693455},
|
||||
{-0.26,-0.31492726918639691,-1.0045934359622013},
|
||||
{-0.24,-0.32607817459151656,-1.0538437511293579},
|
||||
{-0.22,-0.33765083572804089,-1.1063675305908394},
|
||||
{-0.20,-0.34966628059831412,-1.1624421327407444},
|
||||
{-0.18,-0.36214695944532344,-1.2223735894872689},
|
||||
{-0.16,-0.37511686740003325,-1.2865001248178722},
|
||||
{-0.14,-0.38860168003903095,-1.3551961755595476},
|
||||
{-0.12,-0.40262890346626218,-1.4288769964465424},
|
||||
{-0.10,-0.41722804076736686,-1.5080039467563759},
|
||||
{-0.08,-0.43243077695897958,-1.5930905741484096},
|
||||
{-0.06,-0.44827118487571582,-1.6847096335028304},
|
||||
{-0.04,-0.46478595481336826,-1.7835012055832924},
|
||||
{-0.02,-0.48201465118895931,-1.890182113243053},
|
||||
{0.00,-0.5,-2.0055568731891711},
|
||||
{0.02,-0.51878821148283127,-2.1305304708590063},
|
||||
{0.04,-0.53842934310323776,-2.266123307002927},
|
||||
{0.06,-0.55897770888625642,-2.4134887403120628},
|
||||
{0.08,-0.58049234213678902,-2.5739337444577441},
|
||||
{0.10,-0.60303751985624165,-2.7489433156574976},
|
||||
{0.12,-0.62668335867046043,-2.9402094144046673},
|
||||
{0.14,-0.65150649391032422,-3.1496654113984976},
|
||||
{0.16,-0.67759085570376387,-3.3795272436397124},
|
||||
{0.18,-0.70502855864211988,-3.632342787788132},
|
||||
{0.20,-0.73392092489634064,-3.9110513431783414},
|
||||
{0.22,-0.76437966473554453,-4.2190556142460096},
|
||||
{0.24,-0.7965282434431088,-4.5603092252434356},
|
||||
{0.26,-0.83050346989466928,-4.9394236447226962},
|
||||
{0.28,-0.86645734989895151,-5.361799486692064},
|
||||
{0.30,-0.90455925725398401,-5.8337886697644841},
|
||||
{0.32,-0.94499848793067787,-6.3628956591689123},
|
||||
{0.34,-0.98798727865319957,-6.9580293809775497},
|
||||
{0.36,-1.0337643914629273,-7.6298181250132284},
|
||||
{0.38,-1.0825993920687229,-8.3910150599854578},
|
||||
{0.40,-1.1347977838669816,-9.2569895219913114},
|
||||
{0.42,-1.1907072041724314,-10.246454459917672},
|
||||
{0.44,-1.2507249482123559,-11.382048655457382},
|
||||
{0.46,-1.3153071651191928,-12.692443921937331},
|
||||
{0.48,-1.384980176074901,-14.210086389866547},
|
||||
{0.50,-1.4603545088095868,-15.987037203258081},
|
||||
{0.52,-1.5421424407243636,-18.045752498464225},
|
||||
{0.54,-1.6311801184805925,-20.576140424668012},
|
||||
{0.56,-1.728455710203888,-23.218395308861353},
|
||||
{0.58,-1.8351456011744598,-27.764767049040561},
|
||||
{0.60,-1.9526614482240008,-28.111877679505632},
|
||||
{0.62,-2.082712093415199,-47.809694357795301},
|
||||
{0.64,-2.2273861156137764,0.}};
|
||||
static int nzet2=194;
|
||||
static double zet2[194][3]= {{1.88,1.773726366421812,0.},
|
||||
{1.91,1.73817111540579,3.423053940298117},
|
||||
{1.94,1.704997459819862,2.1850871061103181},
|
||||
{1.97,1.673983526341128,2.2347450165483931},
|
||||
{2.00,1.644934066848226,1.9724260665775191},
|
||||
{2.03,1.617676581125839,1.8220425205690245},
|
||||
{2.06,1.592058097974525,1.6660876583057966},
|
||||
{2.09,1.567942487877175,1.532760539302451},
|
||||
{2.12,1.545208207555242,1.4117353539248083},
|
||||
{2.15,1.523746397547932,1.3034334758232795},
|
||||
{2.18,1.50345926998828,1.205747060497862},
|
||||
{2.21,1.484258736211354,1.1175368336960312},
|
||||
{2.24,1.466065233577907,1.0376465612472376},
|
||||
{2.27,1.448806718572342,0.96512777385490744},
|
||||
{2.30,1.432417799315324,0.89914733364618227},
|
||||
{2.33,1.416838985478223,0.83898569101448661},
|
||||
{2.36,1.40201603747119,0.78401543607757262},
|
||||
{2.39,1.387899399906699,0.73368884828547876},
|
||||
{2.42,1.374443706875321,0.68752605820613866},
|
||||
{2.45,1.361607348633357,0.64510551497944379},
|
||||
{2.48,1.349352090988493,0.60605586254545352},
|
||||
{2.51,1.337642740054658,0.57004910836655365},
|
||||
{2.54,1.326446846189378,0.53679482768495013},
|
||||
{2.57,1.315734441872686,0.50603523814671514},
|
||||
{2.60,1.305477809072781,0.47754099830918756},
|
||||
{2.63,1.295651272299547,0.4511076130962095},
|
||||
{2.66,1.286231014096292,0.42655234915648738},
|
||||
{2.69,1.27719491018153,0.40371158023344345},
|
||||
{2.72,1.268522381841718,0.38243849624566623},
|
||||
{2.75,1.260194263504835,0.36260112097532271},
|
||||
{2.78,1.25219268370383,0.34408059237859245},
|
||||
{2.81,1.244500957876449,0.32676966699358972},
|
||||
{2.84,1.237103491650545,0.31057141616997924},
|
||||
{2.87,1.229985693437414,0.29539808681101587},
|
||||
{2.90,1.223133895304355,0.28117010372781304},
|
||||
{2.93,1.216535281225654,0.26781519399862941},
|
||||
{2.96,1.210177821921451,0.25526761693212602},
|
||||
{2.99,1.204050215589313,0.24346748537801852},
|
||||
{3.02,1.198141833915918,0.23236016650498295},
|
||||
{3.05,1.192442672827978,0.22189575163336397},
|
||||
{3.08,1.186943307503945,0.21202858634541449},
|
||||
{3.11,1.181634851222455,0.20271685326768657},
|
||||
{3.14,1.176508917671021,0.19392220095494961},
|
||||
{3.17,1.171557586380136,0.18560941324383382},
|
||||
{3.20,1.166773370984467,0.17774611417765113},
|
||||
{3.23,1.162149190044926,0.17030250422805729},
|
||||
{3.26,1.157678340193674,0.16325112417657045},
|
||||
{3.29,1.153354471389069,0.15656664337357057},
|
||||
{3.32,1.149171564089562,0.15022566965178363},
|
||||
{3.35,1.145123908175106,0.14420657835521838},
|
||||
{3.38,1.141206083461874,0.13848935842254373},
|
||||
{3.41,1.137412941671485,0.13305547357912961},
|
||||
{3.44,1.133739589729561,0.12788773702169684},
|
||||
{3.47,1.130181374280605,0.12297019812194822},
|
||||
{3.50,1.126733867317057,0.11828803987650427},
|
||||
{3.53,1.123392852830047,0.11382748596418971},
|
||||
{3.56,1.12015431439806,0.10957571641301979},
|
||||
{3.59,1.117014423637464,0.10552079099233949},
|
||||
{3.62,1.113969529445854,0.10165157952781294},
|
||||
{3.65,1.111016147975381,0.097957698475736113},
|
||||
{3.68,1.108150953278886,0.094429453090526042},
|
||||
{3.71,1.105370768576718,0.091057784668039821},
|
||||
{3.74,1.102672558096666,0.087834222352239852},
|
||||
{3.77,1.100053419443589,0.084750839084792196},
|
||||
{3.80,1.097510576459004,0.081800211257407274},
|
||||
{3.83,1.095041372534306,0.078975381793226412},
|
||||
{3.86,1.092643264344309,0.076269826250098122},
|
||||
{3.89,1.090313815970592,0.073677421740042881},
|
||||
{3.92,1.08805069338661,0.071192418353802578},
|
||||
{3.95,1.085851659278838,0.068809412910006612},
|
||||
{3.98,1.08371456818028,0.066523324765397865},
|
||||
{4.01,1.081637361894553,0.064329373571598653},
|
||||
{4.04,1.079618065190493,0.06222305872739857},
|
||||
{4.07,1.077654781748776,0.060200140469556893},
|
||||
{4.10,1.075745690343503,0.058256622361922125},
|
||||
{4.13,1.073889041242993,0.056388735160770968},
|
||||
{4.16,1.072083152815211,0.054592921855876442},
|
||||
{4.19,1.070326408324396,0.052865823856315185},
|
||||
{4.22,1.068617252906399,0.051204268171781336},
|
||||
{4.25,1.066954190711215,0.049605255543662775},
|
||||
{4.28,1.065335782201999,0.048065949444368941},
|
||||
{4.31,1.063760641600657,0.04658366583936821},
|
||||
{4.34,1.062227434470793,0.045155863715992485},
|
||||
{4.37,1.060734875429471,0.043780136240254106},
|
||||
{4.40,1.059281725979835,0.042454202569066214},
|
||||
{4.43,1.057866792457208,0.041175900203428546},
|
||||
{4.46,1.05648892408177,0.039943177879255484},
|
||||
{4.49,1.055147011111431,0.038754088939824502},
|
||||
{4.52,1.053839983088912,0.037606785162681358},
|
||||
{4.55,1.052566807177479,0.036499510982131479},
|
||||
{4.58,1.051326486580128,0.03543059812088515},
|
||||
{4.61,1.050118059037379,0.034398460553172631},
|
||||
{4.64,1.048940595399152,0.033401589811060839},
|
||||
{4.67,1.047793198266483,0.032438550586016406},
|
||||
{4.70,1.046675000699131,0.03150797662483848},
|
||||
{4.73,1.045585164985372,0.030608566866521691},
|
||||
{4.76,1.044522881470505,0.029739081858123653},
|
||||
{4.79,1.043487367440835,0.028898340349178375},
|
||||
{4.82,1.042477866060075,0.028085216148126278},
|
||||
{4.85,1.041493645355327,0.027298635128770225},
|
||||
{4.88,1.040533997249946,0.026537572456617805},
|
||||
{4.91,1.039598236640801,0.025801049951032575},
|
||||
{4.94,1.038685700517542,0.025088133644345388},
|
||||
{4.97,1.037795747121678,0.02439793144084252},
|
||||
{5.00,1.03692775514337,0.023729590964713193},
|
||||
{5.03,1.036081122953981,0.023082297496121333},
|
||||
{5.06,1.035255267872542,0.022455272045522098},
|
||||
{5.09,1.034449625464386,0.021847769542722075},
|
||||
{5.12,1.03366364887033,0.021259077120328922},
|
||||
{5.15,1.032896808164854,0.020688512501460612},
|
||||
{5.18,1.032148589741818,0.020135422481716021},
|
||||
{5.21,1.03141849572637,0.019599181489412698},
|
||||
{5.24,1.030706043411723,0.019079190236088738},
|
||||
{5.27,1.030010764719606,0.01857487443023078},
|
||||
{5.30,1.029332205683219,0.018085683577360203},
|
||||
{5.33,1.028669925951618,0.017611089832060688},
|
||||
{5.36,1.028023498314491,0.017150586919361627},
|
||||
{5.39,1.027392508246357,0.016703689117205631},
|
||||
{5.42,1.026776553469275,0.016269930284251142},
|
||||
{5.45,1.026175243533173,0.015848862952556192},
|
||||
{5.48,1.025588199413004,0.015440057457382975},
|
||||
{5.51,1.02501505312192,0.015043101114149762},
|
||||
{5.54,1.024455447339741,0.014657597451369339},
|
||||
{5.57,1.023909035056019,0.01428316546383793},
|
||||
{5.60,1.02337547922703,0.013919438914884102},
|
||||
{5.63,1.022854452446063,0.013566065688307632},
|
||||
{5.66,1.022345636626416,0.013222707131014493},
|
||||
{5.69,1.021848722696524,0.012889037489672295},
|
||||
{5.72,1.021363410306693,0.012564743315914253},
|
||||
{5.75,1.020889407546917,0.012249522946180057},
|
||||
{5.78,1.020426430675303,0.011943085978162262},
|
||||
{5.81,1.019974203856636,0.011645152791711841},
|
||||
{5.84,1.019532458910656,0.011355454098314737},
|
||||
{5.87,1.019100935069616,0.011073730417842379},
|
||||
{5.90,1.01867937874474,0.010799731988834502},
|
||||
{5.93,1.018267543301195,0.010533217160810536},
|
||||
{5.96,1.017865188841217,0.010273956488077638},
|
||||
{5.99,1.017472081995063,0.010021715712312041},
|
||||
{6.02,1.017087995719444,0.0097763175668983369},
|
||||
{6.05,1.016712709103148,0.0095374095008413651},
|
||||
{6.08,1.01634600717954,0.009305329021187813},
|
||||
{6.11,1.015987680745679,0.0090778727223594924},
|
||||
{6.17,1.015295345312643,0.008642775132504197},
|
||||
{6.23,1.014634137971395,0.0082311797290497944},
|
||||
{6.29,1.014002575548372,0.0078407029934879115},
|
||||
{6.35,1.013399251755712,0.007470392233756039},
|
||||
{6.41,1.012822832736173,0.0071190166067899496},
|
||||
{6.47,1.012272052902078,0.006785517079187957},
|
||||
{6.53,1.011745711046087,0.0064688785842991776},
|
||||
{6.59,1.011242666703518,0.0061681576191670367},
|
||||
{6.65,1.010761836747634,0.00588246874707092},
|
||||
{6.71,1.010302192200886,0.0056109826200994037},
|
||||
{6.77,1.009862755246497,0.0053529213704776235},
|
||||
{6.83,1.009442596426046,0.0051075551277877497},
|
||||
{6.89,1.00904083200989,0.0048741986101521152},
|
||||
{6.95,1.008656621528302,0.0046522080444542681},
|
||||
{7.01,1.00828916545217,0.0044409783070365281},
|
||||
{7.07,1.007937703012983,0.004239940300261584},
|
||||
{7.13,1.007601510152612,0.0040485585193739447},
|
||||
{7.19,1.007279897594155,0.0038663288130501468},
|
||||
{7.25,1.006972209025747,0.0036927763089049164},
|
||||
{7.31,1.006677819389865,0.0035274534952119026},
|
||||
{7.37,1.006396133271223,0.0033699384441953448},
|
||||
{7.43,1.006126583376844,0.0032198331650262341},
|
||||
{7.49,1.005868629102374,0.0030767620794814434},
|
||||
{7.55,1.005621755179156,0.0029403706012505774},
|
||||
{7.61,1.005385470396922,0.0028103238230340844},
|
||||
{7.67,1.005159306397399,0.0026863052927333189},
|
||||
{7.73,1.004942816534398,0.0025680158747019366},
|
||||
{7.79,1.004735574796287,0.002455172693089723},
|
||||
{7.85,1.004537174787052,0.0023475081443244251},
|
||||
{7.91,1.004347228762366,0.0022447689799223172},
|
||||
{7.97,1.004165366717389,0.0021467154490857712},
|
||||
{8.03,1.003991235523176,0.002053120498915751},
|
||||
{8.09,1.003824498108849,0.001963769030275603},
|
||||
{8.15,1.003664832686812,0.0018784571981245541},
|
||||
{8.21,1.003511932018526,0.00179699176176089},
|
||||
{8.27,1.003365502718472,0.0017191894751160557},
|
||||
{8.33,1.003225264594127,0.0016448765179545447},
|
||||
{8.39,1.003090950019878,0.001573887948228554},
|
||||
{8.45,1.002962303342977,0.0015060672678566414},
|
||||
{8.51,1.002839080319702,0.0014412656910764718},
|
||||
{8.57,1.002721047580067,0.0013793427004115579},
|
||||
{8.63,1.002607982119472,0.0013201619077992217},
|
||||
{8.69,1.002499670815825,0.0012636045819143255},
|
||||
{8.75,1.00239590997073,0.0012095173500489579},
|
||||
{8.81,1.00229650487343,0.0011579056771549652},
|
||||
{8.87,1.002201269386287,0.0011082102036054395},
|
||||
{8.93,1.002110025550626,0.0010620059777144341},
|
||||
{8.99,1.00202260321186,0.0010129273774552136},
|
||||
{9.05,1.001938839662877,0.00098426748400331446},
|
||||
{9.11,1.001858579304718,0.00088865406038856785},
|
||||
{9.17,1.001781673323641,0.0010517447435007122},
|
||||
{9.20,1.00174443349959,0.}};
|
||||
static int nzet3=54;
|
||||
static double zet3[54][3]= {{0.90,0.009089522010492318,0.},
|
||||
{1.00,0.0072328000492769655,0.049940424539145614},
|
||||
{1.10,0.0057532294481890678,0.02652911791989138},
|
||||
{1.20,0.0045749957805488939,0.024745263849921362},
|
||||
{1.30,0.0036372241858115932,0.018767070422148646},
|
||||
{1.40,0.0028911505710517487,0.015205242447956671},
|
||||
{1.50,0.002297783305219789,0.012035769142756404},
|
||||
{1.60,0.0018259888531061044,0.0095953692119826774},
|
||||
{1.70,0.0014509355190296495,0.0076274248316501999},
|
||||
{1.80,0.0011528349668758395,0.006066600615003909},
|
||||
{1.90,0.00091592859251588747,0.0048226793846485897},
|
||||
{2.00,0.00072767359533962967,0.0038335081566185855},
|
||||
{2.10,0.00057809093856494198,0.0030466922298190919},
|
||||
{2.20,0.00045924394904117114,0.0024211232746551153},
|
||||
{2.30,0.00036482197313021935,0.0019238228392516371},
|
||||
{2.40,0.00028980828964268691,0.0015285608223902069},
|
||||
{2.50,0.00023021545884836824,0.0012144454871157791},
|
||||
{2.60,0.00018287455781888589,0.00096481508804845465},
|
||||
{2.70,0.00014526742069331959,0.00076655250304002586},
|
||||
{2.80,0.00011539316590233686,0.00060870430054149594},
|
||||
{2.90,9.1662038482603366e-05,0.00048450671754362266},
|
||||
{3.00,7.2810999959871433e-05,0.00038132216748494168},
|
||||
{3.20,4.5941763514956608e-05,0.00023870622885571716},
|
||||
{3.40,2.8987742285838751e-05,0.00015113519946172838},
|
||||
{3.60,1.8290207711218413e-05,9.5225971472001988e-05},
|
||||
{3.80,1.154041207421617e-05,6.0121755292975294e-05},
|
||||
{4.00,7.2815360962094121e-06,3.7924956205421726e-05},
|
||||
{4.20,4.5943499611447305e-06,2.3931896326648792e-05},
|
||||
{4.40,2.8988433425797932e-06,1.5099385962944477e-05},
|
||||
{4.60,1.8290482855851182e-06,9.5272940571119381e-06},
|
||||
{4.80,1.1540521610478418e-06,6.0112776772183132e-06},
|
||||
{5.00,7.281579703182653e-07,3.7928853051697142e-06},
|
||||
{5.20,4.5943673213436779e-07,2.3931239839546378e-06},
|
||||
{5.40,2.8988502537869706e-07,1.5100484732457307e-06},
|
||||
{5.60,1.8290510369831115e-07,9.524498843550906e-07},
|
||||
{5.80,1.1540532563983228e-07,6.0217353262003908e-07},
|
||||
{6.00,7.2815840638494783e-08,3.7539994373595081e-07},
|
||||
{6.30,3.6494370945140209e-08,1.8473501971024426e-07},
|
||||
{6.60,1.8290513121224273e-08,9.3500768718982493e-08},
|
||||
{6.90,9.1669717492768626e-09,4.6616335545062541e-08},
|
||||
{7.20,4.5943692309570581e-09,2.342981267593911e-08},
|
||||
{7.50,2.3026392099478531e-09,1.1722580238554488e-08},
|
||||
{7.80,1.1540533768867409e-09,5.8894788997158704e-09},
|
||||
{8.10,5.7839681972433946e-10,2.914789222496127e-09},
|
||||
{8.50,2.3026392143085165e-10,1.1089727185909566e-09},
|
||||
{9.00,7.281584547883154e-11,3.4109996103179389e-10},
|
||||
{9.50,2.3026392147445831e-11,1.1043438017709408e-10},
|
||||
{10.00,7.2815845483192207e-12,3.4234015834048123e-11},
|
||||
{10.50,2.3026392147881895e-12,1.101025086100738e-11},
|
||||
{11.00,7.2815845483628272e-13,3.4321304878213432e-12},
|
||||
{11.50,2.3026392147925503e-13,1.0992966259843476e-12},
|
||||
{12.00,7.2815845483671885e-14,3.4139798491593627e-13},
|
||||
{12.50,2.3026392147929864e-14,1.1891837818809462e-13},
|
||||
{13.00,7.2815845483676246e-15,0.}};
|
||||
static int nzet4= 53;
|
||||
static double zet4[53][3] = {{0.90,-0.009243055244156485,0.},
|
||||
{1.00,-0.0073296843037852326,-0.052814323525202628},
|
||||
{1.10,-0.0058143637920737854,-0.027572963095073559},
|
||||
{1.20,-0.0046135707400949254,-0.025610299934053761},
|
||||
{1.30,-0.0036615640549682469,-0.019257657280021753},
|
||||
{1.40,-0.0029065082749557885,-0.015529614014390146},
|
||||
{1.50,-0.0023074734746324856,-0.012236474475911939},
|
||||
{1.60,-0.001832102981787953,-0.0097230725692241081},
|
||||
{1.70,-0.0014547932913898482,-0.0077077167150478385},
|
||||
{1.80,-0.0011552690638311632,-0.0061173382742368398},
|
||||
{1.90,-0.0009174644067140214,-0.0048546724529303833},
|
||||
{2.00,-0.00072864262972110918,-0.0038536999885796911},
|
||||
{2.10,-0.00057870235837620897,-0.0030594309815580099},
|
||||
{2.20,-0.00045962972904128429,-0.0024291612911735549},
|
||||
{2.30,-0.00036506538392615766,-0.0019288943856263516},
|
||||
{2.40,-0.00028996187150057997,-0.0015317607800506593},
|
||||
{2.50,-0.00023031236246091888,-0.0012164645257209559},
|
||||
{2.60,-0.00018293569986944355,-0.00096608898597696396},
|
||||
{2.70,-0.00014530599872097187,-0.00076735639617336448},
|
||||
{2.80,-0.00011541750699293475,-0.00060921108159023578},
|
||||
{2.90,-9.1677396672773469e-05,-0.00048482812219128224},
|
||||
{3.00,-7.2820690322870533e-05,-0.00038151881179963415},
|
||||
{3.20,-4.5945621317998246e-05,-0.00023878105774549144},
|
||||
{3.40,-2.8989278104899569e-05,-0.00015116582598443502},
|
||||
{3.60,-1.8290819131800366e-05,-9.5238274316694453e-05},
|
||||
{3.80,-1.1540655485134253e-05,-6.0125375713747831e-05},
|
||||
{4.00,-7.2816329998413212e-06,-3.7931397034293605e-05},
|
||||
{4.20,-4.5943885391754506e-06,-2.3915739843136425e-05},
|
||||
{4.40,-2.8988587007704493e-06,-1.5162836932290832e-05},
|
||||
{4.60,-1.8290543997909454e-06,-9.2917430415241644e-06},
|
||||
{4.80,-1.1540545951570238e-06,-6.8908653534506024e-06},
|
||||
{5.00,-7.2815893935458461e-07,-5.1041786939570726e-07},
|
||||
{5.20,-4.59437117914675e-07,-1.4643538323345976e-05},
|
||||
{5.40,-2.8988517896060359e-07,4.4209088789903913e-05},
|
||||
{5.60,-1.8290516484036941e-07,-0.00017157860556134556},
|
||||
{5.80,-1.154053499809241e-07,0.00063618330356635848},
|
||||
{6.10,-5.7839685060107717e-08,-0.0020091375453381894},
|
||||
{6.40,-2.8988510984876702e-08,0.0073984525783966845},
|
||||
{6.70,-1.4528671435507856e-08,-0.02758563219055031},
|
||||
{7.00,-7.2815845968194883e-09,0.10294359533362384},
|
||||
{7.30,-3.6494372283900303e-09,-0.3841889901399097},
|
||||
{7.60,-1.8290513457505598e-09,1.4338122444419159},
|
||||
{8.00,-7.2815845532128537e-10,-4.7302011626775808},
|
||||
{8.50,-2.3026392152775532,-39.381666792126495},
|
||||
{9.00,-7.2815845488521909,98.025521474577729},
|
||||
{9.50,-2.3026392148414865,-113.73104308413625},
|
||||
{10.00,-7.2815845484161246,117.90927483991902},
|
||||
{10.50,-2.30263921479788,-118.91668026291069},
|
||||
{11.00,-7.2815845483725177,118.76807019909455},
|
||||
{11.50,-2.3026392147935195,-117.16622452178022},
|
||||
{12.00,-7.2815845483681567,110.90745187633907},
|
||||
{12.50,-2.3026392147930834,-87.474206971983023},
|
||||
{13.00,-7.2815845483677206,0.}};
|
||||
double hlp;
|
||||
|
||||
if(x>0.96 && x<1.04 && x!=1.0) /*use laurent series -1 and 0 term and spline */
|
||||
{
|
||||
double l,e;
|
||||
l= -log10(fabs(x-1.0));
|
||||
if(l>=6.5) std::cerr <<"Argument of zeta too near to the pole, result will be very imprecise!\n";
|
||||
if(l>12.0) e=0.0;
|
||||
else
|
||||
{
|
||||
if(x>1.0) splintx(zet3,nzet3,l,&e);
|
||||
else splintx(zet4,nzet4,l,&e);
|
||||
}
|
||||
return(c_euler+e+1.0/(x-1.0));
|
||||
}
|
||||
|
||||
if(x==floor(x) && fabs(x)<=(((unsigned int) 1)<<(4*sizeof(int)-1)) )
|
||||
{
|
||||
int n;
|
||||
n=(int)x;
|
||||
if(n==1) laerror("singular point in zeta");
|
||||
if(n==0) return(-0.5);
|
||||
if(n<0)
|
||||
{
|
||||
if(n&1) return(-bernoulli_number(1-n)/(1-n));
|
||||
else return(0.0);
|
||||
}
|
||||
if(n==2||n==4||n==6||n==8) /*precalculated few bernoulli numbers */
|
||||
return(fabs(bernoulli_number(n))*pow(2*M_PI,(double)n)/factorial(n)/2.0);
|
||||
/*otherwise continue like if general real argument*/
|
||||
}
|
||||
if(x<-1.02 || (x>0.58 && x<1.94)) return (2.0*pow(2.0*M_PI,x-1.0)*sin(M_PI*x/2.0)*gamma(1.0-x)*zeta(1.0-x));
|
||||
|
||||
|
||||
/* make different approximations in different intervals - some small non-smoothness is introduced */
|
||||
if(x>=9.1) /* in this region the summation per definition is best and quite exact*/
|
||||
{
|
||||
double s;
|
||||
int i,n;
|
||||
n=(int)ceil(pow(EPS/1000.0,-1.0/x));
|
||||
s=0.0;for(i=n;i>1;i--) s+= pow((double)i,-x);
|
||||
return(1.0+s);
|
||||
}
|
||||
|
||||
if(x>=1.94) { splintx(zet2,nzet2,x,&hlp); return(hlp);}
|
||||
splintx(zet1,nzet1,x,&hlp); return(hlp);
|
||||
}
|
||||
//end zeta
|
||||
|
||||
|
||||
|
||||
|
||||
double gamma(double x)
|
||||
{
|
||||
bool f;
|
||||
double p;
|
||||
|
||||
if(x>MAXFACT) laerror("too big argument in gamma");
|
||||
f= (fabs(x-floor(x+0.5))<EPS);
|
||||
if(x<=0.0 && f) laerror("nonpositive integer argument in gamma");
|
||||
if(f) return factorial(floor(x+0.5)-1);
|
||||
|
||||
if(x<8.0)
|
||||
{
|
||||
p=1.0;
|
||||
while(x<8.0) {p /= x; x += 1.0;}
|
||||
return(p*exp(stirfact(x-1.0)));
|
||||
}
|
||||
return(exp(stirfact(x-1.0)));
|
||||
}
|
||||
|
||||
|
||||
static double factln(double nn)
|
||||
{
|
||||
if(nn<=LIMITFACT+0.5) return factorialln((int)(nn+.5));
|
||||
return stirfact(nn);
|
||||
}
|
||||
|
||||
|
||||
double gammln(double x)
|
||||
{
|
||||
bool f;
|
||||
|
||||
f= (fabs(x-floor(x+0.5))<EPS);
|
||||
if(x<=0.0 && f) laerror("nonpositive integer argument in gamma");
|
||||
if(fabs(x-floor(x+0.5))<EPS) return factln(x-1.0);
|
||||
if(x<8.0)
|
||||
{
|
||||
double p;
|
||||
p=0.0;
|
||||
while(x<8.0) {p -= log(fabs(x)); x += 1.0;}
|
||||
return(p+stirfact(x-1.0));
|
||||
}
|
||||
return(stirfact(x-1.0));
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/* Bernoulli numbers */
|
||||
#define MAXBERN 50
|
||||
double bernoulli_number(int n)
|
||||
{
|
||||
static double bern[MAXBERN]={1.0,0.166666666666666666,-0.03333333333333333333,0.0238095238095238,-0.033333333333333333,0.07575757575757575757};
|
||||
static int top=5;
|
||||
int i;
|
||||
|
||||
if(n<0) laerror("negative index of Bernoulli number");
|
||||
if(n==1) return(-0.5);
|
||||
if(n&1) return(0.0);
|
||||
n /= 2;
|
||||
if( n >= MAXBERN) /* use zeta function for it */
|
||||
{int m;
|
||||
m=2*n;
|
||||
return(zeta((double)m)*(n&1?2.0:-2.0)*factorial(m)/pow(2.0*M_PI,(double)m));
|
||||
}
|
||||
if(n<=top) return(bern[n]);
|
||||
|
||||
/* this is very unexact - recurrent formula with binom. coef.
|
||||
for(i=top+1;i<=n;i++)
|
||||
{
|
||||
int j,d;
|
||||
double s;
|
||||
d=2*i+1;
|
||||
s=0.0;
|
||||
for(j=1;j<i;j++) s += bern[j]*ibinom(d,2*j);
|
||||
bern[i]=0.5-(1.0+s)/d;
|
||||
}
|
||||
*/
|
||||
|
||||
/* precalculate and save for next use, because the calc. is expensive */
|
||||
/* could be calculated just one at need, but let's do all of them */
|
||||
for(i=top+1;i<=n;i++)
|
||||
{int m;
|
||||
m=2*i;
|
||||
bern[i]=zeta((double)m)*(i&1?2.0:-2.0)*factorial(m)/pow(2.0*M_PI,(double)m);
|
||||
}
|
||||
top=n;
|
||||
return(bern[n]);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
}//namespace
|
||||
|
||||
@ -25,10 +25,17 @@
|
||||
namespace LA {
|
||||
|
||||
|
||||
const double c_euler= 0.5772156649015328606065120900824024310422;
|
||||
|
||||
extern unsigned long long factorial(const int n);
|
||||
extern double factorialln(const int n);
|
||||
extern unsigned long long binom(int n, int k);
|
||||
extern unsigned long long longpow(unsigned long long x, int i);
|
||||
extern double zeta(double x);
|
||||
extern double gamma(double x);
|
||||
extern double gammln(double x);
|
||||
extern double bernoulli_number(int n);
|
||||
|
||||
|
||||
|
||||
}//namespace
|
||||
|
||||
Loading…
Reference in New Issue
Block a user