Commit 2e1c0cc0 authored by Lifeng Sun's avatar Lifeng Sun

Imported Upstream version 0.0.20080205

parents
This diff is collapsed.
To install:
in the src directory:
1) Edit the Makefile to reflect your machine. Right now it is set up
for two operating systems: mac os x (Darwin) , and Linux. Set OS_NAME
to reflect your machine. For example if your operating system is Darwin
set the line to:
OS_NAME = Darwin
otherwise set to:
OS_NAME = Linux
2) If you wish to have elliptic curve L-functions you need to
have pari installed on your machine and set the following in
the Makefile:
PREPROCESSOR_DEFINE = -DINCLUDE_PARI
Otherwise (default) just comment out the line:
#PREPROCESSOR_DEFINE = -DINCLUDE_PARI
In the former case (-DINCLUDE_PARI) you also need to make sure that
the following are set to reflect the correct location on your machine
of the pari header file and library. The default in the Makefile
for the location of pari is:
LOCATION_PARI_H = /usr/local/include/pari
LOCATION_PARI_LIBRARY = /usr/local/lib
You will have to modify these lines if pari is installed in
different locations on your machine.
3) at the unix prompt type: make
(If you wish to recompile later, type: make clean
before retyping: make)
This will create the command line application called lcalc
and also a library, libLfunction.a, that can be linked to if
you wish to write your own code. See the README file for usage
and a mini tutorial.
4) If you have root powers, you can do: sudo make install
to copy the application, library, and include files to the
appropriate directories in /usr/local/
This diff is collapsed.
This diff is collapsed.
/*
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.
*/
#ifndef Lcommandline_H
#define Lcommandline_H
//---------------------------------------------------
//
// Command line interface to the L function package
// By Michael Rubinstein
//
//---------------------------------------------------
#include <stdlib.h> //for calling unix commands and things like atoi
#include <string> //for string manipulation
#include "L.h" //for the L_function class
#include "Lcommandline_numbertheory.h" //for number theory functions
#include "Lcommandline_globals.h" //command line global variables
#ifdef INCLUDE_PARI
#include "pari.h" //for pari's elliptic curve functions
#undef init //pari has a '#define init pari_init' which
//causes trouble with the stream.h init.
//pari also causes trouble with things like abs.
//we place the pari include first since otherwise it
//messes up.
#endif //ifdef INCLUDE_PARI
#include "Lcommandline_misc.h"
#include "Lcommandline_elliptic.h" //mainly for initializing an elliptic curve L-function using pari
#include "Lcommandline_twist.h" //for twisting by Dirichlet characters
#include "Lcommandline_values_zeros.h" //for calling value and zero routines
using namespace std;
#endif
/*
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.
*/
#ifndef Lcommandline_elliptic_H
#define Lcommandline_elliptic_H
#include "L.h" //for the L_function class
#include "Lcommandline_numbertheory.h" //for number theory functions
#include "Lcommandline_globals.h" //command line global variables
#include <iostream>
#ifdef INCLUDE_PARI
#include "pari.h" //for pari's elliptic curve functions
#undef init //pari has a '#define init pari_init' which
//causes trouble with the stream.h init.
//pari also causes trouble with things like abs.
//we place the pari include first since otherwise it
//messes up.
void initialize_new_L(char *a1, char *a2, char *a3, char *a4, char *a6, int N_terms);
void data_E(char *a1, char *a2, char *a3, char *a4, char *a6,int N_terms,Double * coeff);
#endif //ifdef INCLUDE_PARI
void verify_rank(int rank);
void compute_rank();
#endif
/*
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.
*/
#ifndef Lcommandline_globals_H
#define Lcommandline_globals_H
#include "L.h"
//-----Global variables--------------------------------
extern int current_L_type; //1,2,3:For the 3 diff types of L-functions: int, Double, Complex
extern L_function<int> int_L,int_L2,int_L3;
extern L_function<Double> Double_L,Double_L2,Double_L3;
extern L_function<Complex> Complex_L,Complex_L2,Complex_L3;
void initialize_commandline_globals();
#endif
/*
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.
*/
#ifndef Lcommandline_misc_H
#define Lcommandline_misc_H
#include <stdlib.h> //for calling unix commands and things like atoi
#include <iostream>
#include <iomanip> //for manipulating output such as setprecision
#include <fstream> //for file input output
#include <strstream> //for ostream:w
#include <iostream> //for ostrstream
#include "L.h" //for the L_function class
#include "Lcommandline_numbertheory.h" //for number theory functions
#include "Lcommandline_globals.h" //command line global variables
void initialize_new_L(char *file_name);
void print_L(int N_terms);
#endif
/*
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.
*/
#ifndef Lcommandline_numbertheory_H
#define Lcommandline_numbertheory_H
#include <stdlib.h> //for things like srand
#include <iostream> //for input and output
#include <iomanip> //for manipulating output such as setprecision
#include <math.h>
//#include <ctime>
#include <limits.h> // for INT_MAX
#include "Lglobals.h"
using namespace std;
long long nextprime(long long n);
bool issquarefree(long long n);
bool isfunddiscriminant(long long n);
long long nextfunddiscriminant(long long n);
int simplified_jacobi(int n,int m);
int simplified_jacobi(long long n,long long m);
int my_kronecker(int n,int m);
int my_kronecker(long long n, long long m);
Double L_1_chi(long long d);
int class_number(long long d);
void ramanujan_tau(Double *c0, int N_terms);
long long gcd(long long a,long long b);
void factor(long long q, long long **factors);
long long power_mod_q(long long a, long long k,long long q);
int prim_root(long long p, int alpha);
void factor(long long q, long long **factors);
int characters();
bool isprime(long long n); //by Kyle Wichert
bool RM(long long a, long long N); //by Kyle Wichert
long long multmodN(long long a, long long b, long long N); //by Kyle Wichert
template <class ttype>
Complex gauss_sum(ttype *chi,long long r,bool cnj=false)
{
Complex SUM=0;
if(cnj==true)
for(int n=1;n<=r; n++) SUM=SUM+conj(chi[n])*exp(n*2*I*Pi/double(r));
else
for(int n=1;n<=r; n++) SUM=SUM+chi[n]*exp(n*2*I*Pi/double(r));
return SUM;
}
#endif
/*
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.
*/
#ifndef Lcommandline_twist_H
#define Lcommandline_twist_H
#include "L.h" //for the L_function class
#include "Lcommandline_numbertheory.h" //for number theory functions
#include "Lcommandline_globals.h" //command line global variables
#include "Lcommandline_elliptic.h"
//===============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);
#endif
/*
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.
*/
#ifndef Lcommandline_values_zeros_H
#define Lcommandline_values_zeros_H
#include "L.h" //for the L_function class
#include "Lcommandline_numbertheory.h" //for number theory functions
#include "Lcommandline_globals.h" //command line global variables
//-----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 L_interpolate(Double x, Double y,Double step_size, int n=1000);
#endif
// When MPFR is enabled and double is passed to a templated function
// The function should use precise(ttype) to make sure calculations run
// within the more precise type
#define precise(T) typename precise<T>::precise_type
template<class T> struct precise { typedef T precise_type; };
template<> struct precise<double> { typedef Double precise_type; };
template<> struct precise<complex<double> > { typedef Complex precise_type; };
#ifdef USE_MPFR
template<> struct precise<long long> { typedef double precise_type; };
#else
template<> struct precise<long long> { typedef long long precise_type; };
#endif
//typedef precise<long long>::precise_type Long;
typedef long long Long;
// To aid conversion to int value
inline double to_double(const Double& x) {
#ifdef USE_MPFR
return x.get_d();
#else
return x;
#endif
}
#ifdef USE_MPFR
inline double 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))
template<class T> inline int sn(T x)
{
if (x>=0) return 1;
else return -1;
}
const bool outputSeries=true; // Whether to output the coefficients or just the answer
// Loop i from m to n
// Useful in tidying up most for loops
#define loop(i,m,n) for(typeof(m) i=(m); i!=(n); i++)
// A class for calculations involving polynomials of small degree
// Not efficient enough for huge polynomials
//
// Polynomial of fixed degree N-1 with coefficients over T
template<class T=Complex> struct smallPoly {
valarray<T> coeffs;
typedef smallPoly<T> poly;
int N;
T zero;
smallPoly(int N=0) {
this->N=N;
coeffs.resize(N);
zero=0;
}
void resize(int N) {
coeffs.resize(N);
loop(i,this->N,N)
coeffs[i]=zero;
this->N=N;
}
// Access a coefficient as a subscript
T& operator[](int n) {
return (n<0 ? zero : coeffs[n]);
}
const T& operator[](int n) const {
return (n<0 ? zero : coeffs[n]);
}
// Multiply two polys together, truncating the result
friend poly operator*(const poly& f, const poly& g) {
poly result(max(f.N,g.N));
loop(i,0,f.N) result[i]=0;
double test;
loop(i,0,f.N) loop(j,0,result.N-i) {
result[i+j]+=f[i]*g[j];
}
return result;
}
// Divide (in place) by (x-1/a) = (1-ax)/(-a) with a!=0
void divideSpecial(T a) {
loop(i,1,N) coeffs[i]+=coeffs[i-1]*a;
coeffs*= -a;
}
// Evaluate the polynomial at x
template<class U> U operator()(U x) {
U sum=coeffs[0];
U cur=1;
loop(i,1,N) {
cur*=x;
sum+=coeffs[i]*x;
}
return sum;
}
// Output to stdout
friend ostream& operator<<(ostream& s, const poly p) {
loop(i,0,p.N) s << (i ? " + " : "") << p[i] << "x^" << i;
return s;
}
// Arithmetic
template<class U> poly& operator*=(const U& x) {
coeffs*=x;
return *this;
}
template<class U> poly& operator/=(const U& x) {
coeffs/=x;
return *this;
}
};
// When MPFR is enabled and double is passed to a templated function
// The function should use precise(ttype) to make sure calculations run
// within the more precise type
#define precise(T) typename precise<T>::precise_type
template<class T> struct precise { typedef T precise_type; };
template<> struct precise<double> { typedef Double precise_type; };
template<> struct precise<complex<double> > { typedef Complex precise_type; };
#ifdef USE_MPFR
template<> struct precise<long long> { typedef double precise_type; };
#else
template<> struct precise<long long> { typedef long long precise_type; };
#endif
//typedef precise<long long>::precise_type Long;
typedef long long Long;
// To aid conversion to int value
inline long double to_double(const Double& x) {
#ifdef USE_MPFR
return x.get_d();
#else
return x;
#endif
}
#ifdef USE_MPFR
inline long double to_double(const double& x) { return x; }
inline long double to_double(const long double& x) { return x; }
#endif
inline long double to_double(const double& x) { return x; }
inline long double to_double(const int& x) { return x; }
inline long double to_double(const long long& x) { return x; }
inline long double to_double(const short& x) { return x; }
inline long 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 long double pow(long double x, double y) {
return pow(x,(long double)y);
}
inline long double pow(double x, long double y) {
return pow((long double)x,y);
}
template<class T> inline int sn(T x)
{
if (x>=0) return 1;
else return -1;
}
const bool outputSeries=true; // Whether to output the coefficients or just the answer
// Loop i from m to n
// Useful in tidying up most for loops
#define loop(i,m,n) for(typeof(m) i=(m); i!=(n); i++)
// A class for calculations involving polynomials of small degree
// Not efficient enough for huge polynomials
//
// Polynomial of fixed degree N-1 with coefficients over T
template<class T=Complex> struct smallPoly {
valarray<T> coeffs;
typedef smallPoly<T> poly;
int N;
T zero;
smallPoly(int N=0) {
this->N=N;
coeffs.resize(N);
zero=0;
}
void resize(int N) {
coeffs.resize(N);
loop(i,this->N,N)
coeffs[i]=zero;
this->N=N;
}
// Access a coefficient as a subscript
T& operator[](int n) {
return (n<0 ? zero : coeffs[n]);
}
const T& operator[](int n) const {
return (n<0 ? zero : coeffs[n]);
}
// Multiply two polys together, truncating the result
friend poly operator*(const poly& f, const poly& g) {
poly result(max(f.N,g.N));
loop(i,0,f.N) result[i]=0;
double test;
loop(i,0,f.N) loop(j,0,result.N-i) {
result[i+j]+=f[i]*g[j];
}
return result;
}
// Divide (in place) by (x-1/a) = (1-ax)/(-a) with a!=0
void divideSpecial(T a) {
loop(i,1,N) coeffs[i]+=coeffs[i-1]*a;
coeffs*= -a;
}
// Evaluate the polynomial at x
template<class U> U operator()(U x) {
U sum=coeffs[0];
U cur=1;
loop(i,1,N) {
cur*=x;
sum+=coeffs[i]*x;
}
return sum;
}
// Output to stdout
friend ostream& operator<<(ostream& s, const poly p) {
loop(i,0,p.N) s << (i ? " + " : "") << p[i] << "x^" << i;
return s;
}
// Arithmetic
template<class U> poly& operator*=(const U& x) {
coeffs*=x;
return *this;
}
template<class U> poly& operator/=(const U& x) {
coeffs/=x;
return *this;
}
};
This diff is collapsed.
/*
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 "Lmisc.h"
template <class ttype>
Complex L_function <ttype>::
partial_dirichlet_series(Complex s, long long N1, long long N2)
{
Complex z=0.;
long long m,n;
if(what_type_L==-1) //i.e. if the Riemann zeta function
for(n=N1;n<=N2;n++) z=z+exp(-s*LOG(n));
else if(what_type_L!=1) //if not periodic
for(n=N1;n<=N2;n++) z=z+dirichlet_coefficient[n]*exp(-s*LOG(n));
else //if periodic
for(n=N1;n<=N2;n++)
{
m=n%period; if(m==0)m=period;
z=z+dirichlet_coefficient[m]*exp(-s*LOG(n));
}
return z;
}
//XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
template <class ttype>
Complex L_function <ttype>::
dirichlet_series(Complex s, long long N=-1)
{
Complex z=0.;
long long m,n;
if(N==-1) N=number_of_dirichlet_coefficients;
if(N>number_of_dirichlet_coefficients&&what_type_L!=-1&&what_type_L!=1)
{
if(print_warning){
print_warning=false;
cout << "WARNING from dirichlet series- we don't have enough Dirichlet coefficients." << endl;
cout << "Will use the maximum possible, though the output ";
cout << "will not necessarily be accurate."