xref: /dports/math/lcalc/lcalc-1.23/src/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 tweak; //used in value_via_Riemann_sum to play with the angle

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

int *prime_table;
int number_primes=0;

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;

bool only_use_dirichlet_series=false; //whether to compute just using the Dirichlet series
int N_use_dirichlet_series; //if so, how many terms in the Dirichlet series to use.

//------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

// Riemann siegel band limited interpolation arrays and variables ----------

const Double sin_cof[]={1,-(Double)1/6,1/120,-(Double)1/5040,(Double)1/362880,-(Double)1/39916800};
const Double sinh_mult_fac=(Double)1/2;
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;

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

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;
Double input_mean_spacing_given=.01; //rough anticipated sampling distance between t's for zeta(1/2+It).
//Should be set if known.

//------------------------------------------------------------------


//-----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,DIGITS);
    tolerance_sqrd=tolerance*tolerance;
    tolerance2=pow(.1,(DIGITS-DIGITS2+1));
    tolerance3=pow(.1,(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
    tweak=1;

    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));
    }

    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;

    extend_prime_table(100);

}

void delete_globals(){

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

    if (prime_table) delete [] prime_table;

}

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;
}

void extend_prime_table(int x)
{
  int *sieve;
  int sqrtx;
  int n, p;
  int number_primes_guess;

  sieve = new int[x+1];

  if (prime_table) delete [] prime_table;
  number_primes = 0;

  // Use explicit prime number theorem to allocate enough memory for primes
  number_primes_guess = max((int)ceil(x/(log(x)-4))+1, 100);

  prime_table = new int[number_primes_guess];

  if (my_verbose > 0)
  {
    cout << "extending prime table to: " << x << "; ";
    cout << "guessed " << number_primes_guess << " primes; ";
  }
  for (n=0;n<=x;n++)
    sieve[n] = 1;

  p = 2;
  sqrtx = (int) sqrt(x);
  while (p <= sqrtx)
  {
    n = 2*p;
    while (n <= x)
    {
      sieve[n] = 0;
      n = n + p;
    }
    do
      p++;
    while (sieve[p] == 0);
  }
  for (n=2;n<=x;n++)
  {
    if (sieve[n] == 1)
    {
      prime_table[number_primes] = n;
      number_primes++;
    }
  }
  delete [] sieve;

  if (my_verbose > 0)
    cout << "found " << number_primes << " primes." << endl;

  return;
}

int get_prime(int j)
{
  while (j >= number_primes)
  {
    extend_prime_table(prime_table[number_primes-1]*2);
  }

  return prime_table[j];
}