liblfunction-dev 1.23+dfsg-5 / usr / include / Lfunction / Lvalue.h

/*

   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.

*/

    template <class ttype>
    Complex L_function <ttype>::
    find_delta (Complex z,Double g)
    {

        //cout << " find delta z g " << z << " " << g<< endl;
        Double sigma=real(z);
        Double t=imag(z); if(t<0) t=-t;
        Double r=abs(z);
        Double theta=atan(t/sigma);
        Double epsilon;

        Double a=-theta;
        Double b=0;
        Double c,f1,f2,f3;
        Double local_tolerance=.01/(t+100); if (local_tolerance<tolerance) local_tolerance=tolerance;

        f1=sigma*log(sigma/(r*cos(theta+a))) - t*a;
        //f2=0;
        if(f1<=DIGITS2*2.3) epsilon=-theta;
        else{
            do{
                c=(a+b)/2;
                f3=sigma*log(sigma/(r*cos(theta+c))) - t*c;
                if(f3>DIGITS2*2.3)a=c;
                else b=c;
                //cout<< "theta+epsilon: " << a+theta<< endl;
            } while(b-a>local_tolerance);
            epsilon=a;
        }
        //if(imag(z)>=0) cout << " returning delta: " << exp(I*(theta+epsilon)*g) << endl;
        //else cout << " returning delta: " << exp(-I*(theta+epsilon)*g) << endl;
        if(imag(z)>=0) return exp(I*(theta+epsilon)*g);
        else return exp(-I*(theta+epsilon)*g);
   }


    //computes (3.2.5) as a Riemann sum using
    //g(w) = exp(A*(w-s)^2) * delta^(-w)

    template <class ttype>
    Complex L_function <ttype>::
    value_via_Riemann_sum(Complex s, const char *return_type)
    {
        int j,k,m,mm,n;
        Complex r,z;
        Complex SUM=0;

        Double tmp;

        Complex L_value=0;
        Double theta;

        Complex *DELTA; // variant on (3.3.10), without the theta

        Double t_0;
        Double c;    //controls speed of convergence but at the expense
                     //of loss of precision.

        //Double v=.5; // the v in (3.2.5)
        Double v=1-real(s); // the v in (3.2.5)
        if(v<.5)v=.5;
        //Double incr; // incrememnt size in the Riemann sum- is now a global variable

        Complex dirichletseries;
        Complex *dirichlet_vector; //used to compute Dirichlet series incrementally
        Complex *dirichlet_vector_copy; //used to compute Dirichlet series incrementally
        Complex *dirichlet_multiplier; // stores powers of n^{-I incr}
        Complex *dirichlet_multiplier_conj; // stores powers of n^{I incr}. Prefer to store since
                                            // this will make vectorization easier.

        //int M; //number of terms to take in the Riemann sum
        //M is no longer used. escape is determined numerically. Is fine to
        //escape this way. Gamma and exp behave predictably

        int N; //the number of terms to take in the Dirichlet series
        Double r1=0,r2=0,r3=0,r4=0,r5=0,mynorm_r,local_average, max_integrand; //used to decide when to truncate the Riemann sum


        //cout << "a= " << a << endl;
        theta=0;
        for(j=1;j<=a;j++)
        {
            theta = theta+gamma[j];
            //cout << "theta = " << theta << endl;
        }

        c=DIGITS2*log(10.); //i.e exp(-c)=10^(-DIGITS2)
                            //this sacrifices roughly at most DIGITS2
                            //out of DIGITS precision.


        //incr is now a global variable
        //incr=2*Pi*v/(log(10.)*DIGITS);

        //M=Int(sqrt(log(10.)*DIGITS/A)/incr); 



        DELTA = new Complex[a+1];

        //if(abs(t_0)<=c/(2*sqrt(A))) tmp= 2*sqrt(A);
        //else tmp=c/t_0;


        Double c1=0;
        for(j=1;j<=a;j++){

            t_0=imag(gamma[j]*s+lambda[j]);

        //cout << "tmp_" << j << " =   " << tmp  << endl;
        //cout << "t_0" << " =   " << t_0  << endl;

            if(abs(t_0)<=2*c/(Pi*a)) tmp= Pi/2;
            else tmp=abs(c/(t_0*a));

            if(t_0>=0)r=1; else r=-1;

            DELTA[j]= exp(I*r*(Pi/2-tweak*tmp)); //tweak defaults to 1. It allows me to globally set a slightly different angle
                                                 //for the purpose of testing precision or looking for L-functions

            c1=c1+gamma[j]*tmp;

            //DELTA[j]=find_delta(s*gamma[j]+lambda[j],gamma[j]);


        }


        for(k=1;k<=number_of_poles;k++){
           z =A*(pole[k]-s)*(pole[k]-s);
           for(j=1;j<=a;j++) z=z-log(DELTA[j])*(gamma[j]*pole[k]+lambda[j]);
           //the 5 below is for kicks. 2.3 would have been fine.
           if(real(z)>-5*DIGITS)
               L_value=L_value+residue[k]*exp(z)/(s-pole[k]);
        }
        //cout << "poles contribute: " << L_value << endl;

        //the rough estimate: G(z,(N DELTA/Q)^2) is, in size,
        //roughly exp(-Re((N*DELTA/Q)^(1/theta))) and we want this
        //to be > 2.3 DIGITS XXXXXXXXX check this

        N=Int(Q*exp(log(2.3*DIGITS*theta/c1)*theta)+10);


        if(N>number_of_dirichlet_coefficients&&what_type_L!=-1&&what_type_L!=1)
        {

            if(print_warning){
                print_warning=false;
                cout << "WARNING from Riemann sum- we don't have enough Dirichlet coefficients." << endl;
                cout << "Will use the maximum possible, though the output ";
                cout << "will not necessarily be accurate." << endl;
            }
            N=number_of_dirichlet_coefficients;
        }

        if(N>number_logs) extend_LG_table(N);

        dirichlet_vector= new Complex[N+1]; //initially stores a(n)/n^{s+v} (or 1-s instead of s)
        dirichlet_vector_copy= new Complex[N+1]; // used to do negative m
        dirichlet_multiplier= new Complex[N+1];
        dirichlet_multiplier_conj= new Complex[N+1];

        #pragma omp parallel for //shared(N,I,incr,dirichlet_multiplier) private(n)
        for(n=1;n<=N;n++){
            dirichlet_multiplier[n]=exp(-I*LOG(n)*incr);
            dirichlet_multiplier_conj[n]=conj(dirichlet_multiplier[n]);

            if(what_type_L==-1)   //i.e. if the Riemann zeta function
                dirichlet_vector[n]=exp(-(s+v)*LOG(n));
            else if(what_type_L!=1) //if not periodic
                dirichlet_vector[n]=dirichlet_coefficient[n]*exp(-(s+v)*LOG(n));
            else //if periodic
            {
                m=n%period; if(m==0) m=period;
                dirichlet_vector[n]=dirichlet_coefficient[m]*exp(-(s+v)*LOG(n));
            }
            dirichlet_vector_copy[n]=dirichlet_vector[n];
        }


        max_n=N;

        if(my_verbose>1)
            cout << "s =  " << s << "  Will use  " << N << " terms of the Dirichlet series" << endl;

        Double log_Q=log(Q);

/* old riemann sum. escape was fixed ahead of time and was not efficient.
        for(m=-M;m<=M;m++){
            r=exp(A*(v+I*incr*m)*(v+I*incr*m)+log_Q*(s+v+I*incr*m));
            for(j=1;j<=a;j++){
                //cout << "gamma[j]*(s+v+I*incr*m)+lambda[j]= " << gamma[j]*(s+v+I*incr*m)+lambda[j] << endl;
                //cout << "DELTA[j] =   " << DELTA[j]  << endl;
                r=r*GAMMA(gamma[j]*(s+v+I*incr*m)+lambda[j],DELTA[j]);
            }
            //cout << "r= " << r << endl;

            //r=r*this->dirichlet_series(s+v+I*incr*m,N,false);

            r=r*this->dirichlet_series(s+v+I*incr*m,N);
            //r=r*this->dirichlet_series(s+v+I*incr*m,N);

            //dirichlet series part needs to be more precise
            SUM=SUM+r/(v+I*incr*m);
            //if(m%100==0)
            cout << "m= " << m << "  SUM1 = " << SUM << endl;
        }
        SUM=SUM*incr/(2*Pi);
        //cout << "m= " << m << "  SUM1 = " << SUM << endl;
*/

        max_integrand=0; //the max of the integrand, without the dirichlet series factor
        mm=0;
        //first do the terms m >=0
        do{
            for(m=mm;m<=mm+99;m++){
                r=exp(A*(v+I*incr*m)*(v+I*incr*m)+log_Q*(s+v+I*incr*m));
                for(j=1;j<=a;j++){
                    r=r*GAMMA(gamma[j]*(s+v+I*incr*m)+lambda[j],DELTA[j]);
                }
                r=r/(v+I*incr*m);

                mynorm_r=my_norm(r);
                if(mynorm_r>max_integrand) max_integrand = mynorm_r;

                r1=r2;r2=r3;r3=r4;r4=r5;r5=mynorm_r;
                local_average=(r1+r2+r3+r4+r5)/5;

                //XXXXXX replaced by vectorized version
                //r=r*this->dirichlet_series(s+v+I*incr*m,N);

                dirichletseries=0;
                for(j=1;j<=N;j++){
                    dirichletseries=dirichletseries+dirichlet_vector[j];
                    dirichlet_vector[j]=dirichlet_vector[j]*dirichlet_multiplier[j];
                }
                r=r*dirichletseries;
                //cout << "1 dirichletseries: " << dirichletseries << endl;
                //cout << "1 thischletseries: " << this->dirichlet_series(s+v+I*incr*m,N) << endl;


                SUM=SUM+r;
                //cout << "m= " << m << "  SUM = " << SUM << "  r= " << r  << "  local average = " << local_average << " max=  "<< max_integrand<< endl;

            }
            mm=m;
        }while(local_average>max_integrand*tolerance_sqrd);

        m=-1;
        //then do the terms negative m
        do{
            r=exp(A*(v+I*incr*m)*(v+I*incr*m)+log_Q*(s+v+I*incr*m));
            for(j=1;j<=a;j++){
                r=r*GAMMA(gamma[j]*(s+v+I*incr*m)+lambda[j],DELTA[j]);
            }
            r=r/(v+I*incr*m);

            mynorm_r=my_norm(r);
            if(mynorm_r>max_integrand) max_integrand = mynorm_r;
            r1=r2;r2=r3;r3=r4;r4=r5;r5=mynorm_r;
            local_average=(r1+r2+r3+r4+r5)/5;

            //r=r*this->dirichlet_series(s+v+I*incr*m,N);

            dirichletseries=0;

            for(j=1;j<=N;j++){
                dirichlet_vector_copy[j]=dirichlet_vector_copy[j]*dirichlet_multiplier_conj[j];
                dirichletseries=dirichletseries+dirichlet_vector_copy[j];
            }
            r=r*dirichletseries;
            //cout << "2 dirichletseries: " << dirichletseries << endl;
            //cout << "2 thischletseries: " << this->dirichlet_series(s+v+I*incr*m,N) << endl;


            SUM=SUM+r;
            //cout << "m= " << m << "  SUM = " << SUM << "  r= " << r  << "  local average = " << local_average << " max=  "<< max_integrand<< endl;
            m--;
        }while(m>-100||local_average>max_integrand*tolerance_sqrd);

        SUM=SUM*incr/(2*Pi);


        //no longer needed
        //r=0;
        //for(j=1;j<=a;j++)r=r+lambda[j];
        //SUM=SUM*exp(log(DELTA)*r);

        //cout << "m= " << m << "  SUM1 = " << SUM << endl;

        L_value=L_value+SUM;


        if(real(s)!=.5){ //do the second sum i.e. for f_2

            v=real(s);

            if(what_type_L==-1)   //i.e. if the Riemann zeta function
            {
                #pragma omp parallel for shared(N,dirichlet_vector,s,v) private(n)
                for(n=1;n<=N;n++) dirichlet_vector[n]=exp(-conj(1-s+v)*LOG(n));
            }
            else if(what_type_L!=1) //if not periodic
            {
                #pragma omp parallel for shared(N,dirichlet_vector,s,v) private(n)
                for(n=1;n<=N;n++) dirichlet_vector[n]=dirichlet_coefficient[n]*exp(-conj(1-s+v)*LOG(n));
            }
            else //if periodic
            {
                for(n=1;n<=N;n++)
                {
                    m=n%period; if(m==0)m=period;
                    dirichlet_vector[n]=dirichlet_coefficient[m]*exp(-conj(1-s+v)*LOG(n));
                }
            }

            for(n=1;n<=N;n++) dirichlet_vector_copy[n]=dirichlet_vector[n];

            SUM=0;


/*
            for(m=-M;m<=M;m++){
                r=exp(A*(v+I*incr*m)*(v+I*incr*m)+log_Q*(1-s+v+I*incr*m));
                for(j=1;j<=a;j++)
                        r=r*GAMMA(gamma[j]*(1-s+v+I*incr*m)+conj(lambda[j]),1/DELTA[j]);
                //r=r*conj(this->dirichlet_series(conj(1-s+v+I*incr*m),N,false));

                r=r*conj(this->dirichlet_series(conj(1-s+v+I*incr*m),N));

                //dirichlet series part needs to be more precise
                SUM=SUM+r/(v+I*incr*m);
                //if(m%100==0)
                //cout << "m= " << m << "  SUM2 = " << SUM << endl;
            }
            SUM=SUM*incr/(2*Pi);
*/

            max_integrand=0; //the max of the integrand, without the dirichlet series factor
            m=0;
            //first do the terms m >=0
            do{
                r=exp(A*(v+I*incr*m)*(v+I*incr*m)+log_Q*(1-s+v+I*incr*m));
                for(j=1;j<=a;j++){
                    r=r*GAMMA(gamma[j]*(1-s+v+I*incr*m)+conj(lambda[j]),1/DELTA[j]);
                }
                r=r/(v+I*incr*m);


                mynorm_r=my_norm(r);
                if(mynorm_r>max_integrand) max_integrand = mynorm_r;
                r1=r2;r2=r3;r3=r4;r4=r5;r5=mynorm_r;
                local_average=(r1+r2+r3+r4+r5)/5;

                //r=r*conj(this->dirichlet_series(conj(1-s+v+I*incr*m),N));
                dirichletseries=0;
                for(j=1;j<=N;j++){
                    dirichletseries=dirichletseries+dirichlet_vector[j];
                }
                for(j=1;j<=N;j++){
                    dirichlet_vector[j]=dirichlet_vector[j]*dirichlet_multiplier_conj[j];
                }
                r=r*conj(dirichletseries);
                //cout << "3 dirichletseries: " << dirichletseries << endl;
                //cout << "3 thischletseries: " << this->dirichlet_series(conj(1-s+v+I*incr*m),N) << endl;

                SUM=SUM+r;
                //cout << "m= " << m << "  SUM = " << SUM << "  r= " << r  << "  local average = " << local_average << " max=  "<< max_integrand<< endl;
                m++;
            }while(m<100||local_average>max_integrand*tolerance_sqrd);

            m=-1;
            //then do the terms negative m
            do{
                r=exp(A*(v+I*incr*m)*(v+I*incr*m)+log_Q*(1-s+v+I*incr*m));
                for(j=1;j<=a;j++){
                    r=r*GAMMA(gamma[j]*(1-s+v+I*incr*m)+conj(lambda[j]),1/DELTA[j]);
                }
                r=r/(v+I*incr*m);

                mynorm_r=my_norm(r);
                if(mynorm_r>max_integrand) max_integrand = mynorm_r;
                r1=r2;r2=r3;r3=r4;r4=r5;r5=mynorm_r;
                local_average=(r1+r2+r3+r4+r5)/5;

                //r=r*conj(this->dirichlet_series(conj(1-s+v+I*incr*m),N));
                dirichletseries=0;
                for(j=1;j<=N;j++){
                    dirichlet_vector_copy[j]=dirichlet_vector_copy[j]*dirichlet_multiplier[j];
                }
                for(j=1;j<=N;j++){
                    dirichletseries=dirichletseries+dirichlet_vector_copy[j];
                }
                r=r*conj(dirichletseries);
                //cout << "4 dirichletseries: " << dirichletseries << endl;
                //cout << "4 thischletseries: " << this->dirichlet_series(conj(1-s+v+I*incr*m),N) << endl;

                SUM=SUM+r;
                //cout << "m= " << m << "  SUM = " << SUM << "  r= " << r  << "  local average = " << local_average << " max=  "<< max_integrand<< endl;
                m--;
            }while(m>-100||local_average>max_integrand*tolerance_sqrd);

            SUM=SUM*incr/(2*Pi);



         }
        else SUM =conj(SUM);

        for(j=1;j<=a;j++){
            r=-gamma[j]-2*real(lambda[j]);
            SUM=SUM*exp(log(DELTA[j])*r);
        }

        //cout << "m= " << m << "  SUM2 = " << SUM << endl;
        //cout << "r= " << r <<  endl;

        L_value=L_value+OMEGA*SUM;

        delete [] dirichlet_vector;
        delete [] dirichlet_vector_copy;
        delete [] dirichlet_multiplier;
        delete [] dirichlet_multiplier_conj;

        //this returns L(s)
        if (!strcmp(return_type,"pure"))
        {
            z=1;
            for(j=1;j<=a;j++)
                z=z*GAMMA(gamma[j]*s+lambda[j],DELTA[j]);
            //cout << "pure  " << L_value*exp(-log(Q)*s)/z << endl;
            delete [] DELTA;
            return L_value*exp(-log(Q)*s)/z;
        }

        //returns L(s) rotated to be real on critical line
        //assumes |OMEGA|=1. Valid assumption since
        //LAMBDA(1/2+it) = OMEGA conj(LAMBDA(1/2+it))
        else if (!strcmp(return_type,"rotated pure"))
        {
            r=1;
            for(j=1;j<=a;j++)
                r=r*GAMMA(gamma[j]*s+lambda[j],DELTA[j]);
            z=0;
            for(j=1;j<=a;j++)
                z=z+log(DELTA[j])*real(gamma[j]*s+lambda[j]);
            //cout << "rotated pure  " <<  L_value*exp(-log(Q)*real(s)-.5*log(OMEGA))*exp(z)/abs(r) << endl;
            delete [] DELTA;
            return L_value*exp(-log(Q)*real(s)-.5*log(OMEGA))*exp(z)/abs(r);
        }

        //else return Lambda(s) OMEGA^(-1/2) delta^(Re(s)).
        //This returns a real number (though, as a Complex)
        //on the critical line assuming |OMEGA|=1. Valid assumption
        //since LAMBDA(1/2+it) = OMEGA conj(LAMBDA(1/2+it))
        else if(!strcmp(return_type,"normalized and real"))
        {
            z=0;
            for(j=1;j<=a;j++)
                z=z+log(DELTA[j])*real(gamma[j]*s+lambda[j]);
            //cout << "normalized and real  " << L_value*exp(z-.5*log(OMEGA)) << endl;
            delete [] DELTA;
            return L_value*exp(z-.5*log(OMEGA));
        }
    }


    // implements (3.3.20) with no precomputations.
    // DIGITS is how much precision we would like.
    // DIGITS2 is how much precision (out of DIGITS)
    // we are willing to sacrifice for the sake of
    template <class ttype>
    Complex L_function <ttype>::
    value_via_gamma_sum(Complex s, const char *return_type)
    {
        Complex L_value=0;
        Double theta=gamma[1];  // equals gamma_1
        Double t_0=imag(s);
        Double c;    //controls speed of convergence but at the expense
                     //of loss of precision.
        Complex DELTA;  //(3.3.10)

        Complex u;
        int k;

        c=DIGITS2*log(10.)/theta; //i.e exp(-c theta)=10^(-DIGITS2)
                                  //this sacrifices roughly at most DIGITS2
                                  //out of DIGITS precision.

        //if(abs(t_0)<=2*c/Pi) DELTA=1;
        //else if(t_0>=0) DELTA = exp(I*theta*(Pi/2-c/t_0));
        //else            DELTA = exp(I*theta*(-Pi/2-c/t_0));
        DELTA=find_delta(s*gamma[1]+lambda[1],gamma[1]);


        u=log(DELTA);
        for(k=1;k<=number_of_poles;k++)
           L_value=L_value+residue[k]*exp(-u*pole[k])/(s-pole[k]);



        u=gamma_sum(s, what_type_L, dirichlet_coefficient,
            number_of_dirichlet_coefficients, gamma[1], lambda[1],
            Q, period, DELTA);


        L_value=L_value+exp(log(DELTA/Q)*lambda[1]/gamma[1])*u;


        if(real(s)!=.5)
            u=gamma_sum(1-conj(s), what_type_L, dirichlet_coefficient,
            number_of_dirichlet_coefficients,gamma[1],lambda[1],Q,period,DELTA);
        u=conj(u);


        L_value=L_value+(OMEGA/DELTA)*exp(-log(DELTA*Q)*conj(lambda[1])/gamma[1])*u;

        //this returns L(s)
        if (!strcmp(return_type,"pure"))
        {
            u=log(DELTA/Q)/gamma[1];
//cout << "returning " << L_value << " divided by  " << (GAMMA(gamma[1]*s+lambda[1],exp(u))*exp(u*lambda[1])) << endl;//XXXXXXXXXXXXXXXXXX
            return L_value/(GAMMA(gamma[1]*s+lambda[1],exp(u))*exp(u*lambda[1]));
        }

        //returns L(s) rotated to be real on critical line
        //assumes |OMEGA|=1.
        else if (!strcmp(return_type,"rotated pure"))
        {
            u=log(DELTA/Q)/gamma[1];
            u=abs(GAMMA(gamma[1]*s+lambda[1],exp(u))*exp(u*lambda[1]));
            //u=GAMMA(gamma s + lambda)*(delta/Q)^(-s)
            return L_value*exp(log(DELTA)*real(s)-.5*log(OMEGA))/u;
        }

        //else return Lambda(s) OMEGA^(-1/2) delta^(Re(s)).
        //This returns a real number (though, as a Complex) 
        //on the critical line assuming |OMEGA|=1
        else if(!strcmp(return_type,"normalized and real"))
            return L_value*exp(log(DELTA)*real(s)-.5*log(OMEGA));

        else // return L(s)
        {
            u=log(DELTA/Q)/gamma[1];
            return L_value/(GAMMA(gamma[1]*s+lambda[1],exp(u))*exp(u*lambda[1]));
        }

    }


    template <class ttype>
    Complex L_function <ttype>::
    value(Complex s, int derivative, const char *return_type)
    {
      Complex L_value;


      if(derivative==0){

          if(my_verbose>1)
              cout << "calling L:  " << s << endl;

           cout << setprecision(DIGITS3);

           if(only_use_dirichlet_series){
               L_value= this->dirichlet_series(s,N_use_dirichlet_series);
               return L_value;
           }


          //if(what_type_L==-1&&real(s)==.5&&abs(imag(s))>500) return Zeta(s,return_type);

          //uses Riemann Siegel. This is good only up to limited precision.
          //last condition in the if takes into account that Riemann Sigel is an asymptotic expansion
          //and that currently we only use the first 5 terms which gives O(t^(-3)) for the remainder
          if(what_type_L==-1&&real(s)==.5&&log(abs(imag(s)))/2.3>DIGITS/3.){

               int success;
               Double error_tol=1E-30;

               if(my_verbose==-33){
                   //cout << " rs blfi " ;
                   L_value= rs(imag(s),error_tol,input_mean_spacing_given,success,return_type);
               }
               else{
                   //cout << " rs " ;
                   L_value = Zeta(s,return_type);
               }

               //1.7725 is Pi^(.5), to account for the Q^\pm s in the approximate functional eqn
               DIGITS3=Int((DIGITS-log(log(1.*max_n*1.7725+3)*abs(imag(s))/6.28+3)/2.3)/pow(2.,abs(global_derivative)))+2;
               cout << setprecision(DIGITS3);
               if (my_verbose>1) cout << "Setting output precision to: " << DIGITS3 << endl;
               tolerance3=pow(.1,(DIGITS3+1));
               //cout << setprecision(DIGITS);
               //cout << s << " riemann siegel (" << DIGITS3 << " ): " <<  L_value << endl;
               return L_value;
          }


         if(a==1){
              //cout << "gamma sum " << endl;
              L_value = this->value_via_gamma_sum(s,return_type);
         }
         else{  //if a>1
              //cout << "riemann sum " << endl;
              L_value = this->value_via_Riemann_sum(s,return_type);
         }

          DIGITS3=Int( (DIGITS-DIGITS2-log(log(1.*max_n*Q+3)*abs(imag(s))/6.28+3)/2.3)/pow(2.,abs(global_derivative)))+2;
          cout << setprecision(DIGITS3);
          if (my_verbose>1) cout << "Setting output precision to: " << DIGITS3 << endl;
          tolerance3=pow(1./10,(DIGITS3+1));


          return L_value;
      }
      else if(derivative>0){
          Double h;
          h=pow(.1,(int)(DIGITS/pow(2.,derivative)));
          return((this->value(s+h,derivative-1,return_type)-this->value(s,derivative-1,return_type))/h);
      }
      else if(derivative==-1){ //simple way to use existing framework. -1 stands for logarithmic derivative
                               //because of this hack I take abs of gloabl_derivative above when determining 
                               //output precision

          L_value=this->value(s,0,return_type);
          return(this->value(s,1,return_type)/L_value); //order, i.e. value then derivative, is important since 
                                                        //derivative sets output to lower precision so should be called 2nd

      }
      else{
          cout << "Error. Specified derivative must be >= -1" << endl;
          exit(1);
      }

    }