/* C program using pari.h -- many examples by Steven Miller, Princeton University Last modified April 2001 */ /* this is a program to investigate possible correction terms to the Sato-Tate distribution for curves with rank r over Q. If r > 1, we expect corrections (on average) to the a_p's of size r/sqrt(p) or maybe r/2sqrt(p). this program inputs an elliptic curve, uses PARI to calculate the coefficients, and then does a very crude test. Namely, it looks to see how many / what percent of the primes give a_p positive, negative and zero. If the Sato-Tate conj is true, we expect an equal percent of a_p positive and negative. */ /* needed input files: eprdat.fil --this is a 100Mb file of primes */ #include "pari.h" /* need this to use PARI functions */ #include #include #include "math.h" void main() { long numprimesless,j, x, p, count; long ltop1; /* these are used in freeing up memory after PARI acts */ long xloop; /* how many primes do */ long printevery; /* how often print results */ FILE *fpr, *fnr; long ap; long numdo = 99900; /* this controls how many a_p's we look at */ long apsigns[4]; double dapsigns[4]; GEN y, a, gap, F, E, b, tempp, temp; /* this is how you declare PARI variables you start with GEN for a general type you then list the variable names when you initialize the PARI variables later in the program, THEN you say what type it is */ numprimesless = 10000000; pari_init(200000000, 2); /* initializes PARI. tells how much memory to set aside, and how many primes to read in. I have it read in no primes, as I prefer to read it in from a C file as needed, saving memory. */ y = cgeti(32); /* initializing PARI variables. this makes it an integer, sets size */ a = cgeti(32); b = cgeti(32); gap = cgeti(32); tempp = cgeti(32); temp = cgetr(MEDDEFAULTPREC); /* initializing PARI variable, makes it REAL */ /* in PARI, one starts with a 5-vector [a1,a2,a3,a4,a6]. Running initell on this, one obtains a 13-vector which represents the elliptic curve with the above coefficients. In PARI, however, you need an extra place in the vectors for saving some system data. Thus, F is defined as a 6 vector, even though we will only use the first 5 components. */ F = cgetg(6, t_VEC); E = cgetg(14, t_VEC); for (x=1; x <= 5; x++) {F[x] = lgeti(MEDDEFAULTPREC);} for (x=1; x <= 13; x++) {E[x] = lgeti(MEDDEFAULTPREC);} /* Note above it is not enough to initialize F by running cgetg(6, t_VEC); we need to initialize each component. Hence the for loop to initialize components 1 to 5. */ gaffsg(0, (GEN) gap); /* assigns the long 0 to the PARI variable gap */ gaffsg(0, (GEN) temp); gaffsg(0, (GEN) tempp); /* lrgprime = 0;*/ xloop = numprimesless; fpr = fopen("eprdat.fil", "r"); /*opens for reading, this file has my list of primes */ xloop = 10; printevery = 1; fnr = fopen("stresult.fil","w"); /*opens for writing*/ fprintf(fnr, "Num\t ap<0\t ap>0\t ap=0\t ap<0\t\t ap>0\t\t ap=0\n"); for (x = 1; x <= 3; ++x) fscanf(fpr,"%ld ",&p); /*curves could be bad mod 2 or 3, starts us later */ /* now we choose the curve to work with */ gaffsg(0,(GEN) F[1]); gaffsg(0,(GEN) F[2]); gaffsg(0,(GEN) F[3]); gaffsg(-4,(GEN) F[4]); gaffsg(16,(GEN) F[5]); E = initell((GEN)F,MEDDEFAULTPREC); /* the above lines initializes E to be the elliptic curve with the above values. we use gaffsg because we are assigning LONG C values to PARI variables */ for (x = 0; x <= 3; ++x) apsigns[x] = 0; /* for (count = 1; count > 0; count++) { fscanf(fpr,"%ld ",&p); if ( (count % 10000) == 0) printf("(approx) numprime = %ld\t prime = %ld\n",count, p); } */ /* the above is to print out the list of primes numprime = 9990000 prime = 179235163 without overflowing */ count = 1; for (count = 1; count <= numdo; count = count + 1) { fscanf(fpr,"%ld ",&p); /* read in the next prime from our file */ ltop1 = avma; /* this is where the PARI memory pointer is currently at. We will run PARI to do some calculations, then reset the pointer here, clearing the calculations we've just done. If we don't need to save the result of calculations, this is the best way to free the memory. */ gaffsg(p, (GEN) tempp); /* this assigns the C LONG p to the PARI variable tempp */ gap = apell((GEN) E,(GEN)tempp); /* apell calculates a_p for the elliptic curve E. Note that we put (GEN) in front of PARI arguments in functions, and not that the two arguments of apell must be PARI variables. That's why we needed to convert the C LONG p to the PARI tempp */ ap = gtolong((GEN) gap); /* this converts the PARI gap to the C LONG ap. it is better to do all calculations in C. Thus, we enter PARI mode for the briefest amount of time and do the one calculation that we must. We then return to C. */ if (ap == 0) apsigns[0] = apsigns[0] + 1; if (ap > 0) apsigns[1] = apsigns[1] + 1; if (ap < 0) apsigns[2] = apsigns[2] + 1; ltop1 = avma; /* this resets the PARI memory pointer to where it was before we calculated a_p. */ if ( (count % 1000) == 0) { for (j = 0; j <= 3; ++j) dapsigns[j] = apsigns[j]*1.0 / count; printf("num ap < 0 is %ld\t, percent is %lf\n", apsigns[2], dapsigns[2]); printf("num ap > 0 is %ld\t, percent is %lf\n", apsigns[1], dapsigns[1]); printf("num ap = 0 is %ld\t, percent is %lf\n\n", apsigns[0], dapsigns[0]); fprintf(fnr,"%ld\t %ld\t %ld\t %ld\t %6.6f\t %6.6f\t %6.6f\n", count, apsigns[2], apsigns[1], apsigns[0], dapsigns[2], dapsigns[1], dapsigns[0]); } } printf("numdo = %ld\n", numdo); printf("num ap < 0 is %ld\t, percent is %lf\n", apsigns[2], 1.0*apsigns[2] / numdo); printf("num ap > 0 is %ld\t, percent is %lf\n", apsigns[1], 1.0*apsigns[1] / numdo); printf("num ap = 0 is %ld\t, percent is %lf\n", apsigns[0], 1.0*apsigns[0] / numdo); fclose(fnr); }