#!perl -w
use strict;
# for faster mult algorithm:
# http://www.dms.auburn.edu/faculty/hankerson/papers/nb.pdf
#
$|=1;

#  degree of a polynomial is the highest exponent occuring in it.
#
#  order of a polynomial over field GF(q)
#  is smallest e for which  (x^e-1) mod P(x) = 0
#
#  a polynomial P(x) of degree n over GF(2) 
#  is primitive if the order = 2^n-1
#
########################################################
package Polynomial;
use strict;

BEGIN {
    no strict 'refs';
    for my $method (qw(add cube divide equals inverse is_zero multiply negative square subtract)) {
        *$method = sub {
            $_[0]->{field}->$method(@_);
        }
    }
}

sub new {
    my ($class, $field, $value)= @_;
    return bless {
        field=>$field,
        value=>$value,
    }, $class;
}
sub copy {
    my ($self)= @_;
    return ref($self)->new($self->{field}, $self->{value});
}

sub asstring {
    my ($self)= @_;
    my @bits= split //,unpack 'b*', pack 'V', $self->{value};
    return join '+',reverse map { "x^$_" } grep $bits[$_],0..$#bits;
}

########################################################
package FieldF2m;
use strict;
use integer;
sub new {
    my ($class, $p, $m, $type)= @_;
    return bless {
        p=>$p,          # polynomial, with highest term x^m
        m=>$m,
        type=>$type,
    }, $class;
}

sub zero {
    my ($self)= @_;
    return $self->value(0);
}

sub value {
    my ($self, $x)= @_;
    return $self->{type}->new($self, $x);
}

sub equals {
    my ($self, $x, $y)= @_;
    
    (undef, $x->{value}) = polydiv($self->{m}, $x->{value}, $self->{p}) if ($x->{value}>$self->{p});
    (undef, $y->{value}) = polydiv($self->{m}, $y->{value}, $self->{p}) if ($y->{value}>$self->{p});

    return $x->{value}==$y->{value};
}
sub is_zero {
    my ($self, $x)= @_;
    return $self->equals($x, $self->zero());
}

sub add {
    my ($self, $x, $y)= @_;
    #printf("F%d: %d + %d = %d\n", $self->{p}, $x->{value}, $y->{value}, ($x->{value}+$y->{value}) % $self->{p} );
    return $self->value( ($x->{value}^$y->{value}) );
}
sub subtract {
    my ($self, $x, $y)= @_;
    #printf("F%d: %d - %d = %d\n", $self->{p}, $x->{value}, $y->{value}, ($x->{value}-$y->{value}) % $self->{p} );
    return $self->value( ($x->{value}^$y->{value}) );
}

sub is_irreducible {
    my ($p, $m)= @_;

    # tests if a prime q exists such that:
    # a[0]..a[m-1] all divisible by q
    # a[0] not divisible by q^2
    # a[m] not divisible by q

    # if m+1 is a prime -> p=2^(m+1)-1  is irreducible

    # product of all prime polynomials in GF(q)[x]
    #   is  x^(q^n)-x

    # todo.
}

# returns ($a,$b)  such that  $a*$x+$b*$y = gcd($x, $y)
# gcd(a, b) = gcd(b, a%b)
sub euclid {
    my ($self, $x, $y)= @_;
    printf("euclid(%b, %b)\n", $x->{value}, $y->{value});
    return ($self->value(1), $self->zero()) if ($y->is_zero());

    my ($q, $r)= polydiv($self->{m}, $x->{value}, $y->{value});
    printf("polydiv: %b = %b * %b + %b\n", $x->{value}, $y->{value}, $q, $r);

    my ($a, $b)= $self->euclid($y, $self->value($r));
    # r=x-y*q
    # a*y+b*r = a*y+b*(x-y*q) = b*x + (a-b*q)*y
    my $c= $a->subtract($b->multiply($self->value($q)));
    my $ggd= $b->multiply($x)->add($c->multiply($y));
    printf("%d * %d + %d * %d = %d\n", $b->{value}, $x->{value}, $c->{value}, $y->{value}, $ggd->{value});
    return ($b, $c);
}



sub find_highestbit {
    my ($x)= @_;
    my $n=0;
    while ($x >>= 1) {
        $n++;
    }
    return $n;
}
# divide x by y, returns x/y, x%y
sub polydiv {
    my ($m, $x, $y)= @_;

    return 1 if ($x==$y);
    my $i= find_highestbit($x);
    printf("highest bit(%b) =%d\n", $x, $i);

    my $z= 0;
    while ($i>=$m) {
        $z<<=1;
        if ($x&(1<<$i)) {
            $x ^= $y<<($i-$m);
            $z|=1;
        }
        $i--;
    }
    return ($z, $x);
}
sub multiply {
    my ($self, $x, $y)= @_;
    #printf("F%d: %d * %d = %d\n", $self->{p}, $x->{value}, (ref $y ? $y->{value} : $y ), ($x->{value}*(ref $y ? $y->{value} : $y )) % $self->{p} );
    if (ref $y) {
        my $prod= 0;
        my $xval= $x->{value};
        my $yval= $y->{value};

        while ($xval>0) {
            if ($xval&1) {
                $prod ^= $yval;
            }

            $yval <<= 1;

            $yval^=$self->{p}  if ($yval & (1<<$self->{m}));

            $xval >>= 1;
        }
        return $self->value($prod);
    }
    elsif ($y&1) {
        return $x;
    }
    else { # $y&1 == 0
        return $self->zero();
    }
}

sub calcparity {
    my ($v)= @_;
    $v ^= $v >> 16;
    $v ^= $v >> 8;
    $v ^= $v >> 4;
    $v &= 0xf;
    return (0x6996 >> $v) & 1;
}

sub inverse {
    my ($self, $x)= @_;

    # todo
    # solve a[i]  from  sum(a[i]*x^i, i=m..0) * sum(x[i]*x^i, i=m..0) = 1
    #
    # sum( c[j]*x^j, j=2m..0 ) = 1
    #
    # c[j] = sum(a[i]*x[k],  i+k=j)
    #
    # ... resolve terms > x^m
    #
    # -> sum( d[i]*x^i, i=m..0 )


    # [ I | P ] * [ a0  0  0   0] * b  = 1
    #             [ a1 a0  0   0]
    #             [ a2 a1 a0   0]
    #             [      ...    ]
    #             [           an]
    #
    #  I = unit matrix
    #  Pk(x) = x^(n+k)= sum(Pki*x^i)
    #  a = vector (a0 .. an) ... number to invert
    #  b = vector (b0 .. bn) ...  value to solve
    #  1 = vector ( 1, 0, 0, ... )

}
sub divide {
    my ($self, $x, $y)= @_;
    return $self->multiply($x, $y->inverse());
}

sub square {
    my ($self, $x)= @_;
    return $self->multiply($x, $x);
}
sub cube {
    my ($self, $x)= @_;
    return $self->multiply($x, $self->square($x));
}

sub test {
    my ($self)= @_;

    printf("pdiv: %b / %b = (%b, %b)\n", 0b1100101, 0b10011, polydiv($self->{m}, 0b1100101, 0b10011));

    {
        my ($x, $y)= ($self->value(0b1100101), $self->value(0b10011));
        my ($a,$b)= $self->euclid($x, $y);
        my $gcd= $a->multiply($x)->add($b->multiply($y));
        printf("euclid: %d * %d + %d * %d = %d\n",
            $a->{value}, $x->{value}, $b->{value}. $y->{value}, $gcd->{value});
    }

    my %vals;
    printf("generators for GF(2)[%d] with f(x)=%d\n", $self->{m}, $self->{p});
    for (my $i=0 ; $i<(1<<$self->{m}) ; $i++) {
        next if (exists $vals{$i});

        my $x= $self->value($i);
        my $j=0;
        my $y;
        my @list;
        my %hist;
        for ($y=$x->copy() ; $j==0 || !$x->equals($y) ; $y=$y->multiply($x)) {
            $vals{$y->{value}}++;
            last if ($hist{$y->{value}});
            $hist{$y->{value}}++;
            push @list, $y->{value};
            $j++;
        }
        printf("%2d:%2d:%s%s\n", $i, $j, join(",",@list), 
            $x->equals($y) ? "" : sprintf(" .. %d", $y->{value}));
    }

    # x^(2^0)-x = 0
    # x^(2^1)-x = x^2-x
    # x^(2^2)-x = x^4-x
    print "x^(2^n)  mod p:  ";
    my $x=$self->value(2);  # x
    for (my $n=0 ; $n<$self->{m} ; $n++)
    {
        printf(" %d", $x->{value});
        $x= $x->square();
    }
    print "\n";
}

# printf("%b / %b = (%b, %b)\n", 0b10010, 0b10011, polydiv(0b10010, 0b10011));
# printf("%b / %b = (%b, %b)\n", 0b100100, 0b10011, polydiv(0b100100, 0b10011));
# printf("%b / %b = (%b, %b)\n", 0b1001000, 0b10011, polydiv(0b1001000, 0b10011));
# 
# printf("%b / %b = (%b, %b)\n", 0b11010, 0b10011, polydiv(0b11010, 0b10011));
# printf("%b / %b = (%b, %b)\n", 0b110100, 0b10011, polydiv(0b110100, 0b10011));
# printf("%b / %b = (%b, %b)\n", 0b1101000, 0b10011, polydiv(0b1101000, 0b10011));
# 
# printf("%b / %b = (%b, %b)\n", 0b10000, 0b10011, polydiv(0b10000, 0b10011));
# 
# 10011/1100101\11
#       10011
#       -----
#        10100
#        10011
#         -----
#          1111

########################################################
package Point;
use strict;
sub new {
    my ($class, $x, $y)= @_;
    return bless {
        x=>$x, y=>$y,
    }, $class;
}
sub copy {
    my ($self)= @_;
    return ref($self)->new($self->{x}, $self->{y});
}
sub distance2 {
    my ($self, $p)= @_;
    return $self->{x}->subtract($p->{x})->square()->add( $self->{y}->subtract($p->{y})->square() );
}
sub equals {
    my ($self, $p)= @_;
    return $self->distance2($p)->is_zero();
}

########################################################
package EllipticCurve;
use strict;
# y^2=x^3+a*x+b
sub new {
    my ($class, $a, $b)= @_;
    return bless {
        a=>$a, b=>$b,
    }, $class;
}

sub is_on_curve {
    my ($self, $pt)= @_;

    return 1 if !defined $pt->{y} ;
    # x^3+a*x^2 + b - y^2 -x*y == 0
    return $pt->{x}->cube()
        ->add($self->{a}->multiply($pt->{x}->square()))
        ->add($self->{b})
        ->subtract($pt->{y}->square())
        ->subtract($pt->{y}->multiply($pt->{x}))->is_zero();
}

sub test {
    my ($self)= @_;

    my $fp= $self->{a}{field};
    my $M= 1<<$fp->{m};

    print "points on curve: ";
    for (my $x=0 ; $x<$M ; $x++) {
    for (my $y=0 ; $y<$M ; $y++) {
        my $p= Point->new($fp->value($x), $fp->value($y));
        printf(" (%d,%d)", $x, $y) if $self->is_on_curve($p);
    }
    }
    print "\n";
}
########################################################
package ECGroup;
use strict;

sub new {
    my ($class, $ec)= @_;

    # 4*a^3 + 27*b^2 == 0
    die "not a group\n" if ( $ec->{a}->cube()->multiply(4)->add($ec->{b}->square()->multiply(27))->is_zero() );
    return bless {
        ec=>$ec,
    }, $class;
}

sub negate {
    my ($self, $pt)= @_;
    die "not on curve\n" unless ($self->{ec}->is_on_curve($pt));

    return Point->new($pt->{x}, defined $pt->{y} ? $pt->{y}->negative() : $pt->{y});
}

sub identity {
    my ($self)= @_;

    return Point->new();
}
sub add {
    my ($self, $p, $q)= @_;

    return $q if (!defined $p->{y});
    return $p if (!defined $q->{y});


    #printf("p=(%g,%g) q=(%g, %g)\n", $p->{x}{value}, $p->{y}{value}, $q->{x}{value}, $q->{y}{value});
    if ($p->equals($q)) {
        if ($p->{y}->is_zero()) {
            return $self->identity;
        }

        # s = y' solved from  d(y^2=x^3+a*x+b)/dx
        my $s= $p->{x}->square()->multiply(3)->add($self->{ec}{a})->divide($p->{y}->multiply(2));
        my $xr= $s->square()->subtract($p->{x}->multiply(2));
        my $yr= $s->multiply($p->{x}->subtract($xr))->subtract($p->{y});

        printf("s=%g xr=%g yr=%g\n", $s->{value}, $xr->{value}, $yr->{value});
        return Point->new($xr, $yr);
    }
    elsif ($p->{x}->equals($q->{x})) {
        return $self->identity;
    }
    else {
        # s = slope of line.
        my $s= $p->{y}->subtract($q->{y})->divide($p->{x}->subtract($q->{x}));
        my $xr= $s->square()->subtract($p->{x}->add($q->{x}));
        #printf("%d %d %d\n", $s->square()->{value}, $p->{x}{value}, $q->{x}{value});
        my $yr= $s->multiply($p->{x}->subtract($xr))->subtract($p->{y});

        #printf("s=%g xr=%g yr=%g\n", $s->{value}, $xr->{value}, $yr->{value});
        return Point->new($xr, $yr);
    }
}
sub mul {
    my ($self, $p, $n)= @_;

    return $self->identity if (!defined $p->{y});

    my $r= $self->identity;
    my $pp= $p->copy();
    while ($n) {
        if ($n&1) {
            $r= $self->add($r, $pp);
        }

        $pp= $self->add($pp, $pp);
        $n >>= 1;
    }
    return $r;
}

########################################################
package main;
use strict;
use warnings;

my $fp= FieldF2m->new(19, 4, "Polynomial");
$fp->test();

my $a= $fp->value(13);
my $b= $fp->value(9);
my $c= $fp->multiply($a, $b);
my $d= $fp->add($a, $b);
my $e= $fp->value(2)->square();
my $f= $e->square();

printf("a=%d b=%d a*b=%d  a+b=%d 2^2=%d 2^2^2=%d\n", $a->{value}, $b->{value}, $c->{value}, $d->{value}, $e->{value}, $f->{value});


my $ec= EllipticCurve->new($fp->value(3), $fp->value(1));
$ec->test();

for (my $m=0 ; $m<6 ; $m++) {
    for (my $p=(1<<$m)+1 ; $p<(2<<$m) ; $p++) {
        my $fp= FieldF2m->new($p, $m, "Polynomial");
        $fp->test();
    }
}