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