To run the code below you need to download and compile Mike Scott's miracl library. Be sure to read the license terms of miracl if you're doing any commercial work, check out the miracl homepage

Here is some code I wrote to handle Points, Divisors and Hyperelliptic Curves. Some examples of the programs contained within are below. Some of the code is experimental, and any corrections/suggestions are welcome. My code is released under the GNU General Public License.

hyper.tar.gz

Recent Updates: Sakai and Sakurai's algorithm to calculate the number of points on the jacobian of a hyperelliptic curve added. It can handle curves of arbitrary genus over fields of characteristic 2.


1. Compiling miracl (on linux)

Download miracl and unzip it in an empty directory like this;

unzip -j -aa -L miracl.zip

Then type;

bash linux

This will create an archive file "miracl.a". To see the large numbers in the point counting example below, instead of in exponential notation, change line 307 of big.cpp from;

if (e>=-2 && e<=2)

to;

if (e>=-7 && e<=7)

You also need to compile other classes;

g++ -c complex.cpp floating.cpp gf2m.cpp big.cpp poly2.cpp poly2xy.cpp poly.cpp polyxy.cpp monty.cpp

Then;

ar r miracl.a complex.o floating.o gf2m.o big.o poly2.o poly2xy.o poly.o polyxy.o monty.o


2. Point Counting over Fields of Chararacteristic 2 - Koblitz algorithm

points.cc is a program I wrote to calculate the number of points on the jacobian of a hyperelliptic curve of genus 2, over fields of characteristic 2. This is the method described by Koblitz.

To compile points.cc on *nix do the following;

g++ -I (miracl_path) points.cc miracl.a -o points

The curve itself is specified by the line; "equation = pow2Y(y,2) + pow2Y(y,1) + pow2X(x,5) + pow2X(x,3) + pow2X(x,1);", where equation is a Poly2XY instance. Simply change this line to count points over different equations. Some other examples are included in comments in the code. The extension field is requested from the user while the program is running

Example (curve and extension field taken from Koblitz - page 149);

./points
Enter field extension(default = 2):
101
Curve: x^5 + x^3 + x + y^2 + y
Working mod: x^2 + x + 1
Points on Curve over field 2^1 = 3
Points on Curve over field 2^2 = 9
Points on the jacobian over F_{2^101}:
6427752177035961102167848369367185711289268433934164747616257

3. Point Counting over Fields of Chararacteristic 2 - Sakai and Sakurai's algorithm

points2.cc is a program I wrote to calculate the number of points on the jacobian of a hyperelliptic curve of arbitrary genus, over fields of characteristic 2. This algorithm was given in the paper "Design of Hyperelliptic Cryptosystems in Small Characteristic and a Software Implementation over F_{2^n} - Sakai and Sakurai".

My program handles curves up to genus 10, over fields of characteristic 2. To handle curves of higher genera, just add in the required irreducible polynomials. To handle fields of characteristic greater than 2, the irreducible polynomials have to be changed.

To compile points2.cc on *nix do the following;

g++ -I (miracl_path) points2.cc miracl.a -o points2

The curve itself is specified by the "equation" variable, where equation is a Poly2XY instance. 5 example curves are included in the code, ranging from genus 2 to genus 6; these curves are taken from the Sakai and Sakurai paper. If changing to a different genus curve, remember to change the "genus" integer to the required genus. The extension field is requested from the user while the program is running

Example (over field F_{2^{89}});

./points2
Enter field extension(default = 1):
89
Curve: x^5 + x^2.y + x.y + x + y^2 + y + 1
Points on Curve over field 2^1: = 1
Points on Curve over field 2^2: = 7
Points on the jacobian over F_{2^89}:
383123885216484912146996836504217327230624063025829938

4. Point Counting over Small Prime Fields using the Hasse-Witt matrix

hasse.cc is a program I wrote to calculate the number of points on the jacobian of a hyperelliptic curve of genus 2, over small prime fields.

This program will handle any prime field that will fit in a c++ "int". There is no point supporting larger prime fields, as the algorithm is too inefficient.

To compile hasse.cc on *nix do the following;

g++ -I (miracl_path) hasse.cc Point.o Divisor.o Cantor.o poly_div.o miracl.a -o hasse

The curve itself is specified by the "f" variable, which is the polynomial f(x) in "y^2 = f(x)". Up to 5 possible group orders may be computed depending on the curve and the field.

Example (over field F_{131});

./hasse
Curve: y^2 = x^5 + 3*x
Prime Field: 131
s1: 0
s2: -62
The order of the jacobian is one of the following:
L: 17100
L: 17362
Checking each possible group order...
u1: x - 41
v1: - 29
The order is: 17100

5. Furukawa, Kawazoe, Takahashi Point Counting method

furukawa.cc is a program I wrote to calculate the number of points on the jacobian of a hyperelliptic curve of genus 2, of the form y^2 = x^5 + ax, a in F_p, p a prime and congruent to 1 (mod 8), and the jacobi symbol (a/p) = -1.

This is the algorithm described in the paper "Counting Points for Hyperelliptic Curves of type $y^2=x^5+ax$ over Finite Prime Fields" by Eisaku Furukawa and Mitsuru Kawazoe and Tetsuya Takahashi. (http://eprint.iacr.org/2002/181.pdf) - p. 11/12

To compile furukawa.cc on *nix do the following;

g++ -I (miracl_path) furukawa.cc Point.o Divisor.o Cantor.o poly_div.o miracl.a -o furukawa

Up to 3 possible group orders may be computed depending on the curve and the field. The correct group order is calculated by generating a random divisor on the curve and multiplying it in turn by each group order, to find out which one gives the identity element of the jacobian.

Example (over field F_{1208925819614629175095961}, y^2 = x^5 + 3x);

./furukawa
p: 1208925819614629175095961
c: 919996162953
d: 425753966374
s1: -1703015865496
t: 241205699450915806567047
The order of the jacobian is: 1461501637332961738997140052922587693332824903682

I have implemented a variant of the above method in the file "furukawa_var.cc", which uses an explicit formula to speed up the computation. This formula was published in the paper "Suitable Curves for Genus-4 HCC over Prime Fields: Point Counting Formulae for Hyperelliptic Curves of type y^2 = x^(2k+1) + ax" by Haneda, Kawazoe and Takahasi.


6. Furukawa, Kawazoe, Takahashi Point Counting method 2

furukawa2.cc is a program I wrote to calculate the number of points on the jacobian of a hyperelliptic curve of genus 2, of the form y^2 = x^5 + a, a in F_p, p a prime and congruent to 1 (mod 5), and p greater than 64.

This is the algorithm described in the paper "Counting Points for Hyperelliptic Curves of type $y^2=x^5+ax$ over Finite Prime Fields" by Eisaku Furukawa and Mitsuru Kawazoe and Tetsuya Takahashi. (http://eprint.iacr.org/2002/181.pdf) - p. 13/14

To compile furukawa2.cc on *nix do the following;

g++ -I (miracl_path) furukawa2.cc miracl.a -o furukawa2

Up to 5 possible group orders may be computed depending on the curve and the field. The correct group order is calculated by generating a random divisor on the curve and multiplying it in turn by each group order, to find out which one gives the identity element of the jacobian.

NOTE: Not included in code bundle (yet).

Example

./furukawa2

7. Point Counting over Prime Fields, Koblitz method

points3.cc is a program I wrote to calculate the number of points on the jacobian of a hyperelliptic curve of genus 2/3, of the form v^2 + v = u^n, over a prime field, where n = 2g + 1 is a prime and p is congruent to 1 (mod n). p must also be a generalized mersenneprime of the form (a^n - 1)/(a - 1).

See Koblitz "Algebraic aspects of Cryptography" p. 150 - 153.

To compile points3.cc on *nix do the following;

g++ -I (miracl_path) points3.cc miracl.a -o points3

This program can also calculate twists of the curve, see Koblitz for details.

Example (a = 100003, p = 100013000640014200121, i = 0, j = 1)

./points3
-100003 1 10000600009 -100003 0
twist: 1
The order of the jacobian is: 10002600296019310791620950350903416514905

8. Point Counting over Prime Fields, Koblitz method (using LLL)

points4.cc is a program I wrote to calculate the number of points on the jacobian of a hyperelliptic curve of genus 2, of the form v^2 + v = u^5, over a prime field, where p is congruent to 1 (mod 5). p does not have to be a generalized mersenne prime as in the previous algorithm. This algorithm makes use of the integral LLL algorithm as detailed elsewhere on this page.

See the Buhler and Koblitz paper, "Lattice basis reduction: Jacobi sums and hyperelliptic cryptosystems".

To compile points4.cc on *nix do the following;

g++ -c -I (miracl_path) lll.cc

g++ -I (miracl_path) points4.cc lll.o miracl.a -o points4

This program can also calculate twists of the curve, see Koblitz for details.

Example (a = 10000600009, p = 100013000640014200121, i = 0, j = 0)

./points4
Transformation matrix (from LLL):
-100002 0 1 -100003
0 0 100004 -100003
-100003 0 1 -100004
0 1 1 -100003
-100003 10000600009 0 1 -100003
twist: 1
The order of the jacobian is: 10002600296019260783120375331503087512625

9. Finding the points on a Curve

curve_points.cc is a program which constructs an object of type Curve over a prime field, given the f and h polynomials. It then finds out the points on that curve.

To compile curve_points.cc on *nix do the following;

g++ -I (miracl_path) -c Point.cc Divisor.cc Curve.cc

g++ -I (miracl_path) curve_points.cc Point.o Divisor.o Curve.o poly_div.o miracl.a -o curve_points

Note that miracl sometimes stores polynomials as minus rather than plus, as we're working with a modulus.

Example curve over prime field F7;

C: y^2 + xy = x^5 - 2*x^4 - x^2 + x - 4
Genus = 2
8 points = (INF)(1,1)(1, - 2)(2,2)(2, - 4)( - 2, - 4)( - 2, - 1)( - 1, - 3)

10. Getting the Divisor of a Polynomial

poly_div.cc is a program I wrote to find out the Divisor of a given polynomial in 2 variables, defined over a hyperelliptic curve. poly_div_test.cc is a test program for it. Note that you need to find the zeros of the polynomial, so this method is impractical if the prime-field over which the curve is defined is large.

To compile poly_div_test.cc on *nix do the following;

g++ -I (miracl_path) -c Point.cc Divisor.cc poly_div.cc

g++ -I (miracl_path) poly_div_test.cc Point.o Divisor.o poly_div.o miracl.a -o poly_div_test

The equation is specified by the line: "equation = powY(u,2) + powX(u,1) * powY(u,1) + 6*powX(u,4) + 6*powX(u,3) + powX(u,2) + 6*powX(u,1);". This can be changed to find the divisors of different polynomials. The hyperelliptic curve also needs to be defined.

Example (Koblitz - page 178 - exercise 4);

equation: - x^4 - x^3 + x^2 + x.y - x + y^2
C: y^2 + yx = x^5 - 2*x^4 - x^2 + x - 4
Divisor: (1,1) + (1, - 2) + 2(2,2) + 2(2, - 4) + 4( - 1, - 3) - 10(INF)
Divisor is not semi-reduced
Divisor is not reduced

11. Calculating the Group Law

Cantor.cc and Cantor.h contain composition and reduction algorithms for adding two divisors on the jacobian of a hyperelliptic curve. It also contains a "doubling" function that will compute kD = D + D + D ... + D (k-times), using the simple "double-and-add" algorithm. Cantor2.cc and Cantor2.h contain the group law over fields of characteristic 2.

CantorTest.cc gives an example of the "doubling" function, where a divisor of weight 1 is multiplied by "1461501637332961738997140052922587693332824903682", which is the order of the curve y^2 = x^5 + 3x over the prime field F_1208925819614629175095961. A group element to the power of the group order should be the identity element of the group, which in the case of divisors, is (1, 0)

To compile CantorTest.cc on *nix do the following;

g++ -I (miracl_path) -c Cantor.cc

g++ -I (miracl_path) CantorTest.cc Cantor.o miracl.a -o CantorTest

Example (output of CantorTest.cc)

./CantorTest
Reduced divisor is:
udash: 1
vdash: 0

12. Calculating a random divisor

Divisor.cc and Divisor.h contain a function called "randomDivisor" which computes a random divisor for a genus 2 curve over a prime field. This divisor can take the form: D = P1 - inf, D = 2*P1 - 2*inf, or D = P1 + P2 - 2*inf

TestRandom.cc gives a test curve of genus 2, over a large prime field, and generates a divisor using the function in Divisor.cc. It then checks the answer by multiplying it by the group order, to get the identity element.

Note: it will generate the same random sequence of divisors every time, so the random functions in the Divisor class need to be seeded to get true randomness

To compile TestRandom.cc on *nix do the following;

g++ -I (miracl_path) -c Point.cc

g++ -I (miracl_path) -c Divisor.cc

g++ -I (miracl_path) -c Cantor.cc

g++ -I (miracl_path) TestRandom.cc Point.o Divisor.o Cantor.o poly_div.o miracl.a -o TestRandom

Example (output of TestRandom.cc)

./TestRandom
u1: x + 167169283870159064269152
v1: 81720149173104813595542
Reduced divisor is:
udash: 1
vdash: 0

13. Integral LLL algorithm

lll.cc contains the integral LLL algorithm, where the basis of a lattice is given by its Gram Matrix. See "A Course in Computational Algebraic Number Theory - Henri Cohen, page 94". This only works for 4x4 gram matrices, as the code is designed for use with genus 2 curves, but I can extend this to arbitrary size gram matrices easily if requested.

llltest.cc constructs a gram matrix, and outputs the transformation matrix after LLL is applied to it (the gram matrix). This example can be checked against LLL algorithms such as in gp-pari.

To compile llltest.cc on *nix do the following;

g++ -I (miracl_path) -c lll.cc

g++ -I (miracl_path) llltest.cc lll.o miracl.a -o llltest

Example (output of llltest.cc)

./llltest
Gram Matrix (Input):
31382404, -81229, -551797, -3845773,
-81229, 214, 1427, 9953,
-551797, 1427, 9706, 67619,
-3845773, 9953, 67619, 471286,
Transformation Matrix (Output):
0 0 0 1
1 -7 0 7
0 1 -7 7
0 0 1 7

14. Reverse Mumford operation

Divisor.cc contains a function called "reverseMumford", which takes in two mumford polynomials, and returns the divisor that they represent. MumfordTest.cc contains a brief program to test this function

To compile MumfordTest.cc on *nix do the following;

g++ -I (miracl_path) -c poly_div.cc Point.cc Divisor.cc

g++ -I (miracl_path) MumfordTest.cc Point.o Divisor.o poly_div.o miracl.a -o MumfordTest

Example (output of MumfordTest.cc)

./MumfordTest
Input Polynomials:
u: x^2 - x - 2
v: - 3*x + 1
Divisor is semi-reduced
Divisor is reduced
Output Divisor:
(2,2) + ( - 1, - 3) - 2(INF)

15. Transforming a Curve

Curve.cc contains a function called "transform", which will transform a curve from the form C:v^2 + h(u)v = f(u) into C:v^2 = f(u). This is done by the variable substitution; v to (v - h(u)/2). See Koblitz, page 156. The characteristic of the field can't be 2! Transform.cc contains a test program to show how to use the function.

To compile Transform.cc on *nix do the following;

g++ -I (miracl_path) -c Point.cc Curve.cc

g++ -I (miracl_path) Transform.cc Point.o Curve.o miracl.a -o Transform

Example (output of Transform.cc: h = 1, f = x^5, p = 100013000640014200121)

./Transform
Old Curve: v^2 + 1v = x^5
New Curve: v^2 = x^5 - 25003250160003550030


Note: All the code on this webpage is released under the General Public License and may be freely altered and distributed under the terms of this license. All code written by me is copyright © 2004 Colm O hEigeartaigh. Please forward any queries to coheigeartaigh AT computing DOT dcu DOT ie