/* C program using pari.h -- many examples by Steven Miller, Princeton University Last modified April 2001 */ #include "pari.h" int main() { long x; double dx; GEN y, z, a, v, b, F, E,c, p, q, r, pr, Er, Fr, fp1, fp2, qr; /* Section 1: by Josko Plazonic initialization procedures GEN stands for generic pari type this section shows how to create GEN variables and goes thru the initialization process */ /*@Ccom take a stack of $10^6$ bytes, no prime table */ pari_init(1000000000, 2); y = cgeti(32); /* need to initialize the variables b4 using this makes them big integers */ z = cgeti(32); a = cgeti(32); b = cgeti(32); c = cgetr(MEDDEFAULTPREC); /* c IS A NEW PARI VARIABLE -- UNLIKE MOST OF THE OTHER PARI VARIABLES, c IS GOING TO BE A REAL NUMBER, HENCE THE R AT THE END INSTEAD OF I: REAL VS INTEGER */ gaffsg(800, y); /* assigns 800 to y */ printf("parametar: "); outbrute(y); /* prints u */ z = phi(y); /* makes z = phi(y), phi a pari function */ printf("\nresult: "); outbrute(z); printf("\n"); for (x=2; x <= 10; ++x) /* note x is a regular C variable */ { gaffsg(x, y); /* assigns x to y */ printf("phi["); outbrute(y); printf("] = "); z = phi(y); /* uses phi, a function from pari */ outbrute(z); printf("\n"); } printf("Made it thru Stage 1 Successfully\n"); /*!!!!!!!!!!!*/ /* */ /* Stage 2 */ /* */ /*!!!!!!!!!!!*/ v = cgetg(3, t_COMPLEX); /* v will be a complex number */ v[1] = lgetr(MEDDEFAULTPREC); /* initializing components */ v[2] = lgetr(MEDDEFAULTPREC); gaffsg(11,b); a = phi(b); /* makes a the pari number phi(11) */ (GEN) v[1] = a; /* assigns a to v[1] */ /* the program works without the (GEN) in front of v[1], but a warning message is generated. the (GEN) should be added to help C correctly cast objects */ printf("v[1] = "); outbrute((GEN) v[1]); printf("\n"); printf("Made it thru Stage 2 Successfully\n"); /*!!!!!!!!!!!*/ /* */ /* Stage 3 */ /* */ /*!!!!!!!!!!!*/ F = cgetg(6, t_VEC); /* makes F a vector */ F[1] = lgeti(MEDDEFAULTPREC); /* initializes components */ F[2] = lgeti(MEDDEFAULTPREC); F[3] = lgeti(MEDDEFAULTPREC); F[4] = lgeti(MEDDEFAULTPREC); F[5] = lgeti(MEDDEFAULTPREC); Fr = cgetg(6, t_VEC); /* makes Fr a vector */ Fr[1] = lgetr(MEDDEFAULTPREC); /* initializes components */ Fr[2] = lgetr(MEDDEFAULTPREC); Fr[3] = lgetr(MEDDEFAULTPREC); Fr[4] = lgetr(MEDDEFAULTPREC); Fr[5] = lgetr(MEDDEFAULTPREC); for (x = 1; x <= 5; ++x) /* assigns values to F's components note we can use the C long variable x to refer to entries of F */ { gaffsg(x*x,(GEN) F[x]); printf("F[%ld] = ",x); outbrute((GEN) F[x]); printf("\n"); } printf("Made it thru Stage 3 Successfully\n"); gaffsg(0,(GEN) F[1]); /* reassigns values of F */ gaffsg(0,(GEN) F[2]); gaffsg(0,(GEN) F[3]); gaffsg(0,(GEN) F[4]); gaffsg(1,(GEN) F[5]); printf("New Values of F[i]\n"); for (x = 1; x <= 5; ++x) /* prints new F values */ { printf("F[%ld] = ",x); outbrute((GEN) F[x]); printf("\n"); } printf("Starting Elliptic Curve Test Stage\n"); /* in pari, elliptic curves are of the form: y^2 + a1 xy + a3 y = x^3 + a2 x^2 + a4 x + a6 represent as a vector [a1,a2,a3,a4,a6] to use elliptic curve functions, need a related vector with 13 components, the first five as above. First we will manually input the data, then we will show how to use the pari built in initialization procedure, where it calculates the other 8 components from a1, ..., a6. */ E = cgetg(14, t_VEC); /* makes E a vector */ for (x = 1; x <= 13; x++) /* itializes components */ { E[x] = lgeti(MEDDEFAULTPREC); } Er = cgetg(14, t_VEC); /* makes Er a vector */ for (x = 1; x <= 13; x++) /* itializes components */ { Er[x] = lgetr(MEDDEFAULTPREC); } gaffsg(0,(GEN) E[1]); gaffsg(0,(GEN) E[2]); gaffsg(0,(GEN) E[3]); gaffsg(0,(GEN) E[4]); gaffsg(1,(GEN) E[5]); gaffsg(0,(GEN) E[6]); gaffsg(0,(GEN) E[7]); gaffsg(4,(GEN) E[8]); gaffsg(0,(GEN) E[9]); gaffsg(0,(GEN) E[10]); gaffsg(-846,(GEN) E[11]); gaffsg(-432,(GEN) E[12]); gaffsg(0,(GEN) E[13]); for (x = 1; x <= 13; ++x) /* prints new E values */ { printf("E[%ld] = ",x); outbrute((GEN) E[x]); printf("\n"); } gaffsg(97,y); a = apell(E,y); printf("a(97) for this curve is "); /*correct ans: 14 */ outbrute(a); printf("\n"); printf("First Section of Elliptic Curve Test Stage Successful\n"); /*second stage of elliptic curve test will try to use the initialization procedures instead of defining everything piecemeal. Below F is a vector with 5 components, and E is a vector with 13. MEDDEFAULTPREC is the precision of storage. The below works fine if all 13 entries are integral */ E = initell(F,MEDDEFAULTPREC); /*this initializes E to be a 13 component vector, and we can now use the pari elliptic curve functions on E */ printf("Printing Check for Second Stage\n"); for (x = 1; x <= 13; ++x) /* prints new E values */ { printf("E[%ld] = ",x); outbrute((GEN) E[x]); printf("\n"); } printf("Second Stage EC Test Successful\n"); /* third stage of EC test: will try to do initialization without having everything integral. one can check these values by running PARI in a dos window */ gaffsg(32,(GEN) F[1]); /* reassigns values of F */ gaffsg(0,(GEN) F[2]); gaffsg(0,(GEN) F[3]); /* c = dbltor(-1.0/4.0); gaffect((GEN) c,(GEN) Fr[4]); c = dbltor(1.0/16.0); gaffect((GEN) c,(GEN) Fr[5]); */ gaffsg(0, (GEN) F[4]); gaffsg(1, (GEN) F[5]); /* this gives us the curve y^2 + 32xy = x^3 + 1. PARI, unfortunately, wants its elliptic curves to have integer coefficients. I constructed this curve to have integer coefficients and a simple rational point, namely (1/4,1/8). */ E = initell(F,MEDDEFAULTPREC); /* again, initializing E */ printf("Printing Check for Third Stage\n"); for (x = 1; x <= 13; ++x) /* prints new E values */ { printf("E[%ld] = ",x); outbrute((GEN) E[x]); printf("\n"); } /*trying to use PARI height functions */ /* IT IS NON-TRIVIAL TO WORK WITH RATIONALS, NOT INTEGERS, AS THE POINTS ON THE ELLIPTIC CURVE. IF YOU STORE AS A DECIMAL, PARI BECOMES QUITE UPSET, AND WILL NOT DO THE COMPUTATION. YOU MUST STORE IT AS A FRACTION, IE, AS THE RATIO OF TWO INTEGERS. TO DO THIS, YOU NEED TO USE THE VARIABLE TYPE t_FRAC */ pr = cgetg(3, t_VEC); /* makes pr a vector of fractions */ pr[1] = cgetg(3, t_FRAC); pr[2] = cgetg(3, t_FRAC); /* It's annoying, but we must do it like this. pr is a vector with two components. each component is a fraction, and each fraction is two components. if the frac is negative, have the numerator less than zero. */ mael2(pr,1,1) = lgeti(MEDDEFAULTPREC); mael2(pr,1,2) = lgeti(MEDDEFAULTPREC); mael2(pr,2,1) = lgeti(MEDDEFAULTPREC); mael2(pr,2,2) = lgeti(MEDDEFAULTPREC); /* mael2 is PARI-speak for *M*ultidimensional *A*rray *EL*ement; the 2 means it is 2-dim. so mael2(pr,1,2) = pr[1][2], the first component of pr, and in the 1st component, take the second entry */ fp1 = cgetg(3, t_FRAC); /* makes p a vector */ fp1[1] = lgeti(MEDDEFAULTPREC); /* initializes components */ fp1[2] = lgeti(MEDDEFAULTPREC); fp2 = cgetg(3, t_FRAC); /* makes p a vector */ fp2[1] = lgeti(MEDDEFAULTPREC); /* initializes components */ fp2[2] = lgeti(MEDDEFAULTPREC); /* These are our two fractions. They'll be nice integers, which we will then put in our vector pr. */ gaffsg(1, (GEN) fp1[1]); gaffsg(4, (GEN) fp1[2]); gaffsg(1, (GEN) fp2[1]); gaffsg(8, (GEN) fp2[2]); /* as we are assigning integers, we use gaffsg; we are now storing the values of the integers in fp1[1], fp1[2], which is the x-coord, and in fp2[1], fp2[2], which is the y-coord. there's probably no need to go thru this intermediate step; we could've probably just put the values straight into the vector of fractions pr. */ printf("hello, good to Point I\n"); /*compiler debug test */ gaffect((GEN) fp1[1], (GEN) mael2(pr,1,1)); gaffect((GEN) fp1[2], (GEN) mael2(pr,1,2)); gaffect((GEN) fp2[1], (GEN) mael2(pr,2,1)); gaffect((GEN) fp2[2], (GEN) mael2(pr,2,2)); /* this assigns our two fracs, x and y, to our vector. remember what mael2 stands for */ printf("\n hello, Take II\n"); /*compiler debug test */ printf("pr[1][1] = "); outbrute(mael2(pr,1,1)); printf(" pr[2][1] = "); outbrute(mael2(pr,2,1)); printf("\n"); printf("pr[1][2] = "); outbrute(mael2(pr,1,2)); printf(" pr[2][2] = "); outbrute(mael2(pr,2,2)); printf("\n"); /* this prints out our vector, to make sure it stored ok */ printf("good - works thru here\n"); /*compiler debug test */ gaffect( (GEN) ghell2(E,pr,MEDDEFAULTPREC), (GEN) c); /* this is the key step -- we use gaffect as we are almost certainly *NOT* going to get an integer. ghell2 is one of the height functions. E is our elliptic curve. pr is a vector of a point that we KNOW is on the curve. it is ESSENTIAL that pr be a fraction and NOT a decimal. we store the result in c */ printf("Height of (1/4,1/8) is "); outbrute(c); printf("\n"); /* we now print out the height of the point (1/4,1/8) which is on the curve y^2 + 32xy = x^3 + 1. the correct answer is 1.6652763.... */ printf("Third Stage EC Test Successful\n"); printf("\nPress Return to Continue\n\n"); scanf("%c", &x); exit(0); }