#!perl -w
use strict;

my $DELTA=0.001;

#
# y^2+a*x*y+b*y = x^3+c*x^2+d*x+e
#
# we can eliminate the  x*y and y terms as follows:
# rewrite the left side as  y*(y+a*x+b)
# then replace y with  y-(a*x+b)/2
# we  get  (y-(a*x+b)/2)*(y+(a*x+b)/2)  = y^2 - {(a*x+b)/2}^2
#
# the x^2 term can be eliminated as follows:
# replace x with  x-c/3
# -> x^3 + (d-c^2/3)*x - c*d/3+2*c^3/27+e
#
# thus the general form of a elliptic curve can be written as:
#
# ec :  y^2 = x^3 + a * x +  b
#       with  4*a^3+27*b^2 != 0  -> no duplicate factors
#
#mupad:
#  plot ( plot::Implicit2d(x^3-3*x+2-y^2, x=-3..3, y=-3..3, Contours=[c $ c=-1..1]) )
#
# 
#
#      
#           /
#          /      curve with b < -sqrt(-4/27*a^3)
#          \          y^2 = x^3 -3*x-3
#           \
#
#      
#           /
#    *     /      curve with b == -sqrt(-4/27*a^3)
#          \          y^2 = x^3 -3*x-2
#           \
#
#   ___
#  /   \    /
# |     \  /      curve with -sqrt(-4/27*a^3) < b < sqrt(-4/27*a^3)
# |     /  \          y^2 = x^3 -3*x
#  \___/    \
#
#   ___
#  /   \  /
# |     \/        curve with b == sqrt(-4/27*a^3)
# |     /\            y^2 = x^3 -3*x+2
#  \___/  \
#
#   ___
#  /   \/
# |               curve with b > sqrt(-4/27*a^3)
# |                   y^2 = x^3 -3*x+3
#  \___/\
#
#
#  derivative of curve at position x:
#
#  y^2 = x^3 + a*x + b   ... diff by x
#  ->  2*y*y' = 3*x^2 + a  ->  y' = (3*x^2+a)/2y
#
# ec-group properties:
#
# negative of a point P=(Px,Py),   -P = (Px, -Py)
# sum of P=(Px,Py)  and Q=(Qx,Qy)
# P+Q = R,   such that  P, Q and -R are on one line.
#           -> P+Q-R = 0
#
# s = (Py-Qy)/(Px-Qx)   ( P!=Q )
# s = (3*Px^2+a)/(2*Py) ( P==Q )
#   Rx = s^2 - (Px+Qx)
#   Ry = -Py + s*(Px-Rx)
#

# identity:  point at infinity,   
#    P+(-P) = 0
#    P+ 0   = P
#  
# doubling P: raaklijn door P, snijdt curve in R.
#    P + P = 2*P = R
#
#
# line through P and Q is defined by:
#   (Px-Qx)*y-(Py-Qy)*x - Px*Qy+Py*Qx = 0
#  

# ... calculating R ( sum of P and Q )
#
# ec : y^2=x^3+a*x+b
# 
# P1=(x1,y1)
# P2=(x2,y2)
# 
# P1 and P2 are on a line
# 
# (y-y1)*(x2-x1)=(x-x1)*(y2-y1)
# => y = (x-x1)*s + y1,   met s= (y2-y1)/(x2-x1)
# 
# substitute y of line in ec
#
# -> x^3+a*x+b - ( (x-x1)*s + y1 )^2
# -> x^3+a*x+b - ( (x-x1)^2*s^2  + 2*(x-x1)*s * y1   + y1^2 )
# 
# -> x^3 - x^2*s^2 + (a + 2*s*(x1*s-y1))*x - (x1*s-y1)^2 + b
# 
# divide out the (x-x1)*(x-x2) known solutions:
#
# x-x1 / x^3 - s^2*x^2 + (a + 2*s^2*x1-2*s*y1)*x - (x1*s-y1)^2 + b \ x^2 +(x1-s^2)*x + (a + 2*s^2*x1-2*s*y1+(x1-s^2)*x1)
#        x^3 - x1*x^2
#       ---------------
#           (x1-s^2)*x^2
#           (x1-s^2)*x^2 -(x1-s^2)*x1*x
#           ----------------------------
#                         (a + 2*s^2*x1-2*s*y1+(x1-s^2)*x1) * x
#                         (a + 2*s^2*x1-2*s*y1+(x1-s^2)*x1) * x -(a + 2*s^2*x1-2*s*y1+(x1-s^2)*x1)*x1
#                         ------------------------------------------------------------------------------
#                                   - (x1*s-y1)^2 + b + (a + 2*s^2*x1-2*s*y1+(x1-s^2)*x1)*x1
# 
# 
# 
# x-x2 / x^2 +(x1-s^2)*x + (a + 2*s^2*x1-2*s*y1+(x1-s^2)*x1) \ x+(x1+x2-s^2)
#        x^2 - x2 * x
#       -------------
#             (x1+x2-s^2)*x
#             (x1+x2-s^2)*x -(x1+x2-s^2)*x2
#             -----------------------------
#                 (a + 2*s^2*x1-2*s*y1+(x1-s^2)*x1) + (x1+x2-s^2)*x2
# .
# -> third solution: x3=s^2-(x1+x2)
# 
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 equals {
    my ($self, $p)= @_;
    return (($self->{x}-$p->{x})**2 + ($self->{y}-$p->{y})**2) < $DELTA;
}
package Line;
use strict;
# a*x+b*y+c=0
sub new {
    my ($class, $a, $b, $c)= @_;
    return bless {
        a=>$a, b=>$b, c=>$c,
    }, $class;
}

# Py=a*Px+b
# Qy=a*Qx+b
# --------- -
# (Py-Qy)/(Px-Qx) = a

# line  a*x+b*y+c=0
# (Py-Qy)*x-(Px-Qx)*y

# P : (Py-Qy)*Px-(Px-Qx)*Py = Py*Px-Qy*Px-Px*Py+Qx*Py = Py*Qx - Px*Qy
# Q : (Py-Qy)*Qx-(Px-Qx)*Qy = Py*Qx-Qy*Qx-Px*Qy+Qx*Qy = Py*Qx - Px*Qy
sub FromPoints {
    my ($class, $p, $q)= @_;
    return $class->new($p->{y}-$q->{y}, $q->{x}-$p->{x}, $p->{x}*$q->{y}-$p->{y}*$q->{x});
}

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} ;
    my $val=  - $pt->{y}**2 + $pt->{x}**3 + $self->{a} * $pt->{x} + $self->{b};
    return 1 if (abs($val)<$DELTA);
    #printf("val=%g\n", $val);
    return 0;
}

package ECGroup;
use strict;

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

    die "not a group\n" if (abs(4*$ec->{a}**3+27*$ec->{b}**2)<$DELTA);
    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} : $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}, $p->{y}, $q->{x}, $q->{y});
    if ($p->equals($q)) {
        if (abs($p->{y})<$DELTA) {
            return $self->identity;
        }

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

        #printf("s=%g xr=%g yr=%g\n", $s, $xr, $yr);
        return Point->new($xr, $yr);
    }
    elsif (abs($p->{x}-$q->{x})<$DELTA) {
        return $self->identity;
    }
    else {
        # s = slope of line.
        my $s= ($p->{y}-$q->{y})/($p->{x}-$q->{x});
        my $xr= $s**2-($p->{x}+$q->{x});
        my $yr= -$p->{y}+$s*($p->{x}-$xr);

        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 {
    # y^2=x^3-7*x
    my $ec= EllipticCurve->new(-7, 0);
    my $g= ECGroup->new($ec);
    my $p= Point->new(-2.34, -sqrt(3.567096));
    my $q= Point->new(-0.1, sqrt(.699));
    printf("p=%d q=%d\n", $ec->is_on_curve($p), $ec->is_on_curve($q));
    my $r= $g->add($p, $q);
    printf("r=(%g,%g)\n", $r->{x}, $r->{y});
}
sub tst2 {
    my $ec= EllipticCurve->new(-3, 5);
    my $g= ECGroup->new($ec);
    my $p= Point->new(2, sqrt(7));
    printf("p=%d\n", $ec->is_on_curve($p));
    my $r= $g->add($p, $p);
    printf("p+p=(%g,%g)\n", $r->{x}, $r->{y});

    $r= $g->mul($p, 7);
    printf("7p=(%g,%g)\n", $r->{x}, $r->{y});
}
sub tst3 {
    my $ec= EllipticCurve->new(-7, -6);
    my $g= ECGroup->new($ec);
}
sub tst4 {
    my $ec= EllipticCurve->new(-3, 2);
    my $g= ECGroup->new($ec);
}


tst1();
tst2();
tst3();
eval { tst4() }; print $@;