/*------------------------------------------------------------------------------- * VERSION: Sept. 23th, 2008. *------------------------------------------------------------------------------ * * Code to compute the Riemann zeta function on the critical line using * band limited interpolation based ideas. * * Note: if run into memory problems, reduce blfi_block_growth * -----------------------------------------------------------------------------*/ //#include //#include //#include //#include #include "L.h" #include "Lriemannsiegel_blfi.h" using namespace std; //typedef complex Complex; //typedef double Double; //typedef int int; const int blfi=1; //-----------------Global constants----------------------------------------------------------------------------- //const Double Pi=3.14159265358979323846264338327950288; //const Complex I=Complex(0.,1.); //const Double bernoulli[] ={1.0, -1./2, 1./6, 0, -1./30, 0, 1./42, 0, -1./30, 0, 5./66, 0, -691./2730, 0, 7./6, 0, -3617./510, 0, 43867./798, 0, -174611./330, 0, 854513./138, 0, -236364091./2730, 0, 8553103./6, 0, -23749461029./870, 0, 8615841276005./14322, 0, -7709321041217./510, 0, 2577687858367./6, 0, -26315271553053477373./1919190, 0, 2929993913841559./6, 0, -261082718496449122051./13530, 0, 1520097643918070802691./1806, 0, -27833269579301024235023./690, 0, 596451111593912163277961./282, 0, -5609403368997817686249127547./46410, 0}; /* -------moved to Lglobals.cc and .h const Double sin_cof[]={1.,-1./6.,1./120.,-1./5040.,1./362880.,-1./39916800.}; const Double sinh_mult_fac=.5; const Double sin_tol=1.e-5; const int sin_terms=2; const Double blfi_block_growth=2.; // how fast blfi blocks grow as we traverse the main sum, keep as is for now const Double beta_fac_mult=6.; // controls the density of blfi sampling and the number of blfi terms needed const Double blfi_fac=.085; // efficiency of the blfi interpolation sum relative to an RS sum of same length const Double pts_array_fac=10000; const int rs_blfi_N=3; //-----------------Global arrays----------------------------------------------------------------------------- Double *klog0; //log(k) at the beginning Double *klog2; //log(k) at the end if needed Double *ksqrt0; // 1/sqrt(k) at the beginning Double *ksqrt2;// 1/sqrt(k) at the end if needed int *num_blocks; // number of blocks int *size_blocks;// size of blocks Double *trig; // stores correction terms Double *zz; // store powers of fractional part Double **klog1; //log(k) in the middle if needed Double **ksqrt1; // 1/sqrt(k) in the middle if needed Double **klog_blfi; //initial term Double **qlog_blfi; //band-width Double **piv_org; //original pivot Double **bbeta; //beta Double **blambda; //lambda Double **bepsilon; //epsilon Double **arg_blfi; //arg_blfi Double **inv_arg_blfi; //inv_arg_blfi Double ***qlog_blfi_dense; // log(1+k/v) terms Double ***qsqrt_blfi_dense; // 1/sqrt(1+k/v) int ***blfi_done_left; //block done or not int ***blfi_done_right; //block done or not Double ***blfi_val_re_left; //real value of block Double ***blfi_val_re_right; //real value of block Double ***blfi_val_im_left; //imag value of block Double ***blfi_val_im_right; //imag value of block //-----------------Global variables----------------------------------------------------------------------------- int length_org=0; // length of the main sum int length_split=0; // length of the portion of the main sum to be evaluated directly int lgdiv=0; // number of divisions of the main sum into intervals of the form [N,2N) int max_pts=0; // max number of interpolation points allowed int range=0; // number of blfi interpolation points needed int blfi_block_size_org=0; // starting length of the blfi block int total_blocks=0; int rs_flag=0; Double bc=0; Double bc2=0.; Double kernel_fac=0.; Double ler=0.; Double mult_fac=0.; Double approx_blfi_mean_spacing=0.; Double interval_length=0.; Double error_tolerance=0.; Double input_mean_spacing=0.; ------------------------------ global variables moved */ //compute sinc(x)=sin(x)/x Double sinc(Double u){ Double ans=1; if(fabs(u) > sin_tol) ans=sin(u)/u; if(fabs(u)<= sin_tol){ Double u2=pow(u,2),temp=u2; for(int j=1;j sin_tol){ Double temp0=exp(u); ans=(temp0-1/temp0)/u; } if(fabs(u)<= sin_tol){ Double u2=pow(u,2),temp=u2; for(int j=1;j=1.) frac=0; z=frac-(Double)1/2; zz[0]=1; for(int i=1;i<51;i++){ zz[i]=z*zz[i-1]; } trig[0]= .3826834323650897717284599840304*zz[0] +1.74896187231008179744118586948533*zz[2] +2.11802520768549637318456427826417*zz[4] -.87072166705114807391892407738239*zz[6] -3.4733112243465167073064116693758*zz[8] -1.66269473089993244964313630119286*zz[10] +1.21673128891923213447689352804417*zz[12] +1.30143041610079757730060538099786*zz[14] +.03051102182736167242108987123981*zz[16] -.37558030515450952427981932122934*zz[18] -.10857844165640659743546975901329*zz[20] +.05183290299954962337576051067322*zz[22] +.02999948061990227592040084956912*zz[24] -.00227593967061256422601994851021*zz[26] -.00438264741658033830594007013585*zz[28] -.00040642301837298469930723272116*zz[30] +.00040060977854221139278910314608*zz[32] +8.97105799138884129783418195378689e-05*zz[34] -2.30256500272391071161029452573899e-05*zz[36] -9.38000660190679248471972940127474e-06*zz[38] +6.32351494760910750424986123959430e-07*zz[40] +6.55102281923150166621223123133411e-07*zz[42] +2.21052374555269725866086890381876e-08*zz[44] -3.32231617644562883503133517017624e-08*zz[46] -3.73491098993365608176460476015222e-09*zz[48] +1.24450670607977391951510000249366e-09*zz[50]; trig[1]= -.05365020525675069405998280791133*zz[1] +.11027818741081482439896362071917*zz[3] +1.23172001543152263131956529162206*zz[5] +1.26349648627994578841755482191213*zz[7] -1.69510899755950301844944739906596*zz[9] -2.99987119676501008895548735894141*zz[11] -.10819944959899208642692257787438*zz[13] +1.94076629462127126879387632539716*zz[15] +.78384235615006865328843457488694*zz[17] -.50548296679003659187902141326173*zz[19] -.3845072349605797405134273885311*zz[21] +.03747264646531532067594447494023*zz[23] +.09092026610973176317258142450576*zz[25] +.01044923755006450921820113972659*zz[27] -.01258297965158341649747892224592*zz[29] -.00339950372115127408505894886137*zz[31] +.00104109505377148912682954240655*zz[33] +.00050109490511184868603556526727*zz[35] -3.95635966900318155954711855696337e-05*zz[37] -4.76245924535718963865409830268035e-05*zz[39] -1.85393553380851322734349064569117e-06*zz[41] +3.19369180800689720404663539343268e-06*zz[43] +4.09078076085060663265089453677018e-07*zz[45] -1.54466243325766321284375723273104e-07*zz[47] -3.46630749176913317222559405934073e-08*zz[49]; trig[2]= .00518854283029316849378458151923*zz[0] +.00123786335522538984133826974438*zz[2] -.18137505725166997411491896409414*zz[4] +.14291492748532126541165603376514*zz[6] +1.33033917666875653250993329998546*zz[8] +.35224723534037336775327655505836*zz[10] -2.4210015958919507237815305433405*zz[12] -1.67607870225381088533346181492372*zz[14] +1.36894167233283721842349153807076*zz[16] +1.55390194302229832214563952655935*zz[18] -.17221642734729980519582586998918*zz[20] -.63590680550454309889704902355845*zz[22] -.09911649873041208105423564341370*zz[24] +.14033480067387008950738254898316*zz[26] +.04782352019827292236438803506512*zz[28] -.01735604064147978079795864709223*zz[30] -.01022501253402859184447660413126*zz[32] +.00092741491597948878994270014371*zz[34] +.00135721943723733853452533619958*zz[36] +6.41369012029388008996238736394533e-05*zz[38] -.00012300805698196629883342322937*zz[40] -1.83135074047892025547675543979621e-05*zz[42] +7.82162860432262730850139938461872e-06*zz[44] +2.00875424847599455034985293919157e-06*zz[46] -3.35327653931857137372749727241453e-07*zz[48] -1.46160209174182309264510097122760e-07*zz[50]; trig[3]= -.00267943218143891380853967145989*zz[1] +.02995372109103514963731329491570*zz[3] -.04257017254182869798501935111688*zz[5] -.28997965779803887506893209478669*zz[7] +.48888319992354459725374746407169*zz[9] +1.23085587639574608119312504336294*zz[11] -.82975607085274087041796910432976*zz[13] -2.24976353666656686652045012659903*zz[15] +.07845139961005471379365473620184*zz[17] +1.74674928008688940039198666645219*zz[19] +.45968080979749935109237306173169*zz[21] -.66193534710397749464339040008983*zz[23] -.31590441036173634578979632973316*zz[25] +.12844792545207495988511847476209*zz[27] +.10073382716626152300969450207513*zz[29] -.00953018384882526775950465984230*zz[31] -.01926442168751408889840098069714*zz[33] -.00124646371587692917124790716458*zz[35] +.00242439696411030857397215245841*zz[37] +.00043764769774185701827561290396*zz[39] -.00020714032687001791275913078304*zz[41] -6.27434450418651556052610958029804e-05*zz[43] +1.15753438145956693483789208989316e-05*zz[45] +5.88385492454037978388597885697078e-06*zz[47] -3.12467740069633622086961449076033e-07*zz[49]; sgn=pow(-(Double)1,length_org-1); res=0; for(int i=0;i=0){ if(blfi_done_right[i][j][n]==1) { if(md==1) res=blfi_val_re_right[i][j][n]; if(md==2) res=blfi_val_im_right[i][j][n]; } if(blfi_done_right[i][j][n]==0){ Double u=(n+piv_org[i][j])*arg_blfi[i][j]; for(int k=0;kmax_pts || abs(shift+range-1) >max_pts){ success=0; return 0; } if(blfi==1){ for(int n=shift-range+1;n=2*range){ denom=blfi_begin+length_blfi; klog_blfi[i][div_blfi]=LOG(denom); qlog_blfi[i][div_blfi]=log(1+(Double)blfi_remain/denom)/2; bbeta[i][div_blfi]=beta_fac_mult*qlog_blfi[i][div_blfi]; blambda[i][div_blfi]=(beta_fac_mult+1)*qlog_blfi[i][div_blfi]/2; bepsilon[i][div_blfi]=(beta_fac_mult-1)*qlog_blfi[i][div_blfi]/2; if(bbeta[i][div_blfi]<= qlog_blfi[i][div_blfi]){ cout<<"Error: choice of beta is producing beta <= tau!"<<"\n"; return; } arg_blfi[i][div_blfi]=Pi/bbeta[i][div_blfi]; inv_arg_blfi[i][div_blfi]=bbeta[i][div_blfi]/Pi; Double temp=bbeta[i][div_blfi]*t/Pi; piv_org[i][div_blfi]=temp-fmod(temp,1); for(int k=0;k=2*range){ denom=blfi_begin+length_blfi; klog_blfi[lgdiv+1][div_blfi]=LOG(denom); qlog_blfi[lgdiv+1][div_blfi]=log(1+(Double)blfi_block_size/denom)/2; bbeta[lgdiv+1][div_blfi]=beta_fac_mult*qlog_blfi[lgdiv+1][div_blfi]; blambda[lgdiv+1][div_blfi]=(beta_fac_mult+1)*qlog_blfi[lgdiv+1][div_blfi]/2; bepsilon[lgdiv+1][div_blfi]=(beta_fac_mult-1)*qlog_blfi[lgdiv+1][div_blfi]/2; if(bbeta[lgdiv+1][div_blfi]<= qlog_blfi[lgdiv+1][div_blfi]){ cout<<"Error: choice of beta is producing beta <= tau !"<<"\n"; return; } arg_blfi[lgdiv+1][div_blfi]=Pi/bbeta[lgdiv+1][div_blfi]; inv_arg_blfi[lgdiv+1][div_blfi]=bbeta[lgdiv+1][div_blfi]/Pi; Double temp=bbeta[lgdiv+1][div_blfi]*t/Pi; piv_org[lgdiv+1][div_blfi]=temp-fmod(temp,1); for(int k=0;k2.)){ cout<<"Error: blfi_block_growth not in [1,2] !"<<"\n"; return 0; } if(beta_fac_mult<=0.){ cout<<"Error: current choice of beta does not work!"<<"\n"; return 0; } if(bc<1.){ cout<<"Error: choose bc >= 1."<<"\n"; return 0; } if(blfi_block_size_org<0){ cout<<"Error: blfi_block_size_org must be >= 0!"<<"\n"; return 0; } if(sin_terms>=5){ cout<<"Error: sin_terms is larger than expected. Please add more terms to sin_cof array!"<<"\n"; return 0; } if(range<1.){ cout<<"Error: range is < 1"<<"\n"; return 0; } return 1; } //function to to collect statistics on various parameter used in blfi void init_blfi_simulate(){ int tail_blfi,length_blfi,div_blfi,blfi_fin; int blfi_block_size_sim=blfi_block_size_org; int blfi_begin_sim=length_split; total_blocks=0; for (int i=1;ierror_tolerance){ bc=bc+1; bc2=pow(bc,2); kernel_fac=bc/sinh(bc); range=(int)((beta_fac_mult/((beta_fac_mult-1)/2))*bc/Pi); blfi_block_size_org=(int)(2*range/blfi_fac); int length_split0=(int)(sqrt(length_split_fac*blfi_block_size_org*length_org*2*Pi/beta_fac_mult)); length_split=min(length_org,length_split0-length_split0%blfi_block_size_org+blfi_block_size_org); lgdiv=(int)(log((Double)length_org/(Double)length_split)/log((Double)2)); mult_fac=kernel_fac*sinh_mult_fac*((beta_fac_mult+1)/2)/beta_fac_mult; approx_blfi_mean_spacing=Pi/(beta_fac_mult*log(1+(Double)blfi_block_size_org/length_split)/2); max_pts=2*(int)((2*interval_length+1)/approx_blfi_mean_spacing+4*range+2); init_blfi_simulate(); error_blocks=2*exp(-bc)*blfi_block_size_org*max((Double) 1,pow(blfi_block_growth/sqrt((Double)2),lgdiv))*sqrt((Double)total_blocks)/(sqrt(length_split)); } if(!check()){ delete [] num_blocks; delete [] size_blocks; return 0; } if(length_split<=0.){ cout<<"Error: length_split must be positive!"<<"\n"; return 0; } init_arrays(md); if(md==1){ Double *klog0_temp=new Double[length_split_prev+1]; Double *ksqrt0_temp=new Double[length_split_prev+1]; for(int i=1;i=2*range){ if(blfi_block_size*blfi_fac<2*range-1){ cout<<"Warning: blfi_block_size "<=2*range){ if(blfi_block_size*blfi_fac<2*range-1){ cout<<"Warning: blfi_block_size "<0.5){ cout<<"Error: mean seperation of input points is too large!"<<"\n"; clean_arrays(0); return 0; } return success; } //function to compute zeta and re-initialize the algorithm if an error is enountered. Complex rs(Double t, Double error_given, Double input_mean_spacing_given, int &success, const char *return_type){ if(!rs_flag) { success=set_up(t,error_given,input_mean_spacing_given); if(!success) return 0; rs_flag=1; } Complex res=0; res=my_zeta(t,success); if(success==0) { clean_arrays(0); success=set_up(t,error_given,input_mean_spacing_given); if(!success) return 0; } if (!strcmp(return_type,"rotated pure")) return res*exp(I*(imag(log_GAMMA(((Double)1/2+I*t)/2)) - (t/2)*log(Pi))); else return res; } /* //the main function int main(){ //------------User input---------------------- Double t=pow(10.,10); Double input_mean_spacing_given=.01; Double error_tol=1E-9; //-------------------------------------------- int success; Complex res; for(int v=0;v<20000;v++){ res=rs(t,error_tol,input_mean_spacing_given,success); if(!success) return 0; cout<