Elliptic curve arithmetic: Difference between revisions
(Updated D entry) |
(Updated D entry) |
||
| Line 127: | Line 127: | ||
double x, y; |
double x, y; |
||
@property static Pt zero() pure nothrow @nogc { |
@property static Pt zero() pure nothrow @nogc @safe { |
||
return Pt(double.infinity, double.infinity); |
return Pt(double.infinity, double.infinity); |
||
} |
} |
||
@property bool isZero() const pure nothrow @nogc { |
@property bool isZero() const pure nothrow @nogc @safe { |
||
return x > 1e20 || x < -1e20; |
return x > 1e20 || x < -1e20; |
||
} |
} |
||
@property static Pt fromY(in double y) nothrow /*pure*/ @nogc { |
@property static Pt fromY(in double y) nothrow /*pure*/ @nogc @safe { |
||
return Pt(cbrt(y ^^ 2 - bCoeff), y); |
return Pt(cbrt(y ^^ 2 - bCoeff), y); |
||
} |
} |
||
@property Pt dbl() const pure nothrow @nogc { |
@property Pt dbl() const pure nothrow @nogc @safe { |
||
if (this.isZero) |
if (this.isZero) |
||
return this; |
return this; |
||
| Line 147: | Line 147: | ||
} |
} |
||
string toString() const |
string toString() const pure /*nothrow*/ @safe { |
||
if (this.isZero) |
if (this.isZero) |
||
return "Zero"; |
return "Zero"; |
||
| Line 154: | Line 154: | ||
} |
} |
||
Pt opUnary(string op)() const pure nothrow @nogc |
Pt opUnary(string op)() const pure nothrow @nogc @safe |
||
if (op == "-") { |
|||
return Pt(this.x, -this.y); |
return Pt(this.x, -this.y); |
||
} |
} |
||
Pt opBinary(string op)(in Pt q) const pure nothrow @nogc |
Pt opBinary(string op)(in Pt q) const pure nothrow @nogc @safe |
||
if (op == "+") { |
if (op == "+") { |
||
if (this.x == q.x && this.y == q.y) |
if (this.x == q.x && this.y == q.y) |
||
| Line 171: | Line 172: | ||
} |
} |
||
Pt opBinary(string op)(in uint n) const pure nothrow @nogc |
Pt opBinary(string op)(in uint n) const pure nothrow @nogc @safe |
||
if (op == "*") { |
if (op == "*") { |
||
auto r = Pt.zero; |
auto r = Pt.zero; |
||
| Line 184: | Line 185: | ||
} |
} |
||
void main() { |
void main() @safe { |
||
immutable a = Pt.fromY(1); |
immutable a = Pt.fromY(1); |
||
immutable b = Pt.fromY(2); |
immutable b = Pt.fromY(2); |
||
Revision as of 16:19, 5 January 2015
Elliptic curves are sometimes used in cryptography as a way to perform digital signatures. The purpose of this task is to implement a simplified (without modular arithmetic) version of the elliptic curve arithmetic which is required by the elliptic curve DSA protocol.
In a nutshell, an elliptic curve is a bi-dimensional curve defined by the following relation between the x and y coordinates of any point on the curve:
a and b are arbitrary parameters that define the specific curve which is used. For this particular task, we'll use the following parameters:
a=0, b=7
The most interesting thing about elliptic curves is the fact that it is possible to define a group structure on it. To do so we define an internal composition rule with an additive notation +, such that for any three distinct points P, Q and R on the curve, whenever these points are aligned, we have:
P + Q + R = 0
Here 0 is the infinity point, for which the x and y values are not defined. It's basically the same kind of point which defines the horizon in projective geometry. We'll also assume here that this infinity point is unique and defines the neutral element of the addition.
This was not the definition of the addition, but only its desired property. For a more accurate definition, we proceed as such:
Given any three aligned points P, Q and R, we define the sum S = P + Q as the point (possibly the infinity point) such that S, R and the infinity point are aligned.
Considering the symmetry of the curve around the x-axis, it's easy to convince oneself that two points S and R can be aligned with the infinity point if and only if S and R are symmetric of one another towards the x-axis (because in that case there is no other candidate than the infinity point to complete the alignment triplet).
S is thus defined as the symmetric of R towards the x axis.
The task consists in defining the addition which, for any two points of the curve, returns the sum of these two points. You will pick two random points on the curve, compute their sum and show that the symmetric of the sum is aligned with the two initial points.
You will use the a and b parameters of secp256k1, i.e. respectively zero and seven.
Hint: You might need to define a "doubling" function, that returns P+P for any given point P.
Extra credit: define the full elliptic curve arithmetic (still not modular, though) by defining a "multiply" function that returns, for any point P and integer n, the point P + P + ... + P (n times).
C
<lang c>#include <stdio.h>
- include <math.h>
- define C 7
typedef struct { double x, y; } pt;
pt zero(void) { return (pt){ INFINITY, INFINITY }; }
// should be INFINITY, but numeric precision is very much in the way int is_zero(pt p) { return p.x > 1e20 || p.x < -1e20; }
pt neg(pt p) { return (pt){ p.x, -p.y }; }
pt dbl(pt p) { if (is_zero(p)) return p;
pt r; double L = (3 * p.x * p.x) / (2 * p.y); r.x = L * L - 2 * p.x; r.y = L * (p.x - r.x) - p.y; return r; }
pt add(pt p, pt q) { if (p.x == q.x && p.y == q.y) return dbl(p); if (is_zero(p)) return q; if (is_zero(q)) return p;
pt r; double L = (q.y - p.y) / (q.x - p.x); r.x = L * L - p.x - q.x; r.y = L * (p.x - r.x) - p.y; return r; }
pt mul(pt p, int n) { int i; pt r = zero();
for (i = 1; i <= n; i <<= 1) { if (i & n) r = add(r, p); p = dbl(p); } return r; }
void show(const char *s, pt p) { printf("%s", s); printf(is_zero(p) ? "Zero\n" : "(%.3f, %.3f)\n", p.x, p.y); }
pt from_y(double y) { pt r; r.x = pow(y * y - C, 1.0/3); r.y = y; return r; }
int main(void) { pt a, b, c, d;
a = from_y(1); b = from_y(2);
show("a = ", a); show("b = ", b); show("c = a + b = ", c = add(a, b)); show("d = -c = ", d = neg(c)); show("c + d = ", add(c, d)); show("a + b + d = ", add(a, add(b, d))); show("a * 12345 = ", mul(a, 12345));
return 0; }</lang>
- Output:
a = (-1.817, 1.000) b = (-1.442, 2.000) c = a + b = (10.375, -33.525) d = -c = (10.375, 33.525) c + d = Zero a + b + d = Zero a * 12345 = (10.759, 35.387)
D
<lang d>import std.stdio, std.math, std.string;
enum bCoeff = 7;
struct Pt {
double x, y;
@property static Pt zero() pure nothrow @nogc @safe {
return Pt(double.infinity, double.infinity);
}
@property bool isZero() const pure nothrow @nogc @safe {
return x > 1e20 || x < -1e20;
}
@property static Pt fromY(in double y) nothrow /*pure*/ @nogc @safe {
return Pt(cbrt(y ^^ 2 - bCoeff), y);
}
@property Pt dbl() const pure nothrow @nogc @safe {
if (this.isZero)
return this;
immutable L = (3 * x * x) / (2 * y);
immutable x2 = L ^^ 2 - 2 * x;
return Pt(x2, L * (x - x2) - y);
}
string toString() const pure /*nothrow*/ @safe {
if (this.isZero)
return "Zero";
else
return format("(%.3f, %.3f)", this.tupleof);
}
Pt opUnary(string op)() const pure nothrow @nogc @safe
if (op == "-") {
return Pt(this.x, -this.y);
}
Pt opBinary(string op)(in Pt q) const pure nothrow @nogc @safe
if (op == "+") {
if (this.x == q.x && this.y == q.y)
return this.dbl;
if (this.isZero)
return q;
if (q.isZero)
return this;
immutable L = (q.y - this.y) / (q.x - this.x);
immutable x = L ^^ 2 - this.x - q.x;
return Pt(x, L * (this.x - x) - this.y);
}
Pt opBinary(string op)(in uint n) const pure nothrow @nogc @safe
if (op == "*") {
auto r = Pt.zero;
Pt p = this;
for (uint i = 1; i <= n; i <<= 1) {
if ((i & n) != 0)
r = r + p;
p = p.dbl;
}
return r;
}
}
void main() @safe {
immutable a = Pt.fromY(1);
immutable b = Pt.fromY(2);
writeln("a = ", a);
writeln("b = ", b);
immutable c = a + b;
writeln("c = a + b = ", c);
immutable d = -c;
writeln("d = -c = ", d);
writeln("c + d = ", c + d);
writeln("a + b + d = ", a + b + d);
writeln("a * 12345 = ", a * 12345);
}</lang>
- Output:
a = (-1.817, 1.000) b = (-1.442, 2.000) c = a + b = (10.375, -33.525) d = -c = (10.375, 33.525) c + d = Zero a + b + d = Zero a * 12345 = (10.759, 35.387)
Go
<lang go>package main
import (
"fmt" "math"
)
const bCoeff = 7
type pt struct{ x, y float64 }
func zero() pt {
return pt{math.Inf(1), math.Inf(1)}
}
func is_zero(p pt) bool {
return p.x > 1e20 || p.x < -1e20
}
func neg(p pt) pt {
return pt{p.x, -p.y}
}
func dbl(p pt) pt {
if is_zero(p) {
return p
}
L := (3 * p.x * p.x) / (2 * p.y)
x := L*L - 2*p.x
return pt{
x: x,
y: L*(p.x-x) - p.y,
}
}
func add(p, q pt) pt {
if p.x == q.x && p.y == q.y {
return dbl(p)
}
if is_zero(p) {
return q
}
if is_zero(q) {
return p
}
L := (q.y - p.y) / (q.x - p.x)
x := L*L - p.x - q.x
return pt{
x: x,
y: L*(p.x-x) - p.y,
}
}
func mul(p pt, n int) pt {
r := zero()
for i := 1; i <= n; i <<= 1 {
if i&n != 0 {
r = add(r, p)
}
p = dbl(p)
}
return r
}
func show(s string, p pt) {
fmt.Printf("%s", s)
if is_zero(p) {
fmt.Println("Zero")
} else {
fmt.Printf("(%.3f, %.3f)\n", p.x, p.y)
}
}
func from_y(y float64) pt {
return pt{
x: math.Cbrt(y*y - bCoeff),
y: y,
}
}
func main() {
a := from_y(1)
b := from_y(2)
show("a = ", a)
show("b = ", b)
c := add(a, b)
show("c = a + b = ", c)
d := neg(c)
show("d = -c = ", d)
show("c + d = ", add(c, d))
show("a + b + d = ", add(a, add(b, d)))
show("a * 12345 = ", mul(a, 12345))
}</lang>
- Output:
a = (-1.817, 1.000) b = (-1.442, 2.000) c = a + b = (10.375, -33.525) d = -c = (10.375, 33.525) c + d = Zero a + b + d = Zero a * 12345 = (10.759, 35.387)
J
Follows the C contribution.
<lang j>C=. 7
zero=: _j_
isZero=: 1e20 < |@{.@+.
neg=: 1 _1&*&.+.
dbl=: 3 :0 if. isZero y do. y return. end. 'px py'=. +. y r j. py-~L*px-r=. (2*px)-~*~L=. (3*px*px)% 2*py )
add=: 4 :0 if. x=y do. dbl x return. end. if. isZero x do. y return. end. if. isZero y do. x return. end. 'px py'=. +. x 'qx qy'=. +. y r j. py-~ L * px-r=. (px+qx)-~*~L=. (qy-py)%qx-px )
mul=: [:{.[:;[:(,dbl)/@]`((add,dbl@])/@])@.[&.>/([:<"0 #:@]),[:< zero,[
from=: j.~[:(**3&%:@|)C-~*~
show =: (,('( ',[,' , ',' )',~])&":/@+.@])`('Zero',~[)@.(isZero@])
task=: 3 :0
a=. from 1
b=. from 2
smoutput 'a = ' show a smoutput 'b = ' show b smoutput 'c = a + b = ' show c =. a add b smoutput 'd = -c = ' show d =. neg c smoutput 'c + d = ' show c add d smoutput 'a + b + d = ' show add/ a, b, d smoutput 'a * 12345 = ' show a mul 12345 )</lang>
- Output:
<lang j> task 'uit' a = ( _1.81712 , 1 ) b = ( _1.44225 , 2 ) c = a + b = ( 10.3754 , _33.5245 ) d = -c = ( 10.3754 , 33.5245 ) c + d = Zero a + b + d = Zero a * 12345 = ( 10.7586 , 35.3875 )</lang>
PARI/GP
The examples were borrowed from C, though the coding is built-in for GP and so not ported. <lang parigp>e=ellinit([0,7]); a=[-6^(1/3),1] b=[-3^(1/3),2] c=elladd(e,a,b) d=ellneg(e,c) elladd(e,c,d) elladd(e,elladd(e,a,b),d) ellmul(e,a,12345)</lang>
- Output:
%1 = [-1.8171205928321396588912117563272605024, 1] %2 = [-1.4422495703074083823216383107801095884, 2] %3 = [10.375375389201411959219947350723254093, -33.524509096269714974732957666465317961] %4 = [10.375375389201411959219947350723254093, 33.524509096269714974732957666465317961] %5 = [0] %6 = [0] %7 = [10.758570529079026647817660298097473136, 35.387434773095871032744887640370612568]
Perl
<lang perl>package EC; our ($a, $b) = (0, 7); {
package EC::Point;
sub new { my $class = shift; bless [ @_ ], $class }
sub zero { bless [], shift }
sub x { shift->[0] }; sub y { shift->[1] };
sub double {
my $self = shift;
return $self unless @$self;
my $L = (3 * $self->x**2) / (2*$self->y);
my $x = $L**2 - 2*$self->x;
bless [ $x, $L * ($self->x - $x) - $self->y ], ref $self;
}
use overload
'==' => sub { my ($p, $q) = @_; $p->x == $q->x and $p->y == $q->y },
'+' => sub {
my ($p, $q) = @_;
return $p->double if $p == $q;
return $p unless @$q;
return $q unless @$p;
my $slope = ($q->y - $p->y) / ($q->x - $p->x);
my $x = $slope**2 - $p->x - $q->x;
bless [ $x, $slope * ($p->x - $x) - $p->y ], ref $p;
},
q{""} => sub {
my $self = shift;
return @$self
? sprintf "EC-point at x=%f, y=%f", @$self
: 'EC point at infinite';
}
}
package Test; my $a = +EC::Point->new(-($EC::b - 1)**(1/3), 1); my $b = +EC::Point->new(-($EC::b - 4)**(1/3), 2); print my $c = $a + $b, "\n"; print "$_\n" for $a, $b, $c; print "check alignement... "; print abs(($b->x - $a->x)*(-$c->y - $a->y) - ($b->y - $a->y)*($c->x - $a->x)) < 0.001
? "ok" : "wrong";
</lang>
- Output:
EC-point at x=-1.817121, y=1.000000 EC-point at x=-1.442250, y=2.000000 EC-point at x=10.375375, y=-33.524509 check alignement... ok
Perl 6
<lang perl6>module EC; our ($A, $B) = (0, 7);
role Horizon { method gist { 'EC Point at horizon' } } class Point {
has ($.x, $.y);
multi method new(
$x, $y where $y**2 ~~ $x**3 + $A*$x + $B
) { self.bless(:$x, :$y) }
multi method new(Horizon $) { self.bless but Horizon }
method gist { "EC Point at x=$.x, y=$.y" }
}
multi prefix:<->(Point $p) { Point.new: x => $p.x, y => -$p.y } multi prefix:<->(Horizon $) { Horizon } multi infix:<->(Point $a, Point $b) { $a + -$b }
multi infix:<+>(Horizon $, Point $p) { $p } multi infix:<+>(Point $p, Horizon) { $p }
multi infix:<*>(Point $u, Int $n) { $n * $u } multi infix:<*>(Int $n, Horizon) { Horizon } multi infix:<*>(0, Point) { Horizon } multi infix:<*>(1, Point $p) { $p } multi infix:<*>(2, Point $p) {
my $l = (3*$p.x**2 + $A) / (2 *$p.y); my $y = $l*($p.x - my $x = $l**2 - 2*$p.x) - $p.y; $p.bless(:$x, :$y);
} multi infix:<*>(Int $n where $n > 2, Point $p) {
2 * ($n div 2 * $p) + $n % 2 * $p;
}
multi infix:<+>(Point $p, Point $q) {
if $p.x ~~ $q.x {
return $p.y ~~ $q.y ?? 2 * $p !! Horizon;
}
else {
my $slope = ($q.y - $p.y) / ($q.x - $p.x);
my $y = $slope*($p.x - my $x = $slope**2 - $p.x - $q.x) - $p.y;
return $p.new(:$x, :$y);
}
}
say my $p = Point.new: x => $_, y => sqrt(abs(1 - $_**3 - $A*$_ - $B)) given 1; say my $q = Point.new: x => $_, y => sqrt(abs(1 - $_**3 - $A*$_ - $B)) given 2; say my $s = $p + $q;
say "checking alignment: ", abs ($p.x - $q.x)*(-$s.y - $q.y) - ($p.y - $q.y)*($s.x - $q.x);</lang>
- Output:
EC Point at x=1, y=2.64575131106459 EC Point at x=2, y=3.74165738677394 EC Point at x=-1.79898987322333, y=0.421678696849803 checking alignment: 8.88178419700125e-16
Racket
<lang racket>
- lang racket
(define a 0) (define b 7) (define (ε? x) (<= (abs x) 1e-14)) (define (== p q) (for/and ([pi p] [qi q]) (ε? (- pi qi)))) (define zero #(0 0)) (define (zero? p) (== p zero)) (define (neg p) (match-define (vector x y) p) (vector x (- y))) (define (⊕ p q)
(cond [(== q (neg p)) zero]
[else
(match-define (vector px py) p)
(match-define (vector qx qy) q)
(define (done λ px py qx)
(define x (- (* λ λ) px qx))
(vector x (- (+ (* λ (- x px)) py))))
(cond [(and (== p q) (ε? py)) zero]
[(or (== p q) (ε? (- px qx)))
(done (/ (+ (* 3 px px) a) (* 2 py)) px py qx)]
[(done (/ (- py qy) (- px qx)) px py qx)])]))
(define (⊗ p n)
(cond [(= n 0) zero]
[(= n 1) p]
[(= n 2) (⊕ p p)]
[(negative? n) (neg (⊗ p (- n)))]
[(even? n) (⊗ (⊗ p (/ n 2)) 2)]
[(odd? n) (⊕ p (⊗ p (- n 1)))]))
</lang> Test: <lang racket> (define (root3 x) (* (sgn x) (expt (abs x) 1/3))) (define (y->point y) (vector (root3 (- (* y y) b)) y)) (define p (y->point 1)) (define q (y->point 2)) (displayln (~a "p = " p)) (displayln (~a "q = " q)) (displayln (~a "p+q = " (⊕ p q))) (displayln (~a "-(p+q) = " (neg (⊕ p q)))) (displayln (~a "(p+q)+(-(p+q)) = " (⊕ (⊕ p q) (neg (⊕ p q))))) (displayln (~a "p+(q+(-(p+q))) = 0 " (zero? (⊕ p (⊕ q (neg (⊕ p q))))))) (displayln (~a "p*12345 " (⊗ p 12345))) </lang> Output: <lang racket> p = #(-1.8171205928321397 1) q = #(-1.4422495703074083 2) p+q = #(10.375375389201409 -33.524509096269696) -(p+q) = #(10.375375389201409 33.524509096269696) (p+q)+(-(p+q)) = #(0 0) p+(q+(-(p+q))) = 0 #t p*12345 #(10.758570529320806 35.387434774282106) </lang>
Tcl
<lang tcl>set C 7 set zero {x inf y inf} proc tcl::mathfunc::cuberoot n {
# General power operator doesn't like negative, but its defined for root3
expr {$n>=0 ? $n**(1./3) : -((-$n)**(1./3))}
} proc iszero p {
dict with p {}
return [expr {$x > 1e20 || $x<-1e20}]
} proc negate p {
dict set p y [expr {-[dict get $p y]}]
} proc double p {
if {[iszero $p]} {return $p}
dict with p {}
set L [expr {(3.0 * $x**2) / (2.0 * $y)}]
set rx [expr {$L**2 - 2.0 * $x}]
set ry [expr {$L * ($x - $rx) - $y}]
return [dict create x $rx y $ry]
} proc add {p q} {
if {[dict get $p x]==[dict get $q x] && [dict get $p y]==[dict get $q y]} {
return [double $p]
}
if {[iszero $p]} {return $q}
if {[iszero $q]} {return $p}
dict with p {}
set L [expr {([dict get $q y]-$y) / ([dict get $q x]-$x)}]
dict set r x [expr {$L**2 - $x - [dict get $q x]}]
dict set r y [expr {$L * ($x - [dict get $r x]) - $y}]
return $r
} proc multiply {p n} {
set r $::zero
for {set i 1} {$i <= $n} {incr i $i} {
if {$i & int($n)} { set r [add $r $p] } set p [double $p]
} return $r
}</lang> Demonstrating: <lang tcl>proc show {s p} {
if {[iszero $p]} {
puts "${s}Zero"
} else {
dict with p {} puts [format "%s(%.3f, %.3f)" $s $x $y]
}
} proc fromY y {
global C
dict set r x [expr {cuberoot($y**2 - $C)}]
dict set r y [expr {double($y)}]
}
set a [fromY 1] set b [fromY 2] show "a = " $a show "b = " $b show "c = a + b = " [set c [add $a $b]] show "d = -c = " [set d [negate $c]] show "c + d = " [add $c $d] show "a + b + d = " [add $a [add $b $d]] show "a * 12345 = " [multiply $a 12345]</lang>
- Output:
a = (-1.817, 1.000) b = (-1.442, 2.000) c = a + b = (10.375, -33.525) d = -c = (10.375, 33.525) c + d = Zero a + b + d = Zero a * 12345 = (10.759, 35.387)