Logo Search packages:      
Sourcecode: lcalc version File versions  Download package

Lriemannsiegel.cc

/*

   Copyright (C) 2001,2002,2003,2004 Michael Rubinstein

   This file is part of the L-function package L.

   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 2
   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.

   Check the License for details. You should have received a copy of it, along
   with the package; see the file 'COPYING'. If not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

*/


#include "Lriemannsiegel.h"           //for global variables


//This implementation of Riemann Siegel for zeta(1/2+it) can be improved.
//To do: Odlyzko Schonhage algorithm XXX

// ZETA FUNCTION
Complex Zeta(Complex s,char *return_type) {

  Complex L_value;

  L_value= siegel(s);

  //this returns Zeta(s)
  if (!strcmp(return_type,"pure"))
      return L_value;

  //returns Zeta(s) rotated to be real on critical line
  else if (!strcmp(return_type,"rotated pure"))
      return L_value*exp(I*(imag(log_GAMMA(s/2)) - (imag(s)/2)*log(Pi)));

  return L_value;

}

// RIEMANN-SIEGEL ZETA EVALUATION
Complex siegel(Complex s, int n) {

  Complex theta, result = 0;
  Complex errorTerms;
  Double tau,p, t;
  int N, i;

  //Complex tmp1,tmp2;
  Double tmp3;

  tau = sqrt(abs(imag(s))/(2*Pi));
  N = Int(floor(tau));
  p = tau - N;

  t = imag(s);

  theta = imag(log_GAMMA(s/2)) - (imag(s)/2)*log(Pi);

  //tmp1=cos(theta);
  //tmp2=sin(theta);

  for (i = N; i > 0; i--){
    tmp3=t*LOG(i);
    result = result+two_inverse_sqrt((double)i)*cos(tmp3-theta);
  }

  max_n=N;


  errorTerms = C(0,p) + C(1,p)*pow(tau,-1) + C(2,p)*pow(tau,-2) + C(3,p)*pow(tau,-3) + C(4,p)*pow(tau,-4) +C(5,p)*pow(tau,-5);

  errorTerms =errorTerms*pow(tau,-.5)*pow(-1., N-1);


  result = result+errorTerms;

  // this is Z.  Now we multiply to get zeta

  result =result* exp(-1.*I*theta);


  return result;

}

Double C(int n,Double p){

    Double z=p-.5;

if (n==0) return
.3826834323650897717284599840304*pow(z,0)
+1.74896187231008179744118586948533*pow(z,2)
+2.11802520768549637318456427826417*pow(z,4)
-.87072166705114807391892407738239*pow(z,6)
-3.4733112243465167073064116693758*pow(z,8)
-1.66269473089993244964313630119286*pow(z,10)
+1.21673128891923213447689352804417*pow(z,12)
+1.30143041610079757730060538099786*pow(z,14)
+.03051102182736167242108987123981*pow(z,16)
-.37558030515450952427981932122934*pow(z,18)
-.10857844165640659743546975901329*pow(z,20)
+.05183290299954962337576051067322*pow(z,22)
+.02999948061990227592040084956912*pow(z,24)
-.00227593967061256422601994851021*pow(z,26)
-.00438264741658033830594007013585*pow(z,28)
-.00040642301837298469930723272116*pow(z,30)
+.00040060977854221139278910314608*pow(z,32)
+8.97105799138884129783418195378689e-05*pow(z,34)
-2.30256500272391071161029452573899e-05*pow(z,36)
-9.38000660190679248471972940127474e-06*pow(z,38)
+6.32351494760910750424986123959430e-07*pow(z,40)
+6.55102281923150166621223123133411e-07*pow(z,42)
+2.21052374555269725866086890381876e-08*pow(z,44)
-3.32231617644562883503133517017624e-08*pow(z,46)
-3.73491098993365608176460476015222e-09*pow(z,48)
+1.24450670607977391951510000249366e-09*pow(z,50);
if (n==1) return
-.05365020525675069405998280791133*pow(z,1)
+.11027818741081482439896362071917*pow(z,3)
+1.23172001543152263131956529162206*pow(z,5)
+1.26349648627994578841755482191213*pow(z,7)
-1.69510899755950301844944739906596*pow(z,9)
-2.99987119676501008895548735894141*pow(z,11)
-.10819944959899208642692257787438*pow(z,13)
+1.94076629462127126879387632539716*pow(z,15)
+.78384235615006865328843457488694*pow(z,17)
-.50548296679003659187902141326173*pow(z,19)
-.3845072349605797405134273885311*pow(z,21)
+.03747264646531532067594447494023*pow(z,23)
+.09092026610973176317258142450576*pow(z,25)
+.01044923755006450921820113972659*pow(z,27)
-.01258297965158341649747892224592*pow(z,29)
-.00339950372115127408505894886137*pow(z,31)
+.00104109505377148912682954240655*pow(z,33)
+.00050109490511184868603556526727*pow(z,35)
-3.95635966900318155954711855696337e-05*pow(z,37)
-4.76245924535718963865409830268035e-05*pow(z,39)
-1.85393553380851322734349064569117e-06*pow(z,41)
+3.19369180800689720404663539343268e-06*pow(z,43)
+4.09078076085060663265089453677018e-07*pow(z,45)
-1.54466243325766321284375723273104e-07*pow(z,47)
-3.46630749176913317222559405934073e-08*pow(z,49);
if (n==2) return
.00518854283029316849378458151923*pow(z,0)
+.00123786335522538984133826974438*pow(z,2)
-.18137505725166997411491896409414*pow(z,4)
+.14291492748532126541165603376514*pow(z,6)
+1.33033917666875653250993329998546*pow(z,8)
+.35224723534037336775327655505836*pow(z,10)
-2.4210015958919507237815305433405*pow(z,12)
-1.67607870225381088533346181492372*pow(z,14)
+1.36894167233283721842349153807076*pow(z,16)
+1.55390194302229832214563952655935*pow(z,18)
-.17221642734729980519582586998918*pow(z,20)
-.63590680550454309889704902355845*pow(z,22)
-.09911649873041208105423564341370*pow(z,24)
+.14033480067387008950738254898316*pow(z,26)
+.04782352019827292236438803506512*pow(z,28)
-.01735604064147978079795864709223*pow(z,30)
-.01022501253402859184447660413126*pow(z,32)
+.00092741491597948878994270014371*pow(z,34)
+.00135721943723733853452533619958*pow(z,36)
+6.41369012029388008996238736394533e-05*pow(z,38)
-.00012300805698196629883342322937*pow(z,40)
-1.83135074047892025547675543979621e-05*pow(z,42)
+7.82162860432262730850139938461872e-06*pow(z,44)
+2.00875424847599455034985293919157e-06*pow(z,46)
-3.35327653931857137372749727241453e-07*pow(z,48)
-1.46160209174182309264510097122760e-07*pow(z,50);
if (n==3) return
-.00267943218143891380853967145989*pow(z,1)
+.02995372109103514963731329491570*pow(z,3)
-.04257017254182869798501935111688*pow(z,5)
-.28997965779803887506893209478669*pow(z,7)
+.48888319992354459725374746407169*pow(z,9)
+1.23085587639574608119312504336294*pow(z,11)
-.82975607085274087041796910432976*pow(z,13)
-2.24976353666656686652045012659903*pow(z,15)
+.07845139961005471379365473620184*pow(z,17)
+1.74674928008688940039198666645219*pow(z,19)
+.45968080979749935109237306173169*pow(z,21)
-.66193534710397749464339040008983*pow(z,23)
-.31590441036173634578979632973316*pow(z,25)
+.12844792545207495988511847476209*pow(z,27)
+.10073382716626152300969450207513*pow(z,29)
-.00953018384882526775950465984230*pow(z,31)
-.01926442168751408889840098069714*pow(z,33)
-.00124646371587692917124790716458*pow(z,35)
+.00242439696411030857397215245841*pow(z,37)
+.00043764769774185701827561290396*pow(z,39)
-.00020714032687001791275913078304*pow(z,41)
-6.27434450418651556052610958029804e-05*pow(z,43)
+1.15753438145956693483789208989316e-05*pow(z,45)
+5.88385492454037978388597885697078e-06*pow(z,47)
-3.12467740069633622086961449076033e-07*pow(z,49);
if (n==4) return
.00046483389361763381853630462560*pow(z,0)
-.00402264294613618830391153989145*pow(z,2)
+.00384717705179612688359130685272*pow(z,4)
+.06581175135809486002088309200741*pow(z,6)
-.19604124343694449117695528448205*pow(z,8)
-.20854053686358853244400012794494*pow(z,10)
+.95077541851417509458477574151058*pow(z,12)
+.53415353129148739760517592459894*pow(z,14)
-1.67634944117634007959116448203404*pow(z,16)
-1.07674715787512899278784663510432*pow(z,18)
+1.23533930165659698528788361189251*pow(z,20)
+1.02578253400572757718348949577914*pow(z,22)
-.40124095793988544378728137523313*pow(z,24)
-.50366639951083034479585257591604*pow(z,26)
+.03573487795502744985807080163387*pow(z,28)
+.14431763086785416624285239495844*pow(z,30)
+.01509152741790346941712677290432*pow(z,32)
-.02609887477919436131761773965448*pow(z,34)
-.00612662837951926174904909908948*pow(z,36)
+.00307750312987084118476787782167*pow(z,38)
+.00115624789340887523161201204220*pow(z,40)
-.00022775966758472127472807733953*pow(z,42)
-.00014189637118181444432681579894*pow(z,44)
+7.46486030795591945312240984450313e-06*pow(z,46)
+1.24797016454091166174449988846871e-05*pow(z,48)
+4.86394518400209461907998084746180e-07*pow(z,50);
if (n==5) return
.00022686811845737363176557957245*pow(z,1)
+.00110812468537183880897586725284*pow(z,3)
-.01621857925555009106408484258686*pow(z,5)
+.05276503405398741662724126665649*pow(z,7)
+.02570880200903323999290010111095*pow(z,9)
-.38058660440806397264435991848146*pow(z,11)
+.22531987892642315322976926989838*pow(z,13)
+1.03445733164952217211304499657389*pow(z,15)
-.55282576970508137898888475296735*pow(z,17)
-1.52877126410780729962736571714169*pow(z,19)
+.32828366427719583672031669394059*pow(z,21)
+1.22911021854008706238425001239677*pow(z,23)
+.04093693938311529830689289790902*pow(z,25)
-.55860404726420193442735876775644*pow(z,27)
-.11241976368059115396788439789609*pow(z,29)
+.15212677711795591829295940144809*pow(z,31)
+.05173718845528038784023625510664*pow(z,33)
-.02561227689700728294043343196050*pow(z,35)
-.01296367251404617794428713962277*pow(z,37)
+.00254555748186116327806192744188*pow(z,39)
+.00211933195108777752885073213414*pow(z,41)
-9.19139194515677754051761292342159e-05*pow(z,43)
-.00024413466533855272657049552509*pow(z,45)
-1.36979826922833871237722416376381e-05*pow(z,47)
+2.06207850332842380642607669485206e-05*pow(z,49);

else return 0;
}


Generated by  Doxygen 1.6.0   Back to index