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:

• A polynomial is a polynomial in t whose coefficients are in {0, 1}\
• When we talk about the sum or product of polynomials we mean sum or product in F_2[t]
• An integer means a non-negative integer
• We convert a polynomial into an integer and vice versa as follows: the polynomial a_0 + a_1 * t + a_2 * t^2 corresponds to the integer a_0 + a_1 * 2 + a_2 * 4 + …
• The only thing we ever do with polynomials is add and take fifth powers
• We perform all operations on the side of the integers and not in terms of polynomials. The reason for this is that I think it will be much faster and consume much less memory.
• A triple is a vector [u1, u2, u3] of integers
• A quadruple is a vector [u0, u1, u2, u3] where u1, u2, u3 are integers and u0 is the integer corresponding to the polynomial f1^5 + f2^5 + f3^5 where f1, f2, f3 are the polynomials corresponding to u1, u2, u3
• A list is an ordered list of quadruples [q1, q2, …]. The ordering is such that q1[1] < q2[1] < … We do not allow consecutive entries to have the same first coordinate.

/* 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));
      );
    );
  );
);
}