Class: EllipticCurve::Math

Inherits:
Object
  • Object
show all
Defined in:
lib/math.rb

Class Method Summary collapse

Class Method Details

._fromJacobian(p, prime) ⇒ Object 



271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
 
# File 'lib/math.rb', line 271

def self._fromJacobian(p, prime)
    # Convert point back from Jacobian coordinates
    #
    # :param p: First Point you want to add
    # :param prime: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    # :return: Point in default coordinates
    if p.y == 0
        return Point.new(0, 0, 0)
    end

    z = self.inv(p.z, prime)
    x = (p.x * z ** 2) % prime
    y = (p.y * z ** 3) % prime

    return Point.new(x, y, 0)
end

._generatorPowersTable(curve) ⇒ Object 



93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
 
# File 'lib/math.rb', line 93

def self._generatorPowersTable(curve)
    # Build [G, 2G, 4G, ..., 2^nBitLength * G] in affine (z=1) form, so each
    # add in multiplyGenerator hits the mixed-add fast path.
    return curve._generatorPowersTable unless curve._generatorPowersTable.nil?
    coeff = curve.a
    prime = curve.p
    current = Point.new(curve.g.x, curve.g.y, 1)
    table = [current]
    # NAF of an nBitLength-bit scalar can be up to nBitLength+1 digits.
    curve.nBitLength.times do
        doubled = self._jacobianDouble(current, coeff, prime)
        if doubled.y == 0
            current = doubled
        else
            zInv = self.inv(doubled.z, prime)
            zInv2 = (zInv * zInv) % prime
            zInv3 = (zInv2 * zInv) % prime
            current = Point.new((doubled.x * zInv2) % prime, (doubled.y * zInv3) % prime, 1)
        end
        table << current
    end
    curve._generatorPowersTable = table
    return table
end

._glvDecompose(k, glv, order) ⇒ Object 



227
228
229
230
231
232
233
234
235
236
237
238
 
# File 'lib/math.rb', line 227

def self._glvDecompose(k, glv, order)
    # Decompose k into (k1, k2) with k == k1 + k2*lambda (mod N) and
    # |k1|, |k2| ~ sqrt(N). Babai rounding against the precomputed basis
    # {(a1, b1), (a2, b2)}; k1 and k2 may be negative.
    a1, b1, a2, b2 = glv[:a1], glv[:b1], glv[:a2], glv[:b2]
    halfN = order / 2
    c1 = (b2 * k + halfN) / order
    c2 = (-b1 * k + halfN) / order
    k1 = k - c1 * a1 - c2 * a2
    k2 = -c1 * b1 - c2 * b2
    return k1, k2
end

._glvMultiplyAndAdd(p1, n1, p2, n2, curve) ⇒ Object 



176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
 
# File 'lib/math.rb', line 176

def self._glvMultiplyAndAdd(p1, n1, p2, n2, curve)
    # Compute n1*p1 + n2*p2 using the GLV endomorphism. Splits each 256-bit
    # scalar into two ~128-bit scalars via k == k1 + k2*lambda (mod N), then
    # runs a 4-scalar simultaneous double-and-add over (p1, phi(p1), p2,
    # phi(p2)) with a 16-entry precomputed table of subset sums. Halves the
    # loop length versus the plain Shamir path.
    glv = curve.glvParams
    order, coeff, prime = curve.n, curve.a, curve.p
    beta = glv[:beta]

    k1, k2 = self._glvDecompose(n1 % order, glv, order)
    k3, k4 = self._glvDecompose(n2 % order, glv, order)

    # Base points (affine, z=1) - phi((x,y)) = (beta*x mod P, y).
    bases = [
        Point.new(p1.x, p1.y, 1),
        Point.new((beta * p1.x) % prime, p1.y, 1),
        Point.new(p2.x, p2.y, 1),
        Point.new((beta * p2.x) % prime, p2.y, 1),
    ]
    scalars = [k1, k2, k3, k4]
    (0...4).each do |i|
        if scalars[i] < 0
            scalars[i] = -scalars[i]
            bases[i] = Point.new(bases[i].x, prime - bases[i].y, 1)
        end
    end

    # Precompute table[idx] = sum of bases[i] selected by bits of idx.
    table = Array.new(16, Point.new(0, 0, 1))
    (1...16).each do |idx|
        low = idx & -idx
        i = low.bit_length - 1
        table[idx] = self._jacobianAdd(table[idx ^ low], bases[i], coeff, prime)
    end

    maxLen = scalars.map(&:bit_length).max
    r = Point.new(0, 0, 1)
    s0, s1, s2, s3 = scalars
    (maxLen - 1).downto(0) do |bit|
        r = self._jacobianDouble(r, coeff, prime)
        idx = ((s0 >> bit) & 1) | (((s1 >> bit) & 1) << 1) \
            | (((s2 >> bit) & 1) << 2) | (((s3 >> bit) & 1) << 3)
        if idx != 0
            r = self._jacobianAdd(r, table[idx], coeff, prime)
        end
    end

    return self._fromJacobian(r, prime)
end

._jacobianAdd(p, q, coeff, prime) ⇒ Object 



318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
 
# File 'lib/math.rb', line 318

def self._jacobianAdd(p, q, coeff, prime)
    # Add two points in elliptic curves
    #
    # :param p: First Point you want to add
    # :param q: Second Point you want to add
    # :param coeff: Coefficient of the first-order term of the equation Y^2 = X^3 + A*X + B (mod p)
    # :param prime: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    # :return: Point that represents the sum of First and Second Point
    if p.y == 0
        return q
    end
    if q.y == 0
        return p
    end

    px, py, pz = p.x, p.y, p.z
    qx, qy, qz = q.x, q.y, q.z

    pz2 = (pz * pz) % prime
    u2 = (qx * pz2) % prime
    s2 = (qy * pz2 * pz) % prime

    if qz == 1
        # Mixed affine+Jacobian add: qz^2=qz^3=1 saves four multiplications.
        u1 = px
        s1 = py
    else
        qz2 = (qz * qz) % prime
        u1 = (px * qz2) % prime
        s1 = (py * qz2 * qz) % prime
    end

    if u1 == u2
        if s1 != s2
            return Point.new(0, 0, 1)
        end
        return self._jacobianDouble(p, coeff, prime)
    end

    h = u2 - u1
    r = s2 - s1
    h2 = (h * h) % prime
    h3 = (h * h2) % prime
    u1h2 = (u1 * h2) % prime
    nx = (r * r - h3 - 2 * u1h2) % prime
    ny = (r * (u1h2 - nx) - s1 * h3) % prime
    nz = qz == 1 ? (h * pz) % prime : (h * pz * qz) % prime

    return Point.new(nx, ny, nz)
end

._jacobianDouble(p, coeff, prime) ⇒ Object 



288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
 
# File 'lib/math.rb', line 288

def self._jacobianDouble(p, coeff, prime)
    # Double a point in elliptic curves
    #
    # :param p: Point you want to double
    # :param coeff: Coefficient of the first-order term of the equation Y^2 = X^3 + A*X + B (mod p)
    # :param prime: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    # :return: Point that represents the sum of First and Second Point
    py = p.y
    if py == 0
        return Point.new(0, 0, 0)
    end

    px, pz = p.x, p.z
    ysq = (py * py) % prime
    s = (4 * px * ysq) % prime
    pz2 = (pz * pz) % prime
    if coeff == 0
        m = (3 * px * px) % prime
    elsif coeff == prime - 3
        m = (3 * (px - pz2) * (px + pz2)) % prime
    else
        m = (3 * px * px + coeff * pz2 * pz2) % prime
    end
    nx = (m * m - 2 * s) % prime
    ny = (m * (s - nx) - 8 * ysq * ysq) % prime
    nz = (2 * py * pz) % prime

    return Point.new(nx, ny, nz)
end

._jacobianMultiply(p, n, order, coeff, prime) ⇒ Object 



369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
 
# File 'lib/math.rb', line 369

def self._jacobianMultiply(p, n, order, coeff, prime)
    # Multiply point and scalar in elliptic curves using Montgomery ladder
    # for constant-time execution.
    #
    # :param p: First Point to multiply
    # :param n: Scalar to multiply
    # :param order: Order of the elliptic curve
    # :param coeff: Coefficient of the first-order term of the equation Y^2 = X^3 + A*X + B (mod p)
    # :param prime: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    # :return: Point that represents the scalar multiplication
    if p.y == 0 or n == 0
        return Point.new(0, 0, 1)
    end

    if n < 0 or n >= order
        n = n % order
    end

    if n == 0
        return Point.new(0, 0, 1)
    end

    # Montgomery ladder: always performs one add and one double per bit
    r0 = Point.new(0, 0, 1)
    r1 = Point.new(p.x, p.y, p.z)

    (n.bit_length - 1).downto(0) do |i|
        if (n >> i) & 1 == 0
            r1 = self._jacobianAdd(r0, r1, coeff, prime)
            r0 = self._jacobianDouble(r0, coeff, prime)
        else
            r0 = self._jacobianAdd(r0, r1, coeff, prime)
            r1 = self._jacobianDouble(r1, coeff, prime)
        end
    end

    return r0
end

._jsfDigits(k0, k1) ⇒ Object 



460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
 
# File 'lib/math.rb', line 460

def self._jsfDigits(k0, k1)
    # Joint Sparse Form of (k0, k1): list of signed-digit pairs (u0, u1) in
    # {-1, 0, 1}, ordered MSB-first. At most one of any two consecutive pairs
    # is non-zero, giving density ~1/2 instead of ~3/4 from raw binary.
    digits = []
    d0 = 0
    d1 = 0
    while k0 + d0 != 0 or k1 + d1 != 0
        a0 = k0 + d0
        a1 = k1 + d1
        if (a0 & 1) != 0
            u0 = (a0 & 3) == 1 ? 1 : -1
            if [3, 5].include?(a0 & 7) and (a1 & 3) == 2
                u0 = -u0
            end
        else
            u0 = 0
        end
        if (a1 & 1) != 0
            u1 = (a1 & 3) == 1 ? 1 : -1
            if [3, 5].include?(a1 & 7) and (a0 & 3) == 2
                u1 = -u1
            end
        else
            u1 = 0
        end
        digits << [u0, u1]
        if 2 * d0 == 1 + u0
            d0 = 1 - d0
        end
        if 2 * d1 == 1 + u1
            d1 = 1 - d1
        end
        k0 >>= 1
        k1 >>= 1
    end
    digits.reverse
end

._shamirMultiply(jp1, n1, jp2, n2, order, coeff, prime) ⇒ Object 



408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
 
# File 'lib/math.rb', line 408

def self._shamirMultiply(jp1, n1, jp2, n2, order, coeff, prime)
    # Compute n1*p1 + n2*p2 using Shamir's trick with Joint Sparse Form
    # (Solinas 2001). JSF picks signed digits in {-1, 0, 1} so at most ~l/2
    # digit pairs are non-zero, versus ~3l/4 for the raw binary form. Not
    # constant-time - use only with public scalars (e.g. verification).
    #
    # :param jp1: First point in Jacobian coordinates
    # :param n1: First scalar
    # :param jp2: Second point in Jacobian coordinates
    # :param n2: Second scalar
    # :param order: Order of the elliptic curve
    # :param coeff: Coefficient of the first-order term of the equation Y^2 = X^3 + A*X + B (mod p)
    # :param prime: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    # :return: Point n1*p1 + n2*p2 in Jacobian coordinates
    if n1 < 0 or n1 >= order
        n1 = n1 % order
    end
    if n2 < 0 or n2 >= order
        n2 = n2 % order
    end

    if n1 == 0 and n2 == 0
        return Point.new(0, 0, 1)
    end

    neg = lambda { |pt| Point.new(pt.x, pt.y == 0 ? 0 : prime - pt.y, pt.z) }

    jp1p2 = self._jacobianAdd(jp1, jp2, coeff, prime)
    jp1mp2 = self._jacobianAdd(jp1, neg.call(jp2), coeff, prime)
    addTable = {
        [1, 0]   => jp1,
        [-1, 0]  => neg.call(jp1),
        [0, 1]   => jp2,
        [0, -1]  => neg.call(jp2),
        [1, 1]   => jp1p2,
        [-1, -1] => neg.call(jp1p2),
        [1, -1]  => jp1mp2,
        [-1, 1]  => neg.call(jp1mp2),
    }

    digits = self._jsfDigits(n1, n2)
    r = Point.new(0, 0, 1)
    digits.each do |u0, u1|
        r = self._jacobianDouble(r, coeff, prime)
        if u0 != 0 or u1 != 0
            r = self._jacobianAdd(r, addTable[[u0, u1]], coeff, prime)
        end
    end

    return r
end

._toJacobian(p) ⇒ Object 



263
264
265
266
267
268
269
 
# File 'lib/math.rb', line 263

def self._toJacobian(p)
    # Convert point to Jacobian coordinates
    #
    # :param p: First Point you want to add
    # :return: Point in Jacobian coordinates
    return Point.new(p.x, p.y, 1)
end

.add(p, q, coeff, prime) ⇒ Object 



132
133
134
135
136
137
138
139
140
141
142
143
 
# File 'lib/math.rb', line 132

def self.add(p, q, coeff, prime)
    # Fast way to add two points in elliptic curves
    #
    # :param p: First Point you want to add
    # :param q: Second Point you want to add
    # :param coeff: Coefficient of the first-order term of the equation Y^2 = X^3 + A*X + B (mod p)
    # :param prime: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    # :return: Point that represents the sum of First and Second Point
    return self._fromJacobian(
        self._jacobianAdd(self._toJacobian(p), self._toJacobian(q), coeff, prime), prime
    )
end

.inv(x, n) ⇒ Object 



240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
 
# File 'lib/math.rb', line 240

def self.inv(x, n)
    # Modular inverse via extended Euclidean algorithm.
    # Roughly 2-3x faster than Fermat's little theorem for 256-bit operands.
    #
    # :param x: Divisor (must be coprime to n)
    # :param n: Mod for division
    # :return: Value representing the modular inverse
    if x % n == 0
        raise ArgumentError, "0 has no modular inverse"
    end

    a = x % n
    b = n
    x0 = 1
    x1 = 0
    while b != 0
        q = a / b
        a, b = b, a - q * b
        x0, x1 = x1, x0 - q * x1
    end
    return x0 % n
end

.modularSquareRoot(value, prime) ⇒ Object 



3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
 
# File 'lib/math.rb', line 3

def self.modularSquareRoot(value, prime)
    # Tonelli-Shanks algorithm for modular square root. Works for all odd primes.
    if value == 0
        return 0
    end
    if prime == 2
        return value % 2
    end

    # Factor out powers of 2: prime - 1 = Q * 2^S
    q = prime - 1
    s = 0
    while q % 2 == 0
        q /= 2
        s += 1
    end

    if s == 1  # prime = 3 (mod 4)
        return value.pow((prime + 1) / 4, prime)
    end

    # Find a quadratic non-residue z
    z = 2
    while z.pow((prime - 1) / 2, prime) != prime - 1
        z += 1
    end

    m = s
    c = z.pow(q, prime)
    t = value.pow(q, prime)
    r = value.pow((q + 1) / 2, prime)

    while true
        if t == 1
            return r
        end

        # Find the least i such that t^(2^i) = 1 (mod prime)
        i = 1
        temp = (t * t) % prime
        while temp != 1
            temp = (temp * temp) % prime
            i += 1
        end

        b = c.pow(1 << (m - i - 1), prime)
        m = i
        c = (b * b) % prime
        t = (t * c) % prime
        r = (r * b) % prime
    end
end

.multiply(p, n, order, coeff, prime) ⇒ Object 



118
119
120
121
122
123
124
125
126
127
128
129
130
 
# File 'lib/math.rb', line 118

def self.multiply(p, n, order, coeff, prime)
    # Fast way to multiply point and scalar in elliptic curves
    #
    # :param p: First Point to multiply
    # :param n: Scalar to multiply
    # :param order: Order of the elliptic curve
    # :param coeff: Coefficient of the first-order term of the equation Y^2 = X^3 + A*X + B (mod p)
    # :param prime: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    # :return: Point that represents the scalar multiplication
    return self._fromJacobian(
        self._jacobianMultiply(self._toJacobian(p), n, order, coeff, prime), prime
    )
end

.multiplyAndAdd(p1, n1, p2, n2, order, coeff, prime, curve = nil) ⇒ Object 



145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
 
# File 'lib/math.rb', line 145

def self.multiplyAndAdd(p1, n1, p2, n2, order, coeff, prime, curve=nil)
    # Compute n1*p1 + n2*p2. If `curve` is given and exposes `glvParams`
    # (e.g. secp256k1), uses the GLV endomorphism to split both scalars
    # into ~128-bit halves and run a 4-scalar simultaneous multi-
    # exponentiation. Otherwise falls back to Shamir's trick with JSF.
    # Not constant-time - use only with public scalars (e.g. verification).
    #
    # :param p1: First point
    # :param n1: First scalar
    # :param p2: Second point
    # :param n2: Second scalar
    # :param order: Order of the elliptic curve (ignored when `curve` is given)
    # :param coeff: Coefficient of the first-order term of the equation Y^2 = X^3 + A*X + B (mod p) (ignored when `curve` is given)
    # :param prime: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p) (ignored when `curve` is given)
    # :param curve: Optional curve object; enables GLV if `curve.glvParams` is set
    # :return: Point n1*p1 + n2*p2
    if not curve.nil?
        order, coeff, prime = curve.n, curve.a, curve.p
        if not curve.glvParams.nil?
            return self._glvMultiplyAndAdd(p1, n1, p2, n2, curve)
        end
    end
    return self._fromJacobian(
        self._shamirMultiply(
            self._toJacobian(p1), n1,
            self._toJacobian(p2), n2,
            order, coeff, prime,
        ), prime,
    )
end

.multiplyGenerator(curve, n) ⇒ Object 



56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
 
# File 'lib/math.rb', line 56

def self.multiplyGenerator(curve, n)
    # Fast scalar multiplication n*G using a precomputed affine table of
    # powers-of-two multiples of G and the width-2 NAF of n. Every non-zero
    # NAF digit triggers one mixed add and zero doublings, trading the ~256
    # doublings of a windowed method for ~86 adds on average - a large net
    # reduction in field multiplications for 256-bit scalars.
    if n < 0 or n >= curve.n
        n = n % curve.n
    end
    if n == 0
        return Point.new(0, 0, 0)
    end

    table = self._generatorPowersTable(curve)
    coeff = curve.a
    prime = curve.p

    r = Point.new(0, 0, 1)
    i = 0
    k = n
    while k > 0
        if (k & 1) != 0
            digit = 2 - (k & 3)  # -1 or +1
            k -= digit
            g = table[i]
            if digit == 1
                r = self._jacobianAdd(r, g, coeff, prime)
            else
                r = self._jacobianAdd(r, Point.new(g.x, prime - g.y, 1), coeff, prime)
            end
        end
        k >>= 1
        i += 1
    end
    return self._fromJacobian(r, prime)
end