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