diff --git a/miscfunc.cc b/miscfunc.cc index 8ee3e8a..9a735c8 100644 --- a/miscfunc.cc +++ b/miscfunc.cc @@ -18,8 +18,9 @@ #include "miscfunc.h" #include "laerror.h" -#include +#include #include +#include 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 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))= 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