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

Lglobals.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 "Lglobals.h"


//-----Global variables---------------------------------------


int my_verbose=0;       // verbosity level: 0 means no verbose

int DIGITS, DIGITS2; // precision and sacrifice
int DIGITS3; // how many digits to output
Double tolerance;
Double tolerance_sqrd;
Double tolerance2;
Double tolerance3;

int global_derivative; //which derivative to compute
int max_n; // largest n used in a dirichlet series while computing an L-value.

Double A; //controls the 'support' of g(w) in the Riemann sum method
Double incr; //the increment in the Riemann sum method

Double *LG;             // lookup table for log(n)
Double *two_inverse_SQUARE_ROOT;    // lookup table for sqrt(n)
int number_sqrts=0;     // how many sqrt(n)'s to store
int number_logs=0;      // how many log(n)'s to store

Double *bernoulli;             // lookup table for bernoulli numbers

bool print_warning=true;

Long my_LLONG_MAX; //equals LLONG_MAX defined in limits.h. For
//reasons doing with compilation bugs, I determine it myself

Double Pi;
Double log_2Pi;
Complex I;

//------Incomplete gamma function global variables--------------

Complex last_z;         // the last z to be considered in inc_GAMMA
Complex last_w;         // the last w to be considered in inc_GAMMA
Complex last_comp_inc_GAMMA; // g(last_z,last_w)

Complex last_z_GAMMA;  //the last z to be considered in GAMMA
Complex last_log_G;    //the last log(GAMMA(z));

Double temme_a[1002],temme_g[501];
//used in Temme's asymptotic expansion of the
//incomplete gamma function
//XXXX might need more terms if I go to higher precision


//-----intializing and cleaning up routines----------------------

//n is the number of log(n)'s to precompute
//call this function after changing the DIGITS global variable
void initialize_globals(int n){


#ifdef USE_MPFR
    if(!DIGITS) DIGITS=40;
    if(!DIGITS2) DIGITS2=2;
    if(!DIGITS3) DIGITS3=DIGITS-DIGITS2;
    //__mpfr_default_fp_bit_precision=Int(DIGITS*log(10.)/log(2.));
    //mpfr_set_default_prec(__mpfr_default_fp_bit_precision);
    mpfr_set_default_prec(Int(DIGITS*log(10.)/log(2.)));
    reset(Pi);
    reset(log_2Pi);
    reset(I);
    reset(tolerance);
    reset(last_z);
    reset(last_w);
    reset(last_comp_inc_GAMMA);
    reset(last_z_GAMMA);
    reset(last_log_G);
    loop(i,0,1002) reset(temme_a[i]);
    loop(i,0,501) reset(temme_g[i]);
    mpfr_const_pi(Pi.get_mpfr_t(), __gmp_default_rounding_mode);
    log_2Pi=log(2*Pi);
#else
    if(!DIGITS){
        typedef std::numeric_limits< Double > tmp_Double;
        DIGITS=tmp_Double::digits10-1;
    }
    if(!DIGITS2) DIGITS2=2;
    if(!DIGITS3) DIGITS3=DIGITS-DIGITS2;
    Pi= 3.14159265358979323846264338328L;
    log_2Pi= 1.837877066409345483560659472811L;
    //Pi= 3.1415926535897932384626433832795028841971693993751L;
    //log_2Pi= 1.8378770664093454835606594728112352797227949472756L;
#endif
    I=Complex(0.,1.);
    tolerance=pow(1./10,DIGITS);
    tolerance_sqrd=tolerance*tolerance;
    tolerance2=pow(1./10,(DIGITS-DIGITS2+1));
    tolerance3=pow(1./10,(DIGITS3+1));

    A=1./(16*Pi*Pi*23*DIGITS/10);
    incr=2*Pi*.5/(log(10.)*DIGITS); //.5 here is current v in Riemann sum method 

    last_z=1.;
    last_w=0.;
    last_comp_inc_GAMMA=0.;
    last_z_GAMMA=1.;
    last_log_G=0.;

    global_derivative=0;
    max_n=1;

    numeric_limits<long long> ll;
    my_LLONG_MAX = ll.max(); //used once, but in the elliptic curve module Lcommandline_ellipitic.cc
    //pari defines max #define max(a,b) ((a)>(b)?(a):(b)) and this causes problems with ll.max
    //so I put it here as a global.
    //also look at: http://www.google.com/answers/threadview?id=334912
    //Re: g++ problems with strtoll (LLONG_MIN, LLONG_MAX, ULLONG_MAX not defined)
    //for an explanation of why we just don't call LLONG_MAX (compiler inconsistencies).
 
    int j,k;
    Double r,x,dsum;

    number_logs=n;
    if(LG) delete[] LG;
    LG = new Double[n+1];
    for(k=1;k<=n;k++) LG[k]=log((Double)k);

    number_sqrts=n;
    if(two_inverse_SQUARE_ROOT) delete[] two_inverse_SQUARE_ROOT;
    two_inverse_SQUARE_ROOT = new Double[n+1];
    for(k=1;k<=n;k++) two_inverse_SQUARE_ROOT[k]=2/sqrt((Double)k);

    if(bernoulli) delete[] bernoulli;
    bernoulli= new Double[DIGITS+1];

    bernoulli[0]=1;
    for(k=1;k<=DIGITS;k++){
        r=k+1; x=0;
        for(j=1;j<=k;j++){
            r=r*(k+1-j)*1./(j+1);
            x=x-bernoulli[k-j]*r;
        }
        bernoulli[k]=x/(k+1);
    }

    //XXXXX this should depend on DIGITS
    temme_a[1]=1;

    temme_a[2]=(Double) 3.;
    temme_a[2]=1./temme_a[2]; //doing temme_a[2]=1./3 causes problems with long doubles or multiprecision. 1.L/3 
                               //could cause trouble for multiprecision.
    temme_a[3]=(Double)36.;
    temme_a[3]=1./temme_a[3];

    for (int i=4;i<=1001;i++) {
      dsum=0;
      for (int j=2;j<=(i-1);j++)
        dsum+=j*temme_a[j]*temme_a[i-j+1];
      temme_a[i]=(temme_a[i-1]-dsum)/((Double) (i+1.0));
    }

    for(int i=1;i<=500;i++)
      temme_g[i]=(1-2*(i%2))*
      dfac(2*i+1)*temme_a[2*i+1];
    temme_g[0]=1;


}

void delete_globals(){

    if(LG) delete [] LG;
    if(two_inverse_SQUARE_ROOT) delete [] two_inverse_SQUARE_ROOT;
    if(bernoulli) delete [] bernoulli;

}

void extend_LG_table(int m){


    int n;

    Double *tmp_LG; //used to copy over the old values
    tmp_LG = new Double[number_logs+1];
    for(n=1;n<=number_logs;n++) tmp_LG[n]=LG[n];

    delete [] LG;
    int new_number_logs=(int)(1.5*m); // extend table by 50 percent

    LG = new Double[new_number_logs+1];
    for(n=1;n<=number_logs;n++) LG[n]=tmp_LG[n];
    for(n=number_logs+1;n<=new_number_logs;n++) LG[n]=log((Double)n);
    number_logs=new_number_logs;

    if(my_verbose>0)
    cout << endl << "extended log table to: " << number_logs << endl;

    delete [] tmp_LG;

}

void extend_sqrt_table(int m){


    int n;

    Double *tmp_sqrt; //used to copy over the old values
    tmp_sqrt = new Double[number_sqrts+1];
    for(n=1;n<=number_sqrts;n++) tmp_sqrt[n]=two_inverse_SQUARE_ROOT[n];

    delete [] two_inverse_SQUARE_ROOT;
    int new_number_sqrts=(int)(1.5*m); // extend table by 50 percent

    two_inverse_SQUARE_ROOT = new Double[new_number_sqrts+1];
    for(n=1;n<=number_sqrts;n++) two_inverse_SQUARE_ROOT[n]=tmp_sqrt[n];
    for(n=number_sqrts+1;n<=new_number_sqrts;n++) two_inverse_SQUARE_ROOT[n]=2/sqrt((Double)n);
    number_sqrts=new_number_sqrts;

    if(my_verbose>0)
    cout << endl << "extended sqrt table to: " << number_sqrts << endl;

    delete [] tmp_sqrt;

}


// a lot of Q time is spend here. This can be precomputed XXXXX
Double dfac(int i)
{
    Double t=1.;
    int n=i;
    if(i==1 || i==0) return 1;
    while(n>0) {
        t=t*n;
        n-=2;
    }
    //return (Double)(double)t;
    return t;
}

Generated by  Doxygen 1.6.0   Back to index