The triples project
On this page we develope code to work with triples as discussed in our meetings.
PART I: low level functions:
Why don't we make a library of low level functions for attacking the problem of finding curves here. Here is how I want to set things up:
/* Take a polynomial and make the corresponding integer */ convert_to_int(f) = { local(d, u); d = poldegree(f, t); u = 0; for(i = 0, d, u = u + polcoeff(f, i, t)*2^i); return(u); }
/* Take an integer and make the corresponding polynomial */ convert_to_pol(u) = { local(i, f); f = 0; i = 0; while(u, f = f + (u % 2)*t^i; u = u >> 1; i = i + 1; ); return(f); }
/* Given two integers find the integer corresponding to the sum of the corresponding polynomials. Note that addition in F_2 is just xor on {0, 1}. Thus this becomes quite easy to do. */ int_sum(u1, u2) = { return(bitxor(u1, u2)); }
/* Given an integer find the integer corresponding to the fifth power of the corresponding polynomial. Here we will use that (a_0 + a_1 t + a_2 t^2 + ...)^5 = (a_0 + a_1 t^4 + a_2 t^8 + ...) * (a_0 + a_1 t + a_2 t^2 + ...) so essentially we only have to shift bits and add */ int_fifth(u) = { local(u1, u2); u1 = 0; u2 = u; while(u, if(u % 2, u1 = int_sum(u1, u2)); u = u >> 1; u2 = u2 << 4; ); return(u1); }
/* Make a random triple of integers which correspond to polynomials of (at most) a given degree. This works because in pari/gp the function random(n) returns a random integer in {0, ..., n - 1}. */ random_triple(d) = { local(n); n = 2^(d + 1); return([random(n), random(n), random(n)]); }
/* Tests that tell us whether our triple is suitable. For now we just test whether the number of polynomials not divisible by t is odd. This might be a good idea because we can always split an actual solution into two parts of three each of which has this property. */ good_triple(v) = { return((v[1] + v[2] + v[3]) % 2); }
/* Find the sum of the fifth power of the triple. */ triple_result(v) = { local(u); u = int_fifth(v[1]); u = int_sum(u, int_fifth(v[2])); u = int_sum(u, int_fifth(v[3])); return(u); }
/* Convert a triple to a quadruple */ make_quadruple(v) = { return([triple_result(v), v[1], v[2], v[3]]); }
PART II: working with the list of quadruples
It turns out that we need a GLOBAL VARIABLE list. This variable will hold our list of quadruples. The reason for making it a global variable is that pari/gp scripting does not allow for pass by reference.
/* Initialize the global variable list if it hasn't been initialized yet. Otherwise set it to the empty list. */ initialize_list() = { if(type(list) == "t_LIST", listkill(list), list=listcreate()); }
The first version of insert_quadruple didn't work because of the issue mentioned above! Here is a version of insert_quadruple that does work.
/* Insert quadruple w into list list at the nth spot */ insert_quadruple(w, n) = { listinsert(list, w, n); }
PART III: Final program
The following code should utilize the functions above to search for curves of degree d. I haven't tried running it yet, because (1) there are probably bugs and (2) what do we do to the list once we've found a curve? This code doesn't insert it into the running list. Also, printing results is terrible, so this will be patched up soon.
/* This is a binary search algorithm that either inserts a quadruple and returns 1 or, if it finds a match, calls go_do_something and returns 0. */ search(listsize,quadr) = { local(value, lower, upper, middle); value = quadr[1]; lower = 1; upper = listsize; middle = 1; while(lower <= upper, middle = floor((upper+lower)/2); if(list[middle][1] < value, lower = middle + 1, if(list[middle][1] > value, upper = middle - 1; , go_do_something(list[middle], quadr); return(0); ); ); ); insert_quadruple(quadr,lower); return(1); }
/* This codes checks the pair of matching triples to see if the solution is "good". If so it prints the six polynomials. This code should be looked at a bit more and it should be hooked up to the code computing splitting types!!! */ go_do_something(v1, v2) = { local(d, ppp, h, A); for(i = 2, 4, for(j = 2, 4, if( v1[i] == v2[j], return))); ppp = [convert_to_pol(v1[2]),convert_to_pol(v1[3]),convert_to_pol(v1[4]),convert_to_pol(v2[2]),convert_to_pol(v2[3]),convert_to_pol(v2[4])]; h = gcd(Mod(1,2)*ppp[1], Mod(1, 2)*ppp[2]); h = gcd(h, Mod(1, 2)*ppp[3]); h = gcd(h, Mod(1, 2)*ppp[4]); h = gcd(h, Mod(1, 2)*ppp[5]); h = gcd(h, Mod(1, 2)*ppp[6]); ppp = lift((Mod(1, 2)*ppp)/h); d = 0; for(i = 1, 6, d = max(d, poldegree(ppp[i], t))); A = matrix(6, d + 1, i, j, polcoeff(ppp[i], j, t)*Mod(1, 2)); if(matrank(A) == 6, print("Here is one: "); print(v1, " ", v2); print(convert_to_pol(v1[2])); print(convert_to_pol(v1[3])); print(convert_to_pol(v1[4])); print(convert_to_pol(v2[2])); print(convert_to_pol(v2[3])); print(convert_to_pol(v2[4])); ); }
/* This runs everything for eternity */ find_curves(d) = { local(i, trip); initialize_list(); i = 0; while(1, trip = random_triple(d); if(good_triple(trip), i = i + search(i, make_quadruple(trip)); if(i % 10000 == 0, print(i)); ); ) }
/* Check the degree of a triple is correct. This assumes the degree is at most d! */ degree_triple(v, d) = { return(bittest(v[1], d) | bittest(v[2], d) | bittest(v[3], d)); }
/* Loop over all triples of a given degree. We may order the triples to reduce the size of the search space a bit. This needs to be checked more carefully! Note that the beginning of the loops are useless, so the program starts out producing no output for quite a while. */ complete_loop(d) = { local(i, trip); trip = vector(3); initialize_list(); i = 0; for(u1 = 0, 2^(d + 1) - 1, for(u2 = u1 + 1, 2^(d + 1) - 1, for(u3 = u2 + 1, 2^(d + 1) - 1, trip[1] = u1; trip[2] = u2; trip[3] = u3; if(good_triple(trip) && degree_triple(trip, d), i = i + search(i, make_quadruple(trip)); if(i % 10000 == 0, print(i)); ); ); ); ); }