Commit 27561f54 authored by Tobias Hansen's avatar Tobias Hansen

Imported Upstream version 1.23

parent 2e1c0cc0
Michael Rubinstein-main author
Ralph Furmaniak- support for gmp
Kyle Wichert- Temme's uniform asymptotics for the incomplete gamma function
Being incorporated into next major version:
Kevin Mcgown- explicit formula test
Sourav Sen Gupta- algorithms from Tim Dokchitser's paper
Ghaith Hiary- Band limited interpolation.
R. Rishikesh- wrapper for sage
......@@ -2,21 +2,17 @@
Dear Colleagues,
Here is the latest version of the L-function c++ class library and, the command line
program lcalc. Two versions, one for linux and one for
mac os x, are available for download. It will probably compile on other flavours
of unix, though you might might need to adjust the makefile.
program lcalc.
To compile:
gunzip it: gunzip Llinux.tar.gz
(or gunzip Ldarwin.tar.gz)
gunzip it: gunzip L-1.22.tar.gz
untar it: tar -xvf Llinux.tar
(or tar -xvf Ldarwin.tar)
untar it: tar -xvf L-1.22.tar
This creates a directory called L
This creates a directory called L-1.22
go to the src directory: cd L/src
go to the src directory: cd L-1.22/src
compile: make
......
......@@ -24,6 +24,10 @@
#ifndef L_H
#define L_H
#ifdef _OPENMP
#include <omp.h>
#endif
#include <iomanip> //for manipulating output such as setprecision
#include <fstream> //for file input output
#include <string.h> //for string functions such as strcpy
......@@ -31,12 +35,15 @@
#include <iostream> //for ostrstream
#include <math.h>
#include<vector>
#include "Lglobals.h" //for global variables
#include "Lmisc.h" //things like sn or LOG
#include "Lgamma.h" //incomplete gamma function code
//#include "Lprecomp.h"
#include "Lriemannsiegel.h" //Riemann Siegel formula
#include "Lriemannsiegel_blfi.h" //Hiary's Riemann Siegel formula using band limited interpolation
//-----THE L Function Class-----------------------------------
......@@ -125,7 +132,7 @@ public:
//-----Constructor: initialize the L-function from given data------------------
L_function (char *NAME, int what_type, int N, ttype *coeff, long long Period,
L_function (const char *NAME, int what_type, int N, ttype *coeff, long long Period,
Double q, Complex w, int A, Double *g, Complex *l,
int n_poles, Complex *p, Complex *r)
{
......@@ -187,7 +194,7 @@ public:
}
//-----Constructor: initialize the L-function from given data------------------
L_function (char *NAME, int what_type, int N, ttype *coeff, long long Period,
L_function (const char *NAME, int what_type, int N, ttype *coeff, long long Period,
Double q, Complex w, int A, Double *g, Complex *l)
{
if(my_verbose>1)
......@@ -381,8 +388,9 @@ public:
//-----addition operator--------------------
//returns the L-function whose basic data (funct eqn, dirichlet coefficients)
//returns the 'L-function' whose basic data (funct eqn, dirichlet coefficients)
//is that of this added to that of L
//not particularly useful
L_function operator + (const L_function &L)
{
L_function L2;
......@@ -423,9 +431,10 @@ public:
return L2;
}
//-----addition operator--------------------
//-----multiplication operator--------------------
//returns the L-function whose basic data (funct eqn, dirichlet coefficients)
//is that of this added to that of L
//is that of this times the scalar t
//not particularly useful
L_function operator * (double t)
{
L_function L2;
......@@ -492,18 +501,27 @@ public:
//#include "Lvalue.h" //value via Riemann sum, via gamma sum, various options for value
Complex find_delta(Complex s,Double g);
Complex value_via_Riemann_sum(Complex s,char *return_type="pure");
Complex value_via_gamma_sum(Complex s, char *return_type="pure");
Complex value(Complex s, int derivative = 0, char *return_type="pure");
Complex value_via_Riemann_sum(Complex s, const char *return_type="pure");
Complex value_via_gamma_sum(Complex s, const char *return_type="pure");
Complex value(Complex s, int derivative = 0, const char *return_type="pure");
//#include "Lfind_zeros.h" //finding zeros routine
Double zeros_zoom_brent(Double L1, Double L2, Double u, Double v);
void find_zeros(Double t1, Double t2, Double step_size, char* filename="cout", char* message_stamp="");
void find_zeros_via_gram(Double t1, Long count=0,Double max_refine=1025,char* filename="cout", char* message_stamp="");
void find_zeros(Double t1, Double t2, Double step_size, const char* filename="cout", const char* message_stamp="");
void find_zeros_v(Double t1, Double t2, Double step_size, vector<Double> &result);//This is the same as above function
void find_zeros_via_gram(Double t1, Long count=0,Double max_refine=1025, const char* filename="cout", const char* message_stamp="");
int compute_rank(bool print_rank=false);
void verify_rank(int rank);
void find_zeros_via_N(Long count=0,bool do_negative=true,Double max_refine=1025, int rank=-1,char* filename="cout", char* message_stamp="");
void find_zeros_elaborate(Double t1, Long count=0,Double max_refine=1025,char* filename="cout", char* message_stamp="");
void find_zeros_via_N(Long count=0,bool do_negative=true,Double max_refine=1025, int rank=-1, bool test_explicit_formula=false, const char* filename="cout", const char* message_stamp="");
void find_zeros_via_N_v(Long count,bool do_negative,Double max_refine, int rank, bool test_explicit_formula, vector<Double> &result);
void find_zeros_elaborate(Double t1, Long count=0,Double max_refine=1025, const char* filename="cout", const char* message_stamp="");
//#include "Ldokchitser.h" //dokchitser algorithm for what he calls phi(t), i.e. inverse
void phi_series(int precision);
//#include "Lexplicit_formula.h"
int dirichlet_coeffs_log_diff(int num_coeff, Complex *c);
int test_explicit_formula(Double A, Double x_0, Double *zero_table, int number_zeroes, Complex *c, int num_coeffs);
};
......@@ -518,5 +536,9 @@ public:
//#include "Ltaylor_series.h" //for computing taylor series for Dirichlet series
#include "Lvalue.h" //value via Riemann sum, via gamma sum, various options for value
#include "Lfind_zeros.h" //finding zeros routine
#include "Ldokchitser.h" //dokchitser algorithm for what he calls phi(t), i.e. inverse
//mellin transform
#include "Lexplicit_formula.h" //for testing zeros with the explicit formula
#endif
......@@ -31,7 +31,7 @@
//===============functions================================================
int quadratic_twists(Long D1, Long D2,Double x,Double y,int count,Double step_size,char *what_to_do,bool do_only_even_twists,int desired_rank=-1);
int all_twists(Long D1, Long D2,Double x,Double y,int count,Double step_size,char *what_to_do,int twist_type,int print_character);
int quadratic_twists(Long D1, Long D2,Double x,Double y,int count,Double step_size, const char *what_to_do,bool do_only_even_twists,bool test_explicit_formula,int desired_rank=-1);
int all_twists(Long D1, Long D2,Double x,Double y,int count,Double step_size, const char *what_to_do,int twist_type,int print_character,bool test_explicit_formula);
#endif
......@@ -30,8 +30,8 @@
//-----functions--------------------------------------------------------------
void compute_values(Double x,Double y,char *file_name="",Double x3=0,Double y3=0,Long count=0);
void compute_zeros(Double x, Double y,Double step_size, Long count=0,int rank=-1);
void compute_values(Double x,Double y,const char *return_type="pure",const char *file_name="",Double x3=0,Double y3=0,Long count=0);
void compute_zeros(Double x, Double y,Double step_size, Long count=0,int rank=-1,bool test_explicit_formula=false);
void L_interpolate(Double x, Double y,Double step_size, int n=1000);
#endif
......
......@@ -15,7 +15,7 @@ typedef long long Long;
// To aid conversion to int value
inline double to_double(const Double& x) {
inline double lcalc_to_double(const Double& x) {
#ifdef USE_MPFR
return x.get_d();
#else
......@@ -23,16 +23,19 @@ inline double to_double(const Double& x) {
#endif
}
#ifdef USE_MPFR
inline double to_double(const double& x) { return x; }
inline double lcalc_to_double(const double& x) { return x; }
#endif
//inline double to_double(const long double& x) { return x; }
inline double to_double(const int& x) { return x; }
inline double to_double(const long long& x) { return x; }
inline double to_double(const short& x) { return x; }
inline double to_double(const char& x) { return x; }
#define Int(x) (int)(to_double(x))
#define Long(x) (Long)(to_double(x))
#define double(x) (double)(to_double(x))
//inline double lcalc_to_double(const long double& x) { return x; }
inline double lcalc_to_double(const int& x) { return x; }
inline double lcalc_to_double(const long long& x) { return x; }
inline double lcalc_to_double(const short& x) { return x; }
inline double lcalc_to_double(const char& x) { return x; }
inline double lcalc_to_double(const long int& x) { return x; }
inline double lcalc_to_double(const unsigned int& x) { return x; }
inline double lcalc_to_double(const long unsigned int& x) { return x; }
#define Int(x) (int)(lcalc_to_double(x))
#define Long(x) (Long)(lcalc_to_double(x))
#define double(x) (double)(lcalc_to_double(x))
template<class T> inline int sn(T x)
{
......
#ifndef Ldokchitser_H
#define Ldokchitser_H
//finding the explicit taylor series for \phi(t) using Dokchitser algo
#define MYDIGITS 5 // estimate of precision ... will be set using the precision variable
void mult_poly_taylor(Complex *, Complex *, Complex *, int );
template <class ttype>
void L_function <ttype>::
phi_series(int precision)
{
cout << "-----------------------------------------------"<< endl << endl;
cout << "phi series for " << name << " L_function" << endl << endl;
int j,k;
// constructing the equivalence classes Lambda[k] for k = 1 to N
int pordtmp[a+1];
Complex diff;
Complex *lambda_k;
int *l;
for (j=1;j<=a;j++)
pordtmp[j]= 1;
for (j=1;j<=a;j++)
for (k=1;k<=a;k++)
if (j != k)
{
diff = 2*(lambda[j] - lambda[k]);
if((imag(diff)==0) && (fmod(real(diff),2) == 0) && (real(diff)<=0))
{
pordtmp[j]+=pordtmp[k];
pordtmp[k]=0;
}
}
Complex temp_lambda_k[a];
int temp_l[a];
j=1;
for (k=1;k<=a;k++)
if (pordtmp[k]!=0)
{
temp_lambda_k[j]= lambda[k];
temp_l[j]=pordtmp[k];
j++;
}
int N = j-1;
lambda_k = new Complex[N+1];
l = new int[N+1];
for (j=1;j<=N;j++)
{
lambda_k[j] = temp_lambda_k[j];
l[j] = temp_l[j];
}
cout << "-----------------------------------------------"<< endl << endl;
cout << "There are "<< N << " equivalence classes Lambda[j]"<<endl;
cout<< "The equivalence classes Lambda[j] for poles are represented by"<<endl;
for (j=1;j<=N;j++)
{
cout << "lambda_k["<<j<<"] = "<< lambda_k[j] << " with order l["<<j<<"] = "<<l[j]<<endl;
}
cout<<endl;
// compute the values m[j] for the respective lambda_k[j]
Complex m[N+1];
for (j=1;j<=N;j++)
m[j] = -2*lambda_k[j] + 2;
// compute sum_{k=1}^a log Gamma((s+m[j]+2*lambda[k])/2) for each j
int n,fact_n;
Complex log_Gamma[N+1][a+1][MYDIGITS+1];
Complex sum_log_Gamma[N+1][MYDIGITS+1];
for (j=1;j<=N;j++)
for (n=0;n<=MYDIGITS;n++)
sum_log_Gamma[j][n] = 0;
for (j=1;j<=N;j++)
{
for (k=1;k<=a;k++)
{
fact_n = 1;
for (n=0;n<=MYDIGITS;n++)
{
if (n!=0)
fact_n = fact_n*n;
log_Gamma[j][k][n] = pow(0.5,n)*(log_GAMMA((m[j]/2 + lambda[k]),n))/fact_n;
sum_log_Gamma[j][n] = sum_log_Gamma[j][n] + log_Gamma[j][k][n];
}
}
}
cout << endl;
// compute the exponential taylor series for gamma = exp(sum_log_Gamma)
Complex exp_sum_log_Gamma[N+1][MYDIGITS+1][MYDIGITS+1]; // symmetric functions
Complex gamma[N+1][MYDIGITS+1]; // gamma(s+m[j]) for j = 1 to N
Complex temp_gamma[MYDIGITS+1];
Complex temp_mult_gamma[MYDIGITS+1];
Complex temp_exp_sum_log_Gamma[MYDIGITS+1];
int fact_n_k;
for (j=1;j<=N;j++)
{
for (n=0;n<=MYDIGITS;n++)
exp_sum_log_Gamma[j][0][n] = 0;
exp_sum_log_Gamma[j][0][0] = exp(sum_log_Gamma[j][0]);
for (k=1;k<=MYDIGITS;k++)
{
fact_n_k = 1;
for (n=0;n<=MYDIGITS;n++)
{
if(n%k == 0)
{
if (n/k != 0)
fact_n_k = fact_n_k*(n/k);
exp_sum_log_Gamma[j][k][n] = pow(sum_log_Gamma[j][k],n/k)/fact_n_k;
}
else
exp_sum_log_Gamma[j][k][n] = 0;
}
}
}
for (j=1;j<=N;j++)
{
for (n=0;n<=MYDIGITS;n++)
temp_mult_gamma[n] = exp_sum_log_Gamma[j][0][n];
for (k=1;k<=MYDIGITS;k++)
{
for (n=0;n<=MYDIGITS;n++)
{
temp_exp_sum_log_Gamma[n] = exp_sum_log_Gamma[j][k][n];
temp_gamma[n] = temp_mult_gamma[n];
}
mult_poly_taylor(temp_gamma,temp_exp_sum_log_Gamma,temp_mult_gamma,MYDIGITS);
}
for (n=0;n<=MYDIGITS;n++)
gamma[j][n] = temp_mult_gamma[n];
}
cout << "-----------------------------------------------"<< endl;
cout << "The gamma(s+m[j]) coefficients are as follows"<< endl<<endl;
for(j=1;j<=N;j++)
{
cout<<"gamma(s+m["<<j<<"]) = "<<"gamma(s+"<<m[j]<<") = "<<gamma[j][0];
for (n=1;n<=MYDIGITS;n++)
cout << " + " << gamma[j][n] <<" s^"<<n;
cout<<endl;
}
cout << "-----------------------------------------------"<< endl;
}
#endif
// This file mainly due to Kevin McGown, with modifications by Michael Rubinstein
#ifndef Lexplicit_formula_H
#define Lexplicit_formula_H
inline Complex xxx_phi(Double A, Double x_0, Complex x); // the phi in the explicit formula
inline Complex xxx_phi_hat(Double A, Double x_0, Complex x); //its fourier transform
/***************************************************************************************************/
inline Complex xxx_phi(Double A, Double x_0, Complex x)
{
return exp(-A*(x-x_0)*(x-x_0));
}
inline Complex xxx_phi_hat(Double A, Double x_0, Complex x)
{
return sqrt(Pi/A) * exp(-2*Pi*I*x_0*x - Pi*Pi*x*x/A);
}
/***************************************************************************************************/
template <class ttype>
int L_function <ttype>::
dirichlet_coeffs_log_diff(int num_coeffs, Complex *c)
{
Complex b[num_coeffs+1];
int j, n, d1, ind;
Complex total, total2, temp;
if (what_type_L != 1 && what_type_L != -1
&& num_coeffs > number_of_dirichlet_coefficients)
{
cout << "Don't have enough Dirichlet coefficients." << endl;
return 1;
}
b[1] = 1;
if(my_verbose!=0) cout << "Computing " << num_coeffs << " Dirichlet coefficients of the logarithmic derivative" << endl;
//XXXXXXXX this should be computed just once, not each time test is called
for (n=2;n<=num_coeffs;n++)
{
total = 0;
total2 = 0;
for (j=1;j<=n/2;j++)
if (n % j == 0)
{
d1 = n/j;
if (what_type_L == -1)
temp = b[j];
else if (what_type_L == 1)
{
ind = d1 % period;
if (ind == 0)
ind = period;
temp = dirichlet_coefficient[ind]*b[j];
}
else
temp = dirichlet_coefficient[d1]*b[j];
total -= temp;
total2 += temp*LOG(d1);
}
b[n] = total;
c[n] = total2;
if(my_verbose>5) cout << "c[" << n << "] = " << c[n] << endl;
}
return 0;
}
/************************************************************************************************/
template <class ttype>
int L_function <ttype>::
test_explicit_formula(Double A, Double x_0, Double *zero_table, int number_zeros, Complex *c, int num_coeffs)
{
Double t, t_begin, t_end, t_step, u;
Double D;
Double total;
Double term1, term2, term3;
int p, n, x, j;
Double temp;
Double LHS, RHS;
//Complex *c;
//int num_coeffs;
int flag;
//num_coeffs = 150; //XXXXXXXXXX should depend on test required
//c = new Complex[num_coeffs+1]; //XXXXXXX move to L.h
//dirichlet_coeffs_log_diff(num_coeffs, c);
//compute the possible contribution from poles
term1 = 0;
for (j=1;j<=number_of_poles;j++)
term1 += real(xxx_phi(A, x_0, (pole[j]-0.5)/I));
// compute the contribution from the Gamma factors (integral of the log diff)
t_step = 0.02;
D = ceil(sqrt(DIGITS*log(10.0)/A)/t_step) * t_step;
t_begin = x_0 - D;
t_end = x_0 + D;
total = 0;
for (t=t_begin;t<=t_end;t=t+t_step)
{
temp = 0;
for (j=1;j<=this->a;j++)
{
temp += 2*real(log_GAMMA(gamma[j]/2 + I*t*gamma[j]+lambda[j], 1)*gamma[j]); //1 here calls the logarithmic derivative
//temp += log_GAMMA(gamma[j]/2 - I*t*gamma[j]+conj(lambda[j]), 1)*gamma[j];
}
total = total + real(xxx_phi(A, x_0, t)) * temp;
//cout << t << " total " << total << endl;
}
term2 = t_step*total + 2*log(Q)*sqrt(Pi/A);
term2 = 1/(2*Pi) * term2;
//compute the contribution from the Dirichlet coefficients
x = num_coeffs;
//extend_prime_table(x);
j = 0;
p = get_prime(j);
term3 = 0;
while (p <= x)
{
n = p;
while (n <= x)
{
temp = 2*real(c[n]*xxx_phi_hat(A, x_0, LOG(n)/(2*Pi)));// + conj(c[n])*xxx_phi_hat(A, x_0, -LOG(n)/(2*Pi));
term3 += two_inverse_sqrt(n)*temp;
n = n * p;
//cout << n << "term3 " << term3 << endl;
}
j = j + 1;
p = get_prime(j);
}
term3 = 1/(4*Pi) * term3;
/*** COMPUTE RHS ***/
RHS = term1 + term2 - term3;
/*** COMPUTE LHS ***/
LHS = 0;
//we will want to truncate this automatically
for (j=0;j<=number_zeros-1;j++)
{
u=zero_table[j]-x_0;
u=u*u*A;
if(u<2.3*DIGITS*2) LHS += exp(-u);
//LHS += exp(-A*(zero_table[j]-x_0)*(zero_table[j]-x_0));
}
/*** Display Results ***/
//if (my_verbose > 3)
if (my_verbose==0)
{
cout << endl << endl;
cout << "*** Testing Explicit Formula for L ***" << endl;
cout << "A = " << A << endl;
cout << "x_0 = " << x_0 << endl;
cout << "D = " << D << endl;
cout << "TERM 1: " << term1 << endl;
cout << "TERM 2: " << term2 << endl;
cout << "TERM 3: " << term3 << endl;
cout << "RHS: " << RHS << endl;
cout << "LHS: " << LHS << endl;
cout << "LHS - RHS: " << LHS-RHS << endl;
cout << "CUTOFF: " << pow(10.0, -DIGITS/3) << endl;
}
//cout << "A = " << A << ", x_0 = " << x_0 << ", ";
cout << " x_0=" << x_0 << ",";
cout << "LHS = " << LHS << ", RHS = " << RHS << ", DIFF = " << abs((LHS - RHS)) << ", ";
// I arbitrarily chose DIGITS/3.
if (abs((LHS - RHS)) < pow(10.0,-DIGITS/3))
{
flag = 0;
//cout << "PASS." << endl;
}
else
{
flag = 1;
//cout << "FAIL!" << endl;
}
return flag;
}
#endif
This diff is collapsed.
......@@ -29,6 +29,7 @@
#include "Lmisc.h"
#include <iomanip> //for manipulating output such as setprecision
#include <iostream> //for ostrstream
#include <cstring>
//using namespace std;
......@@ -43,7 +44,7 @@ Complex GAMMA (ttype z,ttype2 delta);
//computes G(z,w)
template <class ttype>
Complex inc_GAMMA (ttype z,ttype w, char *method="temme", ttype exp_w = 0, bool recycle=false); //computes G(z,w)
Complex inc_GAMMA (ttype z,ttype w, const char *method="temme", ttype exp_w = 0, bool recycle=false); //computes G(z,w)
template <class ttype>
ttype asympt_GAMMA (ttype z,ttype w, ttype exp_w = 0, bool recycle=false); //computes G(z,w) via asymptotic series
......@@ -67,7 +68,7 @@ Complex erfc(Complex z);
Complex erfc2(Complex z);
template <class ttype>
Complex gamma_sum(Complex s, int what_type,ttype *coeff, int N, Double g, Complex l, Double Q, Long Period, Complex delta=1, char *method="temme");
Complex gamma_sum(Complex s, int what_type,ttype *coeff, int N, Double g, Complex l, Double Q, Long Period, Complex delta=1, const char *method="temme");
//computes a sum of G(z,w)'s (as in (3.3.20)).
Complex exp_recycle();
......@@ -229,7 +230,7 @@ Complex GAMMA (ttype z1, ttype2 delta)
//value exp_w which holds exp(-w)
//computes G(z,w), so there's an extra w^(-z) factor.
template <class ttype>
Complex inc_GAMMA (ttype z,ttype w, char *method="temme", ttype exp_w = 0, bool recycle=false)
Complex inc_GAMMA (ttype z,ttype w, const char *method="temme", ttype exp_w = 0, bool recycle=false)
{
Complex G;
......@@ -240,7 +241,7 @@ Complex inc_GAMMA (ttype z,ttype w, char *method="temme", ttype exp_w = 0, bool
if(my_verbose>2) cout << "inc_GAMMA called. G(" << z<< " , " << w << ")" << endl;
//if(my_norm(z)<tolerance_sqrd){
if(my_norm(z)<.1){
if(my_norm(z)<.01){
return cfrac_GAMMA(z,w,exp_w,recycle);
}
......@@ -603,7 +604,7 @@ ttype comp_inc_GAMMA (ttype z,ttype w,ttype exp_w = 0, bool recycle=false) //co
}
template <class ttype>
Complex gamma_sum(Complex s, int what_type, ttype *coeff, int N, Double g, Complex l, Double Q, Long Period, Complex delta=1, char *method="temme")
Complex gamma_sum(Complex s, int what_type, ttype *coeff, int N, Double g, Complex l, Double Q, Long Period, Complex delta=1, const char *method="temme")
{
Complex SUM=0;
......@@ -703,7 +704,7 @@ Complex gamma_sum(Complex s, int what_type, ttype *coeff, int N, Double g, Compl
//which is more time consuming
if(is_z_real&&is_w_real){
G=inc_GAMMA(Double(real(z)),Double(real(w)),method,Double(real(exp_w)),true);
if (my_verbose>2) cout << "GAMMA SUM = " << G << endl;
if (my_verbose>2) cout << "GAMMA SUM with doubles = " << G << endl;
//cout<<"both z and w are real\n";
}
else{
......
......@@ -75,6 +75,9 @@ extern Double Pi;
extern Double log_2Pi;
extern Complex I;
extern bool only_use_dirichlet_series; //whether to compute just using the Dirichlet series
extern int N_use_dirichlet_series; //if so, how many terms in the Dirichlet series to use.
//-----Global variables----------------------------------------
extern int my_verbose; // verbosity level: 0 means no verbose
......@@ -92,6 +95,7 @@ extern int max_n; //the largest n used in a dirichlet series while computing a v
extern Double A; //controls the 'support' of g(w) in Riemann sum method
extern Double incr; //the increment in the Riemann sum method
extern Double tweak; //used in value_via_Riemann_sum to play with the angle
extern Double *LG; // lookup table for log(n)
extern Double *two_inverse_SQUARE_ROOT; // lookup table for sqrt(n)
......@@ -104,6 +108,71 @@ extern bool print_warning;
extern Long my_LLONG_MAX;
extern int *prime_table;
extern int number_primes;
// Riemann Siegel band limited interpolation ----------------------------
extern const Double sin_cof[];//={1.,-1./6.,1./120.,-1./5040.,1./362880.,-1./39916800.};
extern const Double sinh_mult_fac;
extern const Double sin_tol;
extern const int sin_terms;
extern const Double blfi_block_growth; // how fast blfi blocks grow as we traverse the main sum, keep as is for now
extern const Double beta_fac_mult; // controls the density of blfi sampling and the number of blfi terms needed
extern const Double blfi_fac; // efficiency of the blfi interpolation sum relative to an RS sum of same length
extern const Double pts_array_fac;
extern const int rs_blfi_N;
extern Double *klog0; //log(k) at the beginning
extern Double *klog2; //log(k) at the end if needed
extern Double *ksqrt0; // 1/sqrt(k) at the beginning
extern Double *ksqrt2;// 1/sqrt(k) at the end if needed
extern int *num_blocks; // number of blocks
extern int *size_blocks;// size of blocks
extern Double *trig; // stores correction terms
extern Double *zz; // store powers of fractional part
extern Double **klog1; //log(k) in the middle if needed
extern Double **ksqrt1; // 1/sqrt(k) in the middle if needed
extern Double **klog_blfi; //initial term
extern Double **qlog_blfi; //band-width
extern Double **piv_org; //original pivot
extern Double **bbeta; //beta
extern Double **blambda; //lambda
extern Double **bepsilon; //epsilon
extern Double **arg_blfi; //arg_blfi
extern Double **inv_arg_blfi; //inv_arg_blfi
extern Double