liblfunction-dev 1.23+dfsg-6+b1 / usr / include / Lfunction / Lgamma.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.

*/



#ifndef Lgamma_H
#define Lgamma_H

#include "Lglobals.h"
#include "Lmisc.h"
#include <iomanip>          //for manipulating output such as setprecision
#include <iostream>         //for ostrstream
#include <cstring>
//using namespace std;





template <class ttype>
precise(ttype) log_GAMMA (ttype z,int n=0); //return ttype since I also want to
                                         //use it for N(T) which is real.
template <class ttype,class ttype2>
Complex GAMMA (ttype z,ttype2 delta);

//computes G(z,w)
template <class ttype>
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

template <class ttype>
ttype cfrac_GAMMA (ttype z,ttype w, ttype exp_w = 0, bool recycle=false);  //computes G(z,w) via continued fraction

template <class ttype>
ttype comp_inc_GAMMA (ttype z,ttype w,ttype exp_w = 0, bool recycle=false);  //computes g(z,w)

inline Double real(Double z){
    return z;
}
inline Double imag(Double z){
    return 0;
}

//used in Temme's asymptotics of the incomplete gamma function
Complex Q(Complex z, Complex w);
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, const char *method="temme");
//computes a sum of G(z,w)'s (as in (3.3.20)).

Complex exp_recycle();

//compute the nth derivative of log(GAMMA(z))
//n=0 should just give log_GAMMA(z)... thus making log_GAMMA
//code obsolete. But leave log_GAMMA intact anyways.
template <class ttype>
precise(ttype) log_GAMMA (ttype z,int n)
{
    int M;
    precise(ttype) log_G,r,r2,y;

    Double xx=real((Complex)z),yy=imag((Complex)z);
    if(xx<0) xx=-xx;

    Double x;
    int i,m;


    //assume the remainder stopping at the mth term is roughly bernoulli(2m)/(z+M)^(2m).
    //Now bernoulli(2m) is roughly (2m)!/(2Pi)^(2m). So remainder is more or less
    //(2m/(2ePi(z+M))^(2m). Setting 2m = Digits, we see that we need m/(ePi(z+M))
    //to be roughly 1/10 in size, i.e. |z+M| should be roughly 10*m/(ePi)=10*Digits/(2ePi).
    //squaring, gives us how large M should be.

    //n==0 leads to |z+M| >10*2m/(2*Pi*e) with 2m=Digits. However, when we differentiate
    //we end up with extra powers of (z+M) in the denominators, but extra (2m-1)...(2m+n-2)
    //in the numerators. So assume further that |z+M| > 2m+n = Digits+n

    if(n==0){
        //.343 = 100./(4*Pi*Pi*exp(2.))
        if((xx*xx+yy*yy)> .343*DIGITS*DIGITS) M=0;
        else M=Int(ceil(sqrt((DIGITS*DIGITS*.343-yy*yy))-xx+1));
    }
    else{
        if((xx*xx+yy*yy)> .343*(DIGITS+n)*(DIGITS+n)) M=0;
        else M=Int(ceil(sqrt(((DIGITS+n)*(DIGITS+n)*.343-yy*yy))-xx+1));
    }

    if(n==0)
       log_G=log_2Pi/2+(z+M-.5)*log(z+M)-(z+M);
    else if(n==1)
       log_G=log(z+M)-.5/(z+M);
    else{
       r=1;
       for(m=1;m<=n-1;m++){
          r=-r*m/(z+M);
       }
       log_G=log_G-r/(n-1)-.5*r/(z+M);
    }

    r=1;
    for(m=1;m<=n;m++){
       r=-r*m/(z+M);
    }
    r=r/(z+M);

    r2=1/((z+M)*(z+M));
    m=2;
    x=my_norm(r);
    do{
        y=bernoulli[m]*r;
        log_G=log_G+y/(m*(m-1));

        r=r*(m+n-1)*(m+n)*r2/((m-1)*m);
        m=m+2;
    //}while(m<=DIGITS&&my_norm(r)*tolerance_sqrd<x);
    }while(m<=DIGITS);

    for(m=0;m<=M-1;m++){
       if(n==0){
           log_G=log_G-log(z+m); //XXX might be faster to multiply and take one log,
                                 //but careful about overflow errors
       }
       else{
           r=1;
           for(i=1;i<=n;i++){
              r=-r*i/(z+m);
           }
           log_G=log_G+r/n;
       }
    }

    return log_G;
}

// computes GAMMA(z) delta^(-z)
// assumes the ttype is at least real.
// I compute this as exp(log(GAMMA(z))
// since there are 1/2 as many terms and the error analysis
//is easier.
//Temme's book, chapter 3 gives the error analysis for
//log(GAMMA(z)), z>0, and mentions Spira's estimate for complex values of z.
//He also suggests using the reflection formula for negative values.

template <class ttype,class ttype2>
Complex GAMMA (ttype z1, ttype2 delta)
{

    precise(ttype) z=z1;
    Complex log_G;
    Double xx=real((Complex)z),yy=imag((Complex)z);
    if(xx<0) xx=-xx;

    if(z==last_z_GAMMA) log_G=last_log_G;

    else{
        int M;
        precise(ttype) r,r2,y;

        Double x;
        int m;

        //assume the remainder stopping at the mth term is roughly bernoulli(2m)/(z+M)^(2m).
        //Now bernoulli(2m) is roughly (2m)!/(2Pi)^(2m). So remainder is more or less
        //(2m/(2ePi(z+M))^(2m). Setting 2m = Digits, we see that we need m/(ePi(z+M))
        //to be roughly 1/10 in size, i.e. |z+M| should be roughly 10*m/(ePi)=10*Digits/(2ePi).
        //squaring, gives us how large M should be.

        //.343 = 100./(4*Pi*Pi*exp(2.))
        if((xx*xx+yy*yy)> .343*DIGITS*DIGITS) M=0;
        else M=Int(ceil(sqrt((DIGITS*DIGITS*.343-yy*yy))-xx+1));

        log_G=log_2Pi/2+(z+M-.5)*log(z+M)-(z+M);

        r=(z+M)*(z+M);
        r2=(z+M);
        m=2;
        x=my_norm(r2);
        do{
            y=bernoulli[m]/r2;
            log_G=log_G+y/(m*(m-1));
            r2=r2*r;
            m=m+2;
        //}while(m<=DIGITS&&my_norm(r2)*tolerance_sqrd<x);
        }while(m<=DIGITS);

        for(m=0;m<=M-1;m++) log_G=log_G-log(z+m); //XXX could combine into a single log, but careful about overflow XXXX
        last_log_G=log_G;
        last_z_GAMMA=z;
    }


    log_G=log_G-z*log(delta);

    return exp(log_G);

}



//XXXXXXXXXX
//causes problems if z is an integer <=0

//not meant to be used for very negative values of Re(z)
//since a recursion is used to get z positive.
//if recycling is turned on (this is done in gamma_sum, for example), then we use the
//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, const char *method, ttype exp_w, bool recycle)
{

    Complex G;
    Double h,x,y;
    ttype tmp=z+1;


    if(my_verbose>2) cout << "inc_GAMMA called. G(" << z<< " , " << w << ")" << endl;

    //if(my_norm(z)<tolerance_sqrd){
    if(my_norm(z)<.01){
        return cfrac_GAMMA(z,w,exp_w,recycle);
    }

    if(my_norm(z-1)<tolerance_sqrd){
        if (recycle==false) return(exp(-w)/w);
        else return(exp_w/w);
    }


    if(real((Complex)z)<=0){
        //XXXX tmp equals z+1. compiler complains when I try inc_GAMMA(z+1,w)
        //when mpfr is used. Regards z+1 here as a binary expression rather than a
        //Complex or Double
        if (recycle==false) return ((w*inc_GAMMA(tmp,w)-exp(-w))/z);
        else return ((w*inc_GAMMA(tmp,w,method, exp_w,recycle)-exp_w)/z);
    }

    y=my_norm(w);

    if((my_norm(z)>100&&y>my_norm(z*1.01))||!strcmp(method,"continued fraction"))
    {
        return cfrac_GAMMA(z,w,exp_w,recycle);
    }

    x=imag((Complex)z);
    if(x<0)x=-x;
    //x=x-sqrt(x);
    x=.99*x; //temme works better here than x-sqrt(x), at least asymptotically

    if(y<1600||y<x*x)  //XXXX roughly if |w|<|imag(z)| , or if |w| is small
                   //XXXX NOTE:
                   //XXXX one apparent problem with using g(z,w)
                   //XXXX when Re w is large, is that g(z,w) is
                   //XXXX close to GAMMA(z)*w^(-z).
                   //XXXX So subtracting to get G(z,w) loses me precision in
                   //XXXX the 'non leading zero' digits.
                   //XXXX On the other hand... if we count the leading zeros,
                   //XXXX then we have our full precision! And we should count them!
                   //XXXX When we sum the terms in (33) of my Computation methods paper,
                   //XXXX the initial ones start
                   //XXXX off comparatively large... so the leading zeros count in
                   //XXXX the later terms (i.e. even if we had more precision for
                   //XXXX the later G(z,w) terms it would be lost when we summed
                   //XXXX against the earlier terms in (33)).
                   //XXXX note... we can use this to *slightly* speed up
                   //XXXX the computation in that we don't always need to ask g(z,w)
                   //XXXX for too many digits since only the first few will be useful.
                   //XXXX i haven't implemented this, though.
    {
        last_z=z;
        last_w=w;
        last_comp_inc_GAMMA=comp_inc_GAMMA(z,w,exp_w,recycle);
        G=GAMMA(z,w)-last_comp_inc_GAMMA;
        //cout << "series GAMMA("<<z<<","<<w<<")= " << G << endl; //XXXXX
        return G;
    }

    //XXXXX condition here and above should depend on precision
    //if(y>=10&&y<1600)
    //{
        //return cfrac_GAMMA(z,w,exp_w,recycle);
    //}

    //use temme's uniform asymptotics
    if(!strcmp(method,"temme")&&y<my_norm(z*1.2)) 
    {
        G=Q(z,w)*GAMMA(z,w);
        if(my_verbose>3) cout << "temme GAMMA("<<z<<","<<w<<")= " << G << endl; //XXXXX
        return G;
    }

    //asymptotic series... should never be called, but doesn't hurt to have it
    x=abs(z);
    h=(DIGITS+2)*2.3026+1; //XXXX+2 is to be safe
    x=x+h+sqrt(h*h+4*h*x); //XXXX (3.3.48) of my thesis generalized for DIGITS
    if(y>x*x)                //XXXX if abs(w) > (3.3.48)
    {
        return asympt_GAMMA(z,w,exp_w,recycle);
    }


    // this is called if abs(z)<10 and w is around size 40
    //return ((w*inc_GAMMA(z+1,w)-exp(-w))/z);

    if (recycle==false) return ((w*inc_GAMMA(tmp,w)-exp(-w))/z);
    else return ((w*inc_GAMMA(tmp,w,method, exp_w,recycle)-exp_w)/z);


}


template <class ttype>
ttype cfrac_GAMMA (ttype z,ttype w, ttype exp_w, bool recycle)  //computes G(z,w) via continued fraction
{

        ttype G;

        if(my_verbose>3) cout << "called cfrac_GAMMA("<<z<<","<<w<<")"<< endl;
        /*
        //old code commented out. Old code used backward recursion iterated until
        //convergence to within tolerance was satisfied. This requires extra steps
        //compared to forward recursion.

        //Complex tmp=1;
        int m;
        int M=2;
        bool escape=false;


        for(int m=M;m>=1;m--)
            //tmp=w+(m-z)/(1+m/tmp);
            tmp=w+tmp*(m-z)/(tmp+m);

        do{
            M=2*M;
            G=1;

            for(int m=M;m>=1;m--)
                G=w+G*(m-z)/(G+m);
            if(my_norm(1-tmp/G)<tolerance_sqrd)escape=true;
            else tmp=G;
            cout << "cfrac with number terms: "<<M << endl;
        }while(!escape&&M<100000000);
        */

        // newer version uses forward recursion. We also avoid the
        // expense of dividing, keeping track of the numerators
        // and denominators, dividing after convergence is achieved.

        int n;
        ttype P1=1.,P2=w,P3,Q1=0.,Q2=1.,Q3;

        n=0;
        do{
            n++;

            P3=P2+(n-z)*P1;
            Q3=Q2+(n-z)*Q1;
            P1=P2;P2=P3;
            Q1=Q2;Q2=Q3;

            P3=w*P2+n*P1;
            Q3=w*Q2+n*Q1;
            P1=P2;P2=P3;
            Q1=Q2;Q2=Q3;

            //to prevent overflow
            //XXX check this
            if(n%8==0&&(real(P2)>1.e40||real(P2)<-1.e40||imag(P2)>1.e40||imag(P2)<-1.e40)){
           //cout << "before " << P1 ;
                P1=P1*1.e-40;
                P2=P2*1.e-40;
                Q1=Q1*1.e-40;
                Q2=Q2*1.e-40;

          //cout << "  after " << P1 <<endl;
            }

            //cout << "cfrac1 : " << n<< " " <<my_norm(Q2*P1-P2*Q1) << " " << my_norm(Q2*P1*tolerance)<<endl;
            //cout << "cfrac called n equals: " << n << " P2/Q2 " << P2/Q2 << endl;
        } while(n<2||(my_norm(Q2*P1-P2*Q1)>my_norm(Q2*P1*tolerance)&&n<1000000));
        //cout << "cfrac called n equals: " << n << " " << abs(z/w) << endl;

        G=P2/Q2;


        if(n>999999){
             cout << "Continued fraction for G(z,w) failed to converge. z = "
             << z << "  w = " << w << endl;
             exit(1);
        }


        if(recycle==false) G=exp(-w)/G; //XXXXX mpfr should I precise(ttype) w ?
        else G=exp_w/G;

        //cout<< setprecision(30);
        //cout << "cfrac GAMMA("<<z<<","<<w<<")= " << G << endl; //XXXXX
        return G;
}

template <class ttype>
ttype asympt_GAMMA (ttype z,ttype w, ttype exp_w, bool recycle)  //computes G(z,w) via asymptotic series
{

        if(my_verbose>3) cout << "called asympt_GAMMA("<<z<<","<<w<<")"<< endl;
        ttype G=0;
        ttype r=1.;
        int j=0;
        do
        {
            G=G+r;
            r=r*(z-1-j)/w;
            //cout << j << " " << G << " " << r <<" " << tolerance <<endl; //XXXXX
            j++;
        }while(my_norm(r)>tolerance_sqrd);
        if(recycle==false) G=G*exp(-w)/w;
        else G=G*exp_w/w;
        //cout << "asymptotics GAMMA("<<z<<","<<w<<")= " << G << endl; //XXXXX
        return G;
}


template <class ttype>
ttype comp_inc_GAMMA (ttype z,ttype w,ttype exp_w, bool recycle)  //computes g(z,w)
{

    ttype g;
    Double t;
    int m;

    //Complex u=w*w;
    //u=u*u;

    if(my_verbose>3) cout << "called comp_inc_GAMMA("<<z<<","<<w<<")"<< endl;

    t=my_norm(w/z);


    if(t>.9801 || my_norm(w)<.36){

        ttype r=1.;
        m=1;
        g=0;
        do{   //XXXXXXXXXXXXX optimize combine the m++'s and g=g+r and z+m's
            g=g+r;
            r=r*w/(z+m);
            m++;
            g=g+r;
            r=r*w/(z+m);
            m++;
            g=g+r;
            r=r*w/(z+m);
            m++;

            //v=z+m;
            //g=g+r+

        //cout<< "using series for comp inc " << t << " " << abs(w) << " " << m <<" " << g*exp(-w)/z<< endl;

        }while(my_norm(r)>tolerance_sqrd||real((Complex)z)<=-m);
        if(recycle==false) g=g*exp(-w)/z;
        else g=g*exp_w/z;
    }

    else{

        /*-----------------------------------------------------------
        // old iterated backward recursion replaced by forward recursion

        int M=2;
        Complex tmp=1;
        bool escape=false;


        for(int m=M;m>=1;m--){
            if(m%2==0){
               //tmp= z+m-1+(m/2)*w/tmp;
               tmp= z+m-1+m*.5*w/tmp;
            }
            else{
                tmp= z+m-1-(z+(m-1)*.5)*w/tmp;
            }
        }

        do{
            M=2*M;
            g=1;
            for(int m=M;m>=1;m--){
                if(m%2==0){
                    g= z+m-1+m*.5*w/g;
                }
                else{
                    g= z+m-1-(z+(m-1)*.5)*w/g;

                }
            }
            if(abs(1-tmp/g)<tolerance)escape=true;
            else tmp=g;

        }while(!escape&&M<10000000);

        //XXXXXX maybe put an error message + exit here if M is big
        if(M==10000000){
            cout << "continued fraction for complimentary incomplete gamma failed to converge.";
            cout<< endl << "z = " << z << "  w = " << w << endl;
            exit(1);
        }


        */

        //cout <<setprecision(20) << endl;
        //cout << "g1: " << g << endl;

        int n;
        ttype P1=1.,P2=z,P3,Q1=0.,Q2=1.,Q3;
        ttype u=.5*w;
        //ttype t1,t2;



        n=0;
        do{
            n++;
            P3=(z+n)*P2-(z+(n-1)*.5)*w*P1;
            Q3=(z+n)*Q2-(z+(n-1)*.5)*w*Q1;

            //t1=z+n;
            //t2=(z+(n-1)*.5)*w;
            //P3=t1*P2-t2*P1;
            //Q3=t1*Q2-t2*Q1;

            P1=P2;P2=P3;
            Q1=Q2;Q2=Q3;

            n++;
            P3=(z+n)*P2+n*u*P1;
            Q3=(z+n)*Q2+n*u*Q1;
            //t1=t1+1; t2=n*u;
            //P3=t1*P2+t2*P1;
            //Q3=t1*Q2+t2*Q1;

            P1=P2;P2=P3;
            Q1=Q2;Q2=Q3;

//cout << P2/Q2 << " " << norm(Q2*P1-P2*Q1) / norm(Q2*P1*tolerance) <<endl;

            //to prevent overlow
            if(n%8==0&&(real(P2)>1.e50||real(P2)<-1.e50||imag(P2)>1.e50||imag(P2)<-1.e50)){
                P1=P1*1.e-50;
                P2=P2*1.e-50;
                Q1=Q1*1.e-50;
                Q2=Q2*1.e-50;

            }

        //cout<< "using cfrac for comp inc " << t << " " << abs(w) << " " << n <<" " << Q2/P2*exp(-w)<< " GAMMA " << GAMMA(z,w)<< endl;
        } while(n<3||(my_norm(Q2*P1-P2*Q1)>my_norm(Q2*P1*tolerance)&&n<1000000));

        g=P2/Q2;

        //cout<< "using cfrac for comp inc " << t << " " << n << endl;

        if(n>999999){
             cout << "Mofu. Continued fraction for g(z,w) failed to converge. z = "
             << z << "  w = " << w << endl;
             exit(1);
        }



        if(recycle==false) g=exp(-w)/g;
        else g=exp_w/g;

    }

    return g;

}

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, const char *method)
{
    Complex SUM=0;

    Complex z,w;
    Complex G;
    Complex r;
    Complex u;

    Complex e1,e2,e3,exp_w; //used to compute exp(-n w_0) or exp(-n^2 w_0) by repeated multiplication
                                      //saves on taking exponents

    int n=1;
    int n2=1;

    Double x,y,y2,y3=0,MAX=0;
    bool escape=false;
    bool is_z_real=false;
    bool is_w_real=false;


    z=g*s+l;

    if(my_norm(imag(z))<tolerance_sqrd) is_z_real=true;
    if(my_norm(imag(delta))<tolerance_sqrd) is_w_real=true;

    w=delta/Q;
    if(g<.6) w=w*w;     //i.e. if g=1/2

    e1=exp(-w); exp_w=1.;
    e2=e1*e1; e3=1.;

    //for recycling of data
    //two cases: exp(-w), with w= n w_0 or n^2 w_0
    //Denote as e_n. Then

    //exp_w -> exp_w e1 when gamma = 1
    //      = exp_w exp(-(2n-1) w_0) when gamma = 1/2


    //y=abs(z)+abs(real(z))+1;


    if(what_type==-1)   //i.e. if the Riemann zeta function
    do{
        w=Pi*n*n*delta*delta;
        exp_w=exp_w*e1*e3;  //e1 is exp(-(delta/Q)^2) in this case, and e3 is exp(-(2n-2)*(delta/Q)^2), Q =1/sqrt(Pi)
        e3=e3*e2;
        G=inc_GAMMA(z,w,method,exp_w,true);
        SUM=SUM+G;
        n++;

        x=my_norm(SUM);
        if(x>MAX) MAX=x;

        if(real(w-z)>10) //10 is kind of arbitrary . G(z,w) will decay like e(-Re(w)), so 
                         //we'll escape once Re(w) is around log(10)*DIGITS. So, 10 is okay
                         //as a place to start checking
        {
            if (my_norm(G)<MAX*tolerance_sqrd) escape=true;
        }
    }while(!escape);

    //XXXX checking abs(G)>tolerance is not so smart... since two of the ways 
    //XXXX I compute G, series and nielsen, only gives us 10E-15, say, regardless
    //XXXX how small it really is (since we subtract g(z,w) from GAMMA(z)w^(-z)
    //XXXX and these are nearly equal when w is big enough.
    //XXXX Best to compare estimate for largest of terms yet to be added
    //XXXX and see if this is <MAX*tolerance, where MAX is the maximum
    //XXXX attained by the partial sums.

    else
    do
    {


        w=n*delta/Q;
        //i.e. if g=1/2
        if(g<.6) {
            w=w*w;
            exp_w=exp_w*e1*e3;  //e1 is exp(-(delta/Q)^2) in this case, and e3 is exp(-(2n-2)*(delta/Q)^2)
            e3=e3*e2;
        }
        else exp_w=exp_w*e1; //if gamma=1, we need exp_w=exp(-n delta/Q)


        if(l==0){
            u=1;
        }
        else{
            u=exp(LOG(n)*l/g);  //XXX this can be stored if it is called repeatedly
        }


        if(coeff[n2]!=0)
        {
            //if both are real, we should send as Doubles. That way we avoid complex arithmetic
            //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 with doubles = " << G << endl;
                //cout<<"both z and w are real\n";
            }
            else{
                G=inc_GAMMA(z,w,method,exp_w,true);
                if (my_verbose>2) cout << "GAMMA SUM = " << G << endl;
                //cout<<"none are real\n";
            }
            SUM=SUM+G*u*coeff[n2];
        }


        n++; n2++;

        x=my_norm(SUM);
        if(x>MAX) MAX=x;

        if(real(w-z)>10) //10 is kind of arbitrary . G(z,w) will decay like e(-Re(w)), so 
                         //we'll escape once Re(w) is around log(10)*DIGITS. So, 10 is okay
                         //as a place to start checking
        {
            //y3=4*exp(-real(w))/(y2+1);
            //if (my_norm(u)*y3*y3*n*n<MAX*tolerance_sqrd) escape=true;
            if (my_norm(u*G)*n*n<MAX*tolerance_sqrd) escape=true;
        }

        if(n2>Period&&Period>1) n2=(n2-Period);

    }while(n2<=N&&!escape);


    //XXXXnote, we copy the tolerance feature for zeta, but make sure
    //XXXXto also include b(n)<sigma_0(n)< n and n^(Re l/g) factors.
    //XXXX(that is why we have the n*u in the while).
    if(n2>N&&what_type!=-1)
    {
 
        if(print_warning){
            print_warning=false;
            cout << "WARNING from gamma 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;
        }
        //exit(1);

    }


    max_n = n;
    if(my_verbose>0) cout << "s = " << s << "gamma_sum called, number terms used: " << n << endl;

    return SUM;

}

#endif