;; -*- Mode: Irken -*- ;; goal: signatures using nist-p256. ;; bonus: other curves p384, p521, secp256k1 (require "lib/basis.scm") (require "lib/bignum.scm") (require "lib/codecs/hex.scm") (require "lib/time.scm") ;; affine (datatype ecp (:inf) (:point big big) ) ;; homogenous (datatype ecph (:point big big big) ) (define ecp-inf? (ecp:inf) -> #t _ -> #f ) ;; technically X=0 Y!=0 Z=0, but I think in reality ;; you only need to check Z=0? (define ecph-inf? (ecph:point (big:zero) Y (big:zero)) -> #t _ -> #f ) (define point-repr (ecp:inf) -> "Inf" ;;(ecp:point x y) -> (format "(" (big->dec x) ", " (big->dec y) ")") (ecp:point x y) -> (format "(" (big-repr x) ", " (big-repr y) ")") ) ;; point-at-infinity is where x=0 and z=0 and y != 0 (define ecph-repr (ecph:point X Y Z) -> (format "(" (big->dec X) " : " (big->dec Y) " : " (big->dec Z) ")") ) (define (hex->big s) (u256->big (hex->string s))) (define secp256r1 {p = (hex->big "ffffffff00000001000000000000000000000000ffffffffffffffffffffffff") a = (hex->big "ffffffff00000001000000000000000000000000fffffffffffffffffffffffc") b = (hex->big "5ac635d8aa3a93e7b3ebbd55769886bc651d06b0cc53b0f63bce3c3e27d2604b") ;; hey rocky! g = (ecp:point (hex->big "6b17d1f2e12c4247f8bce6e563a440f277037d812deb33a0f4a13945d898c296") (hex->big "4fe342e2fe1a7f9b8ee7eb4a7c0f9e162bce33576b315ececbb6406837bf51f5")) n = (hex->big "ffffffff00000000ffffffffffffffffbce6faada7179e84f3b9cac2fc632551") h = big/1 name = "secp256r1" }) (define secp384r1 {p = (hex->big "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffeffffffff0000000000000000ffffffff") a = (hex->big "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffeffffffff0000000000000000fffffffc") b = (hex->big "b3312fa7e23ee7e4988e056be3f82d19181d9c6efe8141120314088f5013875ac656398d8a2ed19d2a85c8edd3ec2aef") g = (ecp:point (hex->big "aa87ca22be8b05378eb1c71ef320ad746e1d3b628ba79b9859f741e082542a385502f25dbf55296c3a545e3872760ab7") (hex->big "3617de4a96262c6f5d9e98bf9292dc29f8f41dbd289a147ce9da3113b5f0b8c00a60b1ce1d7e819d7a431d7c90ea0e5f")) n = (hex->big "ffffffffffffffffffffffffffffffffffffffffffffffffc7634d81f4372ddf581a0db248b0a77aecec196accc52973") h = big/1 name = "secp384r1" }) (define secp521r1 {p = (hex->big "01ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff") a = (hex->big "01fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffc") b = (hex->big "0051953eb9618e1c9a1f929a21a0b68540eea2da725b99b315f3b8b489918ef109e156193951ec7e937b1652c0bd3bb1bf073573df883d2c34f1ef451fd46b503f00") g = (ecp:point (hex->big "00c6858e06b70404e9cd9e3ecb662395b4429c648139053fb521f828af606b4d3dbaa14b5e77efe75928fe1dc127a2ffa8de3348b3c1856a429bf97e7e31c2e5bd66") (hex->big "011839296a789a3bc0045c8a5fb42c7d1bd998f54449579b446817afbd17273e662c97ee72995ef42640c550b9013fad0761353c7086a272c24088be94769fd16650")) n = (hex->big "01fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffa51868783bf2f966b7fcc0148f709a5d03bb5c9b8899c47aebb6fb71e91386409") h = big/1 name = "secp521r1" }) (define secp256k1 {p = (hex->big "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f") a = big/0 b = big/7 g = (ecp:point (hex->big "79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798") (hex->big "483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8")) n = (hex->big "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141") h = big/1 name = "secp256k1" }) (define demo0 {p = (big (I 97)) a = big/2 b = big/3 g = (ecp:point big/7 big/9) ;; fake n = (big (I 100)) h = big/1 name = "demo0" }) (define demo1 {p = (big (I 43)) a = (big (I -3)) b = (big (I 27)) g = (ecp:point big/1 (big (I 38))) n = (big (I 13)) ;; making this up h = big/1 name = "demo1" }) (define egcd (big:zero) b -> (:tuple b big/0 big/1) a b -> (let (((q r) (big-div b a)) ((g y x) (egcd r a))) (:tuple g (big (- x (* q y))) y)) ) (define (mod-inv a p) (if (big-neg? a) (big-sub p (mod-inv (big-negate a) p)) (let (((g x y) (egcd a p))) (if (not (big= g big/1)) (raise (:NoInverse a)) (big-mod x p))))) ;; using fermat's little theorem - requires prime p ;; [note: a little slower] ;; (define (mod-inv a p) ;; (big-exp-mod a (big-sub p big/2) p) ;; ) (define (make-curve c) (define is-on-curve? (ecp:inf) -> #t (ecp:point x y) -> (big (= (mod (* y y) c.p) (mod (+ (* x x x) (* c.a x) c.b) c.p)))) (define point-negate (ecp:inf) -> (ecp:inf) (ecp:point x y) -> (ecp:point x (big (mod (- y) c.p))) ) ;; affine (define point-add (ecp:inf) b -> b a (ecp:inf) -> a (ecp:point x1 y1) (ecp:point x2 y2) -> (let ((x=? (big= x1 x2))) (if (and x=? (not (big= y1 y2))) (ecp:inf) (let ((m (if x=? (big (* (+ (* (I 3) x1 x1) c.a) (mod-inv (big-lshift y1 1) c.p))) (big (* (- y1 y2) (mod-inv (big-sub x1 x2) c.p)))))) (let ((x3 (big (- (* m m) x1 x2))) (y3 (big (+ y1 (* m (- x3 x1)))))) (ecp:point (big (mod x3 c.p)) (big (mod (- y3) c.p)))) ) )) ) ;; affine (define (scalar-mult k P) (if (or (ecp-inf? P) (big (= big/0 (mod k c.n)))) (ecp:inf) (if (big-neg? k) (scalar-mult (big-negate k) (point-negate P)) (let ((r (ecp:inf)) (addend P)) (while (not (big-zero? k)) (when (big-odd? k) (set! r (point-add r addend))) (set! addend (point-add addend addend)) (set! k (big-rshift k 1))) r )))) (define (to-bits n) (let loop ((r '()) (n n)) (if (big-zero? n) r (loop (list:cons (big-odd? n) r) (big-rshift n 1))))) (define (monty k P) (let loop ((r0 (ecp:inf)) (r1 P) (bits (to-bits k))) (match bits with () -> r0 (#t . tl) -> (loop (point-add r0 r1) (point-add r1 r1) tl) (#f . tl) -> (loop (point-add r0 r0) (point-add r0 r1) tl) ))) (define to-homogenous (ecp:inf) -> (ecph:point big/0 big/1 big/0) (ecp:point x y) -> (ecph:point x y big/1) ) (define to-affine (ecph:point X Y (big:zero)) -> (ecp:inf) (ecph:point X Y Z) -> (let ((Zneg (mod-inv Z c.p))) (ecp:point (big (mod (* X Zneg) c.p)) (big (mod (* Y Zneg) c.p))))) ;; https://www.hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#addition-add-2015-rcb (define (add-2015-rcb X1 Y1 Z1 X2 Y2 Z2) (define (M* a b) (big-mod (big-mul a b) c.p)) (define (M+ a b) (big-mod (big-add a b) c.p)) (define (M- a b) (big-mod (big-sub a b) c.p)) (let ((t0 big/0) (t1 big/0) (t2 big/0) (t3 big/0) (t4 big/0) (t5 big/0) (X3 big/0) (Y3 big/0) (Z3 big/0) (b3 (M* big/3 c.b))) ;; note: should be pre-computed. (set! t0 (M* X1 X2)) ;; t0 = X1*X2 (set! t1 (M* Y1 Y2)) ;; t1 = Y1*Y2 (set! t2 (M* Z1 Z2)) ;; t2 = Z1*Z2 (set! t3 (M+ X1 Y1)) ;; t3 = X1+Y1 (set! t4 (M+ X2 Y2)) ;; t4 = X2+Y2 (set! t3 (M* t3 t4)) ;; t3 = t3*t4 (set! t4 (M+ t0 t1)) ;; t4 = t0+t1 (set! t3 (M- t3 t4)) ;; t3 = t3-t4 (set! t4 (M+ X1 Z1)) ;; t4 = X1+Z1 (set! t5 (M+ X2 Z2)) ;; t5 = X2+Z2 (set! t4 (M* t4 t5)) ;; t4 = t4*t5 (set! t5 (M+ t0 t2)) ;; t5 = t0+t2 (set! t4 (M- t4 t5)) ;; t4 = t4-t5 (set! t5 (M+ Y1 Z1)) ;; t5 = Y1+Z1 (set! X3 (M+ Y2 Z2)) ;; X3 = Y2+Z2 (set! t5 (M* t5 X3)) ;; t5 = t5*X3 (set! X3 (M+ t1 t2)) ;; X3 = t1+t2 (set! t5 (M- t5 X3)) ;; t5 = t5-X3 (set! Z3 (M* c.a t4)) ;; Z3 = a*t4 (set! X3 (M* b3 t2)) ;; X3 = b3*t2 (set! Z3 (M+ X3 Z3)) ;; Z3 = X3+Z3 (set! X3 (M- t1 Z3)) ;; X3 = t1-Z3 (set! Z3 (M+ t1 Z3)) ;; Z3 = t1+Z3 (set! Y3 (M* X3 Z3)) ;; Y3 = X3*Z3 (set! t1 (M+ t0 t0)) ;; t1 = t0+t0 (set! t1 (M+ t1 t0)) ;; t1 = t1+t0 (set! t2 (M* c.a t2)) ;; t2 = a*t2 (set! t4 (M* b3 t4)) ;; t4 = b3*t4 (set! t1 (M+ t1 t2)) ;; t1 = t1+t2 (set! t2 (M- t0 t2)) ;; t2 = t0-t2 (set! t2 (M* c.a t2)) ;; t2 = a*t2 (set! t4 (M+ t4 t2)) ;; t4 = t4+t2 (set! t0 (M* t1 t4)) ;; t0 = t1*t4 (set! Y3 (M+ Y3 t0)) ;; Y3 = Y3+t0 (set! t0 (M* t5 t4)) ;; t0 = t5*t4 (set! X3 (M* t3 X3)) ;; X3 = t3*X3 (set! X3 (M- X3 t0)) ;; X3 = X3-t0 (set! t0 (M* t3 t1)) ;; t0 = t3*t1 (set! Z3 (M* t5 Z3)) ;; Z3 = t5*Z3 (set! Z3 (M+ Z3 t0)) ;; Z3 = Z3+t0 (ecph:point X3 Y3 Z3) )) ;; https://www.hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#doubling-dbl-2015-rcb (define (dbl-2015-rcb X1 Y1 Z1) (define (M* a b) (big-mod (big-mul a b) c.p)) (define (M+ a b) (big-mod (big-add a b) c.p)) (define (M- a b) (big-mod (big-sub a b) c.p)) (let ((t0 big/0) (t1 big/0) (t2 big/0) (t3 big/0) (X3 big/0) (Y3 big/0) (Z3 big/0) (b3 (M* big/3 c.b))) ;; note: should be pre-computed. (set! t0 (M* X1 X1)) ;; t0 = X1^2 (set! t1 (M* Y1 Y1)) ;; t1 = Y1^2 (set! t2 (M* Z1 Z1)) ;; t2 = Z1^2 (set! t3 (M* X1 Y1)) ;; t3 = X1*Y1 (set! t3 (M+ t3 t3)) ;; t3 = t3+t3 (set! Z3 (M* X1 Z1)) ;; Z3 = X1*Z1 (set! Z3 (M+ Z3 Z3)) ;; Z3 = Z3+Z3 (set! X3 (M* c.a Z3)) ;; X3 = a*Z3 (set! Y3 (M* b3 t2)) ;; Y3 = b3*t2 (set! Y3 (M+ X3 Y3)) ;; Y3 = X3+Y3 (set! X3 (M- t1 Y3)) ;; X3 = t1-Y3 (set! Y3 (M+ t1 Y3)) ;; Y3 = t1+Y3 (set! Y3 (M* X3 Y3)) ;; Y3 = X3*Y3 (set! X3 (M* t3 X3)) ;; X3 = t3*X3 (set! Z3 (M* b3 Z3)) ;; Z3 = b3*Z3 (set! t2 (M* c.a t2)) ;; t2 = a*t2 (set! t3 (M- t0 t2)) ;; t3 = t0-t2 (set! t3 (M* c.a t3)) ;; t3 = a*t3 (set! t3 (M+ t3 Z3)) ;; t3 = t3+Z3 (set! Z3 (M+ t0 t0)) ;; Z3 = t0+t0 (set! t0 (M+ Z3 t0)) ;; t0 = Z3+t0 (set! t0 (M+ t0 t2)) ;; t0 = t0+t2 (set! t0 (M* t0 t3)) ;; t0 = t0*t3 (set! Y3 (M+ Y3 t0)) ;; Y3 = Y3+t0 (set! t2 (M* Y1 Z1)) ;; t2 = Y1*Z1 (set! t2 (M+ t2 t2)) ;; t2 = t2+t2 (set! t0 (M* t2 t3)) ;; t0 = t2*t3 (set! X3 (M- X3 t0)) ;; X3 = X3-t0 (set! Z3 (M* t2 t1)) ;; Z3 = t2*t1 (set! Z3 (M+ Z3 Z3)) ;; Z3 = Z3+Z3 (set! Z3 (M+ Z3 Z3)) ;; Z3 = Z3+Z3 (ecph:point X3 Y3 Z3) )) ;; XXX according to the paper (https://eprint.iacr.org/2015/1060.pdf), the formulas ;; are complete and do NOT require checking for zero/equal/etc. (define proj-add ;; (ecph:point _ _ (big:zero)) P2 -> P2 ;; P1 (ecph:point _ _ (big:zero)) -> P1 (ecph:point X1 Y1 Z1) (ecph:point X2 Y2 Z2) -> (add-2015-rcb X1 Y1 Z1 X2 Y2 Z2) ) (define proj-dbl (ecph:point X Y Z) -> (dbl-2015-rcb X Y Z) ) (define proj-scalar-mult k (ecp:point X Y) -> (let ((r (ecph:point big/0 big/1 big/0)) ;; inf (addend (ecph:point X Y big/1))) (while (not (big-zero? k)) (when (big-odd? k) (set! r (proj-add r addend))) (set! addend (proj-dbl addend)) (set! k (big-rshift k 1))) r) k (ecp:inf) -> (ecph:point big/0 big/1 big/0) ;; XXX yes? ) (define proj-monty k (ecp:point X Y) -> (let loop ((r0 (ecph:point big/0 big/1 big/0)) ;; inf (r1 (ecph:point X Y big/1)) (bits (to-bits k))) (match bits with () -> r0 (#t . tl) -> (loop (proj-add r0 r1) (proj-dbl r1) tl) (#f . tl) -> (loop (proj-dbl r0) (proj-add r0 r1) tl) )) k (ecp:inf) -> (ecph:point big/0 big/1 big/0) ;; XXX yes? ) (define (homogenous-add P1 P2) (cond ((ecph-inf? P1) P2) ((ecph-inf? P2) P1) (else (match P1 P2 with (ecph:point X1 Y1 Z1) (ecph:point X2 Y2 Z2) -> (let ((u (big (mod (- (* Y2 Z1) (* Y1 Z2)) c.p))) (v (big (mod (- (* X2 Z1) (* X1 Z2)) c.p))) (u0? (big-zero? u)) (v0? (big-zero? v))) (match u0? v0? with #f #t -> (ecph:point big/0 big/1 big/0) #f #f -> (let (;; v * (Z2 * (Z1 * u^2 - 2 * X1 * v^2) - v^3) (X3 (big (mod (* v (- (* Z2 (- (* Z1 u u) (* (I 2) X1 v v))) (* v v v))) c.p))) ;; Z2 * (3 * X1 * u * v^2 - Y1 * v^3 - Z1 * u^3) + u * v^3 (Y3 (big (mod (+ (* Z2 (- (* (I 3) X1 u v v) (* Y1 v v v) (* Z1 u u u))) (* u v v v)) c.p))) ;; v^3 * Z1 * Z2 (Z3 (big (mod (* v v v Z1 Z2) c.p)))) (ecph:point X3 Y3 Z3)) _ _ -> (let (;; 3 * X1^2 + a * Z1^2 (w (big (mod (+ (* (I 3) X1 X1) (* c.a Z1 Z1)) c.p))) ;; 2 * Y1 * Z1 * (w^2 - 8 * X1 * Y1^2 * Z1) (X3 (big (mod (* (I 2) Y1 Z1 (- (* w w) (* (I 8) X1 Y1 Y1 Z1))) c.p))) ;; 4 * Y1^2 * Z1 * (3 * w * X1 - 2 * Y1^2 * Z1) - w^3 (Y3 (big (mod (- (* (I 4) Y1 Y1 Z1 (- (* (I 3) w X1) (* (I 2) Y1 Y1 Z1))) (* w w w)) c.p))) ;; 8 * (Y1 * Z1)^3 (Z3 (let ((YZ (big (mod (* Y1 Z1) c.p)))) (big (mod (* (I 8) (* YZ YZ YZ)) c.p))))) (printf "w = " (big-repr w) "\n") (ecph:point X3 Y3 Z3)) )) )))) {g = c.g add = point-add neg = point-negate scalarmult = scalar-mult scalarmultp = proj-scalar-mult monty = monty montyp = proj-monty on = is-on-curve? toh = to-homogenous toa = to-affine addh = homogenous-add } ) (let ((c0 (make-curve secp256r1)) (c1 (make-curve demo0))) (printf "on curve? " (bool (c0.on c0.g)) "\n") (let ((P+Q (c1.add (ecp:point (big (I 17)) (big (I 10))) (ecp:point (big (I 95)) (big (I 31)))))) (printf "P+Q = " (point-repr P+Q) "\n")) (let ((P (ecp:point big/3 big/6)) (P50 (c1.scalarmult (big (I 50)) P)) (P33 (c1.scalarmult (big (I 33)) P))) (printf "50 x P = " (point-repr P50) "\n") (printf "33 X P = " (point-repr P33) "\n")) (let ((k (hex->big "cee1c9a2a70d095da48e36c1136b6809114df68f7cfa981c81653b9c27dfdb")) (kP (c0.scalarmult k c0.g))) (printf "affine-scalarmult\n") (printf "kP = " (point-repr kP) "\n") (printf "affine-monty\n") (printf "kP = " (point-repr (c0.monty k c0.g)) "\n") (let ((kP0 (c0.scalarmultp k c0.g))) (printf "proj-scalarmult\n") (printf "kP = " (ecph-repr kP0) "\n") (printf "kP = " (point-repr (c0.toa kP0)) "\n")) (let ((kP1 (c0.montyp k c0.g))) (printf "proj-monty\n") (printf "kP = " (ecph-repr kP1) "\n") (printf "kP = " (point-repr (c0.toa kP1)) "\n")) )) (let ((c0 (make-curve demo1)) (P (ecph:point big/1 (big (I 38)) big/1)) (P2 (c0.addh P P)) (P3 (c0.addh P2 P))) (printf "2P = " (ecph-repr P2) "\n" " = " (point-repr (c0.toa P2)) "\n" "3P = " (ecph-repr P3) "\n" " = " (point-repr (c0.toa P3)) "\n")) (let ((c0 (make-curve demo0))) (let ((P+Q (c0.add (ecp:point (big (I 17)) (big (I 10))) (ecp:point (big (I 95)) (big (I 31)))))) (printf "P+Q = " (point-repr P+Q) "\n") (let ((P (ecph:point (big (I 17)) (big (I 10)) big/1)) (Q (ecph:point (big (I 95)) (big (I 31)) big/1)) (P+Q2 (c0.addh P Q))) (printf "P+Q = " (ecph-repr P+Q2) "\n" " = " (point-repr (c0.toa P+Q2)) "\n") ))) (define (get-nano) (timespec->nanoseconds (clock-gettime CLOCK_REALTIME))) (defmacro timeit (timeit body ...) -> (let (($t0 (get-nano)) ($result (begin body ...)) ($delta (- (get-nano) $t0))) (printf " time: " (int $delta) "\n") $result)) ;; get some timings (let ((c0 (make-curve secp256r1)) (k (hex->big "cee1c9a2a70d095da48e36c1136b6809114df68f7cfa981c81653b9c27dfdb"))) (set-verbose-gc #f) ;; annoys the output (printf "scalarmult: ") (timeit (c0.scalarmult k c0.g)) (printf "monty: ") (timeit (c0.monty k c0.g)) (printf "proj-scalar: ") (timeit (c0.scalarmultp k c0.g)) (printf "proj-monty: ") (timeit (c0.montyp k c0.g)) (set-verbose-gc #t) ) (define (test-curve c k ex ey) (let ((c0 (make-curve c)) (kG0 (c0.scalarmult k c0.g)) (kG1 (c0.toa (c0.montyp k c0.g))) ) (printf "testing curve " c.name "\n") (printf "kG0 = " (point-repr kG0) "\n") (printf "kG1 = " (point-repr kG1) "\n") (printf "equal? " (bool (magic=? kG0 (ecp:point ex ey))) "\n") (printf "equal? " (bool (magic=? kG1 (ecp:point ex ey))) "\n") )) (test-curve secp384r1 (hex->big "099f3c7034d4a2c699884d73a375a67f7624ef7c6b3c0f160647b67414dce655e35b538041e649ee3faef896783ab194") (hex->big "667842d7d180ac2cde6f74f37551f55755c7645c20ef73e31634fe72b4c55ee6de3ac808acb4bdb4c88732aee95f41aa") (hex->big "9482ed1fc0eeb9cafc4984625ccfc23f65032149e0e144ada024181535a0f38eeb9fcff3c2c947dae69b4c634573a81c") ) (test-curve secp521r1 (hex->big "0037ade9319a89f4dabdb3ef411aaccca5123c61acab57b5393dce47608172a095aa85a30fe1c2952c6771d937ba9777f5957b2639bab072462f68c27a57382d4a52") (hex->big "0015417e84dbf28c0ad3c278713349dc7df153c897a1891bd98bab4357c9ecbee1e3bf42e00b8e380aeae57c2d107564941885942af5a7f4601723c4195d176ced3e") (hex->big "017cae20b6641d2eeb695786d8c946146239d099e18e1d5a514c739d7cb4a10ad8a788015ac405d7799dc75e7b7d5b6cf2261a6a7f1507438bf01beb6ca3926f9582") ) (test-curve secp256k1 (hex->big "e9873d79c6d87dc0fb6a5778633389f4453213303da61f20bd67fc233aa33262") (dec->big "40052878126280527701260741223305245603564636128202744842713277751919610658249") (dec->big "112427920116541844817408230468149218341228927370925731589596315545721129686052"))