1 /*
2
3 Copyright (C) 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 Patches borrowed from SageMath.
22
23 */
24
25
26 //---------------------------------------------------
27 //
28 // Command line interface to the L function package
29 // By Michael Rubinstein
30 //
31 //---------------------------------------------------
32
33 #include "Lcommandline.h"
34 #include "cmdline.h"
35
main(int argc,char * argv[])36 int main (int argc, char *argv[])
37 {
38
39
40 int i;
41 Long n;
42 bool do_zeros=false,do_values=false,file_to_open=false,do_twists=false,file_of_s_values=false;
43 bool file2_to_open=false,do_print=false,do_interpolate=false;
44 char data_filename[1000]; //filename of file containing data for L-function.
45 char data_filename2[1000]; //filename of file containing data for L-function.
46 int number_coeff_print=30;
47 bool do_hardy=false;
48 bool do_explicit=false; //whether to test the explicit formula when looking for zeros
49
50 //used for elliptic curves
51 // y^2 + a1 xy + a3 y = x^3 + a2 x^2 + a4 x + a6
52 bool do_elliptic_curve=false;
53 char a1[200];
54 char a2[200];
55 char a3[200];
56 char a4[200];
57 char a6[200];
58 int rank=-1;
59
60 bool do_tau=false;
61 bool check_rank=false;
62 bool do_compute_rank=false;
63 bool limit_the_rank=false;
64 bool do_only_even_twists=false;
65
66 char s_file_name[1000]; //file of s values
67
68 //int derivative=0; //which derivative of L to compute (0 is just the L-value)
69
70 Double x=0,y=0,step_size=1025;
71 Double x2=0,y2=0;
72 Long count=0;
73
74 char twist_type; //type of twist to do: quadratic, primitive, all, nth
75 Long D1,D2; // for twists
76
77 int print_character=0; //1 is all chi(n) , 2 is just chi(-1)
78
79
80 bool do_exit=false;
81 gengetopt_args_info args_info;
82 /* call the cmdline parser */
83 if (cmdline_parser (argc, argv, &args_info) != 0){
84 exit(1) ;
85 }
86
87
88 if(args_info.zeros_given){ /* Whether zeros was given. */
89 do_zeros=true;
90 count=args_info.zeros_arg;
91 }
92 if(args_info.zeros_interval_given){ /* Whether zeros-interval was given. */
93 do_zeros=true;
94 if(!args_info.x_given||!args_info.y_given){
95 cout << "x and y options must be used" << endl;
96 exit(1);
97 }
98 if(!args_info.stepsize_given){
99 cout << "stepsize option must be used" << endl;
100 exit(1);
101 }
102 #ifdef USE_LONG_DOUBLE
103 sscanf(args_info.x_arg,"%Lg",&x);
104 sscanf(args_info.y_arg,"%Lg",&y);
105 #elif USE_MPFR
106 x=args_info.x_arg;
107 y=args_info.y_arg;
108 #else
109 sscanf(args_info.x_arg,"%lg",&x);
110 sscanf(args_info.y_arg,"%lg",&y);
111 #endif
112 if(y<x){
113 cout << "x should be < than y" << endl;
114 exit(1);
115 }
116 #ifdef USE_LONG_DOUBLE
117 sscanf(args_info.stepsize_arg,"%Lg",&step_size);
118 #elif USE_MPFR
119 step_size=args_info.stepsize_arg;
120 #else
121 sscanf(args_info.stepsize_arg,"%lg",&step_size);
122 #endif
123 //step_size=atof(args_info.stepsize_arg);
124 }
125 if(args_info.explicit_given) do_explicit=true; /* Whether to check the explicit formula */
126 if(args_info.value_given){ /* Whether value was given. */
127 do_values=true;
128 if(!args_info.x_given||!args_info.y_given){
129 cout << "x and y options must be used" << endl;
130 exit(1);
131 }
132 #ifdef USE_LONG_DOUBLE
133 sscanf(args_info.x_arg,"%Lg",&x);
134 sscanf(args_info.y_arg,"%Lg",&y);
135 #elif USE_MPFR
136 x=args_info.x_arg;
137 y=args_info.y_arg;
138 #else
139 sscanf(args_info.x_arg,"%lg",&x);
140 sscanf(args_info.y_arg,"%lg",&y);
141 #endif
142 }
143 if(args_info.value_file_given){ /* Whether value-file was given. */
144 do_values=true;
145 file_of_s_values=true;
146 strcpy(s_file_name,args_info.value_file_arg);
147 }
148 if(args_info.value_line_segment_given){/* Whether value-line-segment was given. */
149 do_values=true;
150 if(!args_info.x_given||!args_info.y_given||
151 !args_info.X_given||!args_info.Y_given){
152 cout << "x,y,X,Y options must be used" << endl;
153 do_exit=true;
154 }
155 if(!args_info.number_samples_given){
156 cout << "--number-samples option must be used" << endl;
157 do_exit=true;
158 }
159 if(do_exit) exit(1);
160 #ifdef USE_LONG_DOUBLE
161 sscanf(args_info.x_arg,"%Lg",&x);
162 sscanf(args_info.y_arg,"%Lg",&y);
163 sscanf(args_info.X_arg,"%Lg",&x2);
164 sscanf(args_info.Y_arg,"%Lg",&y2);
165 #elif USE_MPFR
166 x=args_info.x_arg;
167 y=args_info.y_arg;
168 x2=args_info.X_arg;
169 y2=args_info.Y_arg;
170 #else
171 sscanf(args_info.x_arg,"%lg",&x);
172 sscanf(args_info.y_arg,"%lg",&y);
173 sscanf(args_info.X_arg,"%lg",&x2);
174 sscanf(args_info.Y_arg,"%lg",&y2);
175 #endif
176 }
177 if(args_info.number_samples_given){ /* Whether derivative was given. */
178 count=args_info.number_samples_arg;
179 }
180 if(args_info.hardy_given){ /* Whether Hardy Z option was given. */
181 do_hardy=true;
182 }
183
184 if(args_info.elliptic_curve_given){ /* Whether elliptic-curve was given. */
185 do_elliptic_curve=true;
186 if(!args_info.a1_given){ cout << "must specify a1 with --a1" << endl; do_exit=true;}
187 if(!args_info.a2_given){ cout << "must specify a2 with --a2" << endl; do_exit=true;}
188 if(!args_info.a3_given){ cout << "must specify a3 with --a3" << endl; do_exit=true;}
189 if(!args_info.a4_given){ cout << "must specify a4 with --a4" << endl; do_exit=true;}
190 if(!args_info.a6_given){ cout << "must specify a6 with --a6" << endl; do_exit=true;}
191 if(do_exit) exit(1);
192 strcpy(a1,args_info.a1_arg);
193 strcpy(a2,args_info.a2_arg);
194 strcpy(a3,args_info.a3_arg);
195 strcpy(a4,args_info.a4_arg);
196 strcpy(a6,args_info.a6_arg);
197 }
198 if(args_info.file_input_given){ /* Whether file-input was given. */
199 file_to_open=true;
200 strcpy(data_filename,args_info.file_input_arg);
201 }
202 if(args_info.url_given){ /* Whether file-input was given. */
203 char str1[]="wget -O temporary_url_file_lcalc ";
204 if(system(strcat(str1,args_info.url_arg))!=0){
205 cout << "Error retrieving file: " << args_info.url_arg << endl;
206 do_exit=true;
207 }
208 if(do_exit) exit(1);
209 file_to_open=true;
210 strcpy(data_filename,"temporary_url_file_lcalc");
211 }
212 if(args_info.interpolate_given){ /* Whether interpolate was given. */
213 do_interpolate=true;
214 file2_to_open=true;
215 strcpy(data_filename2,args_info.interpolate_arg);
216 }
217 if(args_info.output_character_given){ /* Whether output-character was given. */
218 print_character = atoi(args_info.output_character_arg);
219 }
220 if(args_info.output_data_given){ /* Whether output-data was given. */
221 do_print=true;
222 number_coeff_print=args_info.output_data_arg;
223 }
224 if(args_info.precision_given){ /* Whether precision was given. */
225 DIGITS=args_info.precision_arg;
226 DIGITS3=DIGITS-DIGITS2;
227
228 }
229 if(args_info.sacrifice_given){ /* Whether sacrifice was given. */
230 DIGITS2=args_info.sacrifice_arg;
231 DIGITS3=DIGITS-DIGITS2;
232
233 }
234 if(args_info.rank_compute_given){ /* Whether rank-compute was given. */
235 do_compute_rank=true;
236 }
237 if(args_info.rank_verify_given){ /* Whether rank-verify was given. */
238 check_rank=true;
239 rank=args_info.rank_verify_arg;
240 }
241 if(args_info.rank_limit_given){ /* Whether rank-limit was given. */
242 limit_the_rank=true;
243 rank=args_info.rank_limit_arg;
244 }
245 if(args_info.tau_given){ /* Whether tau was given. */
246 do_tau=true;
247 }
248 if(args_info.twist_quadratic_given){ /* Whether twist-quadratic was given. */
249 do_twists=true;
250 twist_type='q';
251 }
252 if(args_info.twist_quadratic_even_given){ /* Whether twist-quadratic was given. */
253 do_twists=true;
254 do_only_even_twists=true;
255 twist_type='q';
256 }
257 if(args_info.twist_primitive_given){ /* Whether twist-primitive was given. */
258 do_twists=true;
259 twist_type='p';
260 }
261 if(args_info.twist_all_given){ /* Whether twist-all was given. */
262 do_twists=true;
263 twist_type='A';
264 }
265 if(args_info.twist_all_no_conj_pairs_given){ /* Whether twist-all-no-conj-pairs was given. */
266 do_twists=true;
267 twist_type='a';
268 }
269 if(args_info.twist_complex_no_conj_pairs_given){ /* Whether twist-complex-no-conj-pairs was given. */
270 do_twists=true;
271 twist_type='c';
272 }
273 if(args_info.twist_generic_given){ /* Whether twist-complex-no-conj-pairs was given. */
274 do_twists=true;
275 twist_type='g';
276 }
277 if(args_info.degree_given){ /* Whether degree was given. */
278 do_twists=true;
279 twist_type='n';
280 n=args_info.degree_arg;
281 }
282 if(do_twists){
283 if(!args_info.start_given){
284 cout << "must specify starting discriminant/conductor (--finish or -f)" << endl;
285 do_exit=true;
286 }
287 if(!args_info.finish_given){
288 cout << "must specify finishing discriminant/conductor (with --start or -s)" << endl;
289 do_exit=true;
290 }
291 if(do_exit) exit(1);
292 D1=strtoll(args_info.start_arg,0,10); D2=strtoll(args_info.finish_arg,0,10);
293 }
294
295 #ifdef _OPENMP
296 if(args_info.openmp_given) /* whether a number of threads was specified */
297 {
298 /* get the total number of CPUs/cores available for OpenMP */
299 int NCPU = omp_get_num_procs();
300 if (args_info.openmp_arg>NCPU){
301 cout << "ERROR: You only have " << NCPU << " processors available." << endl;
302 cout << "Specify fewer processors" << endl;
303 exit(1);
304 }
305 omp_set_num_threads(args_info.openmp_arg);
306 }
307 else //default is 1
308 omp_set_num_threads(1);
309 #endif
310 #ifndef _OPENMP
311 if(args_info.openmp_given){ /* whether a number of threads was specified */
312 cout << endl;
313 cout << "ERROR: lcalc has not been compiled with openmp enabled." << endl;
314 cout << "You need to uncomment the line: " << endl;
315 cout << "#OPENMP_FLAG = -fopenmp" << endl;
316 cout << "in the Makefile, and recompile by typing:"<< endl;
317 cout << endl;
318 cout << "make clean" << endl;
319 cout << "make" << endl;
320 cout << endl;
321 exit(1);
322 }
323 #endif
324
325 //placed here rather than at top since
326 //precision depends on user input.
327
328 initialize_globals();
329
330
331 //the L-functions used throughout are global and
332 //need to be re-initialized now that we have initialized
333 //Pi in initialize_globals.
334 initialize_commandline_globals();
335 if(args_info.derivative_given){ /* Whether derivative was given. */
336 global_derivative = args_info.derivative_arg;
337 //DIGITS3=(int)(DIGITS3/pow(2.,global_derivative)); // now taken care of in Lvalue.h
338 }
339
340 if(do_zeros){
341 input_mean_spacing_given = (6.28/log(count*1.+100))*2./DIGITS; //only used
342 //for the zeta funciton, and then in the band limited interpolation scheme. Is a
343 //rough estimate for the average distance between sample points.
344 //Dividing by DIGITS/2 more than accounts for searching for and then zooming in
345 //on zeros.
346 //cout << " input_mean_spacing_given " << input_mean_spacing_given << endl;
347
348 }
349
350 if(args_info.use_dirichlet_series_given){ /* Whether to compute just using the Dirichlet series */
351 only_use_dirichlet_series=true;
352 N_use_dirichlet_series=args_info.use_dirichlet_series_arg; //how many terms to use
353 }
354
355
356 if(args_info.verbose_given){ /* Whether vebosity level was specified. */
357 my_verbose=args_info.verbose_arg;
358 }
359
360
361 A=1./(16*Pi*Pi*2.3*DIGITS); //controls the 'support' of g(w)
362
363 if(file2_to_open){
364 initialize_new_L(data_filename2);
365 switch(current_L_type)
366 {
367 case 1:
368 int_L2=int_L;
369 break;
370 case 2:
371 Double_L2=Double_L;
372 break;
373 case 3:
374 Complex_L2=Complex_L;
375 break;
376 }
377
378 }
379
380 if(file_to_open){
381 initialize_new_L(data_filename);
382 if(args_info.url_given){ /* Whether file-input was given. */
383 system("rm temporary_url_file_lcalc");
384 }
385 }
386 else current_L_type=1; //default is zeta, and type 1, i.e. int, is simplest
387
388 if(do_interpolate) //compute the zeros on the critical line
389 //of the L-series that interpolate between f and f2. Must be used in combination with
390 //with the option: -zi x y increment
391 //By interpolate I mean take t*(basic data of L) and (1-t)*(basic data of L2).
392 {
393 L_interpolate(x,y,step_size);
394 return 0;
395 }
396
397
398 Double *coeff;
399
400 if(do_elliptic_curve||do_tau){ //determine how many dirichlet coefficients are needed
401 //and initialize the curve
402
403
404 int N_terms; //number of dirichlet coefficients
405 Double C = DIGITS2*log(10.);
406 Double T;
407 Double tmp_Q; //the maximum normalized conductor (including twists)
408
409 Double t2;
410
411
412 if(do_elliptic_curve){
413 #ifdef INCLUDE_PARI
414
415 t2=.5; //t2=.5 because of the GAMMA(s+1/2)
416
417 pari_init_opts(16000000, 2, INIT_DFTm);
418
419 coeff = new Double[3];
420 //compute the conductor which is copied to coeff[1]
421 data_E(a1,a2,a3,a4,a6,2,coeff);
422
423 tmp_Q=sqrt(coeff[1])/(2*Pi);
424 delete [] coeff;
425 #else
426 cout << "You need to uncomment the line: PARI_DEFINE = -DINCLUDE_PARI" <<endl;
427 cout << "in the Makefile and do: 'make clean', then 'make' if you wish to use" <<endl;
428 cout << "elliptic curve L-functions. Requires that you already have pari installed" <<endl;
429 cout << "on your machine." <<endl;
430 exit(1);
431 #endif //ifdef INCLUDE_PARI
432 }
433
434
435
436 if(do_tau){
437 t2=5.5; //t2=5.5 because of the GAMMA(s+11/2)
438 tmp_Q=1./(2*Pi);
439 }
440
441
442
443 if(do_twists) tmp_Q=tmp_Q*max(abs(1.*D1),abs(1.*D2)); // take into account twists
444
445
446 if(do_zeros)
447 T = max(abs(x),abs(y));
448 else
449 T = max(abs(y),abs(y2));
450
451 if(my_verbose>1) cout << "T = " << T << endl;
452
453 if(count>0&&!args_info.number_samples_given)
454 T=T+(count+100)*Pi/log(T+3);
455
456 if(my_verbose>1) cout << "T = " << T << endl;
457
458 //based on the number of terms used in the gamma_sum routine. possible updates to gamma_sum condition should
459 //thus be reflected here
460 Complex delta=int_L.find_delta(1+I*T+t2,1);
461 N_terms = Int(2.3 * DIGITS*tmp_Q/real(delta));
462 do{
463 N_terms=(int)(N_terms*1.3);
464 if(my_verbose>4) cout << "N_terms to precompute = " << N_terms << endl;
465 }while(N_terms*real(delta)/tmp_Q-log((1.+t2)*N_terms)<2.3*DIGITS);
466 N_terms=(int)(N_terms*1.3+40);
467 if(my_verbose>4) cout << "N_terms to precompute = " << N_terms << endl;
468
469
470
471 #ifdef INCLUDE_PARI
472 if(do_elliptic_curve){
473 // Reallocate PARI stack
474 paristack_setsize((size_t)N_terms*16 + 1000000, 0); //XXXXXXXXX this should depend on whether we're double or long double or mpfr double
475
476 if (my_verbose>0) cout << "Will precompute " << N_terms << " elliptic L-function dirichlet coefficients..." << endl;
477 initialize_new_L(a1,a2,a3,a4,a6,N_terms);
478 }
479 #endif //ifdef INCLUDE_PARI
480
481 if(do_tau){
482
483 tmp_Q=1./(2*Pi);
484
485 coeff = new Double[N_terms+1];
486 ramanujan_tau(coeff,N_terms);
487
488 Double *g;
489 Complex *l;
490 Complex *p;
491 Complex *r;
492
493 g=new Double[2];
494 l=new Complex[2];
495 g[1]=1.;
496 l[1]=5.5;
497
498 p = new Complex[1];
499 r = new Complex[1];
500
501
502
503 if (my_verbose>0) cout << "Will precompute " << N_terms << " Ramanujan tau(n) dirichlet coefficients..." << endl;
504 current_L_type=2; //the normalized dirichlet coeffs are real
505 Double_L=L_function<Double>("Ramanujan Tau",2,N_terms,coeff,0,tmp_Q,1,1,g,l,0,p,r);
506
507 delete [] g;
508 delete [] l;
509 delete [] p;
510 delete [] r;
511 delete [] coeff;
512
513 }
514
515
516
517 }
518
519
520 if(do_print) print_L(number_coeff_print);
521
522 if(check_rank&&rank>=0) verify_rank(rank);
523
524 if(do_compute_rank&&!do_twists) compute_rank();
525
526 if(do_values){
527 if(do_twists){
528 switch(twist_type){
529 case 'q':
530 if(do_compute_rank)
531 quadratic_twists(D1,D2,x,y,0,0,"values and ranks",do_only_even_twists,do_explicit);
532 else if(limit_the_rank)
533 quadratic_twists(D1,D2,x,y,0,0,"values and ranks",do_only_even_twists,rank,do_explicit);
534 else
535 quadratic_twists(D1,D2,x,y,0,0,"values",do_only_even_twists,do_explicit);
536 break;
537 case 'p':
538 if(do_compute_rank)
539 all_twists(D1,D2,x,y,0,0,"values and ranks",0,print_character,do_explicit);
540 else
541 all_twists(D1,D2,x,y,0,0,"values",0,print_character,do_explicit);
542 break;
543 case 'a':
544 if(do_compute_rank)
545 all_twists(D1,D2,x,y,0,0,"values and ranks",1,print_character,do_explicit);
546 else
547 all_twists(D1,D2,x,y,0,0,"values",1,print_character,do_explicit);
548 break;
549 case 'A':
550 if(do_compute_rank)
551 all_twists(D1,D2,x,y,0,0,"values and ranks",2,print_character,do_explicit);
552 else
553 all_twists(D1,D2,x,y,0,0,"values",2,print_character,do_explicit);
554 break;
555 case 'g':
556 if(do_compute_rank)
557 all_twists(D1,D2,x,y,0,0,"values and ranks",-1,print_character,do_explicit);
558 else
559 all_twists(D1,D2,x,y,0,0,"values",-1,print_character,do_explicit);
560 break;
561 case 'c':
562 if(do_compute_rank)
563 all_twists(D1,D2,x,y,0,0,"values and ranks",-2,print_character,do_explicit);
564 else
565 all_twists(D1,D2,x,y,0,0,"values",-2,print_character,do_explicit);
566 break;
567 }
568 }
569 else if(file_of_s_values){
570 if(!do_hardy) compute_values(x,y,"pure",s_file_name);
571 else compute_values(x,y,"rotated pure",s_file_name);
572 }
573 else{
574 input_mean_spacing_given=(y2-y)/count; //used in Riemann Siegel band limited routine
575 //cout << " input_mean_spacing_given " << input_mean_spacing_given << endl;
576 //for zeta, if applicable.
577 if(!do_hardy) compute_values(x,y,"pure","",x2,y2,count);
578 else compute_values(x,y,"rotated pure","",x2,y2,count);
579 }
580
581 }
582 else if(do_zeros){
583 if(do_twists){
584 switch(twist_type){
585 case 'q':
586 if(limit_the_rank)
587 quadratic_twists(D1,D2,x,y,count,step_size,"zeros and ranks",do_only_even_twists,do_explicit,rank);
588 else
589 quadratic_twists(D1,D2,x,y,count,step_size,"zeros",do_only_even_twists,do_explicit);
590 break;
591 case 'p':
592 all_twists(D1,D2,x,y,count,step_size,"zeros",0,print_character,do_explicit);
593 break;
594 case 'a':
595 all_twists(D1,D2,x,y,count,step_size,"zeros",1,print_character,do_explicit);
596 break;
597 case 'A':
598 all_twists(D1,D2,x,y,count,step_size,"zeros",2,print_character,do_explicit);
599 break;
600 case 'g':
601 all_twists(D1,D2,x,y,count,step_size,"zeros",-1,print_character,do_explicit);
602 break;
603 case 'c':
604 all_twists(D1,D2,x,y,count,step_size,"zeros",-2,print_character,do_explicit);
605 break;
606
607 }
608 }
609 else{
610 compute_zeros(x,y,step_size,count,rank,do_explicit);
611 }
612 }
613 else{ //else allow for printing of characters. "print character" is a dummy char str,
614 //so all_twists will print out the character without doing zeros or values
615 if(do_twists){
616 switch(twist_type){
617 case 'p':
618 all_twists(D1,D2,x,y,count,step_size,"print character",0,print_character,do_explicit);
619 break;
620 case 'a':
621 all_twists(D1,D2,x,y,count,step_size,"print character",1,print_character,do_explicit);
622 break;
623 case 'A':
624 all_twists(D1,D2,x,y,count,step_size,"print character",2,print_character,do_explicit);
625 break;
626 case 'g':
627 all_twists(D1,D2,x,y,count,step_size,"print character",-1,print_character,do_explicit);
628 break;
629 case 'c':
630 all_twists(D1,D2,x,y,count,step_size,"print character",-2,print_character,do_explicit);
631 break;
632
633 }
634 }
635
636 }
637
638 delete_globals();
639 cmdline_parser_free (&args_info); /* release allocated memory */
640
641 return 0;
642 }
643