| | 1 | /* |
| | 2 | |
| | 3 | Copyright (C) 2001,2002,2003,2004 Michael Rubinstein |
| | 4 | |
| | 5 | This file is part of the L-function package L. |
| | 6 | |
| | 7 | This program is free software; you can redistribute it and/or |
| | 8 | modify it under the terms of the GNU General Public License |
| | 9 | as published by the Free Software Foundation; either version 2 |
| | 10 | of the License, or (at your option) any later version. |
| | 11 | |
| | 12 | This program is distributed in the hope that it will be useful, |
| | 13 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
| | 14 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| | 15 | GNU General Public License for more details. |
| | 16 | |
| | 17 | Check the License for details. You should have received a copy of it, along |
| | 18 | with the package; see the file 'COPYING'. If not, write to the Free Software |
| | 19 | Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
| | 20 | |
| | 21 | */ |
| | 22 | |
| | 23 | #include "Lcommandline_elliptic.h" |
| | 24 | |
| | 25 | |
| | 26 | //returns 0 if initialization goes smoothly, 1 if there's an error |
| | 27 | // given the elliptic curve |
| | 28 | // y^2 + a1 xy + a3 y = x^3 + a2 x^2 + a4 x + a6 |
| | 29 | // i.e. computes, the sign, the conductor, and first N_terms dirichlet coefficients |
| | 30 | // of the corresponding L-function. The nth dirichlet coefficient is normalized by |
| | 31 | // sqrt(n) |
| | 32 | // |
| | 33 | |
| | 34 | #ifdef INCLUDE_PARI |
| | 35 | void initialize_new_L(char *a1, char *a2, char *a3, char *a4, char *a6, int N_terms) |
| | 36 | { |
| | 37 | |
| | 38 | // basic data for the L-function (see the class L_function for full comments) |
| | 39 | int what_type; |
| | 40 | Long Period; |
| | 41 | Double q; |
| | 42 | Complex w; |
| | 43 | int A; |
| | 44 | Double *g; |
| | 45 | Complex *l; |
| | 46 | int n_poles; |
| | 47 | Complex *p; |
| | 48 | Complex *r; |
| | 49 | |
| | 50 | current_L_type=2; //the normalized dirichlet coeffs are real |
| | 51 | |
| | 52 | what_type=2; //i.e. eliiptic curve L-functions are cusp form L-function |
| | 53 | |
| | 54 | Period=0; |
| | 55 | |
| | 56 | A=1; |
| | 57 | |
| | 58 | //GAMMA factor is GAMMA(s+1/2) |
| | 59 | g=new Double[2]; |
| | 60 | l=new Complex[2]; |
| | 61 | g[1]=1.; |
| | 62 | l[1]=.5; |
| | 63 | |
| | 64 | Double * coeff; |
| | 65 | coeff = new Double[N_terms+1]; |
| | 66 | |
| | 67 | data_E(a1,a2,a3,a4,a6,N_terms,coeff); |
| | 68 | //coeff[n], if n > 1, is the nth dirichlet coefficient nomalized by sqrt(n). |
| | 69 | //coeff[0] comes back with the sign of the functional equation |
| | 70 | //coeff[1] comes back with the conductor of E. We then set it to one. |
| | 71 | |
| | 72 | |
| | 73 | q=sqrt(coeff[1])/(2*Pi); |
| | 74 | coeff[1]=1.; |
| | 75 | |
| | 76 | w=coeff[0]; |
| | 77 | |
| | 78 | n_poles=0; //no poles |
| | 79 | p = new Complex[1]; |
| | 80 | r = new Complex[1]; |
| | 81 | |
| | 82 | Double_L=L_function<Double>("Elliptic curve",what_type,N_terms,coeff,Period,q,w,A,g,l,n_poles,p,r); |
| | 83 | |
| | 84 | delete [] g; |
| | 85 | delete [] l; |
| | 86 | delete [] p; |
| | 87 | delete [] r; |
| | 88 | delete [] coeff; |
| | 89 | |
| | 90 | |
| | 91 | } |
| | 92 | |
| | 93 | |
| | 94 | // SAGE -- used below -- needed for Cygwin. |
| | 95 | #ifndef llrint |
| | 96 | inline long long int llrint (double x) |
| | 97 | { |
| | 98 | long long int llrintres; |
| | 99 | asm |
| | 100 | ("fistpll %0" |
| | 101 | : "=m" (llrintres) : "t" (x) : "st"); |
| | 102 | return llrintres; |
| | 103 | } |
| | 104 | #endif |
| | 105 | |
| | 106 | |
| | 107 | // given the elliptic curve |
| | 108 | // y^2 + a1 xy + a3 y = x^3 + a2 x^2 + a4 x + a6 |
| | 109 | // i.e. computes, the sign, the conductor, and first N_terms dirichlet coefficients |
| | 110 | // of the corresponding L-function. The nth dirichlet coefficient is normalized by |
| | 111 | // sqrt(n) |
| | 112 | void data_E(char *a1, char *a2, char *a3, char *a4, char *a6, int N_terms,Double * coeff) |
| | 113 | { |
| | 114 | |
| | 115 | |
| | 116 | int sign; //sign stores the sign of the functional equation |
| | 117 | Long p; //denotes a prime |
| | 118 | Long m,n; |
| | 119 | Double x,r,tmp,tmp2; |
| | 120 | Long conductor; //the conductor of the elliptic curve |
| | 121 | |
| | 122 | GEN y, F, E, C; |
| | 123 | |
| | 124 | |
| | 125 | y = cgeti(64); |
| | 126 | |
| | 127 | C = cgetg(4, t_VEC); |
| | 128 | |
| | 129 | |
| | 130 | //if a1 etc are integers, we can use gaffsg to |
| | 131 | //assign F[1] etc. However, I am treating a1 etc as character |
| | 132 | //strings to allow for larger integers, and therefore use gaffect |
| | 133 | |
| | 134 | |
| | 135 | F = cgetg(6, t_VEC); |
| | 136 | F[1] = lgeti(BIGDEFAULTPREC); |
| | 137 | F[2] = lgeti(BIGDEFAULTPREC); |
| | 138 | F[3] = lgeti(BIGDEFAULTPREC); |
| | 139 | F[4] = lgeti(BIGDEFAULTPREC); |
| | 140 | F[5] = lgeti(BIGDEFAULTPREC); |
| | 141 | |
| | 142 | //gaffsg(a1,(GEN) F[1]); |
| | 143 | //gaffsg(a2,(GEN) F[2]); |
| | 144 | //gaffsg(a3,(GEN) F[3]); |
| | 145 | //gaffsg(a4,(GEN) F[4]); |
| | 146 | //gaffsg(a6,(GEN) F[5]); |
| | 147 | |
| | 148 | gaffect(strtoGEN(a1), (GEN) F[1]); |
| | 149 | gaffect(strtoGEN(a2), (GEN) F[2]); |
| | 150 | gaffect(strtoGEN(a3), (GEN) F[3]); |
| | 151 | gaffect(strtoGEN(a4), (GEN) F[4]); |
| | 152 | gaffect(strtoGEN(a6), (GEN) F[5]); |
| | 153 | |
| | 154 | E = initell(F,BIGDEFAULTPREC); |
| | 155 | |
| | 156 | C=globalreduction(E); |
| | 157 | |
| | 158 | x=gtodouble((GEN) C[1]); |
| | 159 | |
| | 160 | |
| | 161 | //if(x<1e18) conductor=Long(x+.1); |
| | 162 | if(x<Double(1.*my_LLONG_MAX)) conductor=Long(x+.1); |
| | 163 | |
| | 164 | else{ |
| | 165 | cout << "conductor equals " << x << " and is too large" << endl; |
| | 166 | exit(1); |
| | 167 | } |
| | 168 | |
| | 169 | |
| | 170 | gaffsg(1, (GEN) y); |
| | 171 | sign = ellrootno(E,y); //sign of the functional equation |
| | 172 | |
| | 173 | |
| | 174 | for(n=1;n<=N_terms;n++) coeff[n]=1.; |
| | 175 | |
| | 176 | n=2; |
| | 177 | do{ |
| | 178 | if(isprime(n)){ |
| | 179 | |
| | 180 | p=n; |
| | 181 | gaffsg(p,y); |
| | 182 | coeff[p] = Double(1.*llrint(gtodouble(apell(E,y))))/sqrt(Double(1.*p)); |
| | 183 | //coeff[p] = Double(1.*Long(gtodouble(apell(E,y))+.1))/sqrt(Double(1.*p)); |
| | 184 | |
| | 185 | if(gtolong(gmod((GEN) E[12],(GEN) y))==0) // if p|discriminant, i.e. bad reduction |
| | 186 | { |
| | 187 | |
| | 188 | tmp=coeff[p]; |
| | 189 | r=tmp*tmp; |
| | 190 | x=Double(1.*p)*Double(1.*p); |
| | 191 | m=Long(x+.1); |
| | 192 | if(m<=N_terms) |
| | 193 | do{ |
| | 194 | coeff[m]=coeff[m]*r; |
| | 195 | x=x*p; |
| | 196 | m=Long(x+.1); |
| | 197 | r=r*tmp; |
| | 198 | }while(m<=N_terms); |
| | 199 | |
| | 200 | } |
| | 201 | |
| | 202 | else{ // a(p^(j+1)) = a(p)a(p^j) - p a(p^(j-1)) so normalizng by sqrt |
| | 203 | // gives coeff[p^(j+1)] = coeff[p] coeff[p^j]-coeff[p^(j-1)] |
| | 204 | |
| | 205 | x=Double(1.*p)*Double(1.*p); |
| | 206 | m=Long(x+.1); |
| | 207 | if(m<=N_terms) |
| | 208 | do{ |
| | 209 | tmp2=0; |
| | 210 | coeff[m]=coeff[m]*(coeff[p]*coeff[m/p]-coeff[m/(p*p)]); |
| | 211 | x=x*p; |
| | 212 | m=Long(x+.1); |
| | 213 | r=r*tmp; |
| | 214 | }while(m<=N_terms); |
| | 215 | } |
| | 216 | } |
| | 217 | else{ |
| | 218 | |
| | 219 | p=1; |
| | 220 | do{ |
| | 221 | p++; |
| | 222 | }while(n%p!=0); |
| | 223 | |
| | 224 | m=p; |
| | 225 | do{ |
| | 226 | m=m*p; |
| | 227 | }while(n%m==0&&m<n); |
| | 228 | |
| | 229 | if(n%m!=0) m=m/p; |
| | 230 | |
| | 231 | coeff[n]=coeff[m]*coeff[n/m]; |
| | 232 | } |
| | 233 | //if(n%10000==1) cout << n << endl; |
| | 234 | n++; |
| | 235 | |
| | 236 | }while(n<=N_terms); |
| | 237 | |
| | 238 | coeff[0]=1.*sign; coeff[1]=1.*conductor; //the first two spots are used |
| | 239 | //for the sign and conductor. The rest are |
| | 240 | //used for the Dirichlet coefficients. |
| | 241 | |
| | 242 | } |
| | 243 | |
| | 244 | #endif //ifdef INCLUDE_PARI |
| | 245 | |
| | 246 | void compute_rank(){ |
| | 247 | switch(current_L_type) |
| | 248 | { |
| | 249 | case 1: |
| | 250 | int_L.compute_rank(true); |
| | 251 | break; |
| | 252 | case 2: |
| | 253 | Double_L.compute_rank(true); |
| | 254 | break; |
| | 255 | case 3: |
| | 256 | Complex_L.compute_rank(true); |
| | 257 | break; |
| | 258 | } |
| | 259 | } |
| | 260 | |
| | 261 | void verify_rank(int rank){ |
| | 262 | switch(current_L_type) |
| | 263 | { |
| | 264 | case 1: |
| | 265 | int_L.verify_rank(rank); |
| | 266 | break; |
| | 267 | case 2: |
| | 268 | Double_L.verify_rank(rank); |
| | 269 | break; |
| | 270 | case 3: |
| | 271 | Complex_L.verify_rank(rank); |
| | 272 | break; |
| | 273 | } |
| | 274 | } |
| | 275 | |
