docker: source checking container
[gnuk/gnuk.git] / src / ecc-edwards.c
1 /*                                                    -*- coding: utf-8 -*-
2  * ecc-edwards.c - Elliptic curve computation for
3  *                 the twisted Edwards curve: -x^2 + y^2 = 1 + d*x^2*y^2
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 <stdlib.h>
27 #include <string.h>
28
29 #include "bn.h"
30 #include "mod.h"
31 #include "mod25638.h"
32 #include "sha512.h"
33
34 /*
35  * References:
36  *
37  * [1] Daniel J. Bernstein, Niels Duif, Tanja Lange, Peter Schwabe, Bo-Yin Yang.
38  *     High-speed high-security signatures.
39  *     Journal of Cryptographic Engineering 2 (2012), 77--89.
40  *     http://cr.yp.to/papers.html#ed25519
41  *
42  * [2] Daniel J. Bernstein, Peter Birkner, Marc Joye, Tanja Lange,
43  *     Christiane Peters.
44  *     Twisted Edwards curves.
45  *     Pages 389--405 in Progress in cryptology---AFRICACRYPT 2008.
46  *     http://cr.yp.to/papers.html#twisted
47  */
48
49 /*
50  * IMPLEMENTATION NOTE
51  *
52  * (0) We assume that the processor has no cache, nor branch target
53  *     prediction.  Thus, we don't avoid indexing by secret value.
54  *     We don't avoid conditional jump if both cases have same timing,
55  *     either.
56  *
57  * (1) We use Radix-32 field arithmetic.  It's a representation like
58  *     2^256-38, but it's more redundant.  For example, "1" can be
59  *     represented in three ways in 256-bit: 1, 2^255-18, and
60  *     2^256-37.
61  *
62  * (2) We use fixed base comb multiplication.  Scalar is 252-bit.
63  *     There are various possible choices for 252 = 2 * 2 * 3 * 3 * 7.
64  *     Current choice of total size is 3KB.  We use three tables, and
65  *     a table has 16 points (3 * 1KB).
66  *
67  *     Window size W = 4-bit, E = 21.
68  *                                                       <--21-bit-
69  *                                             <---42-bit----------
70  *     [        ][########][////////][        ][########][////////]
71  *                                   <-------63-bit----------------
72  *                         <-----------84-bit----------------------
73  *               <--------------105-bit----------------------------
74  *
75  *     [        ][########][////////][        ][########][////////]
76  *                                                                 <-126-bit-
77  *                                                       <-147-bit-
78  *                                             <----168-bit--------
79  *
80  *                                   <-------189-bit---------------
81  *                         <----------210-bit----------------------
82  *               <-------------231-bit-----------------------------
83  */
84
85 /*
86  * Identity element: (0,1)
87  * Negation: -(x,y) = (-x,y)
88  *
89  * d: -0x2DFC9311D490018C7338BF8688861767FF8FF5B2BEBE27548A14B235ECA6874A
90  * order:
91  *     0x1000000000000000000000000000000014DEF9DEA2F79CD65812631A5CF5D3ED
92  * Gx: 0x216936D3CD6E53FEC0A4E231FDD6DC5C692CC7609525A7B2C9562D608F25D51A
93  * Gy: 0x6666666666666666666666666666666666666666666666666666666666666658
94  */
95
96 /* d + 2^255 - 19 */
97 static const bn256 coefficient_d[1] = {
98   {{ 0x135978a3, 0x75eb4dca, 0x4141d8ab, 0x00700a4d,
99      0x7779e898, 0x8cc74079, 0x2b6ffe73, 0x52036cee }} };
100
101
102 /**
103  * @brief       Projective Twisted Coordinates
104  */
105 typedef struct
106 {
107   bn256 x[1];
108   bn256 y[1];
109   bn256 z[1];
110 } ptc;
111
112 #include "affine.h"
113
114
115 static int
116 mod25519_is_neg (const bn256 *a)
117 {
118   return (a->word[0] & 1);
119 }
120
121
122 /**
123  * @brief  X = 2 * A
124  *
125  * Compute (X3 : Y3 : Z3) = 2 * (X1 : Y1 : Z1)
126  */
127 static void
128 point_double (ptc *X, const ptc *A)
129 {
130   bn256 b[1], d[1], e[1];
131
132   /* Compute: B = (X1 + Y1)^2 */
133   mod25638_add (b, A->x, A->y);
134   mod25638_sqr (b, b);
135
136   /* Compute: C = X1^2        : E      */
137   mod25638_sqr (e, A->x);
138
139   /* Compute: D = Y1^2             */
140   mod25638_sqr (d, A->y);
141
142   /* E = aC; where a = -1 */
143   /* Compute: D - E = D + C : Y3_tmp */
144   mod25638_add (X->y, e, d);
145
146   /* Compute: -F = -(E + D) = C - D; where a = -1 : E */
147   mod25638_sub (e, e, d);
148
149   /* Compute: H = Z1^2        : D     */
150   mod25638_sqr (d, A->z);
151
152   /* Compute: -J = 2*H - F    : D     */
153   mod25638_add (d, d, d);
154   mod25638_add (d, d, e);
155
156   /* Compute: X3 = (B-C-D)*J = -J*(C+D-B) = -J*(Y3_tmp-B)  */
157   mod25638_sub (X->x, X->y, b);
158   mod25638_mul (X->x, X->x, d);
159
160   /* Compute: Y3 = -F*(D-E) = -F*Y3_tmp            */
161   mod25638_mul (X->y, X->y, e);
162
163   /* Z3 = -F*-J             */
164   mod25638_mul (X->z, e, d);
165 }
166
167
168 /**
169  * @brief       X = A + B
170  *
171  * @param X     Destination PTC
172  * @param A     PTC
173  * @param B     AC
174  *
175  * Compute: (X3 : Y3 : Z3) = (X1 : Y1 : Z1) + (X2 : Y2 : 1)
176  */
177 static void
178 point_add (ptc *X, const ptc *A, const ac *B)
179 {
180   bn256 c[1], d[1], e[1], tmp[1];
181
182   /* Compute: C = X1 * X2 */
183   mod25638_mul (c, A->x, B->x);
184
185   /* Compute: D = Y1 * Y2 */
186   mod25638_mul (d, A->y, B->y);
187
188   /* Compute: E = d * C * D */
189   mod25638_mul (e, c, d);
190   mod25638_mul (e, coefficient_d, e);
191
192   /* Compute: C_1 = C + D */
193   mod25638_add (c, c, d);
194
195   /* Compute: D_1 = Z1^2 : B */
196   mod25638_sqr (d, A->z);
197
198   /* tmp = D_1 - E : F */
199   mod25638_sub (tmp, d, e);
200
201   /* D_2 = D_1 + E : G */
202   mod25638_add (d, d, e);
203
204   /* X3_final = Z1 * tmp * ((X1 + Y1) * (X2 + Y2) - C_1) */
205   mod25638_add (X->x, A->x, A->y);
206   mod25638_add (e, B->x, B->y);
207   mod25638_mul (e, X->x, e);
208   mod25638_sub (e, e, c);
209   mod25638_mul (e, tmp, e);
210   mod25638_mul (X->x, A->z, e);
211
212   /* Y3_final = Z1 * D_2 * C_1 */
213   mod25638_mul (c, d, c);
214   mod25638_mul (X->y, A->z, c);
215
216   /* Z3_final = tmp * D_2 */
217   mod25638_mul (X->z, tmp, d);
218
219   /* A = Z1 */
220   /* B = A^2 */
221   /* C = X1 * X2 */
222   /* D = Y1 * Y2 */
223   /* E = d * C * D */
224   /* F = B - E */
225   /* G = B + E */
226   /* X3 = A * F * ((X1 + Y1) * (X2 + Y2) - C - D) */
227   /* Y3 = A * G * (D - aC); where a = -1 */
228   /* Z3 = F * G */
229 }
230
231
232 /**
233  * @brief       X = convert A
234  *
235  * @param X     Destination AC
236  * @param A     PTC
237  *
238  * (X1:Y1:Z1) represents the affine point (x=X1/Z1, y=Y1/Z1)
239  */
240 static void
241 point_ptc_to_ac (ac *X, const ptc *A)
242 {
243   bn256 z_inv[1];
244
245   /*
246    * A->z may be bigger than p25519, or two times bigger than p25519.
247    * But this is no problem for computation of mod_inv.
248    */
249   mod_inv (z_inv, A->z, p25519);
250
251   mod25638_mul (X->x, A->x, z_inv);
252   mod25519_reduce (X->x);
253   mod25638_mul (X->y, A->y, z_inv);
254   mod25519_reduce (X->y);
255 }
256
257
258 static const ac precomputed_KG[16] = {
259   { {{{ 0, 0, 0, 0, 0, 0, 0, 0 }}},
260     {{{ 1, 0, 0, 0, 0, 0, 0, 0 }}}                         },
261   { {{{ 0x8f25d51a, 0xc9562d60, 0x9525a7b2, 0x692cc760,
262         0xfdd6dc5c, 0xc0a4e231, 0xcd6e53fe, 0x216936d3 }}},
263     {{{ 0x66666658, 0x66666666, 0x66666666, 0x66666666,
264         0x66666666, 0x66666666, 0x66666666, 0x66666666 }}} },
265   { {{{ 0x3713af22, 0xac7137bd, 0xac634604, 0x25ed77a4,
266         0xa815e038, 0xce0d0064, 0xbca90151, 0x041c030f }}},
267     {{{ 0x0780f989, 0xe9b33fcf, 0x3d4445e7, 0xe4e97c2a,
268         0x655e5c16, 0xc67dc71c, 0xee43fb7a, 0x72467625 }}} },
269   { {{{ 0x3ee99893, 0x76a19171, 0x7ba9b065, 0xe647edd9,
270         0x6aeae260, 0x31f39299, 0x5f4a9bb2, 0x6d9e4545 }}},
271     {{{ 0x94cae280, 0xc41433da, 0x79061211, 0x8e842de8,
272         0xa259dc8a, 0xaab95e0b, 0x99013cd0, 0x28bd5fc3 }}} },
273   { {{{ 0x7d23ea24, 0x59e22c56, 0x0460850e, 0x1e745a88,
274         0xda13ef4b, 0x4583ff4c, 0x95083f85, 0x1f13202c }}},
275     {{{ 0x90275f48, 0xad42025c, 0xb55c4778, 0x0085087e,
276         0xfdfd7ffa, 0xf21109e7, 0x6c381b7e, 0x66336d35 }}} },
277   { {{{ 0xd00851f2, 0xaa9476ab, 0x4a61600b, 0xe7838534,
278         0x1a52df87, 0x0de65625, 0xbd675870, 0x5f0dd494 }}},
279     {{{ 0xe23493ba, 0xf20aec1b, 0x3414b0a8, 0x8f7f2741,
280         0xa80e1eb6, 0x497e74bd, 0xe9365b15, 0x1648eaac }}} },
281   { {{{ 0x04ac2b69, 0x5b78dcec, 0x32001a73, 0xecdb66ce,
282         0xb34cf697, 0xb75832f4, 0x3a2bce94, 0x7aaf57c5 }}},
283     {{{ 0x60fdfc6f, 0xb32ed2ce, 0x757924c6, 0x77bf20be,
284         0x48742dd1, 0xaebd15dd, 0x55d38439, 0x6311bb16 }}} },
285   { {{{ 0x42ff5c97, 0x139cdd73, 0xdbd82964, 0xee4c359e,
286         0x70611a3f, 0x91c1cd94, 0x8075dbcb, 0x1d0c34f6 }}},
287     {{{ 0x5f931219, 0x43eaa549, 0xa23d35a6, 0x3737aba7,
288         0x46f167bb, 0x54b1992f, 0xb74a9944, 0x01a11f3c }}} },
289   { {{{ 0xba46b161, 0x67a5310e, 0xd9d67f6c, 0x790f8527,
290         0x2f6cc814, 0x359c5b5f, 0x7786383d, 0x7b6a5565 }}},
291     {{{ 0x663ab0d3, 0xf1431b60, 0x09995826, 0x14a32d8f,
292         0xeddb8571, 0x61d526f6, 0x0eac739a, 0x0cb7acea }}} },
293   { {{{ 0x4a2d009f, 0x5eb1a697, 0xd8df987a, 0xdacb43b4,
294         0x8397f958, 0x4870f214, 0x8a175fbb, 0x5aa0c67c }}},
295     {{{ 0x78887db3, 0x27dbbd4c, 0x64e322ab, 0xe327b707,
296         0x7cbe4e3b, 0x87e293fa, 0xbda72395, 0x17040799 }}} },
297   { {{{ 0x99d1e696, 0xc833a5a2, 0x2d9d5877, 0x969bff8e,
298         0x2216fa67, 0x383a533a, 0x684d3925, 0x338bbe0a }}},
299     {{{ 0xd6cfb491, 0x35b5aae8, 0xaa12f3f8, 0x4a588279,
300         0x2e30380e, 0xa7c2e708, 0x9e4b3d62, 0x69f13e09 }}} },
301   { {{{ 0x27f1cd56, 0xec0dc2ef, 0xdb11cc97, 0x1af11548,
302         0x9ebc7613, 0xb642f86a, 0xcb77c3b9, 0x5ce45e73 }}},
303     {{{ 0x3eddd6de, 0x5d128786, 0x4859eab7, 0x16f9a6b4,
304         0xd8782345, 0x55c53916, 0xdb7b202a, 0x6b1dfa87 }}} },
305   { {{{ 0x19e30528, 0x2461a8ed, 0x665cfb1c, 0xaf756bf9,
306         0x3a6e8673, 0x0fcafd1d, 0x45d10f48, 0x0d264435 }}},
307     {{{ 0x5431db67, 0x543fd4c6, 0x60932432, 0xc153a5b3,
308         0xd2119aa4, 0x41d5b8eb, 0x8b09b6a5, 0x36bd9ab4 }}} },
309   { {{{ 0x21e06738, 0x6d39f935, 0x3765dd86, 0x4e6a7c59,
310         0xa4730880, 0xefc0dd80, 0x4079fe2f, 0x40617e56 }}},
311     {{{ 0x921439b9, 0xbc83cdff, 0x98833c09, 0xd5cccc06,
312         0xda13cdcb, 0xe315c425, 0x67ff5370, 0x37bc6e84 }}} },
313   { {{{ 0xf643b5f5, 0x65e7f028, 0x0ffbf5a8, 0x5b0d4831,
314         0xf4085f62, 0x0f540498, 0x0db7bd1b, 0x6f0bb035 }}},
315     {{{ 0x9733742c, 0x51f65571, 0xf513409f, 0x2fc047a0,
316         0x355facf6, 0x07f45010, 0x3a989a9c, 0x5cd416a9 }}} },
317   { {{{ 0x748f2a67, 0x0bdd7208, 0x415b7f7f, 0x0cf0b80b,
318         0x57aa0119, 0x44afdd5f, 0x430dc946, 0x05d68802 }}},
319     {{{ 0x1a60eeb2, 0x420c46e5, 0x665024f5, 0xc60a9b33,
320         0x48c51347, 0x37520265, 0x00a21bfb, 0x6f4be0af }}} }
321 };
322
323 static const ac precomputed_2E_KG[16] = {
324   { {{{ 0, 0, 0, 0, 0, 0, 0, 0 }}},
325     {{{ 1, 0, 0, 0, 0, 0, 0, 0 }}}                         },
326   { {{{ 0x199c4f7d, 0xec314ac0, 0xb2ebaaf9, 0x66a39c16,
327         0xedd4d15f, 0xab1c92b8, 0x57d9eada, 0x482a4cdf }}},
328     {{{ 0x6e4eb04b, 0xbd513b11, 0x25e4fd6a, 0x3f115fa5,
329         0x14519298, 0x0b3c5fc6, 0x81c2f7a8, 0x7391de43 }}} },
330   { {{{ 0x1254fe02, 0xa57dca18, 0x6da34368, 0xa56a2a14,
331         0x63e7328e, 0x44c6e34f, 0xca63ab3e, 0x3f748617 }}},
332     {{{ 0x7dc1641e, 0x5a13dc52, 0xee4e9ca1, 0x4cbb2899,
333         0x1ba9acee, 0x3938a289, 0x420fc47b, 0x0fed89e6 }}} },
334   { {{{ 0x49cbad08, 0x3c193f32, 0x15e80ef5, 0xdda71ef1,
335         0x9d128c33, 0xda44186c, 0xbf98c24f, 0x54183ede }}},
336     {{{ 0x93d165c1, 0x2cb483f7, 0x177f44aa, 0x51762ace,
337         0xb4ab035d, 0xb3fe651b, 0xa0b0d4e5, 0x426c99c3 }}} },
338   { {{{ 0xef3f3fb1, 0xb3fcf4d8, 0x065060a0, 0x7052292b,
339         0x24240b15, 0x18795ff8, 0x9989ffcc, 0x13aea184 }}},
340     {{{ 0xc2b81f44, 0x1930c101, 0x10600555, 0x672d6ca4,
341         0x1b25e570, 0xfbddbff2, 0x8ca12b70, 0x0884949c }}} },
342   { {{{ 0x00564bbf, 0x9983a033, 0xde61b72d, 0x95587d25,
343         0xeb17ad71, 0xb6719dfb, 0xc0bc3517, 0x46871ad0 }}},
344     {{{ 0xe95a6693, 0xb034fb61, 0x76eabad9, 0x5b0d8d18,
345         0x884785dc, 0xad295dd0, 0x74a1276a, 0x359debad }}} },
346   { {{{ 0xe89fb5ca, 0x2e5a2686, 0x5656c6c5, 0xd3d200ba,
347         0x9c969001, 0xef4c051e, 0x02cb45f4, 0x0d4ea946 }}},
348     {{{ 0x76d6e506, 0xa6f8a422, 0x63209e23, 0x454c768f,
349         0x2b372386, 0x5c12fd04, 0xdbfee11f, 0x1aedbd3e }}} },
350   { {{{ 0x00dbf569, 0x700ab50f, 0xd335b313, 0x9553643c,
351         0xa17dc97e, 0xeea9bddf, 0x3350a2bd, 0x0d12fe3d }}},
352     {{{ 0xa16a3dee, 0xe5ac35fe, 0xf81950c3, 0x4ae4664a,
353         0x3dbbf921, 0x75c63df4, 0x2958a5a6, 0x545b109c }}} },
354   { {{{ 0x0a61b29c, 0xd7a52a98, 0x65aca9ee, 0xe21e0acb,
355         0x5985dcbe, 0x57a69c0f, 0xeb87a534, 0x3c0c1e7b }}},
356     {{{ 0x6384bd2f, 0xf0a0b50d, 0xc6939e4b, 0xff349a34,
357         0x6e2f1973, 0x922c4554, 0xf1347631, 0x74e826b2 }}} },
358   { {{{ 0xa655803c, 0xd7eaa066, 0x38292c5c, 0x09504e76,
359         0x2c874953, 0xe298a02e, 0x8932b73f, 0x225093ed }}},
360     {{{ 0xe69c3efd, 0xf93e2b4d, 0x8a87c799, 0xa2cbd5fc,
361         0x85dba986, 0xdf41da94, 0xccee8edc, 0x36fe85e7 }}} },
362   { {{{ 0x7d742813, 0x78df7dc5, 0x4a193e64, 0x333bcc6d,
363         0x6a966d2d, 0x8242aa25, 0x4cd36d32, 0x03500a94 }}},
364     {{{ 0x580505d7, 0xd5d110fc, 0xfa11e1e9, 0xb2f47e16,
365         0x06eab6b4, 0xd0030f92, 0x62c91d46, 0x2dc80d5f }}} },
366   { {{{ 0x2a75e492, 0x5788b01a, 0xbae31352, 0x992acf54,
367         0x8159db27, 0x4591b980, 0xd3d84740, 0x36c6533c }}},
368     {{{ 0x103883b5, 0xc44c7c00, 0x515d0820, 0x10329423,
369         0x71b9dc16, 0xbd306903, 0xf88f8d32, 0x7edd5a95 }}} },
370   { {{{ 0x005523d7, 0xfd63b1ac, 0xad70dd21, 0x74482e0d,
371         0x02b56105, 0x67c9d9d0, 0x5971b456, 0x4d318012 }}},
372     {{{ 0x841106df, 0xdc9a6f6d, 0xa326987f, 0x7c52ed9d,
373         0x00607ea0, 0x4dbeaa6f, 0x6959e688, 0x115c221d }}} },
374   { {{{ 0xc80f7c16, 0xf8718464, 0xe9930634, 0x05dc8f40,
375         0xc2e9d5f4, 0xefa699bb, 0x021da209, 0x2469e813 }}},
376     {{{ 0xc602a3c4, 0x75c02845, 0x0a200f9d, 0x49d1b2ce,
377         0x2fb3ec8f, 0xd21b75e4, 0xd72a7545, 0x10dd726a }}} },
378   { {{{ 0x63ef1a6c, 0xeda58527, 0x051705e0, 0xb3fc0e72,
379         0x44f1161f, 0xbda6f3ee, 0xf339efe5, 0x7680aebf }}},
380     {{{ 0xb1b070a7, 0xe8d3fd01, 0xdbfbaaa0, 0xc3ff7dbf,
381         0xa320c916, 0xd81ef6f2, 0x62a3b54d, 0x3e22a1fb }}} },
382   { {{{ 0xb1fa18c8, 0xcdbb9187, 0xcb483a17, 0x8ddb5f6b,
383         0xea49af98, 0xc0a880b9, 0xf2dfddd0, 0x53bf600b }}},
384     {{{ 0x9e25b164, 0x4217404c, 0xafb74aa7, 0xfabf06ee,
385         0x2b9f233c, 0xb17712ae, 0xd0eb909e, 0x71f0b344 }}} }
386 };
387
388 static const ac precomputed_4E_KG[16] = {
389   { {{{ 0, 0, 0, 0, 0, 0, 0, 0 }}},
390     {{{ 1, 0, 0, 0, 0, 0, 0, 0 }}}                         },
391   { {{{ 0xe388a820, 0xbb6ec091, 0x5182278a, 0xa928b283,
392         0xa9a6eb83, 0x2259174d, 0x45500054, 0x184b48cb }}},
393     {{{ 0x26e77c33, 0xfe324dba, 0x83faf453, 0x6679a5e3,
394         0x2380ef73, 0xdd60c268, 0x03dc33a9, 0x3ee0e07a }}} },
395   { {{{ 0xce974493, 0x403aff28, 0x9bf6f5c4, 0x84076bf4,
396         0xecd898fb, 0xec57038c, 0xb663ed49, 0x2898ffaa }}},
397     {{{ 0xf335163d, 0xf4b3bc46, 0xfa4fb6c6, 0xe613a0f4,
398         0xb9934557, 0xe759d6bc, 0xab6c9477, 0x094f3b96 }}} },
399   { {{{ 0x6afffe9e, 0x168bb5a0, 0xee748c29, 0x950f7ad7,
400         0xda17203d, 0xa4850a2b, 0x77289e0f, 0x0062f7a7 }}},
401     {{{ 0x4b3829fa, 0x6265d4e9, 0xbdfcd386, 0x4f155ada,
402         0x475795f6, 0x9f38bda4, 0xdece4a4c, 0x560ed4b3 }}} },
403   { {{{ 0x141e648a, 0xdad4570a, 0x019b965c, 0x8bbf674c,
404         0xdb08fe30, 0xd7a8d50d, 0xa2851109, 0x7efb45d3 }}},
405     {{{ 0xd0c28cda, 0x52e818ac, 0xa321d436, 0x792257dd,
406         0x9d71f8b7, 0x867091c6, 0x11a1bf56, 0x0fe1198b }}} },
407   { {{{ 0x06137ab1, 0x4e848339, 0x3e6674cc, 0x5673e864,
408         0x0140502b, 0xad882043, 0x6ea1e46a, 0x34b5c0cb }}},
409     {{{ 0x1d70aa7c, 0x29786814, 0x8cdbb8aa, 0x840ae3f9,
410         0xbd4801fb, 0x78b4d622, 0xcf18ae9a, 0x6cf4e146 }}} },
411   { {{{ 0x36297168, 0x95c270ad, 0x942e7812, 0x2303ce80,
412         0x0205cf0e, 0x71908cc2, 0x32bcd754, 0x0cc15edd }}},
413     {{{ 0x2c7ded86, 0x1db94364, 0xf141b22c, 0xc694e39b,
414         0x5e5a9312, 0xf22f64ef, 0x3c5e6155, 0x649b8859 }}} },
415   { {{{ 0xb6417945, 0x0d5611c6, 0xac306c97, 0x9643fdbf,
416         0x0df500ff, 0xe81faaa4, 0x6f50e615, 0x0792c79b }}},
417     {{{ 0xd2af8c8d, 0xb45bbc49, 0x84f51bfe, 0x16c615ab,
418         0xc1d02d32, 0xdc57c526, 0x3c8aaa55, 0x5fb9a9a6 }}} },
419   { {{{ 0xdee40b98, 0x82faa8db, 0x6d520674, 0xff8a5208,
420         0x446ac562, 0x1f8c510f, 0x2cc6b66e, 0x4676d381 }}},
421     {{{ 0x2e7429f4, 0x8f1aa780, 0x8ed6bdf6, 0x2a95c1bf,
422         0x457fa0eb, 0x051450a0, 0x744c57b1, 0x7d89e2b7 }}} },
423   { {{{ 0x3f95ea15, 0xb6bdacd2, 0x2f1a5d69, 0xc9a9d1b1,
424         0xf4d22d72, 0xd4c2f1a9, 0x4dc516b5, 0x73ecfdf1 }}},
425     {{{ 0x05391e08, 0xa1ce93cd, 0x7b8aac17, 0x98f1e99e,
426         0xa098cbb3, 0x9ba84f2e, 0xf9bdd37a, 0x1425aa8b }}} },
427   { {{{ 0x966abfc0, 0x8a385bf4, 0xf081a640, 0x55e5e8bc,
428         0xee26f5ff, 0x835dff85, 0xe509e1ea, 0x4927e622 }}},
429     {{{ 0x352334b0, 0x164c8dbc, 0xa3fea31f, 0xcac1ad63,
430         0x682fd457, 0x9b87a676, 0x1a53145f, 0x75f382ff }}} },
431   { {{{ 0xc3efcb46, 0x16b944f5, 0x68cb184c, 0x1fb55714,
432         0x9ccf2dc8, 0xf1c2b116, 0x808283d8, 0x7417e00f }}},
433     {{{ 0x930199ba, 0x1ea67a22, 0x718990d8, 0x9fbaf765,
434         0x8f3d5d57, 0x231fc664, 0xe5853194, 0x38141a19 }}} },
435   { {{{ 0x2f81290d, 0xb9f00390, 0x04a9ca6c, 0x44877827,
436         0xe1dbdd65, 0x65d7f9b9, 0xf7c6698a, 0x7133424c }}},
437     {{{ 0xa7cd250f, 0x604cfb3c, 0x5acc18f3, 0x460c3c4b,
438         0xb518e3eb, 0xa53e50e0, 0x98a40196, 0x2b4b9267 }}} },
439   { {{{ 0xc5dbd06c, 0x591b0672, 0xaa1eeb65, 0x10d43dca,
440         0xcd2517af, 0x420cdef8, 0x0b695a8a, 0x513a307e }}},
441     {{{ 0x66503215, 0xee9d6a7b, 0x088fd9a4, 0xdea58720,
442         0x973afe12, 0x8f3cbbea, 0x872f2538, 0x005c2350 }}} },
443   { {{{ 0x35af3291, 0xe5024b70, 0x4f5e669a, 0x1d3eec2d,
444         0x6e79d539, 0xc1f6d766, 0x795b5248, 0x34ec043f }}},
445     {{{ 0x400960b6, 0xb2763511, 0x29e57df0, 0xff7a3d84,
446         0x1666c1f1, 0xaeac7792, 0x66084bc0, 0x72426e97 }}} },
447   { {{{ 0x44f826ca, 0x5b1c3199, 0x790aa408, 0x68b00b73,
448         0x69e9b92b, 0xaf0984b4, 0x3ffe9093, 0x5fe6736f }}},
449     {{{ 0xffd49312, 0xd67f2889, 0x5cb9ed21, 0x3520d747,
450         0x3c65a606, 0x94f893b1, 0x2d65496f, 0x2fee5e8c }}} }
451 };
452
453 /**
454  * @brief       X  = k * G
455  *
456  * @param K     scalar k
457  *
458  * Return -1 on error.
459  * Return 0 on success.
460  */
461 static void
462 compute_kG_25519 (ac *X, const bn256 *K)
463 {
464   ptc Q[1];
465   int i;
466
467   /* identity element */
468   memset (Q, 0, sizeof (ptc));
469   Q->y->word[0] = 1;
470   Q->z->word[0] = 1;
471
472   for (i = 20; i >= 0; i--)
473     {
474       int k0, k1, k2;
475
476       k0 = ((K->word[0] >> i) & 1)
477         | (i < 1 ? ((K->word[1] >> 30) & 2)
478            : (((K->word[2] >> (i-1)) & 1) << 1))
479         | (i < 2 ? ((K->word[3] >> (i+28)) & 4)
480            : (((K->word[4] >> (i-2)) & 1) << 2))
481         | (i < 3 ? ((K->word[5] >> (i+26)) & 8)
482            : (((K->word[6] >> (i-3)) & 1) << 3));
483
484       k1 = (i < 11 ? ((K->word[0] >> (i+21)) & 1)
485             : ((K->word[1] >> (i-11)) & 1))
486         | (i < 12 ? ((K->word[2] >> (i+19)) & 2)
487            : (((K->word[3] >> (i-12)) & 1) << 1))
488         | (i < 13 ? ((K->word[4] >> (i+17)) & 4)
489            : (((K->word[5] >> (i-13)) & 1) << 2))
490         | (i < 14 ? ((K->word[6] >> (i+15)) & 8)
491            : (((K->word[7] >> (i-14)) & 1) << 3));
492
493       k2 = ((K->word[1] >> (i+10)) & 1)
494         | ((K->word[3] >> (i+8)) & 2)
495         | ((K->word[5] >> (i+6)) & 4)
496         | ((K->word[7] >> (i+4)) & 8);
497
498       point_double (Q, Q);
499       point_add (Q, Q, &precomputed_KG[k0]);
500       point_add (Q, Q, &precomputed_2E_KG[k1]);
501       point_add (Q, Q, &precomputed_4E_KG[k2]);
502     }
503
504   point_ptc_to_ac (X, Q);
505 }
506
507
508 #define BN416_WORDS 13
509 #define BN128_WORDS 4
510
511 /* M: The order of the generator G.  */
512 static const bn256 M[1] = {
513   {{  0x5CF5D3ED, 0x5812631A, 0xA2F79CD6, 0x14DEF9DE,
514       0x00000000, 0x00000000, 0x00000000, 0x10000000  }}
515 };
516
517 #define C ((const uint32_t *)M)
518
519 static void
520 bnX_mul_C (uint32_t *r, const uint32_t *q, int q_size)
521 {
522   int i, j, k;
523   int i_beg, i_end;
524   uint32_t r0, r1, r2;
525
526   r0 = r1 = r2 = 0;
527   for (k = 0; k <= q_size + BN128_WORDS - 2; k++)
528     {
529       if (q_size < BN128_WORDS)
530         if (k < q_size)
531           {
532             i_beg = 0;
533             i_end = k;
534           }
535         else
536           {
537             i_beg = k - q_size + 1;
538             i_end = k;
539             if (i_end > BN128_WORDS - 1)
540               i_end = BN128_WORDS - 1;
541           }
542       else
543         if (k < BN128_WORDS)
544           {
545             i_beg = 0;
546             i_end = k;
547           }
548         else
549           {
550             i_beg = k - BN128_WORDS + 1;
551             i_end = k;
552             if (i_end > q_size - 1)
553               i_end = q_size - 1;
554           }
555
556       for (i = i_beg; i <= i_end; i++)
557         {
558           uint64_t uv;
559           uint32_t u, v;
560           uint32_t carry;
561
562           j = k - i;
563           if (q_size < BN128_WORDS)
564             uv = ((uint64_t )q[j])*((uint64_t )C[i]);
565           else
566             uv = ((uint64_t )q[i])*((uint64_t )C[j]);
567           v = uv;
568           u = (uv >> 32);
569           r0 += v;
570           carry = (r0 < v);
571           r1 += carry;
572           carry = (r1 < carry);
573           r1 += u;
574           carry += (r1 < u);
575           r2 += carry;
576         }
577
578       r[k] = r0;
579       r0 = r1;
580       r1 = r2;
581       r2 = 0;
582     }
583
584   r[k] = r0;
585 }
586
587 /**
588  * @brief R = A mod M (using M=2^252+C) (Barret reduction)
589  *
590  * See HAC 14.47.
591  */
592 static void
593 mod_reduce_M (bn256 *R, const bn512 *A)
594 {
595   uint32_t q[BN256_WORDS+1];
596   uint32_t tmp[BN416_WORDS];
597   bn256 r[1];
598   uint32_t carry, next_carry;
599   int i;
600 #define borrow carry
601
602   q[8] = A->word[15]>>28;
603   carry = A->word[15] & 0x0fffffff;
604   for (i = BN256_WORDS - 1; i >= 0; i--)
605     {
606       next_carry = A->word[i+7] & 0x0fffffff;
607       q[i] = (A->word[i+7] >> 28) | (carry << 4);
608       carry = next_carry;
609     }
610   memcpy (R, A, sizeof (bn256));
611   R->word[7] &= 0x0fffffff;
612
613   /* Q_size: 9 */
614   bnX_mul_C (tmp, q, 9); /* TMP = Q*C */
615   /* Q = tmp / 2^252 */
616   carry = tmp[12] & 0x0fffffff;
617   for (i = 4; i >= 0; i--)
618     {
619       next_carry = tmp[i+7] & 0x0fffffff;
620       q[i] = (tmp[i+7] >> 28) | (carry << 4);
621       carry = next_carry;
622     }
623   /* R' = tmp % 2^252 */
624   memcpy (r, tmp, sizeof (bn256));
625   r->word[7] &= 0x0fffffff;
626   /* R -= R' */
627   borrow = bn256_sub (R, R, r);
628   if (borrow)
629     bn256_add (R, R, M);
630   else
631     bn256_add ((bn256 *)tmp, R, M);
632
633   /* Q_size: 5 */
634   bnX_mul_C (tmp, q, 5); /* TMP = Q*C */
635   carry = tmp[8] & 0x0fffffff;
636   q[0] = (tmp[7] >> 28) | (carry << 4);
637   /* R' = tmp % 2^252 */
638   memcpy (r, tmp, sizeof (bn256));
639   r->word[7] &= 0x0fffffff;
640   /* R += R' */
641   bn256_add (R, R, r);
642   borrow = bn256_sub (R, R, M);
643   if (borrow)
644     bn256_add (R, R, M);
645   else
646     bn256_add ((bn256 *)tmp, R, M);
647
648   /* Q_size: 1 */
649   bnX_mul_C (tmp, q, 1); /* TMP = Q*C */
650   /* R' = tmp % 2^252 */
651   memset (((uint8_t *)r)+(sizeof (uint32_t)*5), 0, sizeof (uint32_t)*3);
652   memcpy (r, tmp, sizeof (uint32_t)*5);
653   /* R -= R' */
654   borrow = bn256_sub (R, R, r);
655   if (borrow)
656     bn256_add (R, R, M);
657   else
658     bn256_add ((bn256 *)tmp, R, M);
659 #undef borrow
660 }
661
662
663 int
664 eddsa_sign_25519 (const uint8_t *input, size_t ilen, uint32_t *out,
665                   const bn256 *a, const uint8_t *seed, const bn256 *pk)
666 {
667   bn256 *r, *s;
668   sha512_context ctx;
669   uint8_t hash[64];
670   bn256 tmp[1];
671   ac R[1];
672   uint32_t carry, borrow;
673
674   r = (bn256 *)out;
675   s = (bn256 *)(out+(32/4));
676
677   sha512_start (&ctx);
678   sha512_update (&ctx, seed, sizeof (bn256)); /* It's upper half of the hash */
679   sha512_update (&ctx, input, ilen);
680   sha512_finish (&ctx, hash);
681
682   mod_reduce_M (r, (bn512 *)hash);
683   compute_kG_25519 (R, r);
684
685   /* EdDSA encoding.  */
686   memcpy (tmp, R->y, sizeof (bn256));
687   tmp->word[7] ^= mod25519_is_neg (R->x) * 0x80000000;
688
689   sha512_start (&ctx);
690   sha512_update (&ctx, (uint8_t *)tmp, sizeof (bn256));
691   sha512_update (&ctx, (uint8_t *)pk, sizeof (bn256));
692   sha512_update (&ctx, input, ilen);
693   sha512_finish (&ctx, (uint8_t *)hash);
694
695   mod_reduce_M (s, (bn512 *)hash);
696   bn256_mul ((bn512 *)hash, s, a);
697   mod_reduce_M (s, (bn512 *)hash);
698   carry = bn256_add (s, s, r);
699   borrow = bn256_sub (s, s, M);
700
701   memcpy (r, tmp, sizeof (bn256));
702
703   if ((borrow && !carry))
704     bn256_add (s, s, M);
705   else
706     bn256_add (tmp, s, M);
707
708   return 0;
709 }
710
711 void
712 eddsa_public_key_25519 (bn256 *pk, const bn256 *a)
713 {
714   ac R[1];
715   ptc X[1];
716   bn256 a0[1];
717
718   bn256_shift (a0, a, -3);
719   compute_kG_25519 (R, a0);
720   memcpy (X, R, sizeof (ac));
721   memset (X->z, 0, sizeof (bn256));
722   X->z->word[0] = 1;
723   point_double (X, X);
724   point_double (X, X);
725   point_double (X, X);
726   point_ptc_to_ac (R, X);
727   /* EdDSA encoding.  */
728   memcpy (pk, R->y, sizeof (bn256));
729   pk->word[7] ^= mod25519_is_neg (R->x) * 0x80000000;
730 }
731
732
733 uint8_t *
734 eddsa_compute_public_25519 (const uint8_t *kd)
735 {
736   uint8_t *p0;
737   const bn256 *a = (const bn256 *)kd;
738
739   p0 = (uint8_t *)malloc (sizeof (bn256));
740   if (p0 == NULL)
741     return NULL;
742
743   eddsa_public_key_25519 ((bn256 *)p0, a);
744   return p0;
745 }
746
747
748 #if 0
749 /**
750  * check if P is on the curve.
751  *
752  * Return -1 on error.
753  * Return 0 on success.
754  */
755 static int
756 point_is_on_the_curve (const ac *P)
757 {
758   bn256 s[1], t[1];
759
760   /* Twisted Edwards curve: a*x^2 + y^2 = 1 + d*x^2*y^2 */
761 }
762
763 int
764 compute_kP_25519 (ac *X, const bn256 *K, const ac *P);
765 #endif
766
767 #ifdef PRINT_OUT_TABLE
768 static const ptc G[1] = {{
769   {{{ 0x8f25d51a, 0xc9562d60, 0x9525a7b2, 0x692cc760,
770       0xfdd6dc5c, 0xc0a4e231, 0xcd6e53fe, 0x216936d3 }}},
771   {{{ 0x66666658, 0x66666666, 0x66666666, 0x66666666,
772       0x66666666, 0x66666666, 0x66666666, 0x66666666 }}},
773   {{{ 1, 0, 0, 0, 0, 0, 0, 0 }}},
774 }};
775
776 #include <stdio.h>
777
778 #ifdef TESTING_EDDSA
779 static void
780 print_bn256 (const bn256 *X)
781 {
782   int i;
783
784   for (i = 7; i >= 0; i--)
785     printf ("%08x", X->word[i]);
786   puts ("");
787 }
788 #endif
789
790 #if 0
791 static void
792 print_point (const ac *X)
793 {
794   int i;
795
796 #ifdef PRINT_OUT_TABLE_AS_C
797   fputs ("  { {{{ ", stdout);
798   for (i = 0; i < 4; i++)
799     printf ("0x%08x, ", X->x->word[i]);
800   fputs ("\n        ", stdout);
801   for (; i < 7; i++)
802     printf ("0x%08x, ", X->x->word[i]);
803   printf ("0x%08x }}},\n", X->x->word[i]);
804   fputs ("    {{{ ", stdout);
805   for (i = 0; i < 4; i++)
806     printf ("0x%08x, ", X->y->word[i]);
807   fputs ("\n        ", stdout);
808   for (; i < 7; i++)
809     printf ("0x%08x, ", X->y->word[i]);
810   printf ("0x%08x }}} },\n", X->y->word[i]);
811 #else
812   puts ("--");
813   for (i = 7; i >= 0; i--)
814     printf ("%08x", X->x->word[i]);
815   puts ("");
816   for (i = 7; i >= 0; i--)
817     printf ("%08x", X->y->word[i]);
818   puts ("");
819   puts ("--");
820 #endif
821 }
822
823 static void
824 print_point_ptc (const ptc *X)
825 {
826   int i;
827
828   puts ("---");
829   for (i = 7; i >= 0; i--)
830     printf ("%08x", X->x->word[i]);
831   puts ("");
832   for (i = 7; i >= 0; i--)
833     printf ("%08x", X->y->word[i]);
834   puts ("");
835   for (i = 7; i >= 0; i--)
836     printf ("%08x", X->z->word[i]);
837   puts ("");
838   puts ("---");
839 }
840 #endif
841
842
843 #ifndef TESTING_EDDSA
844 static void power_2 (ac *A, ptc *a, int N)
845 {
846   int i;
847
848   for (i = 0; i < N; i++)
849     ed_double_25638 (a, a);
850   ptc_to_ac_25519 (A, a);
851 }
852
853 static void print_table (ac *a0001, ac *a0010, ac *a0100, ac *a1000)
854 {
855   int i;
856   ptc a[1];
857   ac x[1];
858
859   for (i = 1; i < 16; i++)
860     {
861       /* A := Identity Element  */
862       memset (a, 0, sizeof (ptc));
863       a->y->word[0] = 1;
864       a->z->word[0] = 1;
865
866       if ((i & 1))
867         ed_add_25638 (a, a, a0001);
868       if ((i & 2))
869         ed_add_25638 (a, a, a0010);
870       if ((i & 4))
871         ed_add_25638 (a, a, a0100);
872       if ((i & 8))
873         ed_add_25638 (a, a, a1000);
874
875       ptc_to_ac_25519 (x, a);
876       print_point (x);
877     }
878
879   fputs ("\n", stdout);
880 }
881
882 static void compute_and_print_table (ac *a0001, ac *a0010, ac *a0100, ac *a1000)
883 {
884   ptc a[1];
885
886   memcpy (a, a0001, sizeof (ac));
887   memset (a->z, 0, sizeof (bn256));
888   a->z->word[0] = 1;
889   power_2 (a0010, a, 63);
890   power_2 (a0100, a, 63);
891   power_2 (a1000, a, 63);
892   print_table (a0001, a0010, a0100, a1000);
893 }
894 #endif
895
896 int
897 main (int argc, char *argv[])
898 {
899 #ifdef TESTING_EDDSA
900   uint8_t hash[64];
901   bn256 a[1];
902   uint8_t r_s[64];
903   bn256 pk[1];
904   bn256 *r, *s;
905
906   const bn256 sk[1] = {
907     {{ 0x9db1619d, 0x605afdef, 0xf44a84ba, 0xc42cec92,
908        0x69c54944, 0x1969327b, 0x03ac3b70, 0x607fae1c }} };
909
910   const bn256 r_expected[1] = {
911     {{ 0x004356e5, 0x72ac60c3, 0xcce28690, 0x8a826e80,
912        0x1e7f8784, 0x74d9e5b8, 0x65e073d8, 0x55014922 }} };
913
914   const bn256 s_expected[1] = {
915     {{ 0x1582b85f, 0xac3ba390, 0x70391ec6, 0x6bb4f91c,
916        0xf0f55bd2, 0x24be5b59, 0x43415165, 0x0b107a8e }} };
917
918   r = (bn256 *)r_s;
919   s = (bn256 *)(r_s+32);
920
921   sha512 ((uint8_t *)sk, sizeof (bn256), hash);
922   hash[0] &= 248;
923   hash[31] &= 127;
924   hash[31] |= 64;
925   memcpy (a, hash, sizeof (bn256));
926
927   eddsa_public_key_25519 (pk, a);
928   eddsa_sign_25519 ((const uint8_t *)"", 0, r_s, a, hash+32, pk);
929
930   if (memcmp (r, r_expected, sizeof (bn256)) != 0
931       || memcmp (s, s_expected, sizeof (bn256)) != 0)
932     {
933       print_bn256 (r);
934       print_bn256 (s);
935       return 1;
936     }
937 #else
938   ac a0001[1], a0010[1], a0100[1], a1000[1];
939   ptc a[1];
940
941   memcpy (a, G, sizeof (ptc));
942   ptc_to_ac_25519 (a0001, a);
943   compute_and_print_table (a0001, a0010, a0100, a1000);
944
945   memcpy (a, a0001, sizeof (ac));
946   memset (a->z, 0, sizeof (bn256));
947   a->z->word[0] = 1;
948   power_2 (a0001, a, 21);
949   compute_and_print_table (a0001, a0010, a0100, a1000);
950
951   memcpy (a, a0001, sizeof (ac));
952   memset (a->z, 0, sizeof (bn256));
953   a->z->word[0] = 1;
954   power_2 (a0001, a, 21);
955   compute_and_print_table (a0001, a0010, a0100, a1000);
956 #endif
957
958   return 0;
959 }
960 #endif