Logo Search packages:      
Sourcecode: lcalc version File versions  Download package

Lfind_zeros.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.

*/


    //use Brent's method to locate a zero given a sign change.
    //see the wikipedia article: http://en.wikipedia.org/wiki/Brent's_method
    template <class ttype>
    Double L_function <ttype>::
    zeros_zoom_brent(Double L1, Double L2, Double u, Double v)
    {

        Double a=u,b=v,L_a=L1,L_b=L2,c,L_c,L_s,d;
        bool my_flag=true;
        bool s_bound;
        Double tmp,tmp2;
        Double s,x,y,z;
        //Double tol=tolerance2*sqrt(v*v+1);

        if(my_norm(L_a)<my_norm(L_b)){
            tmp=a;a=b;b=tmp;
            tmp=L_a;L_a=L_b;L_b=tmp;
        }
        c=a; L_c=L_a;

        do{
            //cout << a << " " << b << " " << L_a <<" "  << L_b << endl;
            x=L_a-L_b; y=L_a-L_c; z=L_b-L_c;
            if(y!=0&&z!=0){
                s=a*L_b*L_c/(x*y)-b*L_a*L_c/(x*z)+c*L_a*L_b/(y*z); //inverse quad interpolation
                //cout << "quad     " << s << endl;
            }
            else{
                s=b-L_b*(b-a)/(L_b-L_a); //secant rule
                //cout << "secant   " << s << endl;
            }

            if(a<b){
                tmp=(3*a+b)/4;tmp2=b;
            }
            else{
                tmp2=(3*a+b)/4;tmp=b;
            }
            if ((s<tmp||s>tmp2) || (my_flag &&my_norm(s-b)>= my_norm(b-c) / 2) || (!my_flag && my_norm(s-b)>=my_norm(c-d)/2)){
                s=(a+b)/2; my_flag=true;
                //cout << "midpoint " << s << endl;
            }
            else my_flag=false;

            L_s= real(this->value(.5+I*s,0,"rotated pure"));
            d=c;
            c=b; L_c=L_b;
            if (L_a*L_s<0){
                b=s;L_b=L_s;
            }
            else{
                a=s;
                L_a=L_s;
            }
            if(my_norm(L_a)< my_norm(L_b)){
                tmp=a;a=b;b=tmp;
                tmp=L_a;L_a=L_b;L_b=tmp;
            }

        }while(abs(L_b)>tolerance3 && abs((b-a)/(abs(b)+1))>tolerance2);
        // sqrt(v*v+1) is only good for zeta. Other L-functions have more
        // precision loss.

        //if(abs(L_b)>tol*100)
        //    cout << "Mofu. Big zero  " << b << " " <<L_b<<endl;


        return(b);


    }

    template <class ttype>
    void L_function <ttype>::
    find_zeros(Double t1, Double t2, Double step_size, char* filename, char* message_stamp)
    {
        Double t,x,y;
        Double u,v,tmp2;
        Long count2=0;
        Double previous_zero=t1;

        fstream file;



        if(!strcmp(filename,"cout"))
        {
            cout << setprecision(DIGITS3);
        }
        else
        {
            file.open(filename, ios::out|ios::app);
            file << setprecision(DIGITS3);
        }

        t=t1; u=t;
        x=real(this->value(.5+I*t,0,"rotated pure"));

        if(my_verbose>1)
        cout << "look for sign change " << t << " " << x << endl;


        do{
            t=t+step_size;
            y=real(this->value(.5+I*t,0,"rotated pure"));
            if(my_verbose>1)
            cout << "look for sign change " << t << " " << y << endl;
            v=t;
            if(sn(x)!=sn(y))  //if sign change is found...divide and conquer
                              //until the difference is small.
            {

                tmp2=zeros_zoom_brent(x,y,u,v);
                count2++;

                int tmp_DIGITS = Int(DIGITS3+log(abs(tmp2)+2)/2.3)+1;
                if(!strcmp(filename,"cout"))
                {
                    //cout << tmp2 << endl;
                    cout << message_stamp << " ";
                    if(tmp_DIGITS<DIGITS)
                        cout << setprecision(tmp_DIGITS);
                    else
                        cout << setprecision(DIGITS);
                    cout << tmp2 << " ";
                    cout << setprecision(DIGITS3);
                    if(t1>0)
                    cout << this->N((tmp2+previous_zero)/2)/2 -N(t1)/2-(count2-1.) << endl;
                    else
                    cout << this->N((tmp2+previous_zero)/2)/2 -(count2-1.) << endl;
                    previous_zero=tmp2;
                }
                else
                {
                    //file << tmp2 << endl;
                    file << message_stamp << " ";
                    if(tmp_DIGITS<DIGITS)
                        file << setprecision(tmp_DIGITS);
                    else
                        file << setprecision(DIGITS);
                    file << tmp2 << " ";
                    file << setprecision(DIGITS3);
                    if(t1>0)
                    file << this->N((tmp2+previous_zero)/2)/2 -N(t1)/2-(count2-1.) << endl;
                    else
                    file << this->N((tmp2+previous_zero)/2)/2-(count2-1.) << endl;
                    previous_zero=tmp2;
                }
            }
            u=t;
            x=y;
        }while((t1<t&&t<t2)||(t2<t&&t<t1));
        if(strcmp(filename,"cout")) file.close();

    }


    //find zeros using gram points and Rosser's rule. Not used by package.
    template <class ttype>
    void L_function <ttype>::
    find_zeros_via_gram(Double t1, Long count,Double max_refine,char* filename, char* message_stamp){

        //find_zeros_elaborate(t1,count,max_refine);
 cout <<"GRAM POINT ZERO FINDER"<< endl;
        Long m,n;
        Double r,x;
        int i,j,k;
        int gram_sign; // sign of (-1)^n L(gram_pt_n)

        // compute 100 gram points, check gram's law, and partition into blocks (we might
        // need to stop before the 100th gram point).

        Double gram_pt[100][2];// 0 stores the gram point, 1 stores the corresponding L-value
        x=initialize_gram(t1);
        gram_pt[0][0]=x;
        gram_pt[0][1]=real(this->value(.5+I*x,0,"rotated pure"));
        for(m=1;m<=99;m++){
            x=next_gram(x);
            gram_pt[m][0]=x;
            gram_pt[m][1]=real(this->value(.5+I*x,0,"rotated pure"));
//cout << m << " " << gram_pt[m][1]*pow(-1.,1.*m) << endl;
        }

        r=1;m=n=0;

        for(i=0;i<=99;i++){
            if(r*gram_pt[i][1]>0)m++; else n++;
            r=-r;
        }

        if(m>65) gram_sign=1;
        else if(n>65) gram_sign=-1;
        else{
            cout << "Possible violation of Gram's law" << endl;
for(i=0;i<=99;i++)
cout << i<< " " <<gram_pt[i][0] << " " << gram_pt[i][1]*pow(-1.,i) << " " <<
gram_pt[i][1]*real(this->value(.5-I*gram_pt[i][0],0,"rotated pure"))<< endl;

            exit(1);
        }

        int number_blocks=0; //number of gram blocks
        int gram_block[100][3]; //0 stores the start index, 1 the stop index, 
                                   //2 the number of sign changes in the block.
         
        i=0; r=gram_sign;
        do{
            if(r*gram_pt[i][1]>0){
                j=i;
                do{
                    j++; r=-r;
                }while(r*gram_pt[j][1]<0&&j<99);
                if(r*gram_pt[j][1]>0){
if(j-i>3){ cout<<"is a long block\n";
for(int ii=i;ii<=j;ii++)
cout << ii<< " " <<gram_pt[ii][0] << " " << gram_pt[ii][1]*pow(-1.,ii) << endl;

}

                    number_blocks++;
                    gram_block[number_blocks][0]=i;
                    gram_block[number_blocks][1]=j;
                    gram_block[number_blocks][2]=0;
                }
                i=j;
            }
            else{i++;r=-r;}
      
        }while(i<99);
        
        // check which blocks are good. If any are bad, search in neighboring blocks.

        //if, in the end, any of the front 6 blocks, say, are still bad or if gram_pt_1 > t1
        //then look for the zeros in [t1,gram_pt_1): count how many there ought to be
        //using arg principle, then search for them. Actually, we should use the arg principle
        //to count how many zeros there ought to be from t1 to the first point from where the
        //gram point search kicks in successfully.

        //if all zeros are found dump the first 50, say and repeat the gram block search
        //using the next 50 gram points.

        
        
/*

        int i; Double tmp,tmp2;
        Complex L_1;
        Double y;
         
        Double global_average=0;
        x=.1;i=0;
        do{
          L_1=this->value(x+I*(t1-0.65444));
          //L_1=this->value(x+I*(t1-.1));
          if(L_1!=0) y=imag(log(L_1));


          if(i>0){
             tmp2=y-tmp;
             if(tmp2>5)tmp2=tmp2-2*Pi;
             if(tmp2<-5)tmp2=tmp2+2*Pi;
             global_average=global_average+tmp2;
//cout << "global average: " <<x<<" "<< y << " " <<tmp << " " <<global_average/(2*Pi) << endl;
          }
          tmp=y;
          x=x+.01;
          i++;
        }while(x<=.9);
        Double t2=t1-0.65444;
        //Double t2=t1-0.1;
        do{
          t2=t2+.01;
          L_1=this->value(x+I*t2);
          if(L_1!=0) y=imag(log(L_1));

          tmp2=y-tmp;
          if(tmp2>5)tmp2=tmp2-2*Pi;
          if(tmp2<-5)tmp2=tmp2+2*Pi;
          global_average=global_average+tmp2;
//cout << "global average: " <<x<<" "<< y << " " <<tmp << " " <<global_average/(2*Pi) << endl;
          tmp=y;
        }while(t2<=1);
        do{
          x=x-.01;
          L_1=this->value(x+I*t2);
          if(L_1!=0) y=imag(log(L_1));

          tmp2=y-tmp;
          if(tmp2>5)tmp2=tmp2-2*Pi;
          if(tmp2<-5)tmp2=tmp2+2*Pi;
          global_average=global_average+tmp2;
//cout << "global average: " <<x<<" "<< y << " " <<tmp << " " <<global_average/(2*Pi) << endl;
          tmp=y;
        }while(x>.1);
        do{
          t2=t2-.01;
          L_1=this->value(x+I*t2);
          if(L_1!=0) y=imag(log(L_1));

          tmp2=y-tmp;
          if(tmp2>5)tmp2=tmp2-2*Pi;
          if(tmp2<-5)tmp2=tmp2+2*Pi;
          global_average=global_average+tmp2;
//cout << "global average: " <<x<<" "<< y << " " <<tmp << " " <<global_average/(2*Pi) << endl;
          tmp=y;
        }while(t2>t1-0.65444);
        //}while(t2>t1-.1);
        global_average=global_average/(2*Pi);

cout <<"hihih " << global_average<<endl;
*/

    }

    template <class ttype>
    int L_function <ttype>::
    compute_rank(bool print_rank){

          Complex w,z;
          Double x,y,h1,h2;
          int r;

          z=this->value(.5);
          x=abs(z);
            if(x>.00001){
                //cout << "here 1\n";
                if(print_rank) cout << "analytic rank equals "<< 0 << endl;
                return 0;
          }

            h2=.00001;
          y=abs(this->value(.5+h2));
            if(y>1.e-10){
                //cout << "here 2\n";
              x=abs(this->value(.5+10*h2));
                r = Int(rint(abs(log(x)-log(y))/log(10.)));
              if(print_rank) cout << "analytic rank equals "<< r << endl;
              return r;
          }

            h2=.001;
            y=abs(this->value(.5+h2));
            if(y>1.e-10){
                //cout <<"here 3\n";
                x=abs(this->value(.5+10*h2));
                r = Int(rint(abs(log(x)-log(y))/log(10.)));
                    
                if(print_rank) cout << "analytic rank equals "<< r << endl;
                return r;
            }

          h1=.001;
          do{
                h1=h1*4; 
                x=abs(this->value(.5+h1));
                //cout <<"here 4 " <<x<<endl;
          }while(x<1e-10&&h1<2);
          h2=h1;
          h1=h1*10;
          x=abs(this->value(.5+h1));
          y=abs(this->value(.5+h2));
          //cout << x << " " << y << " "<<(log(x)-log(y))/log(10.)<<endl;
          r = Int(rint(abs(log(x)-log(y))/log(10.)));
          if(print_rank) cout << "analytic rank equals "<< r << endl;
          return (r);
          
    }

    
    template <class ttype>
    void L_function <ttype>::
    verify_rank(int rank){
      int analytic_rank=compute_rank();
      if(rank!=analytic_rank)
      cout<< "given rank "<<rank<< " is different than computed analytic rank "<<
            analytic_rank <<endl;
      //else 
      //cout<< "given rank "<<rank<< " equals analytic rank"<<endl;
    }

    //elaborate zero finder keeps an array of zeros and goes back to find missing zeros
    //if they are detected via N(T) comparison. 
    //Finds zeros of L(1/2+It) or of L(1/2+It)*L(1/2-It) (choice is specified by bool do_negative)
    //Latter must be used if the L-function has complex coefficients.
    template <class ttype>
    void L_function <ttype>::
    find_zeros_via_N(Long count,bool do_negative,Double max_refine, int rank,char* filename, char* message_stamp)
    {


        Double t,x,y;
        Double u,v,tmp=0;
        Long count2=0; //counts how many zeros have been printed
        Double count_all=0; //counts how many zeros have been found
        Double previous_zero=0; 
        Double step_size,refined_step_size;

        fstream file;

        int i,j;
        Double zeros_S[100][3]; // stores list of consecutive zeros: S[n][0]
                                // corresponding S(T) value: S[n][1]
                                // and the sign of the zero (+1 or -1) S[n][2]

        int ii; //counts number of local zeros
        int i1,i2; //used for looping through blocks for missing zeros
        Double local_average;
        int found_missing=0; //number missing zeros found
        int to_find=0; //number of missing zeros to find
        Double x2,y2,u2,v2,tmp2,tmp3,tmp4; //used to search for missing zeros
        Double x3,x3_c,u3;
        Double x_c=1,y_c=1,x2_c=1,y2_c=1; //used to search for zeros below the real axis
        int number_insert; //1 or 2 according to whether we find zeros below and/or above real axis

        bool also_do_end_pt; //when looking for missing zeros between found zeros, 
                                //we need to do one last sign check step in the case 
                                //that we are simultaneously looking for zeros of the 
                                //conjugate as well.

        int analytic_rank; //order of vanishing at the critical point

        for(i=0;i<=99;i++)
        for(j=0;j<=2;j++) 
            zeros_S[i][j]=0;


        if(!strcmp(filename,"cout"))
        {
            cout << setprecision(DIGITS3);
        }
        else
        {
            file.open(filename, ios::out|ios::app);
            file << setprecision(DIGITS3);
        }

        //count multiplicity of zero at 0, then start slightly higher.
        //doesn't matter if I miss a zero inbetween, since this
        //will then be detected by N(T) comparison and searched for.


        ii=0; //will be used to count zeros in the zeros_S array

    if(rank>=0) analytic_rank=rank;
    else analytic_rank=this->compute_rank();

    if(analytic_rank>0){
        count_all=count_all+analytic_rank; 
        if(!do_negative) count_all=count_all*.5;
        ii=analytic_rank;

         for(i=1;i<=ii;i++){
              zeros_S[i][0]=0; 
              zeros_S[i][1]=0;
              zeros_S[i][2]=1; 
         }
     }

    if(analytic_rank==0) t=0;
    else t=exp(log(.00001)/analytic_rank); //ad hoc method that will eventually fail

    u=t;

        x=real(this->value(.5+I*t,0,"rotated pure"));
        if(do_negative) x_c=real(this->value(.5-I*t,0,"rotated pure"));

        if(my_verbose>1)
        cout << "look for sign change " << t << " " << x << endl;


        do{

            step_size=.6/(this->N(t+22)-this->N(t+21)); //i.e. increment by .6 the average
            //the choice of .6 was determined experimentally

            if(do_negative) step_size=step_size/2;



            t=t+step_size;
            y=real(this->value(.5+I*t,0,"rotated pure"));
            if(do_negative) y_c=real(this->value(.5-I*t,0,"rotated pure"));

            if(my_verbose>1)
            cout << "look for sign change " << t << " " << y << " " << y_c<<endl;

            v=t;
            if(sn(x)!=sn(y))  //if sign change is found...divide and conquer
                              //until the difference is small.
            {

                tmp=zeros_zoom_brent(x,y,u,v);
                count_all++; ii++;

                zeros_S[ii][0]=tmp; 
                zeros_S[ii][2]=1; 

                if(do_negative)zeros_S[ii][1]=this->N((tmp+previous_zero)/2) -(count_all-1);
                else zeros_S[ii][1]=this->N((tmp+previous_zero)/2)/2 -(count_all-1);


            }
            if(sn(x_c)!=sn(y_c))  //if sign change is found...divide and conquer
                                  //until the difference is small.
            {

                tmp=-zeros_zoom_brent(y_c,x_c,-v,-u);
                count_all++; ii++;

                zeros_S[ii][0]=tmp;
                zeros_S[ii][2]=-1;
                if(zeros_S[ii][0]<zeros_S[ii-1][0]){

                   tmp2=zeros_S[ii-1][0]; 
                   zeros_S[ii-1][0]= zeros_S[ii][0];
                   zeros_S[ii][0]= tmp2;

                   tmp2=zeros_S[ii-1][2]; 
                   zeros_S[ii-1][2]= zeros_S[ii][2];
                   zeros_S[ii][2]= tmp2;

                   tmp2=zeros_S[ii-1][1]; 
                   zeros_S[ii-1][1]= zeros_S[ii][1]+1;
                   zeros_S[ii][1]= tmp2-1;

                }
                zeros_S[ii][1]=this->N((tmp+previous_zero)/2) -(count_all-1);


            }
                    if(ii>15){
                        local_average=0;
                        for(i=ii;i>=ii-15;i--) local_average=local_average+zeros_S[i][1]/16;

                        if(local_average>.7) to_find=1;
                        if(local_average>1.5) to_find=Int(local_average+.5);
                        if(to_find>0){

                       if(my_verbose>1)
                       cout << "missing zeros detected "<<endl;

                            refined_step_size=step_size/2;
                            i1=ii-7; i2=i1; 



                            found_missing=0;
                            do{

                                if(i2>ii-1+found_missing)i2=ii-1+found_missing; 
                                if(i1<1&&count2>0)i1=1; if(i1<1&&count2==0)i1=0;

                                for(i=i1;i<=i2;i++){


                                    u2=zeros_S[i][0];
                                    if(do_negative){
                                        if(zeros_S[i][2]>0){
                                            x2=real(this->value(.5+I*(u2+refined_step_size),0,"rotated pure"));
                                            x2_c=real(this->value(.5-I*u2,0,"rotated pure"));
                                        }
                                        else{
                                            x2=real(this->value(.5+I*u2,0,"rotated pure"));
                                            x2_c=real(this->value(.5-I*(u2+refined_step_size),0,"rotated pure"));
                                        }


                                    }
                                    else{
                                        u2=zeros_S[i][0]+.5*refined_step_size;
                                        x2=real(this->value(.5+I*u2,0,"rotated pure"));
                                    }


                                    also_do_end_pt=false;

                                    do{
                                        v2=u2+refined_step_size;
                                        if(v2<zeros_S[i+1][0])
                                        {
                                            also_do_end_pt=true;
                                            y2=real(this->value(.5+I*v2,0,"rotated pure"));
                                            if(do_negative) y2_c=real(this->value(.5-I*v2,0,"rotated pure"));
                                            x3=y2;x3_c=y2_c;u3=v2;

                                            if(sn(x2)!=sn(y2)||sn(x2_c)!=sn(y2_c))  //if sign change is found...divide and conquer
                                                                                    //until the difference is small.
                                            {

                                                number_insert=0;
                                                if(sn(x2)!=sn(y2)){
                                                    tmp2=zeros_zoom_brent(x2,y2,u2,v2);
                                                    found_missing++;i2=i2+1;
                                                    count_all++;
                                                    number_insert++;

                                                }
                                                if(sn(x2_c)!=sn(y2_c)){
                                                    tmp3=-zeros_zoom_brent(y2_c,x2_c,-v2,-u2);
                                                    if(sn(x2)!=sn(y2)&&tmp2>tmp3){tmp4=tmp2;tmp2=tmp3;tmp3=tmp4;}
                                                    found_missing++;i2=i2+1;
                                                    count_all++;
                                                    number_insert++;
                                                }
                                                    for(j=ii+found_missing;j>=i+1+number_insert;j--){
                                                        zeros_S[j][0]=zeros_S[j-number_insert][0];
                                                        zeros_S[j][2]=zeros_S[j-number_insert][2];
                                                        zeros_S[j][1]=zeros_S[j-number_insert][1]-number_insert;
                                                    }
                                                    if(number_insert==1&&sn(x2)!=sn(y2)){
                                                        zeros_S[i+1][0]=tmp2;
                                                        zeros_S[i+1][2]=1;
                                                        zeros_S[i+1][1]=zeros_S[i+2][1]-1; //not strictly correct, but close enough 
//XXXXXXXXXXXXX do this more precisely.don't be lazy!
                                                    }
                                                    else if(number_insert==1&&sn(x2_c)!=sn(y2_c)){
                                                        zeros_S[i+1][0]=tmp3;
                                                        zeros_S[i+1][2]=-1;
                                                        zeros_S[i+1][1]=zeros_S[i+2][1]-1; //not strictly correct, but close enough 
//XXXXXXXXXXXXX do this more precisely.don't be lazy!
                                                    }
                                                    else{
                                                        zeros_S[i+2][0]=tmp3;
                                                        zeros_S[i+2][2]=-1;
                                                        zeros_S[i+2][1]=zeros_S[i+3][1]-1; //not strictly correct, but close enough 
                                                        zeros_S[i+1][0]=tmp2;
                                                        zeros_S[i+1][2]=1;
                                                        zeros_S[i+1][1]=zeros_S[i+2][1]-1; //not strictly correct, but close enough 

                                                    }

                                                    i=i-2*number_insert;if(i<1)i=1;
                                            }
                                            //x3=x2;x3_c=x2_c;u3=u2;

                                            x2=y2;x2_c=y2_c;u2=v2;
                                         }

                                    }while(v2<zeros_S[i+1][0]);

                                    //in the do_negative case, there is one final sign change check needed
                                    if(also_do_end_pt&&found_missing==0&&do_negative&&zeros_S[i+1][2]>0){
                                        v2=zeros_S[i+1][0];

                                        y2_c=real(this->value(.5-I*v2,0,"rotated pure"));

                                        if(sn(x3_c)!=sn(y2_c))
                                        {

                                            tmp3=-zeros_zoom_brent(y2_c,x3_c,-v2,-u3);
                                            found_missing++;i2=i2+1;
                                            count_all++;
                                            for(j=ii+found_missing;j>=i+2;j--){
                                                zeros_S[j][0]=zeros_S[j-1][0];
                                                zeros_S[j][2]=zeros_S[j-1][2];
                                                zeros_S[j][1]=zeros_S[j-1][1]-1;
                                            }
                                            zeros_S[i+1][0]=tmp3;
                                            zeros_S[i+1][2]=-1;
                                            zeros_S[i+1][1]=zeros_S[i+2][1]-1;

                                            i=i-2;if(i<1)i=1;

                                        }


                                    } //if(do_negative&&zeros_S[i+1][2]>0)
                                    else if(also_do_end_pt&&found_missing==0&&do_negative&&zeros_S[i+1][2]<0){
                                        v2=zeros_S[i+1][0];
                                        y2=real(this->value(.5+I*v2,0,"rotated pure"));


                                        if(sn(x3)!=sn(y2))
                                        {

                                            tmp2=zeros_zoom_brent(x3,y2,u3,v2);
                                            found_missing++;i2=i2+1;
                                            count_all++;
                                            for(j=ii+found_missing;j>=i+2;j--){
                                                zeros_S[j][0]=zeros_S[j-1][0];
                                                zeros_S[j][1]=zeros_S[j-1][1]-1;
                                                zeros_S[j][2]=zeros_S[j-1][2];
                                            }
                                            zeros_S[i+1][0]=tmp2;
                                            zeros_S[i+1][2]=1;
                                            zeros_S[i+1][1]=zeros_S[i+2][1]-1;

                                            i=i-2;if(i<1)i=1;

                                        }

                                    }
                                }//for i
                                refined_step_size=refined_step_size/2;
                                i2=i2+2;i1=i1-2;
                                if(refined_step_size/step_size< .1){i2=i2+1;i1=i1-8;}

                            }while(found_missing<to_find && refined_step_size > step_size/max_refine);
                            to_find=0;

                            if(refined_step_size <= step_size/max_refine){
                                int tmp_DIGITS = Int(DIGITS3+log(zeros_S[i][0]+2)/2.3)+1;
                                if(!strcmp(filename,"cout"))
                                {
                                    if(tmp_DIGITS<DIGITS)
                                        cout << setprecision(tmp_DIGITS);
                                    else
                                        cout << setprecision(DIGITS);
                                    cout<<"missing zeros detected." << endl;
                                    cout<<"failed to find using refined step sizes of " << refined_step_size << endl;
                                    for(i=1;i<=ii+found_missing;i++)
                                        cout << zeros_S[i][0] << " " << zeros_S[i][1]<< endl;
                                }
                                else{
                                    if(tmp_DIGITS<DIGITS)
                                        file << setprecision(tmp_DIGITS);
                                    else
                                        file << setprecision(DIGITS);
                                    file<<"missing zeros detected." << endl;
                                    file<<"failed to find using refined step sizes of " << refined_step_size << endl;
                                    for(i=1;i<=ii+found_missing;i++)
                                        file << zeros_S[i][0] << " " << zeros_S[i][1]<< endl;
                                }
                                file.close();
                                exit(1);
                            }


                        }

                    }

                    ii=ii+found_missing; found_missing=0;


                    //dump up to 10 zeros
                    ii=ii+found_missing; found_missing=0;
                    if(ii>=50){
                        i=1;
                        do{
                            int tmp_DIGITS = Int(DIGITS3+log(zeros_S[i][0]+2)/2.3)+1;
                            if(!strcmp(filename,"cout"))
                            {

                                cout << message_stamp << " ";
                                if(tmp_DIGITS<DIGITS)
                                    cout << setprecision(tmp_DIGITS);
                                else
                                    cout << setprecision(DIGITS);
                                cout << setprecision(Int(DIGITS3+log(zeros_S[i][0]+2)/2.3));
                                cout << zeros_S[i][0]*zeros_S[i][2]<<endl;

                            }
                            else
                            {
                                file << message_stamp << " ";
                                if(tmp_DIGITS<DIGITS)
                                    file << setprecision(tmp_DIGITS);
                                else
                                    file << setprecision(DIGITS);
                                file << setprecision(Int(DIGITS3+log(zeros_S[i][0]+2)/2.3));
                                file << zeros_S[i][0]*zeros_S[i][2] << endl;

                            }
                            i++; count2++;
                        }while(i<=10&&count2<count);
                        for(i=11;i<=ii;i++){
                            zeros_S[i-10][0]=zeros_S[i][0];
                            zeros_S[i-10][2]=zeros_S[i][2];
                            zeros_S[i-10][1]=zeros_S[i][1];
                        }
                        if(!strcmp(filename,"cout")) cout << setprecision(DIGITS3);
                        else file << setprecision(DIGITS3);
                        ii=ii-10;
                    }
                //}
            previous_zero=tmp;

            u=t;
            x=y;
            x_c=y_c;
        }while(count2<count);
        if(strcmp(filename,"cout")) file.close();

    }

    //elaborate zero finder keeps an array of zeros and goes back to find missing zeros
    //if they are detected via N(T) comparison.
    template <class ttype>
    void L_function <ttype>::
    find_zeros_elaborate(Double t1, Long count,Double max_refine,char* filename, char* message_stamp)
    {
        Double t,x,y;
        Double r;
        Double u,v,L_u,L_v,tmp;
        Long count2=0; //counts how many zeros have been printed
        Long count_all=0; //counts how many zeros have been found
        Double previous_zero=t1; 
        Double step_size,refined_step_size;
  
        fstream file;
    
        int i,j;
        Double zeros_S[100][2]; // stores list of consecutive zeros and S(T) data
        int ii; //counts number of local zeros
        int i1,i2; //used for looping through blocks for missing zeros
        Double average_S;
        Double global_average, local_average;
        int number_big_S; //count how many S(T)'s are suspiciously large
        int found_missing=0; //number missing zeros found
        int to_find=0; //number of missing zeros to find
        Double x2,y2,u2,v2,tmp2; //used to search for missing zeros
       
        Double sum_S=0; //sum of the S(T)'s.
        for(i=0;i<=99;i++)
        for(j=0;j<=1;j++) 
            zeros_S[i][j]=0;

       zeros_S[0][0]=t1; 
    
        if(!strcmp(filename,"cout"))
        {
            cout << setprecision(DIGITS3);
        }
        else
        {
            file.open(filename, ios::out|ios::app);
            file << setprecision(DIGITS3);
        }
    

        global_average=0;

/*

        Complex L_1;
        x=.1;i=0;
        do{
          L_1=this->value(x+I*(t1+.0001));
          if(L_1!=0) y=imag(log(L_1));


          if(i>0){
             tmp2=y-tmp;
             if(tmp2>5)tmp2=tmp2-2*Pi;
             if(tmp2<-5)tmp2=tmp2+2*Pi;
             global_average=global_average+tmp2;
cout << "global average: " <<x<<" "<< y << " " <<tmp << " " <<global_average/(2*Pi) << endl;
          }
          tmp=y;
          x=x+.01;
          i++;
        }while(x<=.9);
     
        global_average=-global_average/(2*Pi);

cout <<"hihih " << global_average<<endl;

*/


       

        //XXXXXXXXXXXXXXXXXx should I start slightly earlier
        t=t1; u=t;
        x=real(this->value(.5+I*t,0,"rotated pure"));

        if(my_verbose>1)
        cout << "look for sign change " << t << " " << x << endl;
    
    
        Double g=0;
        for(i=1;i<=this->a;i++)g=g+gamma[i];

        ii=0; //will be used to count zeros in the zeros_S array
        do{
 
            //step_size=.5/(this->N(t+2)/2-N(t+1)/2); //i.e. increment by half the average
                                                  //the local density

            
            if(count_all<7)
            step_size=g*log((this->Q)*abs(t)/(2*Pi)+4)/(2*Pi*40); 

            else
            step_size=g*log((this->Q)*abs(t)/(2*Pi)+4)/(2*Pi*2); 
            //i.e. increment by half the average
            //XXXXXX I need to take degree into account for the t aspect,
            //i.e. the sum of the gammas.


            t=t+step_size;
            y=real(this->value(.5+I*t,0,"rotated pure"));
            if(my_verbose>1)
            cout << "look for sign change " << t << " " << y << endl;
            v=t;
            if(sn(x)!=sn(y))  //if sign change is found...divide and conquer
                              //until the difference is small.
            {

                tmp=zeros_zoom_brent(x,y,u,v);
//cout << "catching zero: " << tmp << endl;
                count_all++; ii++;

                zeros_S[ii][0]=tmp; 
                //XXXX should not be discontinous
                if(t1>0)
                zeros_S[ii][1]=this->N((tmp+previous_zero)/2)/2 -N(t1)/2-(count_all-1);  
                else
                zeros_S[ii][1]=this->N((tmp+previous_zero)/2)/2 -(count_all-1);
                //zeros_S[ii][1]=this->N((tmp+previous_zero)/2)/2 -N(t1)/2-(count_all-1);  


                //if(ii>=20){
                    //count the number of big S(T)'s. If small then print the first
                    //10 zeros (but not more than count in total).
                    //Otherwise go back to look for missing zeros 
                    //searching with step sizes which are half as big. look for
                    //a few consecutive blocks where S(T) is large, and look for missing
                    //zeros in a neighborhood of the transition. If none are found
                    //halve the step size and search again. Do this a few times. If still 
                    //missing then search all i blocks with very small steps.
                    //If still not found then print a message and keep going.
                    //average_S=0;
                    //for(i=1;i<=ii;i++) average_S=average_S+zeros_S[i][1]/ii;


                    if(ii==7) for(i=ii;i>=ii-6;i--) global_average=global_average+zeros_S[i][1]/7;

                    if(ii>6){
                        local_average=0;
                        for(i=ii;i>=ii-6;i--) local_average=local_average+zeros_S[i][1]/7;
                        
                        if(local_average-global_average>.7) to_find=1;
                        if(local_average-global_average>1.5) to_find=Int(local_average-global_average+.5);
                        if(to_find>0){

//cout << "will try to find " << to_find << " zeros" << endl;
//cout <<"local vs global "<< local_average << " " << global_average << endl;
                            refined_step_size=step_size/2;
                            i1=ii-7; i2=i1; 


//for(i=1;i<=ii;i++) cout << "missing zeros detected: " << zeros_S[i][0] << " " << zeros_S[i][1]<< endl;

                            found_missing=0;
                            do{
    
                                if(i2>ii-1+found_missing)i2=ii-1+found_missing; 
                                if(i1<1&&count2>0)i1=1; if(i1<1&&count2==0)i1=0;
    
                                for(i=i1;i<=i2;i++){
                                    //search for missing zeros in block i with refined_step_size
                                    //each time a zero is found update all relevant info
                                    //cout << "searching block " << i << " stepsize " << refined_step_size << endl;
                                    u2=zeros_S[i][0]+refined_step_size/2;
                                    x2=real(this->value(.5+I*u2,0,"rotated pure"));
                                    do{
                                        v2=u2+refined_step_size;
    //cout <<"u2 v2: " << u2 << " " <<v2 << endl;
                                        if(v2<zeros_S[i+1][0])
                                        {
                                            y2=real(this->value(.5+I*v2,0,"rotated pure"));
                                            if(sn(x2)!=sn(y2))  //if sign change is found...divide and conquer
                                                              //until the difference is small.
                                            {
                                
    //cout <<"zoom in on: " << u2 << " " <<v2 << endl;
                                                tmp2=zeros_zoom_brent(x2,y2,u2,v2);
    //cout << " catching missing zero: " << tmp2 << endl;
                                                found_missing++;i2=i2+1;
                            //cout << "found a missing zero: " << tmp2 << " "<<endl;
                                                count_all++; 
                                                for(j=ii+found_missing;j>=i+2;j--){
                                                    zeros_S[j][0]=zeros_S[j-1][0];
                                                    zeros_S[j][1]=zeros_S[j-1][1]-1;
                                                }
                                                zeros_S[i+1][0]=tmp2;
                                                zeros_S[i+1][1]=zeros_S[i+2][1]-1; //not strictly correct, but close enough 
                                                i=i-2;if(i<1)i=1;
    //cout << "well blow my lips off: "<< zeros_S[i+1][0] << " " << zeros_S[i+1][1]<<endl;
    //cout << "well blow my lips off: "<< zeros_S[i+2][0] << " " << zeros_S[i+1][1]<<endl;
                                
                                            }
                                            x2=y2;u2=v2;
                                         }
                            
                                    }while(v2<zeros_S[i+1][0]);
                                
                                }
    
                                refined_step_size=refined_step_size/2;
                                i2=i2+2;i1=i1-2;
                                if(refined_step_size/step_size< .1){i2=i2+1;i1=i1-8;}
        
    
                            }while(found_missing<to_find && refined_step_size > step_size/max_refine);
                            to_find=0;
    
    
    
                            if(refined_step_size <= step_size/max_refine){
                                if(!strcmp(filename,"cout"))
                                {
                                    cout<<"missing zeros detected." << endl;
                                    cout<<"failed to find using refined step sizes of " << refined_step_size << endl;
                                    for(i=1;i<=ii+found_missing;i++)
                                        cout << zeros_S[i][0] << " " << zeros_S[i][1]<< endl;
                                }
                                else{
                                    file<<"missing zeros detected." << endl;
                                    file<<"failed to find using refined step sizes of " << refined_step_size << endl;
                                    for(i=1;i<=ii+found_missing;i++)
                                        file << zeros_S[i][0] << " " << zeros_S[i][1]<< endl;
                                }
                                file.close();
                                exit(1);
                            }


                        }

                    }

                    ii=ii+found_missing; found_missing=0;



                    number_big_S=0;
                    for(i=1;i<=ii;i++) if(zeros_S[i][1]>.9+global_average)number_big_S=number_big_S+1;
                    if(number_big_S>17) //i.e. if missing zero detected
                    {
//cout << "missing zero detected" << endl;
                   
                        //for(i=1;i<=ii;i++)
                            //cout << i << " : " <<zeros_S[i][0] << " " << zeros_S[i][1]<< endl;
                        refined_step_size=step_size/2;
                        i=ii;
                        //find transition zone 
                        do{
                            i--;
                        }while(zeros_S[i][1]>1&&i>=1);
                        i1=i; i2=i;

                        found_missing=0;
                        do{

                            if(i2>ii-1+found_missing)i2=ii-1+found_missing; 
                            if(i1<1&&count2>0)i1=1; if(i1<1&&count2==0)i1=0;
                            for(i=i1;i<=i2;i++){
                                //search for missing zeros in block i with refined_step_size
                                //each time a zero is found update all relevant info
                                //cout << "searching block " << i << " stepsize " << refined_step_size << endl;
                                u2=zeros_S[i][0]+refined_step_size/2;
                                x2=real(this->value(.5+I*u2,0,"rotated pure"));
                                do{
                                    v2=u2+refined_step_size;
//cout <<"u2 v2: " << u2 << " " <<v2 << endl;
                                    if(v2<zeros_S[i+1][0])
                                    {
                                        y2=real(this->value(.5+I*v2,0,"rotated pure"));
                                        if(sn(x2)!=sn(y2))  //if sign change is found...divide and conquer
                                                          //until the difference is small.
                                        {
                            
//cout <<"zoom in on: " << u2 << " " <<v2 << endl;
                                            tmp2=zeros_zoom_brent(x2,y2,u2,v2);
//cout << " catching missing zero: " << tmp2 << endl;
                                            found_missing++;i2=i2+1;
                        //cout << "found a missing zero: " << tmp2 << " "<<endl;
                                            count_all++; 
                                            for(j=ii+found_missing;j>=i+2;j--){
                                                zeros_S[j][0]=zeros_S[j-1][0];
                                                zeros_S[j][1]=zeros_S[j-1][1]-1;
                                            }
                                            zeros_S[i+1][0]=tmp2;
                                            zeros_S[i+1][1]=zeros_S[i+2][1]-1; //not strictly correct, but close enough 
                                            i=i-2;if(i<1)i=1;
//cout << "well blow my lips off: "<< zeros_S[i+1][0] << " " << zeros_S[i+1][1]<<endl;
//cout << "well blow my lips off: "<< zeros_S[i+2][0] << " " << zeros_S[i+1][1]<<endl;
                            
                                        }
                                        x2=y2;u2=v2;
                                     }
                        
                                }while(v2<zeros_S[i+1][0]);
                            
                            }
                            //average_S=0;
                            //for(i=1;i<=ii+found_missing;i++) average_S=average_S+zeros_S[i][1]/ii;
                            number_big_S=0;
                            for(i=1;i<=ii+found_missing;i++) if(zeros_S[i][1]>.9+global_average)number_big_S=number_big_S+1;

                            refined_step_size=refined_step_size/2;
                            i2=i2+2;i1=i1-3;
                            if(refined_step_size/step_size< .24){i2=i2+1;i1=i1-8;}
    

                        }while(number_big_S>17 && refined_step_size > step_size/max_refine);
                        if(refined_step_size <= step_size/max_refine){
                            if(!strcmp(filename,"cout"))
                            {
                                cout<<"missing zeros detected." << endl;
                                cout<<"failed to find using refined step sizes of " << refined_step_size << endl;
                                for(i=1;i<=ii+found_missing;i++)
                                    cout << zeros_S[i][0] << " " << zeros_S[i][1]<< endl;
                            }
                            else{
                                file<<"missing zeros detected." << endl;
                                file<<"failed to find using refined step sizes of " << refined_step_size << endl;
                                for(i=1;i<=ii+found_missing;i++)
                                    file << zeros_S[i][0] << " " << zeros_S[i][1]<< endl;
                            }
                            file.close();
                            exit(1);
                        }
        
                    }
                    //dump up to 10 zeros
                    ii=ii+found_missing; found_missing=0;
                    if(ii>=50){
                        i=1;
                        do{
                            if(!strcmp(filename,"cout"))
                            {
                                sum_S=sum_S+zeros_S[i][1];
                                cout << message_stamp << " ";
                                cout << setprecision(DIGITS3);
                                cout << zeros_S[i][0] << " ";
                                cout << setprecision(4) << zeros_S[i][1]<<" " << sum_S/(count2+1.)<< " "<<global_average<<" ";
                                cout << real(this->value(.5+I*zeros_S[i][0],0,"rotated pure")) << endl;
                                //cout << zeros_S[i][0] << endl;
  
                            }
                            else
                            {
                                file << message_stamp << " ";
                                //file << zeros_S[i][0] << " " << zeros_S[i][1]<< endl;
                                file << zeros_S[i][0] << endl;
                            }
                            i++; count2++;
                        }while(i<=10&&count2<count);
                        for(i=11;i<=ii;i++){
                            zeros_S[i-10][0]=zeros_S[i][0];
                            zeros_S[i-10][1]=zeros_S[i][1];
    //cout << "buffered: " <<zeros_S[i-10][0] << " " << zeros_S[i-10][1]<< endl;
                        }
                        if(!strcmp(filename,"cout")) cout << setprecision(DIGITS3);
                        else file << setprecision(DIGITS3);
                        ii=ii-10;
                    }
                //}
                previous_zero=tmp;

            }
            u=t;
            x=y;
        }while(count2<count);
        if(strcmp(filename,"cout")) file.close(); 

    }


Generated by  Doxygen 1.6.0   Back to index