Lcommandline_values_zeros.cc 5.86 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
/*

   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.

*/



25
#include "L.h"
26 27 28 29
#include "Lcommandline_values_zeros.h"


//-----functions--------------------------------------------------------------
30
void compute_values(Double x,Double y,const char *return_type,const char *file_name,Double x3,Double y3,Long count)
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
{

    fstream file;

    Complex s,u;
    Double x2,y2;

    //cout << setprecision(DIGITS3);

    if(strcmp(file_name,"")) //i.e. if a file_name has been provided...
    {
        file.open(file_name, ios::in);
        if(!file.is_open())
        {
            cout << endl << "FAILED TO FIND THE FILE: " << file_name<< endl;
            exit(1);
        }

        while (!file.eof())
        {
            file >> x2;
            file >> y2;
            s=Complex(x2,y2);
            if(!file.eof())
            switch(current_L_type)
            {
                case 1:
58
                    u=int_L.value(s,global_derivative,return_type); 
59 60 61 62 63 64
                    //global_derivative specifies which derivative to compute. Default is 0 which
                    //is just the plain L-value
                    cout << real(u) << " " << imag(u) << endl;
                    //cout << imag(s) << " " << abs(u)<< endl;
                    break;
                case 2:
65
                    u=Double_L.value(s,global_derivative,return_type);
66 67 68
                    cout << real(u) << " " << imag(u) << endl;
                    break;
                case 3:
69
                    u=Complex_L.value(s,global_derivative,return_type);
70 71 72 73 74 75 76 77 78 79 80
                    cout << real(u) << " " << imag(u) << endl;
                    break;
            }
        }
        file.close();
    }

    else  // else use the x and y provided
    {
        Long n=0;
        s=Complex(x,y);
81 82
        //#pragma omp sections
        {
83 84 85 86 87 88 89 90
        do{
            n++;
            cout << setprecision(DIGITS);
            if(count>0)cout << real(s) << " " << imag(s) <<" ";
            cout << setprecision(DIGITS3);
            switch(current_L_type)
            {
                case 1:
91
                    u=int_L.value(s,global_derivative,return_type);
92 93 94 95
                    cout << real(u) << " " << imag(u) << endl;
                    //cout << imag(s) << " " << abs(u)<< endl;
                    break;
                case 2:
96
                    u=Double_L.value(s,global_derivative,return_type);
97 98 99
                    cout << real(u) << " " << imag(u) << endl;
                    break;
                case 3:
100
                    u=Complex_L.value(s,global_derivative,return_type);
101 102 103 104 105
                    cout << real(u) << " " << imag(u) << endl;
                    break;
            }
            s=s+Complex(x3-x,y3-y)/(double)count;
        } while(n<count);
106
        }
107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
    }


}


//------------------------------------------------
//  Zero finding code

//looks for zeros on the critical line starting at 1/2+It1 by
//advancing in steps of size step_size (which can be negative)
//and looking for sign changes. If count==0 we look for
//zeros between t1=x and t2=y. If count>0, we look for the first
//count zeros (starting at t1 not yet implemented) using steps of no smaller
//than step_size.
//if count >0 then t2 is not used, so it's value is irrelevant.

124
void compute_zeros(Double x, Double y,Double step_size, Long count,int rank, bool test_explicit_formula)
125 126 127 128 129 130
{

    if(current_L_type==1){
        if(count==0)
        int_L.find_zeros(x,y,step_size);
        //else int_L.find_zeros_via_gram(x,count,step_size);
131
        else int_L.find_zeros_via_N(count,false,step_size,rank,test_explicit_formula);
132 133 134 135 136 137
        //else int_L.find_zeros_elaborate(x,count,step_size);
    }
    if(current_L_type==2){
        if(count==0)
        Double_L.find_zeros(x,y,step_size);
        //else Double_L.find_zeros_via_gram(x,count,step_size);
138
        else Double_L.find_zeros_via_N(count,false,step_size,rank,test_explicit_formula);
139 140 141 142 143 144
        //else Double_L.find_zeros_elaborate(x,count,step_size);
    }
    if(current_L_type==3){
        if(count==0)
        Complex_L.find_zeros(x,y,step_size);
        //Complex_L.find_zeros_via_gram(x,count,step_size);
145
        Complex_L.find_zeros_via_N(count,true,step_size,rank,test_explicit_formula);
146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
    }

}


//find the zeros of the function obtained by interpolating L and L2. Assumes
//they are of the same quasi degree and have same number of dirichlet 
//coefficients. The interpolated L-function has all its basic data as t times the
//data of the first L-function and (1-t) times the data of the second L-function.
//does not work with integer L_functions since interpolation gives at least doubles
void L_interpolate(Double x, Double y,Double step_size, int n)
{

    double t;
    char message_stamp[300];
    ostrstream os2(message_stamp,300);


    t=0;

    do{

        os2.seekp(0);
        os2 << t << ends;

        if(current_L_type==2){
            Double_L3=Double_L*(1-t)+Double_L2*t;
            Double_L3.find_zeros(x,y,step_size,"cout",message_stamp);
        }
        if(current_L_type==3){
            Complex_L3=Complex_L*(1-t)+Complex_L2*t;
            Complex_L3.find_zeros(x,y,step_size,"cout",message_stamp);
        }
        t=t+1./n;
    }while(t<=1);

}