#!perl -w
use strict;
$|=1;

########################################################
package Value;
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;
}


########################################################
package FieldFp;
use strict;
use integer;
sub new {
    my ($class, $p)= @_;
    return bless {
        p=>$p,
    }, $class;
}

sub zero {
    my ($self)= @_;
    return Value->new($self, 0);
}
sub one {
    my ($self)= @_;
    return Value->new($self, 1);
}
sub value {
    my ($self, $x)= @_;
    return Value->new($self, $x);
}

sub equals {
    my ($self, $x, $y)= @_;
    return (($x->{value}-$y->{value})%$self->{p})==0;
}
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}) % $self->{p} );
}
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}) % $self->{p} );
}
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} );
    return $self->value(($x->{value} * (ref $y ? $y->{value} : $y )) % $self->{p});
}

# returns ($a,$b)  such that  $a*$x+$b*$y = gcd($x, $y)
# gcd(a, b) = gcd(b, a%b)
sub euclid {
    my ($x, $y)= @_;
    return (1, 0) if ($y==0);

    my ($q, $r)= (int($x/$y), $x%$y);

    my ($a, $b)= euclid($y, $r);
    # r=x-y*q
    # a*y+b*r = a*y+b*(x-y*q) = b*x + (a-b*q)*y
    #printf("%d * %d + %d * %d = %d\n", $b, $x, $a-$b*$q, $y, $b*$x+($a-$b*$q)*$y);
    return ($b, $a-$b*$q);
}

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

    my ($a, $b)= main::euclid($x->{value}, $self->{p});
    #printf("F%d: %d^-1 = %d\n", $self->{p}, $x->{value}, $a);
    return $self->value($a);
}
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));
}
########################################################
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 + b - y^2 == 0
    return $pt->{x}->cube()->add($self->{a}->multiply($pt->{x}))->add($self->{b})->subtract($pt->{y}->square())->is_zero();
}

########################################################
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;
sub tst1 {
    for my $p (2,3,5,7,11,13,17,19,23) {
        my $fp= FieldFp->new($p);
        for my $a (0..$p-1) {
            for my $b (0..$p-1) {
                my $ec= EllipticCurve->new($fp->value($a), $fp->value($b));
                my $g= eval { ECGroup->new($ec) }; if ($@) { printf("%d.%d %s\n", $a, $b, $@); next; }

                printf("F%d ec(%d, %d)\n", $p, $a, $b);
                my $count= 0;
                for my $y (0..$p-1) {
                    printf("%2d:", $y);
                    for my $x (0..$p-1) {
                        #printf("testing %d,%d in F%d(%d,%d)\n", $x, $y, $p, $a, $b);
                        my $pt= Point->new($fp->value($x), $fp->value($y));
                        if ($ec->is_on_curve($pt)) {
                            #printf("%d.%d is on %d.%d\n", $x, $y, $a, $b);
                            print "*";
                            $count++;
                        }
                        else {
                            print " ";
                        }
                    }
                    print "\n";
                }
                printf("F%d ec(%d, %d) has %d points\n", $p, $a, $b, $count+1);
            }
        }
    }
}
sub tst3 {
    my $p = 17;
    my $fp= FieldFp->new($p);
    my $a= 10;
    my $b= 5;
    my $ec= EllipticCurve->new($fp->value($a), $fp->value($b));
    my $g= eval { ECGroup->new($ec) }; if ($@) { printf("%d.%d %s\n", $a, $b, $@); }
}
sub tst4 {
    my $p = 17;
    my $fp= FieldFp->new($p);
    my $a= 1;
    my $b= 7;
    my $ec= EllipticCurve->new($fp->value($a), $fp->value($b));
    my $g= eval { ECGroup->new($ec) }; if ($@) { printf("%d.%d %s\n", $a, $b, $@); }

    my $P= Point->new($fp->value(2), $fp->value(0));
    my $Q= Point->new($fp->value(6), $fp->value(3));
    printf("P=%d Q=%d\n", $ec->is_on_curve($P), $ec->is_on_curve($Q));

}
sub tst5 {
    my $p = 17;
    my $fp= FieldFp->new($p);
    my $a= 1;
    my $b= 7;
    my $ec= EllipticCurve->new($fp->value($a), $fp->value($b));
    my $g= eval { ECGroup->new($ec) }; if ($@) { printf("%d.%d %s\n", $a, $b, $@); }

    my $P= Point->new($fp->value(2), $fp->value(0));
    my $Q= Point->new($fp->value(1), $fp->value(3));

    my $R= $g->add($P,$Q);

    $P= Point->new($fp->value(1), $fp->value(3));
    my $PP= $g->add($P,$P);

    printf("P+Q=(%d,%d)\n", $R->{x}{value}, $R->{y}{value});
    printf("2P=(%d,%d)\n", $PP->{x}{value}, $PP->{y}{value});
}

tst1();
#tst3();
#tst4();
#tst5();