/*
    LA: linear algebra C++ interface library
    Copyright (C) 2024 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include <iostream>
#include <stdint.h>
#include <math.h>

#include "miscfunc.h"
#include "laerror.h"
#include "polynomial.h"

namespace LA {


#define MAXFACT 20
unsigned long long factorial(const int n)
{
static int ntop=20;
static  unsigned long long a[MAXFACT+1]={1,1,2,6,24,120,720,5040,
40320,
362880,
3628800,
39916800ULL,
479001600ULL,
6227020800ULL,
87178291200ULL,
1307674368000ULL,
20922789888000ULL,
355687428096000ULL,
6402373705728000ULL,
121645100408832000ULL,
2432902008176640000ULL};

int j;
if (n < 0) laerror("negative argument of factorial");
if (n > MAXFACT) laerror("overflow in factorial");
        while (ntop<n) {
                j=ntop++;
                a[ntop]=a[j]*ntop;
                }
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= std::log(z);
if(z>=2.0) /*does not have sense for smaller anyway and could make troubles*/
        {
        int i,ii,m;
        m=(int)std::ceil(0.5-std::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=std::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]= std::log(factorial(n)));

                {
                int j;
                if(!a[top]) a[top]= std::log(factorial(top));
                while (top<n) {
                        j=top++;
                        a[top]=a[j]+std::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)))
unsigned long long binom(int n, int k)
{
unsigned long long p,value;
int d;
static unsigned long long ibitab[ibidxmaxeven(MAXBINOM)]= /*only nontrivial are stored,
                                            zero initialization by compiler assumed*/
{
6,
10,
15,20,
21,35,
28,56,70
};


if(k>n||k<0) return(0);
if(k>n/2) k=n-k;
if(k==0) return(1);
if(k==1) return(n);
int ind=0;
if(n<=MAXBINOM)
        {
        ind=ibidx(n,k);
        if (ibitab[ind]) return ibitab[ind];
        }
/* nonrecurent method used anyway */
d=n-k;
p=1;
for(;n>d;n--) p *= n;
value=p/factorial(k);
if(n<=MAXBINOM) ibitab[ind]=value;
return value;
}
#undef ibidx
#undef ibidxmaxeven
#undef MAXBINOM


unsigned long long longpow(unsigned long long x, int i)
{
if(i<0) return 0;
unsigned long long y=1;
do
        {
        if(i&1) y *= x;
	i >>= 1;
        if(i) x *= x;
        }
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= -std::log10(std::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==std::floor(x) && std::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(std::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)std::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= (std::fabs(x-std::floor(x+0.5))<EPS);
if(x<=0.0 && f) laerror("nonpositive integer argument in gamma");
if(f) return factorial(std::floor(x+0.5)-1);

if(x<8.0)
        {
        p=1.0;
        while(x<8.0) {p /= x; x += 1.0;}
        return(p* std::exp(stirfact(x-1.0)));
        }
return( std::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= (std::fabs(x-std::floor(x+0.5))<EPS);
if(x<=0.0 && f) laerror("nonpositive integer argument in gamma");
if(std::fabs(x-std::floor(x+0.5))<EPS) return factln(x-1.0);
if(x<8.0)
        {
        double p;
        p=0.0;  
        while(x<8.0) {p -= std::log(std::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]);
}



//there is a closed form s(d,n) = binom(n+d-1,d)
//we could also pre-generate these polynomials for small d
//but that would be inefficient, better is to buffer all values
//we use recurrence s(d,n) = s(d,n-1) + s(d-1,n)
//we cannot also use the above efficient binom() since here we need laerge n while d is small
//
//enlarge storage size if needed
//
#define SIMPLICIAL_MAXD 8
#define SIMPLICIAL_MAXN 1024

unsigned long long simplicial(int d, int n)
{
static unsigned long long stored[SIMPLICIAL_MAXD-2][SIMPLICIAL_MAXN]={{0,1,3}};
static int maxd=2;
static int maxn=2;

if(n<=0) return 0;
if(d==0||n==1) return 1;
if(d==1) return n;

//resize precomputed part as needed
if(n>maxn)
	{
	if(n>=SIMPLICIAL_MAXN) laerror("storage overflow in simplicial");
	int newdim=2*n; 
	if(newdim<2*maxn) newdim=2*maxn;
	if(newdim>=SIMPLICIAL_MAXN) newdim=SIMPLICIAL_MAXN-1;
	for(int i=maxn+1; i<=newdim;++i) stored[0][i] = i*(i+1)/2;
	for(int dd=1; dd<=maxd-2; ++dd)
		{
		for(int i=maxn+1; i<=newdim;++i) stored[dd][i] =  stored[dd][i-1] + stored[dd-1][i];
		}
	maxn=newdim;
	}

if(d>maxd)
	{
	if(d>=SIMPLICIAL_MAXD) laerror("storage overflow in simplicial");
	for(int dd=maxd+1; dd<=d; ++dd)
		{
		stored[dd-2][0]=0;
		stored[dd-2][1]=1;
		for(int i=2; i<=maxn; ++i) stored[dd-2][i] = stored[dd-3][i] + stored[dd-2][i-1];
		}
	maxd=d;
	}

return stored[d-2][n];
}


//find largest n such that simplicial(d,n)<=s
int inverse_simplicial(int d,  unsigned long long s)
{
if(s==0) return 0;
if(d==0 || s==1) return 1;
if(d==1) return (int)s;
if(d>=SIMPLICIAL_MAXD) laerror("storage overflow in inverse_simplicial");

static int maxd=0;
static Polynomial<double> polynomials[SIMPLICIAL_MAXD];
static Polynomial<double> polynomials_der[SIMPLICIAL_MAXD];
for(int i=maxd+1; i<=d; ++i) //generate as needed
	{
	NRVec<double> roots(i);
	for(int j=0; j<i;++j) roots[j] = -j;
	polynomials[i-1] = polyfromroots(roots) * (1./factorial(i));
	polynomials_der[i-1] = polynomials[i-1].derivative(1);
	maxd=i;
	}

//solve polynomial equation by newton
double x=std::pow(factorial(d)*(double)s,1./d); //good init guess
double dx;
do 
	{
	double residual = value(polynomials[d-1],x)-s;
	dx = residual/value(polynomials_der[d-1],x);
	x -= dx;
	}
while(std::fabs(dx)>.5); //so usually 2 iterations are enough
int n=  std::floor(x);
while(simplicial(d,n)>s) --n;
return n;
}

#undef SIMPLICIAL_MAXD
#undef SIMPLICIAL_MAXN


}//namespace