1<?php
2
3/**
4 * Curves over y^2 = x^3 + a*x + b
5 *
6 * These are curves used in SEC 2 over prime fields: http://www.secg.org/SEC2-Ver-1.0.pdf
7 * The curve is a weierstrass curve with a[1], a[3] and a[2] set to 0.
8 *
9 * Uses Jacobian Coordinates for speed if able:
10 *
11 * https://en.wikipedia.org/wiki/Jacobian_curve
12 * https://en.wikibooks.org/wiki/Cryptography/Prime_Curve/Jacobian_Coordinates
13 *
14 * PHP version 5 and 7
15 *
16 * @category  Crypt
17 * @package   EC
18 * @author    Jim Wigginton <terrafrost@php.net>
19 * @copyright 2017 Jim Wigginton
20 * @license   http://www.opensource.org/licenses/mit-license.html  MIT License
21 * @link      http://pear.php.net/package/Math_BigInteger
22 */
23
24namespace phpseclib3\Crypt\EC\BaseCurves;
25
26use phpseclib3\Common\Functions\Strings;
27use phpseclib3\Math\BigInteger;
28use phpseclib3\Math\Common\FiniteField\Integer;
29use phpseclib3\Math\PrimeField;
30use phpseclib3\Math\PrimeField\Integer as PrimeInteger;
31
32/**
33 * Curves over y^2 = x^3 + a*x + b
34 *
35 * @package Prime
36 * @author  Jim Wigginton <terrafrost@php.net>
37 * @access  public
38 */
39class Prime extends Base
40{
41    /**
42     * Prime Field Integer factory
43     *
44     * @var \phpseclib3\Math\PrimeFields
45     */
46    protected $factory;
47
48    /**
49     * Cofficient for x^1
50     *
51     * @var object
52     */
53    protected $a;
54
55    /**
56     * Cofficient for x^0
57     *
58     * @var object
59     */
60    protected $b;
61
62    /**
63     * Base Point
64     *
65     * @var object
66     */
67    protected $p;
68
69    /**
70     * The number one over the specified finite field
71     *
72     * @var object
73     */
74    protected $one;
75
76    /**
77     * The number two over the specified finite field
78     *
79     * @var object
80     */
81    protected $two;
82
83    /**
84     * The number three over the specified finite field
85     *
86     * @var object
87     */
88    protected $three;
89
90    /**
91     * The number four over the specified finite field
92     *
93     * @var object
94     */
95    protected $four;
96
97    /**
98     * The number eight over the specified finite field
99     *
100     * @var object
101     */
102    protected $eight;
103
104    /**
105     * The modulo
106     *
107     * @var BigInteger
108     */
109    protected $modulo;
110
111    /**
112     * The Order
113     *
114     * @var BigInteger
115     */
116    protected $order;
117
118    /**
119     * Sets the modulo
120     */
121    public function setModulo(BigInteger $modulo)
122    {
123        $this->modulo = $modulo;
124        $this->factory = new PrimeField($modulo);
125        $this->two = $this->factory->newInteger(new BigInteger(2));
126        $this->three = $this->factory->newInteger(new BigInteger(3));
127        // used by jacobian coordinates
128        $this->one = $this->factory->newInteger(new BigInteger(1));
129        $this->four = $this->factory->newInteger(new BigInteger(4));
130        $this->eight = $this->factory->newInteger(new BigInteger(8));
131    }
132
133    /**
134     * Set coefficients a and b
135     */
136    public function setCoefficients(BigInteger $a, BigInteger $b)
137    {
138        if (!isset($this->factory)) {
139            throw new \RuntimeException('setModulo needs to be called before this method');
140        }
141        $this->a = $this->factory->newInteger($a);
142        $this->b = $this->factory->newInteger($b);
143    }
144
145    /**
146     * Set x and y coordinates for the base point
147     *
148     * @param BigInteger|PrimeInteger $x
149     * @param BigInteger|PrimeInteger $y
150     * @return PrimeInteger[]
151     */
152    public function setBasePoint($x, $y)
153    {
154        switch (true) {
155            case !$x instanceof BigInteger && !$x instanceof PrimeInteger:
156                throw new \UnexpectedValueException('Argument 1 passed to Prime::setBasePoint() must be an instance of either BigInteger or PrimeField\Integer');
157            case !$y instanceof BigInteger && !$y instanceof PrimeInteger:
158                throw new \UnexpectedValueException('Argument 2 passed to Prime::setBasePoint() must be an instance of either BigInteger or PrimeField\Integer');
159        }
160        if (!isset($this->factory)) {
161            throw new \RuntimeException('setModulo needs to be called before this method');
162        }
163        $this->p = [
164            $x instanceof BigInteger ? $this->factory->newInteger($x) : $x,
165            $y instanceof BigInteger ? $this->factory->newInteger($y) : $y
166        ];
167    }
168
169    /**
170     * Retrieve the base point as an array
171     *
172     * @return array
173     */
174    public function getBasePoint()
175    {
176        if (!isset($this->factory)) {
177            throw new \RuntimeException('setModulo needs to be called before this method');
178        }
179        /*
180        if (!isset($this->p)) {
181            throw new \RuntimeException('setBasePoint needs to be called before this method');
182        }
183        */
184        return $this->p;
185    }
186
187    /**
188     * Adds two "fresh" jacobian form on the curve
189     *
190     * @return FiniteField[]
191     */
192    protected function jacobianAddPointMixedXY(array $p, array $q)
193    {
194        list($u1, $s1) = $p;
195        list($u2, $s2) = $q;
196        if ($u1->equals($u2)) {
197            if (!$s1->equals($s2)) {
198                return [];
199            } else {
200                return $this->doublePoint($p);
201            }
202        }
203        $h = $u2->subtract($u1);
204        $r = $s2->subtract($s1);
205        $h2 = $h->multiply($h);
206        $h3 = $h2->multiply($h);
207        $v = $u1->multiply($h2);
208        $x3 = $r->multiply($r)->subtract($h3)->subtract($v->multiply($this->two));
209        $y3 = $r->multiply(
210            $v->subtract($x3)
211        )->subtract(
212            $s1->multiply($h3)
213        );
214        return [$x3, $y3, $h];
215    }
216
217    /**
218     * Adds one "fresh" jacobian form on the curve
219     *
220     * The second parameter should be the "fresh" one
221     *
222     * @return FiniteField[]
223     */
224    protected function jacobianAddPointMixedX(array $p, array $q)
225    {
226        list($u1, $s1, $z1) = $p;
227        list($x2, $y2) = $q;
228
229        $z12 = $z1->multiply($z1);
230
231        $u2 = $x2->multiply($z12);
232        $s2 = $y2->multiply($z12->multiply($z1));
233        if ($u1->equals($u2)) {
234            if (!$s1->equals($s2)) {
235                return [];
236            } else {
237                return $this->doublePoint($p);
238            }
239        }
240        $h = $u2->subtract($u1);
241        $r = $s2->subtract($s1);
242        $h2 = $h->multiply($h);
243        $h3 = $h2->multiply($h);
244        $v = $u1->multiply($h2);
245        $x3 = $r->multiply($r)->subtract($h3)->subtract($v->multiply($this->two));
246        $y3 = $r->multiply(
247            $v->subtract($x3)
248        )->subtract(
249            $s1->multiply($h3)
250        );
251        $z3 = $h->multiply($z1);
252        return [$x3, $y3, $z3];
253    }
254
255    /**
256     * Adds two jacobian coordinates on the curve
257     *
258     * @return FiniteField[]
259     */
260    protected function jacobianAddPoint(array $p, array $q)
261    {
262        list($x1, $y1, $z1) = $p;
263        list($x2, $y2, $z2) = $q;
264
265        $z12 = $z1->multiply($z1);
266        $z22 = $z2->multiply($z2);
267
268        $u1 = $x1->multiply($z22);
269        $u2 = $x2->multiply($z12);
270        $s1 = $y1->multiply($z22->multiply($z2));
271        $s2 = $y2->multiply($z12->multiply($z1));
272        if ($u1->equals($u2)) {
273            if (!$s1->equals($s2)) {
274                return [];
275            } else {
276                return $this->doublePoint($p);
277            }
278        }
279        $h = $u2->subtract($u1);
280        $r = $s2->subtract($s1);
281        $h2 = $h->multiply($h);
282        $h3 = $h2->multiply($h);
283        $v = $u1->multiply($h2);
284        $x3 = $r->multiply($r)->subtract($h3)->subtract($v->multiply($this->two));
285        $y3 = $r->multiply(
286            $v->subtract($x3)
287        )->subtract(
288            $s1->multiply($h3)
289        );
290        $z3 = $h->multiply($z1)->multiply($z2);
291        return [$x3, $y3, $z3];
292    }
293
294    /**
295     * Adds two points on the curve
296     *
297     * @return FiniteField[]
298     */
299    public function addPoint(array $p, array $q)
300    {
301        if (!isset($this->factory)) {
302            throw new \RuntimeException('setModulo needs to be called before this method');
303        }
304
305        if (!count($p) || !count($q)) {
306            if (count($q)) {
307                return $q;
308            }
309            if (count($p)) {
310                return $p;
311            }
312            return [];
313        }
314
315        // use jacobian coordinates
316        if (isset($p[2]) && isset($q[2])) {
317            if (isset($p['fresh']) && isset($q['fresh'])) {
318                return $this->jacobianAddPointMixedXY($p, $q);
319            }
320            if (isset($p['fresh'])) {
321                return $this->jacobianAddPointMixedX($q, $p);
322            }
323            if (isset($q['fresh'])) {
324                return $this->jacobianAddPointMixedX($p, $q);
325            }
326            return $this->jacobianAddPoint($p, $q);
327        }
328
329        if (isset($p[2]) || isset($q[2])) {
330            throw new \RuntimeException('Affine coordinates need to be manually converted to Jacobi coordinates or vice versa');
331        }
332
333        if ($p[0]->equals($q[0])) {
334            if (!$p[1]->equals($q[1])) {
335                return [];
336            } else { // eg. doublePoint
337                list($numerator, $denominator) = $this->doublePointHelper($p);
338            }
339        } else {
340            $numerator = $q[1]->subtract($p[1]);
341            $denominator = $q[0]->subtract($p[0]);
342        }
343        $slope = $numerator->divide($denominator);
344        $x = $slope->multiply($slope)->subtract($p[0])->subtract($q[0]);
345        $y = $slope->multiply($p[0]->subtract($x))->subtract($p[1]);
346
347        return [$x, $y];
348    }
349
350    /**
351     * Returns the numerator and denominator of the slope
352     *
353     * @return FiniteField[]
354     */
355    protected function doublePointHelper(array $p)
356    {
357        $numerator = $this->three->multiply($p[0])->multiply($p[0])->add($this->a);
358        $denominator = $this->two->multiply($p[1]);
359        return [$numerator, $denominator];
360    }
361
362    /**
363     * Doubles a jacobian coordinate on the curve
364     *
365     * @return FiniteField[]
366     */
367    protected function jacobianDoublePoint(array $p)
368    {
369        list($x, $y, $z) = $p;
370        $x2 = $x->multiply($x);
371        $y2 = $y->multiply($y);
372        $z2 = $z->multiply($z);
373        $s = $this->four->multiply($x)->multiply($y2);
374        $m1 = $this->three->multiply($x2);
375        $m2 = $this->a->multiply($z2->multiply($z2));
376        $m = $m1->add($m2);
377        $x1 = $m->multiply($m)->subtract($this->two->multiply($s));
378        $y1 = $m->multiply($s->subtract($x1))->subtract(
379            $this->eight->multiply($y2->multiply($y2))
380        );
381        $z1 = $this->two->multiply($y)->multiply($z);
382        return [$x1, $y1, $z1];
383    }
384
385    /**
386     * Doubles a "fresh" jacobian coordinate on the curve
387     *
388     * @return FiniteField[]
389     */
390    protected function jacobianDoublePointMixed(array $p)
391    {
392        list($x, $y) = $p;
393        $x2 = $x->multiply($x);
394        $y2 = $y->multiply($y);
395        $s = $this->four->multiply($x)->multiply($y2);
396        $m1 = $this->three->multiply($x2);
397        $m = $m1->add($this->a);
398        $x1 = $m->multiply($m)->subtract($this->two->multiply($s));
399        $y1 = $m->multiply($s->subtract($x1))->subtract(
400            $this->eight->multiply($y2->multiply($y2))
401        );
402        $z1 = $this->two->multiply($y);
403        return [$x1, $y1, $z1];
404    }
405
406    /**
407     * Doubles a point on a curve
408     *
409     * @return FiniteField[]
410     */
411    public function doublePoint(array $p)
412    {
413        if (!isset($this->factory)) {
414            throw new \RuntimeException('setModulo needs to be called before this method');
415        }
416
417        if (!count($p)) {
418            return [];
419        }
420
421        // use jacobian coordinates
422        if (isset($p[2])) {
423            if (isset($p['fresh'])) {
424                return $this->jacobianDoublePointMixed($p);
425            }
426            return $this->jacobianDoublePoint($p);
427        }
428
429        list($numerator, $denominator) = $this->doublePointHelper($p);
430
431        $slope = $numerator->divide($denominator);
432
433        $x = $slope->multiply($slope)->subtract($p[0])->subtract($p[0]);
434        $y = $slope->multiply($p[0]->subtract($x))->subtract($p[1]);
435
436        return [$x, $y];
437    }
438
439    /**
440     * Returns the X coordinate and the derived Y coordinate
441     *
442     * @return array
443     */
444    public function derivePoint($m)
445    {
446        $y = ord(Strings::shift($m));
447        $x = new BigInteger($m, 256);
448        $xp = $this->convertInteger($x);
449        switch ($y) {
450            case 2:
451                $ypn = false;
452                break;
453            case 3:
454                $ypn = true;
455                break;
456            default:
457                throw new \RuntimeException('Coordinate not in recognized format');
458        }
459        $temp = $xp->multiply($this->a);
460        $temp = $xp->multiply($xp)->multiply($xp)->add($temp);
461        $temp = $temp->add($this->b);
462        $b = $temp->squareRoot();
463        if (!$b) {
464            throw new \RuntimeException('Unable to derive Y coordinate');
465        }
466        $bn = $b->isOdd();
467        $yp = $ypn == $bn ? $b : $b->negate();
468        return [$xp, $yp];
469    }
470
471    /**
472     * Tests whether or not the x / y values satisfy the equation
473     *
474     * @return boolean
475     */
476    public function verifyPoint(array $p)
477    {
478        list($x, $y) = $p;
479        $lhs = $y->multiply($y);
480        $temp = $x->multiply($this->a);
481        $temp = $x->multiply($x)->multiply($x)->add($temp);
482        $rhs = $temp->add($this->b);
483
484        return $lhs->equals($rhs);
485    }
486
487    /**
488     * Returns the modulo
489     *
490     * @return \phpseclib3\Math\BigInteger
491     */
492    public function getModulo()
493    {
494        return $this->modulo;
495    }
496
497    /**
498     * Returns the a coefficient
499     *
500     * @return \phpseclib3\Math\PrimeField\Integer
501     */
502    public function getA()
503    {
504        return $this->a;
505    }
506
507    /**
508     * Returns the a coefficient
509     *
510     * @return \phpseclib3\Math\PrimeField\Integer
511     */
512    public function getB()
513    {
514        return $this->b;
515    }
516
517    /**
518     * Multiply and Add Points
519     *
520     * Adapted from https://git.io/vxPUH
521     *
522     * @return int[]
523     */
524    public function multiplyAddPoints(array $points, array $scalars)
525    {
526        $length = count($points);
527
528        foreach ($points as &$point) {
529            $point = $this->convertToInternal($point);
530        }
531
532        $wnd = [$this->getNAFPoints($points[0], 7)];
533        $wndWidth = [isset($points[0]['nafwidth']) ? $points[0]['nafwidth'] : 7];
534        for ($i = 1; $i < $length; $i++) {
535            $wnd[] = $this->getNAFPoints($points[$i], 1);
536            $wndWidth[] = isset($points[$i]['nafwidth']) ? $points[$i]['nafwidth'] : 1;
537        }
538
539        $naf = [];
540
541        // comb all window NAFs
542
543        $max = 0;
544        for ($i = $length - 1; $i >= 1; $i -= 2) {
545            $a = $i - 1;
546            $b = $i;
547            if ($wndWidth[$a] != 1 || $wndWidth[$b] != 1) {
548                $naf[$a] = $scalars[$a]->getNAF($wndWidth[$a]);
549                $naf[$b] = $scalars[$b]->getNAF($wndWidth[$b]);
550                $max = max(count($naf[$a]), count($naf[$b]), $max);
551                continue;
552            }
553
554            $comb = [
555                $points[$a], // 1
556                null,        // 3
557                null,        // 5
558                $points[$b]  // 7
559            ];
560
561            $comb[1] = $this->addPoint($points[$a], $points[$b]);
562            $comb[2] = $this->addPoint($points[$a], $this->negatePoint($points[$b]));
563
564            $index = [
565                -3, /* -1 -1 */
566                -1, /* -1  0 */
567                -5, /* -1  1 */
568                -7, /*  0 -1 */
569                 0, /*  0 -1 */
570                 7, /*  0  1 */
571                 5, /*  1 -1 */
572                 1, /*  1  0 */
573                 3  /*  1  1 */
574            ];
575
576            $jsf = self::getJSFPoints($scalars[$a], $scalars[$b]);
577
578            $max = max(count($jsf[0]), $max);
579            if ($max > 0) {
580                $naf[$a] = array_fill(0, $max, 0);
581                $naf[$b] = array_fill(0, $max, 0);
582            } else {
583                $naf[$a] = [];
584                $naf[$b] = [];
585            }
586
587            for ($j = 0; $j < $max; $j++) {
588                $ja = isset($jsf[0][$j]) ? $jsf[0][$j] : 0;
589                $jb = isset($jsf[1][$j]) ? $jsf[1][$j] : 0;
590
591                $naf[$a][$j] = $index[3 * ($ja + 1) + $jb + 1];
592                $naf[$b][$j] = 0;
593                $wnd[$a] = $comb;
594            }
595        }
596
597        $acc = [];
598        $temp = [0, 0, 0, 0];
599        for ($i = $max; $i >= 0; $i--) {
600            $k = 0;
601            while ($i >= 0) {
602                $zero = true;
603                for ($j = 0; $j < $length; $j++) {
604                    $temp[$j] = isset($naf[$j][$i]) ? $naf[$j][$i] : 0;
605                    if ($temp[$j] != 0) {
606                        $zero = false;
607                    }
608                }
609                if (!$zero) {
610                    break;
611                }
612                $k++;
613                $i--;
614            }
615
616            if ($i >= 0) {
617                $k++;
618            }
619            while ($k--) {
620                $acc = $this->doublePoint($acc);
621            }
622
623            if ($i < 0) {
624                break;
625            }
626
627            for ($j = 0; $j < $length; $j++) {
628                $z = $temp[$j];
629                $p = null;
630                if ($z == 0) {
631                    continue;
632                }
633                $p = $z > 0 ?
634                    $wnd[$j][($z - 1) >> 1] :
635                    $this->negatePoint($wnd[$j][(-$z - 1) >> 1]);
636                $acc = $this->addPoint($acc, $p);
637            }
638        }
639
640        return $this->convertToAffine($acc);
641    }
642
643    /**
644     * Precomputes NAF points
645     *
646     * Adapted from https://git.io/vxY1f
647     *
648     * @return int[]
649     */
650    private function getNAFPoints($point, $wnd)
651    {
652        if (isset($point['naf'])) {
653            return $point['naf'];
654        }
655
656        $res = [$point];
657        $max = (1 << $wnd) - 1;
658        $dbl = $max == 1 ? null : $this->doublePoint($point);
659        for ($i = 1; $i < $max; $i++) {
660            $res[] = $this->addPoint($res[$i - 1], $dbl);
661        }
662
663        $point['naf'] = $res;
664
665        /*
666        $str = '';
667        foreach ($res as $re) {
668            $re[0] = bin2hex($re[0]->toBytes());
669            $re[1] = bin2hex($re[1]->toBytes());
670            $str.= "            ['$re[0]', '$re[1]'],\r\n";
671        }
672        file_put_contents('temp.txt', $str);
673        exit;
674        */
675
676        return $res;
677    }
678
679    /**
680     * Precomputes points in Joint Sparse Form
681     *
682     * Adapted from https://git.io/vxrpD
683     *
684     * @return int[]
685     */
686    private static function getJSFPoints(Integer $k1, Integer $k2)
687    {
688        static $three;
689        if (!isset($three)) {
690            $three = new BigInteger(3);
691        }
692
693        $jsf = [[], []];
694        $k1 = $k1->toBigInteger();
695        $k2 = $k2->toBigInteger();
696        $d1 = 0;
697        $d2 = 0;
698
699        while ($k1->compare(new BigInteger(-$d1)) > 0 || $k2->compare(new BigInteger(-$d2)) > 0) {
700            // first phase
701            $m14 = $k1->testBit(0) + 2 * $k1->testBit(1);
702            $m14 += $d1;
703            $m14 &= 3;
704
705            $m24 = $k2->testBit(0) + 2 * $k2->testBit(1);
706            $m24 += $d2;
707            $m24 &= 3;
708
709            if ($m14 == 3) {
710                $m14 = -1;
711            }
712            if ($m24 == 3) {
713                $m24 = -1;
714            }
715
716            $u1 = 0;
717            if ($m14 & 1) { // if $m14 is odd
718                $m8 = $k1->testBit(0) + 2 * $k1->testBit(1) + 4 * $k1->testBit(2);
719                $m8 += $d1;
720                $m8 &= 7;
721                $u1 = ($m8 == 3 || $m8 == 5) && $m24 == 2 ? -$m14 : $m14;
722            }
723            $jsf[0][] = $u1;
724
725            $u2 = 0;
726            if ($m24 & 1) { // if $m24 is odd
727                $m8 = $k2->testBit(0) + 2 * $k2->testBit(1) + 4 * $k2->testBit(2);
728                $m8 += $d2;
729                $m8 &= 7;
730                $u2 = ($m8 == 3 || $m8 == 5) && $m14 == 2 ? -$m24 : $m24;
731            }
732            $jsf[1][] = $u2;
733
734            // second phase
735            if (2 * $d1 == $u1 + 1) {
736                $d1 = 1 - $d1;
737            }
738            if (2 * $d2 == $u2 + 1) {
739                $d2 = 1 - $d2;
740            }
741            $k1 = $k1->bitwise_rightShift(1);
742            $k2 = $k2->bitwise_rightShift(1);
743        }
744
745        return $jsf;
746    }
747
748    /**
749     * Returns the affine point
750     *
751     * A Jacobian Coordinate is of the form (x, y, z).
752     * To convert a Jacobian Coordinate to an Affine Point
753     * you do (x / z^2, y / z^3)
754     *
755     * @return \phpseclib3\Math\PrimeField\Integer[]
756     */
757    public function convertToAffine(array $p)
758    {
759        if (!isset($p[2])) {
760            return $p;
761        }
762        list($x, $y, $z) = $p;
763        $z = $this->one->divide($z);
764        $z2 = $z->multiply($z);
765        return [
766            $x->multiply($z2),
767            $y->multiply($z2)->multiply($z)
768        ];
769    }
770
771    /**
772     * Converts an affine point to a jacobian coordinate
773     *
774     * @return \phpseclib3\Math\PrimeField\Integer[]
775     */
776    public function convertToInternal(array $p)
777    {
778        if (isset($p[2])) {
779            return $p;
780        }
781
782        $p[2] = clone $this->one;
783        $p['fresh'] = true;
784        return $p;
785    }
786}
787