804419db9c418848196b2b4288d798e3af3fb737
[gnuk/gnuk.git] / src / ecc-mont.c
1 /*                                                    -*- coding: utf-8 -*-
2  * ecc-mont.c - Elliptic curve computation for
3  *              the Montgomery curve: y^2 = x^3 + 486662*x^2 + x.
4  *
5  * Copyright (C) 2014 Free Software Initiative of Japan
6  * Author: NIIBE Yutaka <gniibe@fsij.org>
7  *
8  * This file is a part of Gnuk, a GnuPG USB Token implementation.
9  *
10  * Gnuk is free software: you can redistribute it and/or modify it
11  * under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * Gnuk is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
22  *
23  */
24
25 #include <stdint.h>
26 #include <string.h>
27 #include "bn.h"
28 #include "mod25638.h"
29 #include "mod.h"
30
31 /*
32  * References:
33  *
34  * [1] D. J. Bernstein. Curve25519: new Diffie-Hellman speed records.
35  *     Proceedings of PKC 2006, to appear. 
36  *     http://cr.yp.to/papers.html#curve25519. Date: 2006.02.09.
37  *
38  * [2] D. J. Bernstein. Can we avoid tests for zero in fast
39  *     elliptic-curve arithmetic?
40  *     http://cr.yp.to/papers.html#curvezero. Date: 2006.07.26.
41  *
42  */
43
44 /*
45  * IMPLEMENTATION NOTE
46  *
47  * (0) We assume that the processor has no cache, nor branch target
48  *     prediction.  Thus, we don't avoid indexing by secret value. 
49  *     We don't avoid conditional jump if both cases have same timing,
50  *     either.
51  *
52  * (1) We use Radix-32 field arithmetic.  It's a representation like
53  *     2^256-38, but it's more redundant.  For example, "1" can be
54  *     represented in three ways in 256-bit: 1, 2^255-18, and
55  *     2^256-37.
56  *
57  * (2) We use Montgomery double-and-add.
58  *
59  */
60
61 #ifndef BN256_C_IMPLEMENTATION
62 #define ASM_IMPLEMENTATION 1
63 #endif
64 /*
65  *
66  * 121665 = 0x1db41
67  *            1 1101 1011 0100 0001
68  */
69 static void
70 mod25638_mul_121665 (bn256 *x, const bn256 *a)
71 {
72 #if ASM_IMPLEMENTATION
73 #include "muladd_256.h"
74   const uint32_t *s;
75   uint32_t *d;
76   uint32_t w;
77   uint32_t c;
78
79   s = a->word;
80   d = x->word;
81   w = 121665;
82   MULADD_256_ASM (s, d, w, c);
83 #else
84   uint32_t c, c1;
85   bn256 m[1];
86
87   c = c1 = bn256_shift (m, a, 6); c += bn256_add (x, a, m);
88   c1 <<= 2; c1 |= bn256_shift (m, m, 2); c = c + c1 + bn256_add (x, x, m);
89   c1 <<= 1; c1 |= bn256_shift (m, m, 1); c = c + c1 + bn256_add (x, x, m);
90   c1 <<= 2; c1 |= bn256_shift (m, m, 2); c = c + c1 + bn256_add (x, x, m);
91   c1 <<= 1; c1 |= bn256_shift (m, m, 1); c = c + c1 + bn256_add (x, x, m);
92   c1 <<= 2; c1 |= bn256_shift (m, m, 2); c = c + c1 + bn256_add (x, x, m);
93   c1 <<= 1; c1 |= bn256_shift (m, m, 1); c = c + c1 + bn256_add (x, x, m);
94   c1 <<= 1; c1 |= bn256_shift (m, m, 1); c = c + c1 + bn256_add (x, x, m);
95 #endif
96   c = bn256_add_uint (x, x, c*38);
97   x->word[0] += c * 38;
98 }
99
100
101 typedef struct
102 {
103   bn256 x[1];
104   bn256 z[1];
105 } pt;
106
107
108 /**
109  * @brief  Process Montgomery double-and-add
110  *
111  * With Q0, Q1, DIF (= Q0 - Q1), compute PRD = 2Q0, SUM = Q0 + Q1
112  * Q0 and Q1 are clobbered.
113  *
114  */
115 static void
116 mont_d_and_a (pt *prd, pt *sum, pt *q0, pt *q1, const bn256 *dif_x)
117 {
118                                         mod25638_add (sum->x, q1->x, q1->z);
119                                         mod25638_sub (q1->z, q1->x, q1->z);
120   mod25638_add (prd->x, q0->x, q0->z);
121   mod25638_sub (q0->z, q0->x, q0->z);
122                                         mod25638_mul (q1->x, q0->z, sum->x);
123                                         mod25638_mul (q1->z, prd->x, q1->z);
124   mod25638_sqr (q0->x, prd->x);
125   mod25638_sqr (q0->z, q0->z);
126                                         mod25638_add (sum->x, q1->x, q1->z);
127                                         mod25638_sub (q1->z, q1->x, q1->z);
128   mod25638_mul (prd->x, q0->x, q0->z);
129   mod25638_sub (q0->z, q0->x, q0->z);
130                                         mod25638_sqr (sum->x, sum->x);
131                                         mod25638_sqr (sum->z, q1->z);
132   mod25638_mul_121665 (prd->z, q0->z);
133                                         mod25638_mul (sum->z, sum->z, dif_x);
134   mod25638_add (prd->z, q0->x, prd->z);
135   mod25638_mul (prd->z, prd->z, q0->z);
136 }
137
138
139 /**
140  * @brief       RES  = x-coordinate of [n]Q
141  *
142  * @param N     Scalar N (three least significant bits are 000)
143  * @param Q_X   x-coordinate of Q
144  *
145  */
146 void
147 compute_nQ (bn256 *res, const bn256 *n, const bn256 *q_x)
148 {
149   int i, j;
150   pt p0[1], p1[1], p0_[1], p1_[1];
151
152   /* P0 = O = (1:0)  */
153   memset (p0->x, 0, sizeof (bn256));
154   p0->x->word[0] = 1;
155   memset (p0->z, 0, sizeof (bn256));
156
157   /* P1 = (X:1) */
158   memcpy (p1->x, q_x, sizeof (bn256));
159   memset (p1->z, 0, sizeof (bn256));
160   p1->z->word[0] = 1;
161
162   for (i = 0; i < 8; i++)
163     {
164       uint32_t u = n->word[7-i];
165
166       for (j = 0; j < 16; j++)
167         {
168           pt *q0, *q1;
169           pt *sum_n, *prd_n;
170
171           if ((u & 0x80000000))
172             q0 = p1,  q1 = p0,  sum_n = p0_, prd_n = p1_;
173           else
174             q0 = p0,  q1 = p1,  sum_n = p1_, prd_n = p0_;
175           mont_d_and_a (prd_n, sum_n, q0, q1, q_x);
176
177           if ((u & 0x40000000))
178             q0 = p1_, q1 = p0_, sum_n = p0,  prd_n = p1;
179           else
180             q0 = p0_, q1 = p1_, sum_n = p1,  prd_n = p0;
181           mont_d_and_a (prd_n, sum_n, q0, q1, q_x);
182
183           u <<= 2;
184         }
185     }
186
187   /* We know the LSB of N is always 0.  Thus, result is always in P0.  */
188   /*
189    * p0->z may be zero here, but our mod_inv doesn't raise error for 0,
190    * but returns 0 (like the implementation of z^(p-2)), thus, RES will
191    * be 0 in that case, which is correct value.
192    */
193   mod_inv (res, p0->z, p25519);
194   mod25638_mul (res, res, p0->x);
195   mod25519_reduce (res);
196 }