1 /* 2 * Copyright 2013 The Android Open Source Project 3 * 4 * Redistribution and use in source and binary forms, with or without 5 * modification, are permitted provided that the following conditions are met: 6 * * Redistributions of source code must retain the above copyright 7 * notice, this list of conditions and the following disclaimer. 8 * * Redistributions in binary form must reproduce the above copyright 9 * notice, this list of conditions and the following disclaimer in the 10 * documentation and/or other materials provided with the distribution. 11 * * Neither the name of Google Inc. nor the names of its contributors may 12 * be used to endorse or promote products derived from this software 13 * without specific prior written permission. 14 * 15 * THIS SOFTWARE IS PROVIDED BY Google Inc. ``AS IS'' AND ANY EXPRESS OR 16 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 17 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO 18 * EVENT SHALL Google Inc. BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 19 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 20 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 21 * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 22 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 23 * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF 24 * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25 */ 26 27 // This is an implementation of the P256 elliptic curve group. It's written to 28 // be portable 32-bit, although it's still constant-time. 29 // 30 // WARNING: Implementing these functions in a constant-time manner is far from 31 // obvious. Be careful when touching this code. 32 // 33 // See http://www.imperialviolet.org/2010/12/04/ecc.html ([1]) for background. 34 35 #include <stdint.h> 36 #include <stdio.h> 37 38 #include <string.h> 39 #include <stdlib.h> 40 41 #include "constrainedcrypto/p256.h" 42 43 typedef uint8_t u8; 44 typedef uint32_t u32; 45 typedef int32_t s32; 46 typedef uint64_t u64; 47 48 /* Our field elements are represented as nine 32-bit limbs. 49 * 50 * The value of an felem (field element) is: 51 * x[0] + (x[1] * 2**29) + (x[2] * 2**57) + ... + (x[8] * 2**228) 52 * 53 * That is, each limb is alternately 29 or 28-bits wide in little-endian 54 * order. 55 * 56 * This means that an felem hits 2**257, rather than 2**256 as we would like. A 57 * 28, 29, ... pattern would cause us to hit 2**256, but that causes problems 58 * when multiplying as terms end up one bit short of a limb which would require 59 * much bit-shifting to correct. 60 * 61 * Finally, the values stored in an felem are in Montgomery form. So the value 62 * |y| is stored as (y*R) mod p, where p is the P-256 prime and R is 2**257. 63 */ 64 typedef u32 limb; 65 #define NLIMBS 9 66 typedef limb felem[NLIMBS]; 67 68 static const limb kBottom28Bits = 0xfffffff; 69 static const limb kBottom29Bits = 0x1fffffff; 70 71 /* kOne is the number 1 as an felem. It's 2**257 mod p split up into 29 and 72 * 28-bit words. */ 73 static const felem kOne = { 74 2, 0, 0, 0xffff800, 75 0x1fffffff, 0xfffffff, 0x1fbfffff, 0x1ffffff, 76 0 77 }; 78 static const felem kZero = {0}; 79 static const felem kP = { 80 0x1fffffff, 0xfffffff, 0x1fffffff, 0x3ff, 81 0, 0, 0x200000, 0xf000000, 82 0xfffffff 83 }; 84 static const felem k2P = { 85 0x1ffffffe, 0xfffffff, 0x1fffffff, 0x7ff, 86 0, 0, 0x400000, 0xe000000, 87 0x1fffffff 88 }; 89 /* kPrecomputed contains precomputed values to aid the calculation of scalar 90 * multiples of the base point, G. It's actually two, equal length, tables 91 * concatenated. 92 * 93 * The first table contains (x,y) felem pairs for 16 multiples of the base 94 * point, G. 95 * 96 * Index | Index (binary) | Value 97 * 0 | 0000 | 0G (all zeros, omitted) 98 * 1 | 0001 | G 99 * 2 | 0010 | 2**64G 100 * 3 | 0011 | 2**64G + G 101 * 4 | 0100 | 2**128G 102 * 5 | 0101 | 2**128G + G 103 * 6 | 0110 | 2**128G + 2**64G 104 * 7 | 0111 | 2**128G + 2**64G + G 105 * 8 | 1000 | 2**192G 106 * 9 | 1001 | 2**192G + G 107 * 10 | 1010 | 2**192G + 2**64G 108 * 11 | 1011 | 2**192G + 2**64G + G 109 * 12 | 1100 | 2**192G + 2**128G 110 * 13 | 1101 | 2**192G + 2**128G + G 111 * 14 | 1110 | 2**192G + 2**128G + 2**64G 112 * 15 | 1111 | 2**192G + 2**128G + 2**64G + G 113 * 114 * The second table follows the same style, but the terms are 2**32G, 115 * 2**96G, 2**160G, 2**224G. 116 * 117 * This is ~2KB of data. */ 118 static const limb kPrecomputed[NLIMBS * 2 * 15 * 2] = { 119 0x11522878, 0xe730d41, 0xdb60179, 0x4afe2ff, 0x12883add, 0xcaddd88, 0x119e7edc, 0xd4a6eab, 0x3120bee, 120 0x1d2aac15, 0xf25357c, 0x19e45cdd, 0x5c721d0, 0x1992c5a5, 0xa237487, 0x154ba21, 0x14b10bb, 0xae3fe3, 121 0xd41a576, 0x922fc51, 0x234994f, 0x60b60d3, 0x164586ae, 0xce95f18, 0x1fe49073, 0x3fa36cc, 0x5ebcd2c, 122 0xb402f2f, 0x15c70bf, 0x1561925c, 0x5a26704, 0xda91e90, 0xcdc1c7f, 0x1ea12446, 0xe1ade1e, 0xec91f22, 123 0x26f7778, 0x566847e, 0xa0bec9e, 0x234f453, 0x1a31f21a, 0xd85e75c, 0x56c7109, 0xa267a00, 0xb57c050, 124 0x98fb57, 0xaa837cc, 0x60c0792, 0xcfa5e19, 0x61bab9e, 0x589e39b, 0xa324c5, 0x7d6dee7, 0x2976e4b, 125 0x1fc4124a, 0xa8c244b, 0x1ce86762, 0xcd61c7e, 0x1831c8e0, 0x75774e1, 0x1d96a5a9, 0x843a649, 0xc3ab0fa, 126 0x6e2e7d5, 0x7673a2a, 0x178b65e8, 0x4003e9b, 0x1a1f11c2, 0x7816ea, 0xf643e11, 0x58c43df, 0xf423fc2, 127 0x19633ffa, 0x891f2b2, 0x123c231c, 0x46add8c, 0x54700dd, 0x59e2b17, 0x172db40f, 0x83e277d, 0xb0dd609, 128 0xfd1da12, 0x35c6e52, 0x19ede20c, 0xd19e0c0, 0x97d0f40, 0xb015b19, 0x449e3f5, 0xe10c9e, 0x33ab581, 129 0x56a67ab, 0x577734d, 0x1dddc062, 0xc57b10d, 0x149b39d, 0x26a9e7b, 0xc35df9f, 0x48764cd, 0x76dbcca, 130 0xca4b366, 0xe9303ab, 0x1a7480e7, 0x57e9e81, 0x1e13eb50, 0xf466cf3, 0x6f16b20, 0x4ba3173, 0xc168c33, 131 0x15cb5439, 0x6a38e11, 0x73658bd, 0xb29564f, 0x3f6dc5b, 0x53b97e, 0x1322c4c0, 0x65dd7ff, 0x3a1e4f6, 132 0x14e614aa, 0x9246317, 0x1bc83aca, 0xad97eed, 0xd38ce4a, 0xf82b006, 0x341f077, 0xa6add89, 0x4894acd, 133 0x9f162d5, 0xf8410ef, 0x1b266a56, 0xd7f223, 0x3e0cb92, 0xe39b672, 0x6a2901a, 0x69a8556, 0x7e7c0, 134 0x9b7d8d3, 0x309a80, 0x1ad05f7f, 0xc2fb5dd, 0xcbfd41d, 0x9ceb638, 0x1051825c, 0xda0cf5b, 0x812e881, 135 0x6f35669, 0x6a56f2c, 0x1df8d184, 0x345820, 0x1477d477, 0x1645db1, 0xbe80c51, 0xc22be3e, 0xe35e65a, 136 0x1aeb7aa0, 0xc375315, 0xf67bc99, 0x7fdd7b9, 0x191fc1be, 0x61235d, 0x2c184e9, 0x1c5a839, 0x47a1e26, 137 0xb7cb456, 0x93e225d, 0x14f3c6ed, 0xccc1ac9, 0x17fe37f3, 0x4988989, 0x1a90c502, 0x2f32042, 0xa17769b, 138 0xafd8c7c, 0x8191c6e, 0x1dcdb237, 0x16200c0, 0x107b32a1, 0x66c08db, 0x10d06a02, 0x3fc93, 0x5620023, 139 0x16722b27, 0x68b5c59, 0x270fcfc, 0xfad0ecc, 0xe5de1c2, 0xeab466b, 0x2fc513c, 0x407f75c, 0xbaab133, 140 0x9705fe9, 0xb88b8e7, 0x734c993, 0x1e1ff8f, 0x19156970, 0xabd0f00, 0x10469ea7, 0x3293ac0, 0xcdc98aa, 141 0x1d843fd, 0xe14bfe8, 0x15be825f, 0x8b5212, 0xeb3fb67, 0x81cbd29, 0xbc62f16, 0x2b6fcc7, 0xf5a4e29, 142 0x13560b66, 0xc0b6ac2, 0x51ae690, 0xd41e271, 0xf3e9bd4, 0x1d70aab, 0x1029f72, 0x73e1c35, 0xee70fbc, 143 0xad81baf, 0x9ecc49a, 0x86c741e, 0xfe6be30, 0x176752e7, 0x23d416, 0x1f83de85, 0x27de188, 0x66f70b8, 144 0x181cd51f, 0x96b6e4c, 0x188f2335, 0xa5df759, 0x17a77eb6, 0xfeb0e73, 0x154ae914, 0x2f3ec51, 0x3826b59, 145 0xb91f17d, 0x1c72949, 0x1362bf0a, 0xe23fddf, 0xa5614b0, 0xf7d8f, 0x79061, 0x823d9d2, 0x8213f39, 146 0x1128ae0b, 0xd095d05, 0xb85c0c2, 0x1ecb2ef, 0x24ddc84, 0xe35e901, 0x18411a4a, 0xf5ddc3d, 0x3786689, 147 0x52260e8, 0x5ae3564, 0x542b10d, 0x8d93a45, 0x19952aa4, 0x996cc41, 0x1051a729, 0x4be3499, 0x52b23aa, 148 0x109f307e, 0x6f5b6bb, 0x1f84e1e7, 0x77a0cfa, 0x10c4df3f, 0x25a02ea, 0xb048035, 0xe31de66, 0xc6ecaa3, 149 0x28ea335, 0x2886024, 0x1372f020, 0xf55d35, 0x15e4684c, 0xf2a9e17, 0x1a4a7529, 0xcb7beb1, 0xb2a78a1, 150 0x1ab21f1f, 0x6361ccf, 0x6c9179d, 0xb135627, 0x1267b974, 0x4408bad, 0x1cbff658, 0xe3d6511, 0xc7d76f, 151 0x1cc7a69, 0xe7ee31b, 0x54fab4f, 0x2b914f, 0x1ad27a30, 0xcd3579e, 0xc50124c, 0x50daa90, 0xb13f72, 152 0xb06aa75, 0x70f5cc6, 0x1649e5aa, 0x84a5312, 0x329043c, 0x41c4011, 0x13d32411, 0xb04a838, 0xd760d2d, 153 0x1713b532, 0xbaa0c03, 0x84022ab, 0x6bcf5c1, 0x2f45379, 0x18ae070, 0x18c9e11e, 0x20bca9a, 0x66f496b, 154 0x3eef294, 0x67500d2, 0xd7f613c, 0x2dbbeb, 0xb741038, 0xe04133f, 0x1582968d, 0xbe985f7, 0x1acbc1a, 155 0x1a6a939f, 0x33e50f6, 0xd665ed4, 0xb4b7bd6, 0x1e5a3799, 0x6b33847, 0x17fa56ff, 0x65ef930, 0x21dc4a, 156 0x2b37659, 0x450fe17, 0xb357b65, 0xdf5efac, 0x15397bef, 0x9d35a7f, 0x112ac15f, 0x624e62e, 0xa90ae2f, 157 0x107eecd2, 0x1f69bbe, 0x77d6bce, 0x5741394, 0x13c684fc, 0x950c910, 0x725522b, 0xdc78583, 0x40eeabb, 158 0x1fde328a, 0xbd61d96, 0xd28c387, 0x9e77d89, 0x12550c40, 0x759cb7d, 0x367ef34, 0xae2a960, 0x91b8bdc, 159 0x93462a9, 0xf469ef, 0xb2e9aef, 0xd2ca771, 0x54e1f42, 0x7aaa49, 0x6316abb, 0x2413c8e, 0x5425bf9, 160 0x1bed3e3a, 0xf272274, 0x1f5e7326, 0x6416517, 0xea27072, 0x9cedea7, 0x6e7633, 0x7c91952, 0xd806dce, 161 0x8e2a7e1, 0xe421e1a, 0x418c9e1, 0x1dbc890, 0x1b395c36, 0xa1dc175, 0x1dc4ef73, 0x8956f34, 0xe4b5cf2, 162 0x1b0d3a18, 0x3194a36, 0x6c2641f, 0xe44124c, 0xa2f4eaa, 0xa8c25ba, 0xf927ed7, 0x627b614, 0x7371cca, 163 0xba16694, 0x417bc03, 0x7c0a7e3, 0x9c35c19, 0x1168a205, 0x8b6b00d, 0x10e3edc9, 0x9c19bf2, 0x5882229, 164 0x1b2b4162, 0xa5cef1a, 0x1543622b, 0x9bd433e, 0x364e04d, 0x7480792, 0x5c9b5b3, 0xe85ff25, 0x408ef57, 165 0x1814cfa4, 0x121b41b, 0xd248a0f, 0x3b05222, 0x39bb16a, 0xc75966d, 0xa038113, 0xa4a1769, 0x11fbc6c, 166 0x917e50e, 0xeec3da8, 0x169d6eac, 0x10c1699, 0xa416153, 0xf724912, 0x15cd60b7, 0x4acbad9, 0x5efc5fa, 167 0xf150ed7, 0x122b51, 0x1104b40a, 0xcb7f442, 0xfbb28ff, 0x6ac53ca, 0x196142cc, 0x7bf0fa9, 0x957651, 168 0x4e0f215, 0xed439f8, 0x3f46bd5, 0x5ace82f, 0x110916b6, 0x6db078, 0xffd7d57, 0xf2ecaac, 0xca86dec, 169 0x15d6b2da, 0x965ecc9, 0x1c92b4c2, 0x1f3811, 0x1cb080f5, 0x2d8b804, 0x19d1c12d, 0xf20bd46, 0x1951fa7, 170 0xa3656c3, 0x523a425, 0xfcd0692, 0xd44ddc8, 0x131f0f5b, 0xaf80e4a, 0xcd9fc74, 0x99bb618, 0x2db944c, 171 0xa673090, 0x1c210e1, 0x178c8d23, 0x1474383, 0x10b8743d, 0x985a55b, 0x2e74779, 0x576138, 0x9587927, 172 0x133130fa, 0xbe05516, 0x9f4d619, 0xbb62570, 0x99ec591, 0xd9468fe, 0x1d07782d, 0xfc72e0b, 0x701b298, 173 0x1863863b, 0x85954b8, 0x121a0c36, 0x9e7fedf, 0xf64b429, 0x9b9d71e, 0x14e2f5d8, 0xf858d3a, 0x942eea8, 174 0xda5b765, 0x6edafff, 0xa9d18cc, 0xc65e4ba, 0x1c747e86, 0xe4ea915, 0x1981d7a1, 0x8395659, 0x52ed4e2, 175 0x87d43b7, 0x37ab11b, 0x19d292ce, 0xf8d4692, 0x18c3053f, 0x8863e13, 0x4c146c0, 0x6bdf55a, 0x4e4457d, 176 0x16152289, 0xac78ec2, 0x1a59c5a2, 0x2028b97, 0x71c2d01, 0x295851f, 0x404747b, 0x878558d, 0x7d29aa4, 177 0x13d8341f, 0x8daefd7, 0x139c972d, 0x6b7ea75, 0xd4a9dde, 0xff163d8, 0x81d55d7, 0xa5bef68, 0xb7b30d8, 178 0xbe73d6f, 0xaa88141, 0xd976c81, 0x7e7a9cc, 0x18beb771, 0xd773cbd, 0x13f51951, 0x9d0c177, 0x1c49a78, 179 }; 180 181 182 /* Field element operations: */ 183 184 /* NON_ZERO_TO_ALL_ONES returns: 185 * 0xffffffff for 0 < x <= 2**31 186 * 0 for x == 0 or x > 2**31. 187 * 188 * x must be a u32 or an equivalent type such as limb. */ 189 #define NON_ZERO_TO_ALL_ONES(x) ((((u32)(x) - 1) >> 31) - 1) 190 191 /* felem_reduce_carry adds a multiple of p in order to cancel |carry|, 192 * which is a term at 2**257. 193 * 194 * On entry: carry < 2**3, inout[0,2,...] < 2**29, inout[1,3,...] < 2**28. 195 * On exit: inout[0,2,..] < 2**30, inout[1,3,...] < 2**29. */ 196 static void felem_reduce_carry(felem inout, limb carry) { 197 const u32 carry_mask = NON_ZERO_TO_ALL_ONES(carry); 198 199 inout[0] += carry << 1; 200 inout[3] += 0x10000000 & carry_mask; 201 /* carry < 2**3 thus (carry << 11) < 2**14 and we added 2**28 in the 202 * previous line therefore this doesn't underflow. */ 203 inout[3] -= carry << 11; 204 inout[4] += (0x20000000 - 1) & carry_mask; 205 inout[5] += (0x10000000 - 1) & carry_mask; 206 inout[6] += (0x20000000 - 1) & carry_mask; 207 inout[6] -= carry << 22; 208 /* This may underflow if carry is non-zero but, if so, we'll fix it in the 209 * next line. */ 210 inout[7] -= 1 & carry_mask; 211 inout[7] += carry << 25; 212 } 213 214 /* felem_sum sets out = in+in2. 215 * 216 * On entry, in[i]+in2[i] must not overflow a 32-bit word. 217 * On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29 */ 218 static void felem_sum(felem out, const felem in, const felem in2) { 219 limb carry = 0; 220 unsigned i; 221 222 for (i = 0;; i++) { 223 out[i] = in[i] + in2[i]; 224 out[i] += carry; 225 carry = out[i] >> 29; 226 out[i] &= kBottom29Bits; 227 228 i++; 229 if (i == NLIMBS) 230 break; 231 232 out[i] = in[i] + in2[i]; 233 out[i] += carry; 234 carry = out[i] >> 28; 235 out[i] &= kBottom28Bits; 236 } 237 238 felem_reduce_carry(out, carry); 239 } 240 241 #define two31m3 (((limb)1) << 31) - (((limb)1) << 3) 242 #define two30m2 (((limb)1) << 30) - (((limb)1) << 2) 243 #define two30p13m2 (((limb)1) << 30) + (((limb)1) << 13) - (((limb)1) << 2) 244 #define two31m2 (((limb)1) << 31) - (((limb)1) << 2) 245 #define two31p24m2 (((limb)1) << 31) + (((limb)1) << 24) - (((limb)1) << 2) 246 #define two30m27m2 (((limb)1) << 30) - (((limb)1) << 27) - (((limb)1) << 2) 247 248 /* zero31 is 0 mod p. */ 249 static const felem zero31 = { two31m3, two30m2, two31m2, two30p13m2, two31m2, two30m2, two31p24m2, two30m27m2, two31m2 }; 250 251 /* felem_diff sets out = in-in2. 252 * 253 * On entry: in[0,2,...] < 2**30, in[1,3,...] < 2**29 and 254 * in2[0,2,...] < 2**30, in2[1,3,...] < 2**29. 255 * On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29. */ 256 static void felem_diff(felem out, const felem in, const felem in2) { 257 limb carry = 0; 258 unsigned i; 259 260 for (i = 0;; i++) { 261 out[i] = in[i] - in2[i]; 262 out[i] += zero31[i]; 263 out[i] += carry; 264 carry = out[i] >> 29; 265 out[i] &= kBottom29Bits; 266 267 i++; 268 if (i == NLIMBS) 269 break; 270 271 out[i] = in[i] - in2[i]; 272 out[i] += zero31[i]; 273 out[i] += carry; 274 carry = out[i] >> 28; 275 out[i] &= kBottom28Bits; 276 } 277 278 felem_reduce_carry(out, carry); 279 } 280 281 /* felem_reduce_degree sets out = tmp/R mod p where tmp contains 64-bit words 282 * with the same 29,28,... bit positions as an felem. 283 * 284 * The values in felems are in Montgomery form: x*R mod p where R = 2**257. 285 * Since we just multiplied two Montgomery values together, the result is 286 * x*y*R*R mod p. We wish to divide by R in order for the result also to be 287 * in Montgomery form. 288 * 289 * On entry: tmp[i] < 2**64 290 * On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29 */ 291 static void felem_reduce_degree(felem out, u64 tmp[17]) { 292 /* The following table may be helpful when reading this code: 293 * 294 * Limb number: 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10... 295 * Width (bits): 29| 28| 29| 28| 29| 28| 29| 28| 29| 28| 29 296 * Start bit: 0 | 29| 57| 86|114|143|171|200|228|257|285 297 * (odd phase): 0 | 28| 57| 85|114|142|171|199|228|256|285 */ 298 limb tmp2[18], carry, x, xMask; 299 unsigned i; 300 301 /* tmp contains 64-bit words with the same 29,28,29-bit positions as an 302 * felem. So the top of an element of tmp might overlap with another 303 * element two positions down. The following loop eliminates this 304 * overlap. */ 305 tmp2[0] = (limb)(tmp[0] & kBottom29Bits); 306 307 /* In the following we use "(limb) tmp[x]" and "(limb) (tmp[x]>>32)" to try 308 * and hint to the compiler that it can do a single-word shift by selecting 309 * the right register rather than doing a double-word shift and truncating 310 * afterwards. */ 311 tmp2[1] = ((limb) tmp[0]) >> 29; 312 tmp2[1] |= (((limb)(tmp[0] >> 32)) << 3) & kBottom28Bits; 313 tmp2[1] += ((limb) tmp[1]) & kBottom28Bits; 314 carry = tmp2[1] >> 28; 315 tmp2[1] &= kBottom28Bits; 316 317 for (i = 2; i < 17; i++) { 318 tmp2[i] = ((limb)(tmp[i - 2] >> 32)) >> 25; 319 tmp2[i] += ((limb)(tmp[i - 1])) >> 28; 320 tmp2[i] += (((limb)(tmp[i - 1] >> 32)) << 4) & kBottom29Bits; 321 tmp2[i] += ((limb) tmp[i]) & kBottom29Bits; 322 tmp2[i] += carry; 323 carry = tmp2[i] >> 29; 324 tmp2[i] &= kBottom29Bits; 325 326 i++; 327 if (i == 17) 328 break; 329 tmp2[i] = ((limb)(tmp[i - 2] >> 32)) >> 25; 330 tmp2[i] += ((limb)(tmp[i - 1])) >> 29; 331 tmp2[i] += (((limb)(tmp[i - 1] >> 32)) << 3) & kBottom28Bits; 332 tmp2[i] += ((limb) tmp[i]) & kBottom28Bits; 333 tmp2[i] += carry; 334 carry = tmp2[i] >> 28; 335 tmp2[i] &= kBottom28Bits; 336 } 337 338 tmp2[17] = ((limb)(tmp[15] >> 32)) >> 25; 339 tmp2[17] += ((limb)(tmp[16])) >> 29; 340 tmp2[17] += (((limb)(tmp[16] >> 32)) << 3); 341 tmp2[17] += carry; 342 343 /* Montgomery elimination of terms. 344 * 345 * Since R is 2**257, we can divide by R with a bitwise shift if we can 346 * ensure that the right-most 257 bits are all zero. We can make that true by 347 * adding multiplies of p without affecting the value. 348 * 349 * So we eliminate limbs from right to left. Since the bottom 29 bits of p 350 * are all ones, then by adding tmp2[0]*p to tmp2 we'll make tmp2[0] == 0. 351 * We can do that for 8 further limbs and then right shift to eliminate the 352 * extra factor of R. */ 353 for (i = 0;; i += 2) { 354 tmp2[i + 1] += tmp2[i] >> 29; 355 x = tmp2[i] & kBottom29Bits; 356 xMask = NON_ZERO_TO_ALL_ONES(x); 357 tmp2[i] = 0; 358 359 /* The bounds calculations for this loop are tricky. Each iteration of 360 * the loop eliminates two words by adding values to words to their 361 * right. 362 * 363 * The following table contains the amounts added to each word (as an 364 * offset from the value of i at the top of the loop). The amounts are 365 * accounted for from the first and second half of the loop separately 366 * and are written as, for example, 28 to mean a value <2**28. 367 * 368 * Word: 3 4 5 6 7 8 9 10 369 * Added in top half: 28 11 29 21 29 28 370 * 28 29 371 * 29 372 * Added in bottom half: 29 10 28 21 28 28 373 * 29 374 * 375 * The value that is currently offset 7 will be offset 5 for the next 376 * iteration and then offset 3 for the iteration after that. Therefore 377 * the total value added will be the values added at 7, 5 and 3. 378 * 379 * The following table accumulates these values. The sums at the bottom 380 * are written as, for example, 29+28, to mean a value < 2**29+2**28. 381 * 382 * Word: 3 4 5 6 7 8 9 10 11 12 13 383 * 28 11 10 29 21 29 28 28 28 28 28 384 * 29 28 11 28 29 28 29 28 29 28 385 * 29 28 21 21 29 21 29 21 386 * 10 29 28 21 28 21 28 387 * 28 29 28 29 28 29 28 388 * 11 10 29 10 29 10 389 * 29 28 11 28 11 390 * 29 29 391 * -------------------------------------------- 392 * 30+ 31+ 30+ 31+ 30+ 393 * 28+ 29+ 28+ 29+ 21+ 394 * 21+ 28+ 21+ 28+ 10 395 * 10 21+ 10 21+ 396 * 11 11 397 * 398 * So the greatest amount is added to tmp2[10] and tmp2[12]. If 399 * tmp2[10/12] has an initial value of <2**29, then the maximum value 400 * will be < 2**31 + 2**30 + 2**28 + 2**21 + 2**11, which is < 2**32, 401 * as required. */ 402 tmp2[i + 3] += (x << 10) & kBottom28Bits; 403 tmp2[i + 4] += (x >> 18); 404 405 tmp2[i + 6] += (x << 21) & kBottom29Bits; 406 tmp2[i + 7] += x >> 8; 407 408 /* At position 200, which is the starting bit position for word 7, we 409 * have a factor of 0xf000000 = 2**28 - 2**24. */ 410 tmp2[i + 7] += 0x10000000 & xMask; 411 /* Word 7 is 28 bits wide, so the 2**28 term exactly hits word 8. */ 412 tmp2[i + 8] += (x - 1) & xMask; 413 tmp2[i + 7] -= (x << 24) & kBottom28Bits; 414 tmp2[i + 8] -= x >> 4; 415 416 tmp2[i + 8] += 0x20000000 & xMask; 417 tmp2[i + 8] -= x; 418 tmp2[i + 8] += (x << 28) & kBottom29Bits; 419 tmp2[i + 9] += ((x >> 1) - 1) & xMask; 420 421 if (i+1 == NLIMBS) 422 break; 423 tmp2[i + 2] += tmp2[i + 1] >> 28; 424 x = tmp2[i + 1] & kBottom28Bits; 425 xMask = NON_ZERO_TO_ALL_ONES(x); 426 tmp2[i + 1] = 0; 427 428 tmp2[i + 4] += (x << 11) & kBottom29Bits; 429 tmp2[i + 5] += (x >> 18); 430 431 tmp2[i + 7] += (x << 21) & kBottom28Bits; 432 tmp2[i + 8] += x >> 7; 433 434 /* At position 199, which is the starting bit of the 8th word when 435 * dealing with a context starting on an odd word, we have a factor of 436 * 0x1e000000 = 2**29 - 2**25. Since we have not updated i, the 8th 437 * word from i+1 is i+8. */ 438 tmp2[i + 8] += 0x20000000 & xMask; 439 tmp2[i + 9] += (x - 1) & xMask; 440 tmp2[i + 8] -= (x << 25) & kBottom29Bits; 441 tmp2[i + 9] -= x >> 4; 442 443 tmp2[i + 9] += 0x10000000 & xMask; 444 tmp2[i + 9] -= x; 445 tmp2[i + 10] += (x - 1) & xMask; 446 } 447 448 /* We merge the right shift with a carry chain. The words above 2**257 have 449 * widths of 28,29,... which we need to correct when copying them down. */ 450 carry = 0; 451 for (i = 0; i < 8; i++) { 452 /* The maximum value of tmp2[i + 9] occurs on the first iteration and 453 * is < 2**30+2**29+2**28. Adding 2**29 (from tmp2[i + 10]) is 454 * therefore safe. */ 455 out[i] = tmp2[i + 9]; 456 out[i] += carry; 457 out[i] += (tmp2[i + 10] << 28) & kBottom29Bits; 458 carry = out[i] >> 29; 459 out[i] &= kBottom29Bits; 460 461 i++; 462 out[i] = tmp2[i + 9] >> 1; 463 out[i] += carry; 464 carry = out[i] >> 28; 465 out[i] &= kBottom28Bits; 466 } 467 468 out[8] = tmp2[17]; 469 out[8] += carry; 470 carry = out[8] >> 29; 471 out[8] &= kBottom29Bits; 472 473 felem_reduce_carry(out, carry); 474 } 475 476 /* felem_square sets out=in*in. 477 * 478 * On entry: in[0,2,...] < 2**30, in[1,3,...] < 2**29. 479 * On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29. */ 480 static void felem_square(felem out, const felem in) { 481 u64 tmp[17]; 482 483 tmp[0] = ((u64) in[0]) * in[0]; 484 tmp[1] = ((u64) in[0]) * (in[1] << 1); 485 tmp[2] = ((u64) in[0]) * (in[2] << 1) + 486 ((u64) in[1]) * (in[1] << 1); 487 tmp[3] = ((u64) in[0]) * (in[3] << 1) + 488 ((u64) in[1]) * (in[2] << 1); 489 tmp[4] = ((u64) in[0]) * (in[4] << 1) + 490 ((u64) in[1]) * (in[3] << 2) + ((u64) in[2]) * in[2]; 491 tmp[5] = ((u64) in[0]) * (in[5] << 1) + ((u64) in[1]) * 492 (in[4] << 1) + ((u64) in[2]) * (in[3] << 1); 493 tmp[6] = ((u64) in[0]) * (in[6] << 1) + ((u64) in[1]) * 494 (in[5] << 2) + ((u64) in[2]) * (in[4] << 1) + 495 ((u64) in[3]) * (in[3] << 1); 496 tmp[7] = ((u64) in[0]) * (in[7] << 1) + ((u64) in[1]) * 497 (in[6] << 1) + ((u64) in[2]) * (in[5] << 1) + 498 ((u64) in[3]) * (in[4] << 1); 499 /* tmp[8] has the greatest value of 2**61 + 2**60 + 2**61 + 2**60 + 2**60, 500 * which is < 2**64 as required. */ 501 tmp[8] = ((u64) in[0]) * (in[8] << 1) + ((u64) in[1]) * 502 (in[7] << 2) + ((u64) in[2]) * (in[6] << 1) + 503 ((u64) in[3]) * (in[5] << 2) + ((u64) in[4]) * in[4]; 504 tmp[9] = ((u64) in[1]) * (in[8] << 1) + ((u64) in[2]) * 505 (in[7] << 1) + ((u64) in[3]) * (in[6] << 1) + 506 ((u64) in[4]) * (in[5] << 1); 507 tmp[10] = ((u64) in[2]) * (in[8] << 1) + ((u64) in[3]) * 508 (in[7] << 2) + ((u64) in[4]) * (in[6] << 1) + 509 ((u64) in[5]) * (in[5] << 1); 510 tmp[11] = ((u64) in[3]) * (in[8] << 1) + ((u64) in[4]) * 511 (in[7] << 1) + ((u64) in[5]) * (in[6] << 1); 512 tmp[12] = ((u64) in[4]) * (in[8] << 1) + 513 ((u64) in[5]) * (in[7] << 2) + ((u64) in[6]) * in[6]; 514 tmp[13] = ((u64) in[5]) * (in[8] << 1) + 515 ((u64) in[6]) * (in[7] << 1); 516 tmp[14] = ((u64) in[6]) * (in[8] << 1) + 517 ((u64) in[7]) * (in[7] << 1); 518 tmp[15] = ((u64) in[7]) * (in[8] << 1); 519 tmp[16] = ((u64) in[8]) * in[8]; 520 521 felem_reduce_degree(out, tmp); 522 } 523 524 /* felem_mul sets out=in*in2. 525 * 526 * On entry: in[0,2,...] < 2**30, in[1,3,...] < 2**29 and 527 * in2[0,2,...] < 2**30, in2[1,3,...] < 2**29. 528 * On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29. */ 529 static void felem_mul(felem out, const felem in, const felem in2) { 530 u64 tmp[17]; 531 532 tmp[0] = ((u64) in[0]) * in2[0]; 533 tmp[1] = ((u64) in[0]) * (in2[1] << 0) + 534 ((u64) in[1]) * (in2[0] << 0); 535 tmp[2] = ((u64) in[0]) * (in2[2] << 0) + ((u64) in[1]) * 536 (in2[1] << 1) + ((u64) in[2]) * (in2[0] << 0); 537 tmp[3] = ((u64) in[0]) * (in2[3] << 0) + ((u64) in[1]) * 538 (in2[2] << 0) + ((u64) in[2]) * (in2[1] << 0) + 539 ((u64) in[3]) * (in2[0] << 0); 540 tmp[4] = ((u64) in[0]) * (in2[4] << 0) + ((u64) in[1]) * 541 (in2[3] << 1) + ((u64) in[2]) * (in2[2] << 0) + 542 ((u64) in[3]) * (in2[1] << 1) + 543 ((u64) in[4]) * (in2[0] << 0); 544 tmp[5] = ((u64) in[0]) * (in2[5] << 0) + ((u64) in[1]) * 545 (in2[4] << 0) + ((u64) in[2]) * (in2[3] << 0) + 546 ((u64) in[3]) * (in2[2] << 0) + ((u64) in[4]) * 547 (in2[1] << 0) + ((u64) in[5]) * (in2[0] << 0); 548 tmp[6] = ((u64) in[0]) * (in2[6] << 0) + ((u64) in[1]) * 549 (in2[5] << 1) + ((u64) in[2]) * (in2[4] << 0) + 550 ((u64) in[3]) * (in2[3] << 1) + ((u64) in[4]) * 551 (in2[2] << 0) + ((u64) in[5]) * (in2[1] << 1) + 552 ((u64) in[6]) * (in2[0] << 0); 553 tmp[7] = ((u64) in[0]) * (in2[7] << 0) + ((u64) in[1]) * 554 (in2[6] << 0) + ((u64) in[2]) * (in2[5] << 0) + 555 ((u64) in[3]) * (in2[4] << 0) + ((u64) in[4]) * 556 (in2[3] << 0) + ((u64) in[5]) * (in2[2] << 0) + 557 ((u64) in[6]) * (in2[1] << 0) + 558 ((u64) in[7]) * (in2[0] << 0); 559 /* tmp[8] has the greatest value but doesn't overflow. See logic in 560 * felem_square. */ 561 tmp[8] = ((u64) in[0]) * (in2[8] << 0) + ((u64) in[1]) * 562 (in2[7] << 1) + ((u64) in[2]) * (in2[6] << 0) + 563 ((u64) in[3]) * (in2[5] << 1) + ((u64) in[4]) * 564 (in2[4] << 0) + ((u64) in[5]) * (in2[3] << 1) + 565 ((u64) in[6]) * (in2[2] << 0) + ((u64) in[7]) * 566 (in2[1] << 1) + ((u64) in[8]) * (in2[0] << 0); 567 tmp[9] = ((u64) in[1]) * (in2[8] << 0) + ((u64) in[2]) * 568 (in2[7] << 0) + ((u64) in[3]) * (in2[6] << 0) + 569 ((u64) in[4]) * (in2[5] << 0) + ((u64) in[5]) * 570 (in2[4] << 0) + ((u64) in[6]) * (in2[3] << 0) + 571 ((u64) in[7]) * (in2[2] << 0) + 572 ((u64) in[8]) * (in2[1] << 0); 573 tmp[10] = ((u64) in[2]) * (in2[8] << 0) + ((u64) in[3]) * 574 (in2[7] << 1) + ((u64) in[4]) * (in2[6] << 0) + 575 ((u64) in[5]) * (in2[5] << 1) + ((u64) in[6]) * 576 (in2[4] << 0) + ((u64) in[7]) * (in2[3] << 1) + 577 ((u64) in[8]) * (in2[2] << 0); 578 tmp[11] = ((u64) in[3]) * (in2[8] << 0) + ((u64) in[4]) * 579 (in2[7] << 0) + ((u64) in[5]) * (in2[6] << 0) + 580 ((u64) in[6]) * (in2[5] << 0) + ((u64) in[7]) * 581 (in2[4] << 0) + ((u64) in[8]) * (in2[3] << 0); 582 tmp[12] = ((u64) in[4]) * (in2[8] << 0) + ((u64) in[5]) * 583 (in2[7] << 1) + ((u64) in[6]) * (in2[6] << 0) + 584 ((u64) in[7]) * (in2[5] << 1) + 585 ((u64) in[8]) * (in2[4] << 0); 586 tmp[13] = ((u64) in[5]) * (in2[8] << 0) + ((u64) in[6]) * 587 (in2[7] << 0) + ((u64) in[7]) * (in2[6] << 0) + 588 ((u64) in[8]) * (in2[5] << 0); 589 tmp[14] = ((u64) in[6]) * (in2[8] << 0) + ((u64) in[7]) * 590 (in2[7] << 1) + ((u64) in[8]) * (in2[6] << 0); 591 tmp[15] = ((u64) in[7]) * (in2[8] << 0) + 592 ((u64) in[8]) * (in2[7] << 0); 593 tmp[16] = ((u64) in[8]) * (in2[8] << 0); 594 595 felem_reduce_degree(out, tmp); 596 } 597 598 static void felem_assign(felem out, const felem in) { 599 memcpy(out, in, sizeof(felem)); 600 } 601 602 /* felem_inv calculates |out| = |in|^{-1} 603 * 604 * Based on Fermat's Little Theorem: 605 * a^p = a (mod p) 606 * a^{p-1} = 1 (mod p) 607 * a^{p-2} = a^{-1} (mod p) 608 */ 609 static void felem_inv(felem out, const felem in) { 610 felem ftmp, ftmp2; 611 /* each e_I will hold |in|^{2^I - 1} */ 612 felem e2, e4, e8, e16, e32, e64; 613 unsigned i; 614 615 felem_square(ftmp, in); /* 2^1 */ 616 felem_mul(ftmp, in, ftmp); /* 2^2 - 2^0 */ 617 felem_assign(e2, ftmp); 618 felem_square(ftmp, ftmp); /* 2^3 - 2^1 */ 619 felem_square(ftmp, ftmp); /* 2^4 - 2^2 */ 620 felem_mul(ftmp, ftmp, e2); /* 2^4 - 2^0 */ 621 felem_assign(e4, ftmp); 622 felem_square(ftmp, ftmp); /* 2^5 - 2^1 */ 623 felem_square(ftmp, ftmp); /* 2^6 - 2^2 */ 624 felem_square(ftmp, ftmp); /* 2^7 - 2^3 */ 625 felem_square(ftmp, ftmp); /* 2^8 - 2^4 */ 626 felem_mul(ftmp, ftmp, e4); /* 2^8 - 2^0 */ 627 felem_assign(e8, ftmp); 628 for (i = 0; i < 8; i++) { 629 felem_square(ftmp, ftmp); 630 } /* 2^16 - 2^8 */ 631 felem_mul(ftmp, ftmp, e8); /* 2^16 - 2^0 */ 632 felem_assign(e16, ftmp); 633 for (i = 0; i < 16; i++) { 634 felem_square(ftmp, ftmp); 635 } /* 2^32 - 2^16 */ 636 felem_mul(ftmp, ftmp, e16); /* 2^32 - 2^0 */ 637 felem_assign(e32, ftmp); 638 for (i = 0; i < 32; i++) { 639 felem_square(ftmp, ftmp); 640 } /* 2^64 - 2^32 */ 641 felem_assign(e64, ftmp); 642 felem_mul(ftmp, ftmp, in); /* 2^64 - 2^32 + 2^0 */ 643 for (i = 0; i < 192; i++) { 644 felem_square(ftmp, ftmp); 645 } /* 2^256 - 2^224 + 2^192 */ 646 647 felem_mul(ftmp2, e64, e32); /* 2^64 - 2^0 */ 648 for (i = 0; i < 16; i++) { 649 felem_square(ftmp2, ftmp2); 650 } /* 2^80 - 2^16 */ 651 felem_mul(ftmp2, ftmp2, e16); /* 2^80 - 2^0 */ 652 for (i = 0; i < 8; i++) { 653 felem_square(ftmp2, ftmp2); 654 } /* 2^88 - 2^8 */ 655 felem_mul(ftmp2, ftmp2, e8); /* 2^88 - 2^0 */ 656 for (i = 0; i < 4; i++) { 657 felem_square(ftmp2, ftmp2); 658 } /* 2^92 - 2^4 */ 659 felem_mul(ftmp2, ftmp2, e4); /* 2^92 - 2^0 */ 660 felem_square(ftmp2, ftmp2); /* 2^93 - 2^1 */ 661 felem_square(ftmp2, ftmp2); /* 2^94 - 2^2 */ 662 felem_mul(ftmp2, ftmp2, e2); /* 2^94 - 2^0 */ 663 felem_square(ftmp2, ftmp2); /* 2^95 - 2^1 */ 664 felem_square(ftmp2, ftmp2); /* 2^96 - 2^2 */ 665 felem_mul(ftmp2, ftmp2, in); /* 2^96 - 3 */ 666 667 felem_mul(out, ftmp2, ftmp); /* 2^256 - 2^224 + 2^192 + 2^96 - 3 */ 668 } 669 670 /* felem_scalar_3 sets out=3*out. 671 * 672 * On entry: out[0,2,...] < 2**30, out[1,3,...] < 2**29. 673 * On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29. */ 674 static void felem_scalar_3(felem out) { 675 limb carry = 0; 676 unsigned i; 677 678 for (i = 0;; i++) { 679 out[i] *= 3; 680 out[i] += carry; 681 carry = out[i] >> 29; 682 out[i] &= kBottom29Bits; 683 684 i++; 685 if (i == NLIMBS) 686 break; 687 688 out[i] *= 3; 689 out[i] += carry; 690 carry = out[i] >> 28; 691 out[i] &= kBottom28Bits; 692 } 693 694 felem_reduce_carry(out, carry); 695 } 696 697 /* felem_scalar_4 sets out=4*out. 698 * 699 * On entry: out[0,2,...] < 2**30, out[1,3,...] < 2**29. 700 * On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29. */ 701 static void felem_scalar_4(felem out) { 702 limb carry = 0, next_carry; 703 unsigned i; 704 705 for (i = 0;; i++) { 706 next_carry = out[i] >> 27; 707 out[i] <<= 2; 708 out[i] &= kBottom29Bits; 709 out[i] += carry; 710 carry = next_carry + (out[i] >> 29); 711 out[i] &= kBottom29Bits; 712 713 i++; 714 if (i == NLIMBS) 715 break; 716 717 next_carry = out[i] >> 26; 718 out[i] <<= 2; 719 out[i] &= kBottom28Bits; 720 out[i] += carry; 721 carry = next_carry + (out[i] >> 28); 722 out[i] &= kBottom28Bits; 723 } 724 725 felem_reduce_carry(out, carry); 726 } 727 728 /* felem_scalar_8 sets out=8*out. 729 * 730 * On entry: out[0,2,...] < 2**30, out[1,3,...] < 2**29. 731 * On exit: out[0,2,...] < 2**30, out[1,3,...] < 2**29. */ 732 static void felem_scalar_8(felem out) { 733 limb carry = 0, next_carry; 734 unsigned i; 735 736 for (i = 0;; i++) { 737 next_carry = out[i] >> 26; 738 out[i] <<= 3; 739 out[i] &= kBottom29Bits; 740 out[i] += carry; 741 carry = next_carry + (out[i] >> 29); 742 out[i] &= kBottom29Bits; 743 744 i++; 745 if (i == NLIMBS) 746 break; 747 748 next_carry = out[i] >> 25; 749 out[i] <<= 3; 750 out[i] &= kBottom28Bits; 751 out[i] += carry; 752 carry = next_carry + (out[i] >> 28); 753 out[i] &= kBottom28Bits; 754 } 755 756 felem_reduce_carry(out, carry); 757 } 758 759 /* felem_is_zero_vartime returns 1 iff |in| == 0. It takes a variable amount of 760 * time depending on the value of |in|. */ 761 static char felem_is_zero_vartime(const felem in) { 762 limb carry; 763 int i; 764 limb tmp[NLIMBS]; 765 766 felem_assign(tmp, in); 767 768 /* First, reduce tmp to a minimal form. */ 769 do { 770 carry = 0; 771 for (i = 0;; i++) { 772 tmp[i] += carry; 773 carry = tmp[i] >> 29; 774 tmp[i] &= kBottom29Bits; 775 776 i++; 777 if (i == NLIMBS) 778 break; 779 780 tmp[i] += carry; 781 carry = tmp[i] >> 28; 782 tmp[i] &= kBottom28Bits; 783 } 784 785 felem_reduce_carry(tmp, carry); 786 } while (carry); 787 788 /* tmp < 2**257, so the only possible zero values are 0, p and 2p. */ 789 return memcmp(tmp, kZero, sizeof(tmp)) == 0 || 790 memcmp(tmp, kP, sizeof(tmp)) == 0 || 791 memcmp(tmp, k2P, sizeof(tmp)) == 0; 792 } 793 794 795 /* Group operations: 796 * 797 * Elements of the elliptic curve group are represented in Jacobian 798 * coordinates: (x, y, z). An affine point (x', y') is x'=x/z**2, y'=y/z**3 in 799 * Jacobian form. */ 800 801 /* point_double sets {x_out,y_out,z_out} = 2*{x,y,z}. 802 * 803 * See http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-dbl-2009-l */ 804 static void point_double(felem x_out, felem y_out, felem z_out, const felem x, 805 const felem y, const felem z) { 806 felem delta, gamma, alpha, beta, tmp, tmp2; 807 808 felem_square(delta, z); 809 felem_square(gamma, y); 810 felem_mul(beta, x, gamma); 811 812 felem_sum(tmp, x, delta); 813 felem_diff(tmp2, x, delta); 814 felem_mul(alpha, tmp, tmp2); 815 felem_scalar_3(alpha); 816 817 felem_sum(tmp, y, z); 818 felem_square(tmp, tmp); 819 felem_diff(tmp, tmp, gamma); 820 felem_diff(z_out, tmp, delta); 821 822 felem_scalar_4(beta); 823 felem_square(x_out, alpha); 824 felem_diff(x_out, x_out, beta); 825 felem_diff(x_out, x_out, beta); 826 827 felem_diff(tmp, beta, x_out); 828 felem_mul(tmp, alpha, tmp); 829 felem_square(tmp2, gamma); 830 felem_scalar_8(tmp2); 831 felem_diff(y_out, tmp, tmp2); 832 } 833 834 /* point_add_mixed sets {x_out,y_out,z_out} = {x1,y1,z1} + {x2,y2,1}. 835 * (i.e. the second point is affine.) 836 * 837 * See http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-add-2007-bl 838 * 839 * Note that this function does not handle P+P, infinity+P nor P+infinity 840 * correctly. */ 841 static void point_add_mixed(felem x_out, felem y_out, felem z_out, 842 const felem x1, const felem y1, const felem z1, 843 const felem x2, const felem y2) { 844 felem z1z1, z1z1z1, s2, u2, h, i, j, r, rr, v, tmp; 845 846 felem_square(z1z1, z1); 847 felem_sum(tmp, z1, z1); 848 849 felem_mul(u2, x2, z1z1); 850 felem_mul(z1z1z1, z1, z1z1); 851 felem_mul(s2, y2, z1z1z1); 852 felem_diff(h, u2, x1); 853 felem_sum(i, h, h); 854 felem_square(i, i); 855 felem_mul(j, h, i); 856 felem_diff(r, s2, y1); 857 felem_sum(r, r, r); 858 felem_mul(v, x1, i); 859 860 felem_mul(z_out, tmp, h); 861 felem_square(rr, r); 862 felem_diff(x_out, rr, j); 863 felem_diff(x_out, x_out, v); 864 felem_diff(x_out, x_out, v); 865 866 felem_diff(tmp, v, x_out); 867 felem_mul(y_out, tmp, r); 868 felem_mul(tmp, y1, j); 869 felem_diff(y_out, y_out, tmp); 870 felem_diff(y_out, y_out, tmp); 871 } 872 873 /* point_add sets {x_out,y_out,z_out} = {x1,y1,z1} + {x2,y2,z2}. 874 * 875 * See http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-add-2007-bl 876 * 877 * Note that this function does not handle P+P, infinity+P nor P+infinity 878 * correctly. */ 879 static void point_add(felem x_out, felem y_out, felem z_out, const felem x1, 880 const felem y1, const felem z1, const felem x2, 881 const felem y2, const felem z2) { 882 felem z1z1, z1z1z1, z2z2, z2z2z2, s1, s2, u1, u2, h, i, j, r, rr, v, tmp; 883 884 felem_square(z1z1, z1); 885 felem_square(z2z2, z2); 886 felem_mul(u1, x1, z2z2); 887 888 felem_sum(tmp, z1, z2); 889 felem_square(tmp, tmp); 890 felem_diff(tmp, tmp, z1z1); 891 felem_diff(tmp, tmp, z2z2); 892 893 felem_mul(z2z2z2, z2, z2z2); 894 felem_mul(s1, y1, z2z2z2); 895 896 felem_mul(u2, x2, z1z1); 897 felem_mul(z1z1z1, z1, z1z1); 898 felem_mul(s2, y2, z1z1z1); 899 felem_diff(h, u2, u1); 900 felem_sum(i, h, h); 901 felem_square(i, i); 902 felem_mul(j, h, i); 903 felem_diff(r, s2, s1); 904 felem_sum(r, r, r); 905 felem_mul(v, u1, i); 906 907 felem_mul(z_out, tmp, h); 908 felem_square(rr, r); 909 felem_diff(x_out, rr, j); 910 felem_diff(x_out, x_out, v); 911 felem_diff(x_out, x_out, v); 912 913 felem_diff(tmp, v, x_out); 914 felem_mul(y_out, tmp, r); 915 felem_mul(tmp, s1, j); 916 felem_diff(y_out, y_out, tmp); 917 felem_diff(y_out, y_out, tmp); 918 } 919 920 /* point_add_or_double_vartime sets {x_out,y_out,z_out} = {x1,y1,z1} + 921 * {x2,y2,z2}. 922 * 923 * See http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-add-2007-bl 924 * 925 * This function handles the case where {x1,y1,z1}={x2,y2,z2}. */ 926 static void point_add_or_double_vartime( 927 felem x_out, felem y_out, felem z_out, const felem x1, const felem y1, 928 const felem z1, const felem x2, const felem y2, const felem z2) { 929 felem z1z1, z1z1z1, z2z2, z2z2z2, s1, s2, u1, u2, h, i, j, r, rr, v, tmp; 930 char x_equal, y_equal; 931 932 felem_square(z1z1, z1); 933 felem_square(z2z2, z2); 934 felem_mul(u1, x1, z2z2); 935 936 felem_sum(tmp, z1, z2); 937 felem_square(tmp, tmp); 938 felem_diff(tmp, tmp, z1z1); 939 felem_diff(tmp, tmp, z2z2); 940 941 felem_mul(z2z2z2, z2, z2z2); 942 felem_mul(s1, y1, z2z2z2); 943 944 felem_mul(u2, x2, z1z1); 945 felem_mul(z1z1z1, z1, z1z1); 946 felem_mul(s2, y2, z1z1z1); 947 felem_diff(h, u2, u1); 948 x_equal = felem_is_zero_vartime(h); 949 felem_sum(i, h, h); 950 felem_square(i, i); 951 felem_mul(j, h, i); 952 felem_diff(r, s2, s1); 953 y_equal = felem_is_zero_vartime(r); 954 if (x_equal && y_equal) { 955 point_double(x_out, y_out, z_out, x1, y1, z1); 956 return; 957 } 958 felem_sum(r, r, r); 959 felem_mul(v, u1, i); 960 961 felem_mul(z_out, tmp, h); 962 felem_square(rr, r); 963 felem_diff(x_out, rr, j); 964 felem_diff(x_out, x_out, v); 965 felem_diff(x_out, x_out, v); 966 967 felem_diff(tmp, v, x_out); 968 felem_mul(y_out, tmp, r); 969 felem_mul(tmp, s1, j); 970 felem_diff(y_out, y_out, tmp); 971 felem_diff(y_out, y_out, tmp); 972 } 973 974 /* copy_conditional sets out=in if mask = 0xffffffff in constant time. 975 * 976 * On entry: mask is either 0 or 0xffffffff. */ 977 static void copy_conditional(felem out, const felem in, limb mask) { 978 int i; 979 980 for (i = 0; i < NLIMBS; i++) { 981 const limb tmp = mask & (in[i] ^ out[i]); 982 out[i] ^= tmp; 983 } 984 } 985 986 /* select_affine_point sets {out_x,out_y} to the index'th entry of table. 987 * On entry: index < 16, table[0] must be zero. */ 988 static void select_affine_point(felem out_x, felem out_y, const limb* table, 989 limb index) { 990 limb i, j; 991 992 memset(out_x, 0, sizeof(felem)); 993 memset(out_y, 0, sizeof(felem)); 994 995 for (i = 1; i < 16; i++) { 996 limb mask = i ^ index; 997 mask |= mask >> 2; 998 mask |= mask >> 1; 999 mask &= 1; 1000 mask--; 1001 for (j = 0; j < NLIMBS; j++, table++) { 1002 out_x[j] |= *table & mask; 1003 } 1004 for (j = 0; j < NLIMBS; j++, table++) { 1005 out_y[j] |= *table & mask; 1006 } 1007 } 1008 } 1009 1010 /* select_jacobian_point sets {out_x,out_y,out_z} to the index'th entry of 1011 * table. On entry: index < 16, table[0] must be zero. */ 1012 static void select_jacobian_point(felem out_x, felem out_y, felem out_z, 1013 const limb* table, limb index) { 1014 limb i, j; 1015 1016 memset(out_x, 0, sizeof(felem)); 1017 memset(out_y, 0, sizeof(felem)); 1018 memset(out_z, 0, sizeof(felem)); 1019 1020 /* The implicit value at index 0 is all zero. We don't need to perform that 1021 * iteration of the loop because we already set out_* to zero. */ 1022 table += 3 * NLIMBS; 1023 1024 // Hit all entries to obscure cache profiling. 1025 for (i = 1; i < 16; i++) { 1026 limb mask = i ^ index; 1027 mask |= mask >> 2; 1028 mask |= mask >> 1; 1029 mask &= 1; 1030 mask--; 1031 for (j = 0; j < NLIMBS; j++, table++) { 1032 out_x[j] |= *table & mask; 1033 } 1034 for (j = 0; j < NLIMBS; j++, table++) { 1035 out_y[j] |= *table & mask; 1036 } 1037 for (j = 0; j < NLIMBS; j++, table++) { 1038 out_z[j] |= *table & mask; 1039 } 1040 } 1041 } 1042 1043 /* scalar_base_mult sets {nx,ny,nz} = scalar*G where scalar is a little-endian 1044 * number. Note that the value of scalar must be less than the order of the 1045 * group. */ 1046 static void scalar_base_mult(felem nx, felem ny, felem nz, 1047 const p256_int* scalar) { 1048 int i, j; 1049 limb n_is_infinity_mask = -1, p_is_noninfinite_mask, mask; 1050 u32 table_offset; 1051 1052 felem px, py; 1053 felem tx, ty, tz; 1054 1055 memset(nx, 0, sizeof(felem)); 1056 memset(ny, 0, sizeof(felem)); 1057 memset(nz, 0, sizeof(felem)); 1058 1059 /* The loop adds bits at positions 0, 64, 128 and 192, followed by 1060 * positions 32,96,160 and 224 and does this 32 times. */ 1061 for (i = 0; i < 32; i++) { 1062 if (i) { 1063 point_double(nx, ny, nz, nx, ny, nz); 1064 } 1065 table_offset = 0; 1066 for (j = 0; j <= 32; j += 32) { 1067 char bit0 = p256_get_bit(scalar, 31 - i + j); 1068 char bit1 = p256_get_bit(scalar, 95 - i + j); 1069 char bit2 = p256_get_bit(scalar, 159 - i + j); 1070 char bit3 = p256_get_bit(scalar, 223 - i + j); 1071 limb index = bit0 | (bit1 << 1) | (bit2 << 2) | (bit3 << 3); 1072 1073 select_affine_point(px, py, kPrecomputed + table_offset, index); 1074 table_offset += 30 * NLIMBS; 1075 1076 /* Since scalar is less than the order of the group, we know that 1077 * {nx,ny,nz} != {px,py,1}, unless both are zero, which we handle 1078 * below. */ 1079 point_add_mixed(tx, ty, tz, nx, ny, nz, px, py); 1080 /* The result of point_add_mixed is incorrect if {nx,ny,nz} is zero 1081 * (a.k.a. the point at infinity). We handle that situation by 1082 * copying the point from the table. */ 1083 copy_conditional(nx, px, n_is_infinity_mask); 1084 copy_conditional(ny, py, n_is_infinity_mask); 1085 copy_conditional(nz, kOne, n_is_infinity_mask); 1086 1087 /* Equally, the result is also wrong if the point from the table is 1088 * zero, which happens when the index is zero. We handle that by 1089 * only copying from {tx,ty,tz} to {nx,ny,nz} if index != 0. */ 1090 p_is_noninfinite_mask = NON_ZERO_TO_ALL_ONES(index); 1091 mask = p_is_noninfinite_mask & ~n_is_infinity_mask; 1092 copy_conditional(nx, tx, mask); 1093 copy_conditional(ny, ty, mask); 1094 copy_conditional(nz, tz, mask); 1095 /* If p was not zero, then n is now non-zero. */ 1096 n_is_infinity_mask &= ~p_is_noninfinite_mask; 1097 } 1098 } 1099 } 1100 1101 /* point_to_affine converts a Jacobian point to an affine point. If the input 1102 * is the point at infinity then it returns (0, 0) in constant time. */ 1103 static void point_to_affine(felem x_out, felem y_out, const felem nx, 1104 const felem ny, const felem nz) { 1105 felem z_inv, z_inv_sq; 1106 felem_inv(z_inv, nz); 1107 felem_square(z_inv_sq, z_inv); 1108 felem_mul(x_out, nx, z_inv_sq); 1109 felem_mul(z_inv, z_inv, z_inv_sq); 1110 felem_mul(y_out, ny, z_inv); 1111 } 1112 1113 /* scalar_base_mult sets {nx,ny,nz} = scalar*{x,y}. */ 1114 static void scalar_mult(felem nx, felem ny, felem nz, const felem x, 1115 const felem y, const p256_int* scalar) { 1116 int i; 1117 felem px, py, pz, tx, ty, tz; 1118 felem precomp[16][3]; 1119 limb n_is_infinity_mask, index, p_is_noninfinite_mask, mask; 1120 1121 /* We precompute 0,1,2,... times {x,y}. */ 1122 memset(precomp, 0, sizeof(felem) * 3); 1123 memcpy(&precomp[1][0], x, sizeof(felem)); 1124 memcpy(&precomp[1][1], y, sizeof(felem)); 1125 memcpy(&precomp[1][2], kOne, sizeof(felem)); 1126 1127 for (i = 2; i < 16; i += 2) { 1128 point_double(precomp[i][0], precomp[i][1], precomp[i][2], 1129 precomp[i / 2][0], precomp[i / 2][1], precomp[i / 2][2]); 1130 1131 point_add_mixed(precomp[i + 1][0], precomp[i + 1][1], precomp[i + 1][2], 1132 precomp[i][0], precomp[i][1], precomp[i][2], x, y); 1133 } 1134 1135 memset(nx, 0, sizeof(felem)); 1136 memset(ny, 0, sizeof(felem)); 1137 memset(nz, 0, sizeof(felem)); 1138 n_is_infinity_mask = -1; 1139 1140 /* We add in a window of four bits each iteration and do this 64 times. */ 1141 for (i = 0; i < 256; i += 4) { 1142 if (i) { 1143 point_double(nx, ny, nz, nx, ny, nz); 1144 point_double(nx, ny, nz, nx, ny, nz); 1145 point_double(nx, ny, nz, nx, ny, nz); 1146 point_double(nx, ny, nz, nx, ny, nz); 1147 } 1148 1149 index = (p256_get_bit(scalar, 255 - i - 0) << 3) | 1150 (p256_get_bit(scalar, 255 - i - 1) << 2) | 1151 (p256_get_bit(scalar, 255 - i - 2) << 1) | 1152 p256_get_bit(scalar, 255 - i - 3); 1153 1154 /* See the comments in scalar_base_mult about handling infinities. */ 1155 select_jacobian_point(px, py, pz, precomp[0][0], index); 1156 point_add(tx, ty, tz, nx, ny, nz, px, py, pz); 1157 copy_conditional(nx, px, n_is_infinity_mask); 1158 copy_conditional(ny, py, n_is_infinity_mask); 1159 copy_conditional(nz, pz, n_is_infinity_mask); 1160 1161 p_is_noninfinite_mask = NON_ZERO_TO_ALL_ONES(index); 1162 mask = p_is_noninfinite_mask & ~n_is_infinity_mask; 1163 1164 copy_conditional(nx, tx, mask); 1165 copy_conditional(ny, ty, mask); 1166 copy_conditional(nz, tz, mask); 1167 n_is_infinity_mask &= ~p_is_noninfinite_mask; 1168 } 1169 } 1170 1171 #define kRDigits {2, 0, 0, 0xfffffffe, 0xffffffff, 0xffffffff, 0xfffffffd, 1} // 2^257 mod p256.p 1172 1173 #define kRInvDigits {0x80000000, 1, 0xffffffff, 0, 0x80000001, 0xfffffffe, 1, 0x7fffffff} // 1 / 2^257 mod p256.p 1174 1175 static const p256_int kR = { kRDigits }; 1176 static const p256_int kRInv = { kRInvDigits }; 1177 1178 /* to_montgomery sets out = R*in. */ 1179 static void to_montgomery(felem out, const p256_int* in) { 1180 p256_int in_shifted; 1181 int i; 1182 1183 p256_init(&in_shifted); 1184 p256_modmul(&SECP256r1_p, in, 0, &kR, &in_shifted); 1185 1186 for (i = 0; i < NLIMBS; i++) { 1187 if ((i & 1) == 0) { 1188 out[i] = P256_DIGIT(&in_shifted, 0) & kBottom29Bits; 1189 p256_shr(&in_shifted, 29, &in_shifted); 1190 } else { 1191 out[i] = P256_DIGIT(&in_shifted, 0) & kBottom28Bits; 1192 p256_shr(&in_shifted, 28, &in_shifted); 1193 } 1194 } 1195 1196 p256_clear(&in_shifted); 1197 } 1198 1199 /* from_montgomery sets out=in/R. */ 1200 static void from_montgomery(p256_int* out, const felem in) { 1201 p256_int result, tmp; 1202 int i, top; 1203 1204 p256_init(&result); 1205 p256_init(&tmp); 1206 1207 p256_add_d(&tmp, in[NLIMBS - 1], &result); 1208 for (i = NLIMBS - 2; i >= 0; i--) { 1209 if ((i & 1) == 0) { 1210 top = p256_shl(&result, 29, &tmp); 1211 } else { 1212 top = p256_shl(&result, 28, &tmp); 1213 } 1214 top |= p256_add_d(&tmp, in[i], &result); 1215 } 1216 1217 p256_modmul(&SECP256r1_p, &kRInv, top, &result, out); 1218 1219 p256_clear(&result); 1220 p256_clear(&tmp); 1221 } 1222 1223 /* p256_base_point_mul sets {out_x,out_y} = nG, where n is < the 1224 * order of the group. */ 1225 void p256_base_point_mul(const p256_int* n, p256_int* out_x, p256_int* out_y) { 1226 felem x, y, z; 1227 1228 scalar_base_mult(x, y, z, n); 1229 1230 { 1231 felem x_affine, y_affine; 1232 1233 point_to_affine(x_affine, y_affine, x, y, z); 1234 from_montgomery(out_x, x_affine); 1235 from_montgomery(out_y, y_affine); 1236 } 1237 } 1238 1239 /* p256_points_mul_vartime sets {out_x,out_y} = n1*G + n2*{in_x,in_y}, where 1240 * n1 and n2 are < the order of the group. 1241 * 1242 * As indicated by the name, this function operates in variable time. This 1243 * is safe because it's used for signature validation which doesn't deal 1244 * with secrets. */ 1245 void p256_points_mul_vartime( 1246 const p256_int* n1, const p256_int* n2, const p256_int* in_x, 1247 const p256_int* in_y, p256_int* out_x, p256_int* out_y) { 1248 felem x1, y1, z1, x2, y2, z2, px, py; 1249 1250 /* If both scalars are zero, then the result is the point at infinity. */ 1251 if (p256_is_zero(n1) != 0 && p256_is_zero(n2) != 0) { 1252 p256_clear(out_x); 1253 p256_clear(out_y); 1254 return; 1255 } 1256 1257 to_montgomery(px, in_x); 1258 to_montgomery(py, in_y); 1259 scalar_base_mult(x1, y1, z1, n1); 1260 scalar_mult(x2, y2, z2, px, py, n2); 1261 1262 if (p256_is_zero(n2) != 0) { 1263 /* If n2 == 0, then {x2,y2,z2} is zero and the result is just 1264 * {x1,y1,z1}. */ 1265 } else if (p256_is_zero(n1) != 0) { 1266 /* If n1 == 0, then {x1,y1,z1} is zero and the result is just 1267 * {x2,y2,z2}. */ 1268 memcpy(x1, x2, sizeof(x2)); 1269 memcpy(y1, y2, sizeof(y2)); 1270 memcpy(z1, z2, sizeof(z2)); 1271 } else { 1272 /* This function handles the case where {x1,y1,z1} == {x2,y2,z2}. */ 1273 point_add_or_double_vartime(x1, y1, z1, x1, y1, z1, x2, y2, z2); 1274 } 1275 1276 point_to_affine(px, py, x1, y1, z1); 1277 from_montgomery(out_x, px); 1278 from_montgomery(out_y, py); 1279 } 1280
