1 #include "sympow.h"
2 #define DEBUG (FALSE || GLOBAL_DEBUG)
3
num_terms(int s,int sl)4 static llint num_terms(int s,int sl)
5 {QD x={0.0,0.0,0.0,0.0},R={1.0,0.0,0.0,0.0}; int i,which; double STEP,TAR,u;
6 which=whi[NUM_SUMS]; x[0]=DECAY[s][0];
7 TAR=Pow(10,(double) (-wprec[NUM_SUMS]+sl));
8 while (1)
9 {if (x[0]>=EXPAND0_LIM[which]) get_wt_large(which,x,R,1,NUM_SUMS);
10 if (R[0]>TAR) x[0]=x[0]+x[0]; else break;}
11 STEP=x[0]*0.25; R[0]=1.0; x[0]-=STEP; STEP*=0.5;
12 for (i=0;i<10;i++)
13 {if (x[0]>=EXPAND0_LIM[which]) get_wt_large(which,x,R,1,NUM_SUMS);
14 if (R[0]<TAR) x[0]-=STEP; else x[0]+=STEP; STEP*=0.5;}
15 u=x[0]/DECAY[s][0]; if (u>1e15) return (((llint) 1)<<50);
16 return (llint) (x[0]/DECAY[s][0]);}
17
prepare_decay_hecke(int s,int W,int sl)18 static llint prepare_decay_hecke(int s,int W,int sl)
19 {QD o; if (DEBUG) printf("pdh %i %i %i\n",s,W,sl);
20 if (!GLOBAL) return 0; QD_sqr(W,QD_twopi,o);
21 QD_div(W,o,COND[s],o); QD_sqrt(W,o,DECAY[s]);
22 if (VERBOSE>=2) {printf("DCY %i ",s); QD_output(W,wprec[NUM_SUMS],DECAY[s]);}
23 return num_terms(s,sl);}
24
prepare_decay(int s,int W,int sl)25 llint prepare_decay(int s,int W,int sl)
26 {QD o; if (DEBUG) printf("prepare_decay %i %i %i\n",s,W,sl);
27 if (!GLOBAL) return 0; QD_powi(W,QD_twopi,s,o);
28 if (s&1) QD_mul(W,o,QD_twopi,o); else QD_mul(W,o,QD_pi,o);
29 QD_div(W,o,COND[s],o); QD_sqrt(W,o,DECAY[s]);
30 if (VERBOSE>=2) {printf("DCY %i ",s); QD_output(W,wprec[NUM_SUMS],DECAY[s]);}
31 return num_terms(s,sl);}
32
prepare_sympow_hecke(int sp,int dv,int sl)33 static llint prepare_sympow_hecke(int sp,int dv,int sl)
34 {llint NT=-1; ASSERT(CM_CASE!=0);
35 if (DEBUG) printf("psh %i %i %i\n",sp,dv,sl);
36 if (sp&1)
37 {whi[NUM_SUMS]=wlo[NUM_SUMS]=WHICH;
38 load_files_hecke(WHICH,sp,1+sp/2,dv); WHICH++;}
39 else
40 {wlo[NUM_SUMS]=WHICH; WHICH++; whi[NUM_SUMS]=WHICH;
41 load_files_hecke(wlo[NUM_SUMS],sp,sp/2,dv);
42 load_files_hecke(whi[NUM_SUMS],sp,1+sp/2,dv); WHICH++;}
43 if (COND[sp][0]==0.0) compute_conductor_hecke(sp);
44 NT=prepare_decay_hecke(sp,w[NUM_SUMS],sl);
45 if (VERBOSE) printf("NT %id%i: %lli\n",sp,dv,NT);
46 NUM_SUMS++; return NT;}
47
prepare_sympow_check_hecke(int sp,int dv,int prec,int nw,int sl,int bk)48 static llint prepare_sympow_check_hecke
49 (int sp,int dv,int prec,int nw,int sl,int bk)
50 {int i,W; llint NT; QD wig;
51 if (DEBUG) printf("psch %i %i %i %i %i %i\n",sp,dv,prec,nw,sl,bk);
52 wprec[NUM_SUMS]=prec; W=w[NUM_SUMS]=1+((prec-1)>>4); NUM_WIGS[NUM_SUMS]=nw;
53 QD_copy(W,QD_one,WIGGLE[NUM_SUMS]); QD_copy(W,QD_one,WIGSQI[NUM_SUMS]);
54 BLOCH_KATO[NUM_SUMS]=bk; NT=prepare_sympow_hecke(sp,dv,sl);
55 if (!GLOBAL) return NT;
56 for (i=0;i<nw;i++)
57 {w[NUM_SUMS]=W; wprec[NUM_SUMS]=prec; QD_copy(W,QD_one,wig); wig[0]=0.5;
58 QD_powi(1,wig,3+i,wig); QD_add(W,QD_one,wig,wig); SLOPPY[NUM_SUMS]=sl;
59 if (VERBOSE>=2) {printf("wiggle %i ",i); QD_output(W,16*W,wig);}
60 whi[NUM_SUMS]=whi[NUM_SUMS-1]; wlo[NUM_SUMS]=wlo[NUM_SUMS-1];
61 BLOCH_KATO[NUM_SUMS]=BLOCH_KATO[NUM_SUMS-1];
62 QD_copy(W,wig,WIGGLE[NUM_SUMS]); QD_div(W,QD_one,wig,WIGSQI[NUM_SUMS]);
63 QD_sqr(W,WIGSQI[NUM_SUMS],WIGSQI[NUM_SUMS]); NUM_WIGS[NUM_SUMS]=nw;
64 if (VERBOSE>=2) {printf("wigsqi %i ",i); QD_output(W,16*W,WIGSQI[NUM_SUMS]);}
65 NUM_SUMS++;} return NT;}
66
prepare_sympow(int sp,int dv,int sl)67 static llint prepare_sympow(int sp,int dv,int sl)
68 {llint NT=-1; if (DEBUG) printf("prepare_sympow %i %i %i\n",sp,dv,sl);
69 if (sp&1)
70 {whi[NUM_SUMS]=wlo[NUM_SUMS]=WHICH;
71 load_files(WHICH,sp,(sp+1)>>1,dv); WHICH++;}
72 else
73 {wlo[NUM_SUMS]=WHICH; WHICH++; whi[NUM_SUMS]=WHICH;
74 load_files(wlo[NUM_SUMS],sp,sp>>1,dv);
75 load_files(whi[NUM_SUMS],sp,1+(sp>>1),dv); WHICH++;}
76 if (COND[sp][0]==0.0) compute_conductor(sp);
77 NT=prepare_decay(sp,w[NUM_SUMS],sl);
78 if (VERBOSE) printf("NT %id%i: %lli\n",sp,dv,NT);
79 NUM_SUMS++; return NT;}
80
prepare_sympow_check(int sp,int dv,int prec,int nw,int sl,int bk)81 static llint prepare_sympow_check(int sp,int dv,int prec,int nw,int sl,int bk)
82 {int i,W; llint NT; QD wig;
83 if (DEBUG) printf("psc %i %i %i %i %i %i\n",sp,dv,prec,nw,sl,bk);
84 if (HECKE) return prepare_sympow_check_hecke(sp,dv,prec,nw,sl,bk);
85 wprec[NUM_SUMS]=prec; W=w[NUM_SUMS]=1+((prec-1)>>4); NUM_WIGS[NUM_SUMS]=nw;
86 QD_copy(W,QD_one,WIGGLE[NUM_SUMS]); QD_copy(W,QD_one,WIGSQI[NUM_SUMS]);
87 BLOCH_KATO[NUM_SUMS]=bk;
88 NT=prepare_sympow(sp,dv,sl); if (!GLOBAL) return NT;
89 for (i=0;i<nw;i++)
90 {w[NUM_SUMS]=W; wprec[NUM_SUMS]=prec; QD_copy(W,QD_one,wig); wig[0]=0.5;
91 QD_powi(1,wig,3+i,wig); QD_add(W,QD_one,wig,wig); SLOPPY[NUM_SUMS]=sl;
92 if (VERBOSE>=2) {printf("wiggle %i ",i); QD_output(W,16*W,wig);}
93 whi[NUM_SUMS]=whi[NUM_SUMS-1]; wlo[NUM_SUMS]=wlo[NUM_SUMS-1];
94 BLOCH_KATO[NUM_SUMS]=BLOCH_KATO[NUM_SUMS-1];
95 QD_copy(W,wig,WIGGLE[NUM_SUMS]); QD_div(W,QD_one,wig,WIGSQI[NUM_SUMS]);
96 QD_sqr(W,WIGSQI[NUM_SUMS],WIGSQI[NUM_SUMS]); NUM_WIGS[NUM_SUMS]=nw;
97 if (VERBOSE>=2) {printf("wigsqi %i ",i); QD_output(W,16*W,WIGSQI[NUM_SUMS]);}
98 NUM_SUMS++;} return NT;}
99
process_string(char * IN,llint UB)100 llint process_string(char *IN,llint UB)
101 {llint NT=0,nt,CP; char LAST[256],*NEXT,*CURR,*c,*n,COMMA[2]=",\0";
102 int i,prec,sp,lim=0,nwigs,dvs=0,sl,bk,DO;
103 c=CURR=malloc(256); n=NEXT=malloc(256); strcpy(CURR,IN);
104 do
105 {strcpy(LAST,CURR); ASSERT(ISA_NUMBER(CURR[0]));
106 if (!ISA_NUMBER(CURR[1])) sp=(int) (CURR[0]-'0');
107 else {sp=(int) ((CURR[0]-'0')*10+(CURR[1]-'0')); CURR++;}
108 nwigs=1; CURR++; CP=NT; sl=0; bk=0; DO=FALSE;
109 if (CURR[0]=='b') {bk=1; CURR++;}
110 if (CURR[0]=='w') {nwigs=CURR[1]-'0'; CURR+=2;}
111 if (CURR[0]=='s') {sl=CURR[1]-'0'; CURR+=2;}
112 if (CURR[0]=='P') {DO=TRUE;} else ASSERT(CURR[0]=='p');
113 if (CURR[1]=='z')
114 {if (!RERUN) {prec=6; ZEROCHECK=1;} else {prec=12; ZEROCHECK=0;} CURR+=2;}
115 else
116 {if (ISA_NUMBER(CURR[2])) {prec=10*(CURR[1]-'0')+(CURR[2]-'0'); CURR+=3;}
117 else {prec=(CURR[1]-'0'); CURR+=2;}}
118 if ((sp&1)==0)
119 {if ((CURR[0]=='d') || (CURR[0]=='D')) errorit("Derivative with even power");
120 dvs=1; nt=prepare_sympow_check(sp,0,prec,nwigs,sl,bk); if (nt>NT) NT=nt;}
121 else
122 {if (CURR[0]=='D')
123 {lim=CURR[1]-'0'; CURR+=2; dvs=1;
124 nt=prepare_sympow_check(sp,lim,prec,nwigs,sl,bk); if (nt>NT) NT=nt;}
125 else
126 {ASSERT(CURR[0]=='d'); lim=CURR[1]-'0'; CURR+=2; dvs=(lim+1);
127 for (i=0;i<=lim;i++)
128 {nt=prepare_sympow_check(sp,i,prec,nwigs,sl,bk); if (nt>NT) NT=nt;}}}
129 if ((NT>UB) && (DO)) UB=NT;
130 if ((NT>UB) || (((sp&1)==0) && NO_QT))
131 {for (i=0;((LAST[i]!=0) && (LAST[i]!=','));i++); LAST[i]=0;
132 if (VERBOSE)
133 {if (NT>UB) printf("Ignoring %s due to bound\n",LAST);
134 else printf("Ignoring %s since curve is non qtwist-minimal\n",LAST);}
135 NT=CP; NUM_SUMS-=dvs*(1+nwigs);}
136 else CP=NT;
137 NEXT=strstr(CURR,COMMA); CURR=NEXT+1; } while(NEXT);
138 free(c); free(n); return(NT);}
139
preparation(int K,char * IN,llint UB)140 llint preparation(int K,char *IN,llint UB)
141 {int i; llint NT;
142 POWSER=malloc(K*sizeof(QD**)); TABLE=malloc(K*sizeof(QD**));
143 SYMPOW=malloc(K*sizeof(int)); evalpt=malloc(K*sizeof(int));
144 BLOCH_KATO=malloc(K*sizeof(int)); MESH_COUNT=malloc(K*sizeof(int));
145 whi=malloc(K*sizeof(int)); wlo=malloc(K*sizeof(int));
146 derivative=malloc(K*sizeof(int)); HALF_ZERO=malloc(K*sizeof(QD));
147 WIGGLE=malloc(K*sizeof(QD)); WIGSQI=malloc(K*sizeof(QD));
148 EXPAND0_LIM=malloc(K*sizeof(QD)); STEP_SIZE=malloc(K*sizeof(QD));
149 TOO_BIG=malloc(K*sizeof(QD)); NUM_LOGS=malloc(K*sizeof(QD));
150 TACKON=malloc(K*sizeof(int)); TACKS=malloc(K*sizeof(QD*));
151 w=malloc(K*sizeof(int)); wprec=malloc(K*sizeof(int));
152 DECAY=malloc(K*sizeof(QD)); NUM_WIGS=malloc(K*sizeof(QD));
153 for (i=0;i<K;i++) QD_copy(wmax,QD_zero,COND[i]); WHICH=0; NUM_SUMS=0;
154 NT=process_string(IN,UB);
155 if (MODDEG) {if (MD_SPEED>0.0) NT=(llint) ((double) NT/MD_SPEED);
156 else NT=postpare_moddeg();}
157 printf("Maximal number of terms is %lli\n",NT);
158 if (NT<MIN_TERMS) NT=MIN_TERMS; return(NT);}
159