source: sage/libs/ec/equation.c @ 0:039f6310c6fe

Revision 0:039f6310c6fe, 12.4 KB checked in by tornaria@…, 7 years ago (diff)

[project @ original sage-0.10.12]
Line 
1#include <stdio.h>
2#include <math.h>
3#include <stdlib.h>
4#include "ellcurv.h"
5
6GEN avecfromc4c6(GEN c4,GEN c6)
7{GEN a1,a2,a3,a4,a6,b2,AI,BI,DI,thevec,z; int T=0; int g;
8
9 a1=gmod(c4,gdeux); a2=gmod(gneg(gadd(a1,c6)),stoi(3));
10 if (gequal(a2,gdeux)) a2=gneg(gun);
11 T=0; DI=gcd0(gdiv(gsub(gmul(c4,gsqr(c4)),gsqr(c6)),stoi(1728)),gsqr(c6),0);
12 g=ggval(DI,gdeux)/12;
13 AI=gmul2n(c4,-g*4); BI=gmul2n(c6,-g*6);
14 if ((gequal(gmod(BI,gsqr(gdeux)),stoi(3))==1) &&
15     (gequal(gmod(AI,gdeux),gun)==1))
16   T=itos(gmod(gmul2n(gadd(BI,gun),-2),gsqr(gdeux)));
17 z=gmod(BI,gmul2n(gun,5));
18 if ((gequal(gmod(AI,gmul2n(gun,4)),gzero)==1) &&
19     (gequal(z,gzero) || (gequal(z,gmul2n(gun,3))))==1)
20   T=itos(gmod(gmul2n(BI,-3),gdeux));
21 if (T&1) a3=gsub(gun,gmod(gmul(a1,a2),gdeux));
22 else a3=gmod(gmul(a1,a2),gdeux);
23 b2=gadd(a1,gmul2n(a2,2));
24 a4=gdiv(gsub(gsqr(b2),gadd(c4,gmul(stoi(24),gmul(a1,a3)))),stoi(48));
25 a6=gdiv(gsub(gmul(stoi(36),gmul(b2,gadd(gmul(a1,a3),gmul2n(a4,1)))),
26              gadd(gmul(b2,gsqr(b2)),gadd(c6,gmul(stoi(216),a3)))),
27         stoi(864));
28 thevec=cgetg(6,t_VEC);
29 thevec[1]=(long) a1; thevec[2]=(long) a2; thevec[3]=(long) a3;
30 thevec[4]=(long) a4; thevec[5]=(long) a6;
31 //output(a1); output(a2); output(a3); output(a4); output(a6);
32 return(ellinit0(thevec,1,ELLACC));
33}
34
35
36GEN mineqfromc4c6(GEN c4,GEN c6)
37{GEN AI,BI,F,disc,g,u,p,z; int i,j,number,d;
38
39 u=gcd0(denom(c4),denom(c6),0); c4=gmul(c4,gsqr(gsqr(u)));
40 c6=gmul(c6,gmul(gsqr(u),gsqr(gsqr(u))));
41 disc=gdiv(gsub(gmul(c4,gsqr(c4)),gsqr(c6)),stoi(1728));
42 g=gabs(gcd0(gsqr(c6),disc,0),-1); F=factor(gmul(stoi(6),g));
43   ((GEN)F[2])[1] = (long) gsub((GEN) ((GEN) F[2])[1],gun);
44   ((GEN)F[2])[2] = (long) gsub((GEN) ((GEN) F[2])[2],gun); 
45
46 number=itos((GEN) matsize(F)[1]); u=gun;
47 for(i=1;i<=number;i++)
48 {p=(GEN) ((GEN) F[1])[i];
49  d=itos(gfloor(gdiv((GEN) ((GEN) F[2])[i],stoi(12))));
50  if (gequal(p,gdeux))
51  {AI=gmul2n(c4,-4*d); BI=gmul2n(c6,-6*d);
52   z=gmod(BI,gmul2n(gun,5));
53   if ((!gequal(gmod(AI,gdeux),gun) || (!gequal(gmod(BI,gsqr(gdeux)),stoi(3))))
54       && ((!gequal(gmod(AI,gmul2n(gun,4)),gzero)) ||
55           (!gequal(z,gzero) && (!gequal(z,gmul2n(gun,3)))))==1) d--;
56  }
57  if (gequal(p,stoi(3))) {if (ggval(c6,stoi(3))==6*d+2) d--;}
58  if (d==-1)
59  {c4=gmul(c4,gsqr(gsqr(p))); c6=gmul(c6,gmul(gsqr(p),gsqr(gsqr(p))));}
60  for (j=0;j<d;j++)
61  {c4=gdiv(c4,gsqr(gsqr(p))); c6=gdiv(c6,gmul(gsqr(gsqr(p)),gsqr(p)));}
62 }
63 //printf("OUT\n"); output(c4); output(c6);
64 return(avecfromc4c6(c4,c6));
65}
66
67
68GEN qtwist(GEN ellc,GEN p)
69{GEN c4; GEN c6;
70
71 if (gequal(p,gun)) return(ellc);
72 c4=gmul((GEN) ellc[10],gsqr(p));
73 c6=gmul(p,gmul((GEN) ellc[11],(GEN) gsqr(p)));
74 return(mineqfromc4c6(c4,c6));
75}
76
77/* #define DBG(arg) printf("retminimaltwist: %s\n",arg); */
78#define DBG(arg)
79
80/* minimize twist at primes bigger than 3, then does 3, then does 2. */
81GEN retminimaltwist(GEN EC,GEN *TWEF)
82{
83  GEN ELLAP,TEMP,currp,INVOL,OUTVOL;
84  int number,i,kodsym,twocond,locv,toptwo;
85
86 number = itos((GEN) matsize(CF)[1]); *TWEF=gun;
87
88 DBG("1");
89
90 for(i=1;i<=number;i++)
91   {
92     DBG("2");
93     currp=(GEN) ((GEN) CF[1])[i];
94     if (mpcmp(currp,gdeux)==0)
95       {twocond=itos((GEN) ((GEN) CF[2])[i]);
96        if ((twocond==4) || (twocond==6))
97          {
98            toptwo=twocond; INVOL=myvol(EC);
99            if (twocond==6) EC=qtwist(EC,gdeux);
100            twocond=itos(gel(elllocalred(EC,gdeux),1));
101            if (twocond==4) EC=qtwist(EC,gneg(gun));
102            toptwo-=itos(gel(elllocalred(EC,gdeux),1));
103            OUTVOL=myvol(EC); *TWEF=gmul2n(volratio(OUTVOL,INVOL),toptwo);
104            twocond=itos(gel(elllocalred(EC,gdeux),1));
105            if (twocond==0) {
106              ELLAP=ellap0(EC,gdeux,0); ELLAP=gsub(gmul(ELLAP,ELLAP),gdeux);
107              TEMP=gsub(gmul(gdeux,gsqr(gdeux)),gmul(gdeux,ELLAP));
108              TEMP=gsub(gadd(TEMP,ELLAP),gun); *TWEF=gmul2n(gmul(*TWEF,TEMP),-3);
109            }
110            if (twocond==1)
111              {TEMP=gmul2n(gadd(gdeux,gun),-2); *TWEF=gmul(*TWEF,TEMP);}
112          }
113        if (twocond>=7)
114          {if ((ggval(gel(EC,10),gdeux)>=6) &&
115               (ggval((GEN) (EC)[11],gdeux)>=6) &&
116               (ggval((GEN) (EC)[12],gdeux)>=6))
117             {EC=qtwist(EC,gdeux); *TWEF=gmul(*TWEF,gdeux);}
118         }
119        if ((twocond>=5) && (gsigne((GEN) EC[12])==-1)) EC=qtwist(EC,gneg(gun));
120      }
121     else {
122       kodsym=itos(gel(elllocalred(EC,currp),2));
123       if ((kodsym==-1) || (kodsym<=-5)) {
124         if (mod4(currp)==3) currp=gneg(currp); 
125         EC=qtwist(EC,currp);
126         if (gsigne(currp)==-1) currp=gneg(currp);
127         locv=itos(gel(elllocalred(EC,currp),1));
128         if (locv==0) {
129           ELLAP=ellap0(EC,currp,0); 
130           ELLAP=gsub(gmul(ELLAP,ELLAP),currp);
131           TEMP=gsub(gmul(currp,gsqr(currp)),gmul(currp,ELLAP));
132           TEMP=gsub(gadd(TEMP,ELLAP),gun);
133           *TWEF=gmul(*TWEF,TEMP); // (p+1+ap)(p+1-ap)(p-1) when good red
134         }
135         if (locv==1) {
136           TEMP=gsub(gsqr(currp),gun); *TWEF=gmul(*TWEF,TEMP);
137         }
138       }
139       else if (ggval((GEN) EC[10],currp)>=2) {
140         if (ggval((GEN) EC[11],currp)>=3)
141           {if (!gequal(currp,stoi(3)))
142              {if (mod4(currp)==3) currp=gneg(currp); EC=qtwist(EC,currp);
143               if (gsigne(currp)==-1) currp=gneg(currp);
144               *TWEF=gmul(*TWEF,currp);
145               }
146           else
147             {if ((ggval((GEN) EC[11],currp)!=5) &&
148                  (ggval((GEN) EC[12],currp)>=6))
149                {EC=qtwist(EC,gneg(currp)); *TWEF=gmul(*TWEF,currp);}
150            }
151          }
152       }
153     }
154   }
155 return(EC);
156}
157
158void minimaltwist()
159{GEN ELLAP,TEMP,currp,INVOL,OUTVOL;
160 int number,i,kodsym,twocond,locv,toptwo;
161
162 number=itos((GEN) matsize(CF)[1]); TWCOND=gcopy(COND);
163 TWCURVE=gcopy(CURVE); TWPROD=gun; TWEFF=gun; TWCF=gcopy(CF);
164
165 for(i=1;i<=number;i++)
166 {currp=(GEN) ((GEN) CF[1])[i];
167  if (mpcmp(currp,gdeux)==0)
168  {twocond=itos((GEN) ((GEN) CF[2])[i]);
169   if ((twocond==4) || (twocond==6))
170   {toptwo=twocond; INVOL=myvol(TWCURVE);
171    if (twocond==6)
172    {TWCURVE=qtwist(TWCURVE,gdeux); TWPROD=gdeux;}
173    twocond=itos(gel(elllocalred(TWCURVE,currp),1));
174    if (twocond==4)
175    {TWCURVE=qtwist(TWCURVE,gneg(gun)); TWPROD=gneg(TWPROD);}
176    twocond=itos(gel(elllocalred(TWCURVE,currp),1));
177    TWCOND=gmul2n(TWCOND,twocond-toptwo); OUTVOL=myvol(TWCURVE);
178    TWEFF=gmul2n(volratio(OUTVOL,INVOL),toptwo-twocond);
179    TWPROD=gmul(TWPROD,gsqr(gdeux));
180    if(twocond==0)
181    {ELLAP=ellap0(TWCURVE,gdeux,0); ELLAP=gsub(gmul(ELLAP,ELLAP),gdeux);
182     TEMP=gsub(gmul(gdeux,gsqr(gdeux)),gmul(gdeux,ELLAP));
183     TEMP=gsub(gadd(TEMP,ELLAP),gun); TWEFF=gmul2n(gmul(TWEFF,TEMP),-3);
184    }
185    if (twocond==1)
186    {TEMP=gmul2n(gadd(gdeux,gun),-2); TWEFF=gmul(TWEFF,TEMP);}
187   }
188   if (twocond>=7)
189   {if ((ggval((GEN) TWCURVE[10],gdeux)>=6) &&
190        (ggval((GEN) TWCURVE[11],gdeux)>=6) &&
191        (ggval((GEN) TWCURVE[12],gdeux)>=6))
192    {TWCURVE=qtwist(TWCURVE,gdeux); TWPROD=gmul(TWPROD,gdeux);
193     TWEFF=gmul(TWEFF,gdeux);
194    }
195   }
196   if ((twocond>=5) && (gsigne((GEN) (TWCURVE)[11])==-1))
197   {TWCURVE=qtwist(TWCURVE,gneg(gun)); TWPROD=gneg(TWPROD);}
198   ((GEN) TWCF[2])[i]=(long) stoi(twocond);
199  }
200  else
201  {kodsym=itos(gel(elllocalred(TWCURVE,currp),2));
202   if ((kodsym==-1) || (kodsym<=-5))
203   {if (mod4(currp)==3) currp=gneg(currp);
204    TWCURVE=qtwist(TWCURVE,currp); TWPROD=gmul(TWPROD,currp);
205    if (gsigne(currp)==-1) currp=gneg(currp);
206    locv=itos(gel(elllocalred(TWCURVE,currp),1));
207    if (locv==0)
208    {ELLAP=ellap0(TWCURVE,currp,0); ELLAP=gsub(gmul(ELLAP,ELLAP),currp);
209     TEMP=gsub(gmul(currp,gsqr(currp)),gmul(currp,ELLAP));
210     TEMP=gsub(gadd(TEMP,ELLAP),gun); TWEFF=gmul(TWEFF,TEMP);
211     TWCOND=gdiv(TWCOND,gsqr(currp)); // (p+1+ap)(p+1-ap)(p-1) when good red
212    }
213    if (locv==1)
214    {TEMP=gsub(gsqr(currp),gun); TWEFF=gmul(TWEFF,TEMP);
215     TWCOND=gdiv(TWCOND,currp);
216    }
217    ((GEN) TWCF[2])[i]=(long) stoi(locv);   
218   }
219   else if (ggval((GEN) TWCURVE[10],currp)>=2)
220   {if (ggval((GEN) TWCURVE[11],currp)>=3)
221    {if (!gequal(currp,stoi(3)))
222     {if (mod4(currp)==3) currp=gneg(currp);
223      TWCURVE=qtwist(TWCURVE,currp); TWPROD=gmul(TWPROD,currp);
224      if (gsigne(currp)==-1) currp=gneg(currp); TWEFF=gmul(TWEFF,currp);
225     }
226     else
227     {if ((ggval((GEN) TWCURVE[11],currp)!=5) &&
228          (ggval((GEN) TWCURVE[12],currp)>=6))
229      {TWCURVE=qtwist(TWCURVE,gneg(currp));
230       TWPROD=gmul(TWPROD,gneg(currp)); TWEFF=gmul(TWEFF,currp);
231      }
232     }
233    }
234   }
235  }
236 }
237
238 if (gequal(gmod(TWPROD,gsqr(gdeux)),gdeux)) TWPROD=gmul2n(TWPROD,2);
239 if (gequal(TWPROD,gneg(gun))) TWPROD=gneg(gsqr(gdeux));
240 ISSQFREE=1; ISPRIME=0;
241 for (i=1;i<=itos((GEN) matsize(TWCF)[1]);i++)
242 {locv=itos((GEN) ((GEN) CF[2])[i]); if (locv>1) ISSQFREE=0; ISPRIME+=locv;}
243 if (ISPRIME>1) ISPRIME=0;
244
245 if (0 && VERBOSE)
246 {if (gequal(TWPROD,gun)==1) return;
247  printf("Minimal Rational Quadratic Twist is ");
248  printcurven(TWCURVE); printf("of conductor "); output((GEN) (TWCOND));
249  printf("Product of Twisted Primes is "); output(TWPROD);
250  printf("Twists give a factor of "); output(TWEFF);
251 }
252}
253
254GEN minimalequation()
255{int i,number; GEN p,l;
256
257 CURVE=mineqfromc4c6((GEN) CURVE[10],(GEN) CURVE[11]);
258 CF=factor(gabs((GEN) CURVE[12],-1)); number=itos((GEN) matsize(CF)[1]);
259 for(i=1;i<=number;i++)
260 {p=(GEN) ((GEN) CF[1])[i];
261  ROOTNO*=ellrootno(CURVE,p); l=elllocalred(CURVE,p);
262  COND=gmul(COND,gpow(p,(GEN) l[1],-1)); TAMA=gmul(TAMA,(GEN) l[4]);
263  ((GEN) CF[2])[i]=lcopy((GEN) l[1]);
264 }
265 return(CURVE);
266}
267
268GEN gettama(GEN C)
269{int i,n; GEN p,l,T;
270
271 T=gun; n=itos((GEN) matsize(CF)[1]);
272 for(i=1;i<=n;i++)
273 {p=(GEN) ((GEN) CF[1])[i]; l=elllocalred(C,p); T=gmul(T,(GEN) l[4]);}
274 return(T);
275}
276
277void symsq()
278{int n,i,myp,mypow; GEN p,pow;
279
280 SSCOND=gcopy(TWCOND); EXPOS=gun; EXNEG=gun; EXEFF=gun;
281 n=itos((GEN) matsize(TWCF)[1]);
282 for(i=1;i<=n;i++)
283 {p=(GEN) ((GEN) TWCF[1])[i];  pow=(GEN) ((GEN) TWCF[2])[i];
284  if (itos(pow)>=2)
285  {myp=itos(p);
286   if (myp==2)
287   {mypow=itos(pow);
288    if (mypow==2)
289    {SSCOND=gmul2n(SSCOND,-1); EXPOS=gmul2n(EXPOS,1);
290     EXEFF=gmul2n(gmul(EXEFF,stoi(3)),-1);
291    }
292    if (mypow==3) {SSCOND=gmul2n(SSCOND,-1); /* Non-Exotic */}
293    if (mypow==5) {SSCOND=gmul2n(SSCOND,-2); /* Non-Exotic */}
294    if (mypow==7) {SSCOND=gmul2n(SSCOND,-3); /* Non-Exotic */}
295    if (mypow==8)
296    {if (ggval((GEN) TWCURVE[11],p)>=9)
297     {SSCOND=gmul2n(SSCOND,-4); /* Non-Exotic */}
298     else
299     {SSCOND=gmul2n(SSCOND,-5);
300      if (itos((GEN) gmod((GEN) TWCURVE[10],stoi(128)))==32)
301      {EXPOS=gmul2n(EXPOS,1); EXEFF=gmul2n(gmul(EXEFF,stoi(3)),-1);}
302      if (itos((GEN) gmod((GEN) TWCURVE[10],stoi(128)))==96)
303      {EXNEG=gmul2n(EXNEG,1); EXEFF=gmul2n(EXEFF,-1);}
304     }
305    }
306   }
307   else if (myp==3)
308   {mypow=itos(pow);
309    if (mypow==2)
310    {SSCOND=gdiv(SSCOND,stoi(3)); EXPOS=gmul(EXPOS,stoi(3));
311     EXEFF=gdiv(gmul2n(EXEFF,2),stoi(3));
312    }
313    if (mypow==3) {SSCOND=gdiv(SSCOND,stoi(3)); /* Non-Exotic */}
314    if (mypow==4)
315    {SSCOND=gdiv(SSCOND,stoi(9));
316     if (itos((GEN) gmod((GEN) TWCURVE[11],stoi(243)))==0)
317     {if (itos((GEN) gmod((GEN) TWCURVE[10],stoi(81)))==27)
318      {EXNEG=gmul(EXNEG,stoi(3)); EXEFF=gdiv(gmul2n(EXEFF,1),stoi(3));}
319      if (itos((GEN) gmod((GEN) TWCURVE[10],stoi(81)))==54)
320      {EXPOS=gmul(EXPOS,stoi(3)); EXEFF=gdiv(gmul2n(EXEFF,2),stoi(3));}
321     }
322     if (itos((GEN) gmod((GEN) TWCURVE[10],stoi(27)))==9)
323     {if (itos((GEN) gmod((GEN) TWCURVE[11],stoi(243)))==54)
324      {EXPOS=gmul(EXPOS,stoi(3)); EXEFF=gdiv(gmul2n(EXEFF,2),stoi(3));}
325      if (itos((GEN) gmod((GEN) TWCURVE[11],stoi(243)))==108)
326      {EXNEG=gmul(EXNEG,stoi(3)); EXEFF=gdiv(gmul2n(EXEFF,1),stoi(3));}
327      if (itos((GEN) gmod((GEN) TWCURVE[11],stoi(243)))==135)
328      {EXNEG=gmul(EXNEG,stoi(3)); EXEFF=gdiv(gmul2n(EXEFF,1),stoi(3));}
329      if (itos((GEN) gmod((GEN) TWCURVE[11],stoi(243)))==189)
330      {EXPOS=gmul(EXPOS,stoi(3)); EXEFF=gdiv(gmul2n(EXEFF,2),stoi(3));}
331     }
332    }
333    if (mypow==5) {SSCOND=gdiv(SSCOND,stoi(9)); /* Non-Exotic */}
334   }
335   else
336   {SSCOND=gdiv(SSCOND,p);
337    if ((myp%12)==1)
338    {EXNEG=gmul(EXNEG,p); EXEFF=gdiv(gmul(EXEFF,gsub(p,gun)),p);}
339    if ((myp%12)==5)
340    {if ((ggval((GEN) TWCURVE[10],p)<2) &&
341         (ggval((GEN) TWCURVE[11],p)>=2))
342     {EXNEG=gmul(EXNEG,p); EXEFF=gdiv(gmul(EXEFF,gsub(p,gun)),p);}
343     else
344     {EXPOS=gmul(EXPOS,p); EXEFF=gdiv(gmul(EXEFF,gadd(p,gun)),p);}
345    }
346    if ((myp%12)==7)
347    {if ((ggval((GEN) TWCURVE[11],p)<2) ||
348         (ggval((GEN) TWCURVE[10],p)>=2))
349     {EXNEG=gmul(EXNEG,p); EXEFF=gdiv(gmul(EXEFF,gsub(p,gun)),p);}
350     else
351     {EXPOS=gmul(EXPOS,p); EXEFF=gdiv(gmul(EXEFF,gadd(p,gun)),p);}
352    }
353    if ((myp%12)==11)
354    {EXPOS=gmul(EXPOS,p); EXEFF=gdiv(gmul(EXEFF,gadd(p,gun)),p);}
355   }
356  }
357 }
358
359 if (0 && (VERBOSE) && (PRINT))
360 {if (gequal(SSCOND,TWCOND)==1) return;
361  printf("Symmetric Square Conductor is "); output(SSCOND);
362  if (gequal(EXPOS,gun)!=1)
363  {printf("Plus product is "); output(EXPOS);}
364  if (gequal(EXNEG,gun)!=1)
365  {printf("Minus product is ");output(EXNEG);}
366  printf("Even valuation gives a factor of "); output(EXEFF);
367 }
368}
Note: See TracBrowser for help on using the repository browser.