1 #include "sympow.h"
2 #define DEBUG (FALSE || GLOBAL_DEBUG)
3
generic_plusminus(llint p,int m,int A,int B,QDpoly * v)4 static void generic_plusminus(llint p,int m,int A,int B,QDpoly *v)
5 {QD P,t; QDpoly v1,v2,w;
6 if (DEBUG) printf("generic_plusminus %lli %i %i %i\n",p,m,A,B);
7 QD_copy(wmax,QD_zero,P); P[0]=(double) p;
8 if (m&1)
9 {initQDpoly(&w,2); ASSERT(A==B);
10 QD_copy(wmax,QD_one,w.coeff[0]);
11 QD_powi(wmax,P,m,t); QD_copy(wmax,t,w.coeff[2]);
12 QDpoly_pow(w,A,v,-1); delQDpoly(&w); return;}
13 initQDpoly(&w,1); QD_copy(wmax,QD_one,w.coeff[0]);
14 QD_powi(wmax,P,m/2,t); QD_copy(wmax,t,w.coeff[1]);
15 QDpoly_pow(w,A,&v1,-1); QD_neg(wmax,w.coeff[1],w.coeff[1]);
16 QDpoly_pow(w,B,&v2,-1); delQDpoly(&w); QDpoly_mul(v1,v2,v,-1);
17 delQDpoly(&v1); delQDpoly(&v2);}
18
cyclic_nonabelian(llint p,int m,int eps,QDpoly * v)19 static void cyclic_nonabelian(llint p,int m,int eps,QDpoly *v)
20 {int A,B;
21 if (DEBUG) printf("cyclic_nonabelian %lli %i %i\n",p,m,eps);
22 if (m&1) {ASSERT((eps&1)==0); A=eps/2; B=eps/2;}
23 else if (m&2) {A=(eps+1)/2; B=A-1;} else {A=(eps-1)/2; B=A+1;}
24 return generic_plusminus(p,m,A,B,v);}
25
cyclic_abelian_ap(llint p,int ap,int m,int d,QDpoly * v)26 static void cyclic_abelian_ap(llint p,int ap,int m,int d,QDpoly *v)
27 {QDpoly qp,lp,A; QD P,t,u; int i;
28 initQDpoly(&A,0); QD_copy(wmax,QD_one,A.coeff[0]);
29 initQDpoly(&qp,2); QD_copy(wmax,QD_zero,P); P[0]=(double) p;
30 QD_powi(wmax,P,m,t); QD_copy(wmax,t,qp.coeff[2]);
31 QD_copy(wmax,QD_one,qp.coeff[0]);
32 for (i=0;i<(m+1)/2;i++)
33 {if (((2*i-m)%d)==0)
34 {power_trace_from_ap(wmax,p,ap,m-2*i,t);
35 QD_neg(wmax,t,t); QD_powi(wmax,P,i,u); QD_mul(wmax,t,u,qp.coeff[1]);
36 QDpoly_mul(A,qp,v,-1); delQDpoly(&A); A=*v;}}
37 delQDpoly(&qp);
38 if ((m&1)==0)
39 {initQDpoly(&lp,1); QD_copy(wmax,QD_one,lp.coeff[0]);
40 QD_powi(wmax,P,m/2,t); QD_neg(wmax,t,t); QD_copy(wmax,t,lp.coeff[1]);
41 QDpoly_mul(A,lp,v,-1); delQDpoly(&A); delQDpoly(&lp);}}
42
hecke_good(llint p,int ap,int m,QDpoly * v)43 static void hecke_good(llint p,int ap,int m,QDpoly *v)
44 {QD T; if (m==0) errorit("m=0 in hecke_good?!");
45 initQDpoly(v,2); QD_copy(wmax,QD_one,(*v).coeff[0]);
46 if (ap!=0)
47 {QD_copy(wmax,QD_zero,T); T[0]=(double) p; QD_powi(wmax,T,m,T);
48 power_trace_from_ap(wmax,p,ap,m,(*v).coeff[1]);
49 QD_neg(wmax,(*v).coeff[1],(*v).coeff[1]); QD_copy(wmax,T,(*v).coeff[2]);}
50 else
51 {QD_copy(wmax,QD_zero,T); T[0]=(double) -p; QD_powi(wmax,T,m,T);
52 QD_neg(wmax,T,(*v).coeff[2]); QD_copy(wmax,QD_zero,(*v).coeff[1]);}}
53
cyclic_abelian(llint p,int m,int d,QDpoly * v,int bpt)54 static void cyclic_abelian(llint p,int m,int d,QDpoly *v,int bpt)
55 {double P=(double) p; int v4,v6,ap,i,c4,c6; QD Et4,Et6;
56 if (DEBUG) printf("cyclic abelian p:%lli sp:%i Cd:%i bpt:%i\n",p,m,d,bpt);
57 if ((d>2) && (p<=3)) {cyclic_abelian_ap(p,p,m,d,v); return;}
58 if ((d>1) && (p>3))
59 {QD_copy(wmax,QD_zero,Et4); QD_copy(wmax,QD_zero,Et6);
60 v4=QD_valuation(Ec4,P); v6=QD_valuation(Ec6,P);
61 if (3*v4>=2*v6)
62 {QD_copy(wmax,Ec6,Et6); for (i=1;i<=v6;i++) QD_div1(wmax,Et6,P,Et6);}
63 if (3*v4<=2*v6)
64 {QD_copy(wmax,Ec4,Et4); for (i=1;i<=v4;i++) QD_div1(wmax,Et4,P,Et4);}}
65 if ((d==2) && (p==2)) /* need to worry here */
66 {if (bpt==16) {QD_div1(wmax,Ec4,16.0,Et4); QD_div1(wmax,Ec6,64.0,Et6);}
67 else
68 {QD_div1(wmax,Ec4,4.0,Et4); QD_div1(wmax,Ec6,8.0,Et6);
69 c4=QD_modi(Et4,1<<29); c6=QD_modi(Et6,1<<29);
70 if ((((c4&31)==16) && ((c6&255)==192)) ||
71 (((c4&255)==0) && ((c6&2047)==512)))
72 {QD_div1(wmax,Et4,16.0,Et4); QD_div1(wmax,Et6,64.0,Et6);}
73 else if ((((c4&31)==16) && ((c6&255)==64)) ||
74 (((c4&255)==0) && ((c6&2047)==1536)))
75 {QD_div1(wmax,Et4,16.0,Et4); QD_div1(wmax,Et6,-64.0,Et6);}
76 }
77 }
78 if ((d==2) && (p==3)) /* think this is OK */
79 {QD_div1(wmax,Ec4,9.0,Et4); QD_div1(wmax,Ec6,-27.0,Et6);}
80 if (d==1) ap=(int) ec_do(p); else ap=(int) ec_ap(Et4,Et6,p);
81 if ((HECKE) && (d==1)) return hecke_good(p,ap,m,v);
82 cyclic_abelian_ap(p,ap,m,d,v);}
83
euler_factor_good(llint p,int m,QDpoly * v)84 static void euler_factor_good(llint p,int m,QDpoly *v)
85 {cyclic_abelian(p,m,1,v,0);}
86
euler_factor_bad(llint p,int bpt,int m,QDpoly * v)87 void euler_factor_bad(llint p,int bpt,int m,QDpoly *v)
88 {int ap,eps,lc; if (DEBUG) printf("euler_factor_bad %lli %i %i\n",p,bpt,m);
89 lc=tame_local_conductor(bpt,m); eps=m+1-lc;
90 if (lc==m+1) {initQDpoly(v,0); QD_copy(wmax,QD_one,(*v).coeff[0]); return;}
91 if ((bpt==1) || (bpt>=29))
92 {ap=ec_ap(Ec4,Ec6,p); initQDpoly(v,1); QD_copy(wmax,QD_one,(*v).coeff[0]);
93 if ((ap>=0) || !(m&1)) QD_neg(wmax,(*v).coeff[0],(*v).coeff[1]);
94 else QD_copy(wmax,QD_one,(*v).coeff[1]); return;}
95 if ((bpt>=2) && (bpt<=6))
96 {if ((p%bpt)!=1) cyclic_nonabelian(p,m,eps,v);
97 else cyclic_abelian(p,m,bpt,v,bpt); return;}
98 if (((bpt>=8) && (bpt<=10)) || ((bpt>=24) && (bpt<=27)))
99 {if (((m&7)==2) || ((m&7)==4)) generic_plusminus(2,m,eps/2,eps/2,v);
100 if ((m&7)==6) generic_plusminus(2,m,(eps+1)/2,(eps-1)/2,v);
101 if ((m&7)==0) generic_plusminus(2,m,(eps-1)/2,(eps+1)/2,v); return;
102 }
103 if (bpt==22) {cyclic_abelian(2,m,4,v,bpt); return;}
104 if (bpt==23) {cyclic_nonabelian(p,m,eps,v); return;}
105 if ((bpt==12) || (bpt==13))
106 {if ((eps&1)==0) generic_plusminus(3,m,eps/2,eps/2,v);
107 else if ((m&3)==0) generic_plusminus(3,m,(eps-1)/2,(eps+1)/2,v);
108 else if ((m&3)==2) generic_plusminus(3,m,(eps+1)/2,(eps-1)/2,v);}
109 if ((bpt==16) || (bpt==17)) cyclic_abelian(2,m,2,v,bpt);
110 if (bpt==14) cyclic_abelian(3,m,3,v,bpt);
111 if (bpt==15) cyclic_abelian(3,m,6,v,bpt);
112 if ((bpt==18) || (bpt==19)) cyclic_nonabelian(p,m,eps,v);
113 if ((bpt==20) || (bpt==21)) cyclic_nonabelian(p,m,eps,v); return;}
114
deflate(int CM)115 static int deflate(int CM)
116 {if ((CM==-27) || (CM== -12)) return 3; if (CM==-28) return 7;
117 if (CM==-16) return 4; return -CM;}
118
euler_factor_hecke_bad(llint p,int bpt,int m,QDpoly * v)119 void euler_factor_hecke_bad(llint p,int bpt,int m,QDpoly *v)
120 {QDpoly ef1,ef2,poly,R; int i,k; QD P,T; int A[8]={0,1,0,-1,0,-1,0,1};
121 if (DEBUG) printf("euler_factor_hecke_bad p:%lli bpt:%i sp:%i\n",p,bpt,m);
122 QD_copy(wmax,QD_zero,P); P[0]=(double) p; euler_factor_bad(p,bpt,m,&ef1);
123 if (m>2) euler_factor_bad(p,bpt,m-2,&ef2);
124 else {initQDpoly(&ef2,0); QD_copy(wmax,QD_one,ef2.coeff[0]);}
125 if ((m&1)==0)
126 {initQDpoly(&poly,1); QD_copy(wmax,QD_one,poly.coeff[0]);
127 k=-deflate(CM_CASE); if (p>2) k=kronll((llint) k,p); else k=A[k&7];
128 if ((m&3)==2)
129 {QD_powi(wmax,P,m/2-1,T); if (k==1) QD_neg(wmax,T,T);
130 if (k==0) poly.deg=0; else QD_copy(wmax,T,poly.coeff[1]);
131 QDpoly_mul(ef2,poly,&R,-1); delQDpoly(&ef2); delQDpoly(&poly); ef2=R;}
132 if ((m&3)==0)
133 {QD_powi(wmax,P,m/2,T); if (k==1) QD_neg(wmax,T,T);
134 if (k==0) poly.deg=0; else QD_copy(wmax,T,poly.coeff[1]);
135 QDpoly_mul(ef1,poly,&R,-1); delQDpoly(&ef1); delQDpoly(&poly); ef1=R;}
136 initQDpoly(&poly,1); QD_copy(wmax,QD_one,poly.coeff[0]);
137 if ((m&3)==0)
138 {QD_powi(wmax,P,m/2-1,T); QD_neg(wmax,T,poly.coeff[1]);
139 QDpoly_mul(ef2,poly,&R,-1); delQDpoly(&ef2); delQDpoly(&poly); ef2=R;}
140 if ((m&3)==2)
141 {if (m<=2) delQDpoly(&poly);
142 else
143 {QD_powi(wmax,P,m/2,T); QD_neg(wmax,T,poly.coeff[1]);
144 QDpoly_mul(ef1,poly,&R,-1); delQDpoly(&ef1); delQDpoly(&poly); ef1=R;}}}
145 QD_copy(wmax,P,T);
146 for (i=1;i<=ef2.deg;i++)
147 {QD_mul(wmax,ef2.coeff[i],T,ef2.coeff[i]); QD_mul(wmax,T,P,T);}
148 QDpoly_inv(ef2,ef1.deg-ef2.deg,&R); QDpoly_mul(ef1,R,v,ef1.deg-ef2.deg);
149 delQDpoly(&R); delQDpoly(&ef1); delQDpoly(&ef2);}
150
euler_factor(llint p,int m,QDpoly * v)151 void euler_factor(llint p,int m,QDpoly *v)
152 {int i;
153 if ((COND0%p)!=0) {euler_factor_good(p,m,v);}
154 else
155 {for (i=0;p!=badprimes[i];i++);
156 if (HECKE) euler_factor_hecke_bad(p,badprimetype[i],m,v);
157 else euler_factor_bad(p,badprimetype[i],m,v);}
158 QDpoly_intround(v);}
159
localinfo(llint p,int sp)160 static void localinfo(llint p,int sp)
161 {QDpoly v; printf("Euler factor at %lli is ",p);
162 euler_factor(p,sp,&v); QDpoly_intout(v); delQDpoly(&v);
163 if ((!HECKE) && (p<=3))
164 {printf("Tame cond exponent at %lli is %i\n",p,get_tame_conductor(p,sp));
165 printf("Wild cond exponent at %lli is %i\n",p,get_wild_conductor(p,sp));}
166 else printf("Conductor exponent at %lli is %i\n",p,get_conductor(p,sp));}
167
primetest(llint x)168 static int primetest(llint x)
169 {int i; llint p; if (x==1) return FALSE;
170 for(i=0;PRIMES[i]!=0;i++)
171 {p=PRIMES[i]; if (x==p) return 1; if ((x%p)==0) return 0;} return 1;}
172
localinfos(char * p,char * sp)173 void localinfos(char *p,char *sp)
174 {int k=0; llint prl=0,prh=0,pr; int ml=0,mh=0,m;
175 while ((p[k]!='-') && (p[k]!=0))
176 {prl*=10; ASSERT(ISA_NUMBER(p[k])); prl+=p[k]-'0'; k++;}
177 if (p[k]=='-') k++; else prh=prl;
178 while (p[k]!=0)
179 {prh*=10; ASSERT(ISA_NUMBER(p[k])); prh+=p[k]-'0'; k++;}
180 k=0; while ((sp[k]!='-') && (sp[k]!=0))
181 {ml*=10; ASSERT(ISA_NUMBER(sp[k])); ml+=sp[k]-'0'; k++;}
182 if (sp[k]=='-') k++; else mh=ml;
183 while (sp[k]!=0)
184 {mh*=10; ASSERT(ISA_NUMBER(sp[k])); mh+=sp[k]-'0'; k++;}
185 if ((ml<=0) || (mh>99)) errorit("Symmetric power range invalid");
186 for (m=ml;m<=mh;m++)
187 {if (HECKE) printf("Hecke "); printf("Symmetric power %i\n",m);
188 for (pr=prl;pr<=prh;pr++) if (primetest(pr)) localinfo(pr,m);}}
189