f25519.c

   1 /* Arithmetic mod p = 2^255-19
   2  * Daniel Beer <dlbeer@gmail.com>, 5 Jan 2014
   3  *
   4  * This file is in the public domain.
   5  */
   6
   7 #include "f25519.h"
   8
   9 const uint8_t f25519_one[F25519_SIZE] = {1};
  10
  11 void f25519_load(uint8_t *x, uint32_t c)
  12 {
  13         int i;
  14
  15         for (i = 0; i < sizeof(c); i++) {
  16                 x[i] = c;
  17                 c >>= 8;
  18         }
  19
  20         for (; i < F25519_SIZE; i++)
  21                 x[i] = 0;
  22 }
  23
  24 void f25519_normalize(uint8_t *x)
  25 {
  26         uint8_t minusp[F25519_SIZE];
  27         uint16_t c;
  28         int i;
  29
  30         /* Reduce using 2^255 = 19 mod p */
  31         c = (x[31] >> 7) * 19;
  32         x[31] &= 127;
  33
  34         for (i = 0; i < F25519_SIZE; i++) {
  35                 c += x[i];
  36                 x[i] = c;
  37                 c >>= 8;
  38         }
  39
  40         /* The number is now less than 2^255 + 18, and therefore less than
  41          * 2p. Try subtracting p, and conditionally load the subtracted
  42          * value if underflow did not occur.
  43          */
  44         c = 19;
  45
  46         for (i = 0; i + 1 < F25519_SIZE; i++) {
  47                 c += x[i];
  48                 minusp[i] = c;
  49                 c >>= 8;
  50         }
  51
  52         c += ((uint16_t)x[i]) - 128;
  53         minusp[31] = c;
  54
  55         /* Load x-p if no underflow */
  56         f25519_select(x, minusp, x, (c >> 15) & 1);
  57 }
  58
  59 uint8_t f25519_eq(const uint8_t *x, const uint8_t *y)
  60 {
  61         uint8_t sum = 0;
  62         int i;
  63
  64         for (i = 0; i < F25519_SIZE; i++)
  65                 sum |= x[i] ^ y[i];
  66
  67         sum |= (sum >> 4);
  68         sum |= (sum >> 2);
  69         sum |= (sum >> 1);
  70
  71         return (sum ^ 1) & 1;
  72 }
  73
  74 void f25519_select(uint8_t *dst,
  75                    const uint8_t *zero, const uint8_t *one,
  76                    uint8_t condition)
  77 {
  78         const uint8_t mask = -condition;
  79         int i;
  80
  81         for (i = 0; i < F25519_SIZE; i++)
  82                 dst[i] = zero[i] ^ (mask & (one[i] ^ zero[i]));
  83 }
  84
  85 void f25519_add(uint8_t *r, const uint8_t *a, const uint8_t *b)
  86 {
  87         uint16_t c = 0;
  88         int i;
  89
  90         /* Add */
  91         for (i = 0; i < F25519_SIZE; i++) {
  92                 c >>= 8;
  93                 c += ((uint16_t)a[i]) + ((uint16_t)b[i]);
  94                 r[i] = c;
  95         }
  96
  97         /* Reduce with 2^255 = 19 mod p */
  98         r[31] &= 127;
  99         c = (c >> 7) * 19;
 100
 101         for (i = 0; i < F25519_SIZE; i++) {
 102                 c += r[i];
 103                 r[i] = c;
 104                 c >>= 8;
 105         }
 106 }
 107
 108 void f25519_sub(uint8_t *r, const uint8_t *a, const uint8_t *b)
 109 {
 110         uint32_t c = 0;
 111         int i;
 112
 113         /* Calculate a + 2p - b, to avoid underflow */
 114         c = 218;
 115         for (i = 0; i + 1 < F25519_SIZE; i++) {
 116                 c += 65280 + ((uint32_t)a[i]) - ((uint32_t)b[i]);
 117                 r[i] = c;
 118                 c >>= 8;
 119         }
 120
 121         c += ((uint32_t)a[31]) - ((uint32_t)b[31]);
 122         r[31] = c & 127;
 123         c = (c >> 7) * 19;
 124
 125         for (i = 0; i < F25519_SIZE; i++) {
 126                 c += r[i];
 127                 r[i] = c;
 128                 c >>= 8;
 129         }
 130 }
 131
 132 void f25519_neg(uint8_t *r, const uint8_t *a)
 133 {
 134         uint32_t c = 0;
 135         int i;
 136
 137         /* Calculate 2p - a, to avoid underflow */
 138         c = 218;
 139         for (i = 0; i + 1 < F25519_SIZE; i++) {
 140                 c += 65280 - ((uint32_t)a[i]);
 141                 r[i] = c;
 142                 c >>= 8;
 143         }
 144
 145         c -= ((uint32_t)a[31]);
 146         r[31] = c & 127;
 147         c = (c >> 7) * 19;
 148
 149         for (i = 0; i < F25519_SIZE; i++) {
 150                 c += r[i];
 151                 r[i] = c;
 152                 c >>= 8;
 153         }
 154 }
 155
 156 void f25519_mul__distinct(uint8_t *r, const uint8_t *a, const uint8_t *b)
 157 {
 158         uint32_t c = 0;
 159         int i;
 160
 161         for (i = 0; i < F25519_SIZE; i++) {
 162                 int j;
 163
 164                 c >>= 8;
 165                 for (j = 0; j <= i; j++)
 166                         c += ((uint32_t)a[j]) * ((uint32_t)b[i - j]);
 167
 168                 for (; j < F25519_SIZE; j++)
 169                         c += ((uint32_t)a[j]) *
 170                              ((uint32_t)b[i + F25519_SIZE - j]) * 38;
 171
 172                 r[i] = c;
 173         }
 174
 175         r[31] &= 127;
 176         c = (c >> 7) * 19;
 177
 178         for (i = 0; i < F25519_SIZE; i++) {
 179                 c += r[i];
 180                 r[i] = c;
 181                 c >>= 8;
 182         }
 183 }
 184
 185 static void f25519_mul_c(uint8_t *r, const uint8_t *a, uint32_t b)
 186 {
 187         uint32_t c = 0;
 188         int i;
 189
 190         for (i = 0; i < F25519_SIZE; i++) {
 191                 c >>= 8;
 192                 c += b * ((uint32_t)a[i]);
 193                 r[i] = c;
 194         }
 195
 196         r[31] &= 127;
 197         c >>= 7;
 198         c *= 19;
 199
 200         for (i = 0; i < F25519_SIZE; i++) {
 201                 c += r[i];
 202                 r[i] = c;
 203                 c >>= 8;
 204         }
 205 }
 206
 207 void f25519_inv__distinct(uint8_t *r, const uint8_t *x)
 208 {
 209         uint8_t s[F25519_SIZE];
 210         int i;
 211
 212         /* This is a prime field, so by Fermat's little theorem:
 213          *
 214          *     x^(p-1) = 1 mod p
 215          *
 216          * Therefore, raise to (p-2) = 2^255-21 to get a multiplicative
 217          * inverse.
 218          *
 219          * This is a 255-bit binary number with the digits:
 220          *
 221          *     11111111... 01011
 222          *
 223          * We compute the result by the usual binary chain, but
 224          * alternate between keeping the accumulator in r and s, so as
 225          * to avoid copying temporaries.
 226          */
 227
 228         /* 1 1 */
 229         f25519_mul__distinct(s, x, x);
 230         f25519_mul__distinct(r, s, x);
 231
 232         /* 1 x 248 */
 233         for (i = 0; i < 248; i++) {
 234                 f25519_mul__distinct(s, r, r);
 235                 f25519_mul__distinct(r, s, x);
 236         }
 237
 238         /* 0 */
 239         f25519_mul__distinct(s, r, r);
 240
 241         /* 1 */
 242         f25519_mul__distinct(r, s, s);
 243         f25519_mul__distinct(s, r, x);
 244
 245         /* 0 */
 246         f25519_mul__distinct(r, s, s);
 247
 248         /* 1 */
 249         f25519_mul__distinct(s, r, r);
 250         f25519_mul__distinct(r, s, x);
 251
 252         /* 1 */
 253         f25519_mul__distinct(s, r, r);
 254         f25519_mul__distinct(r, s, x);
 255 }
 256
 257 /* Raise x to the power of (p-5)/8 = 2^252-3, using s for temporary
 258  * storage.
 259  */
 260 static void exp2523(uint8_t *r, const uint8_t *x, uint8_t *s)
 261 {
 262         int i;
 263
 264         /* This number is a 252-bit number with the binary expansion:
 265          *
 266          *     111111... 01
 267          */
 268
 269         /* 1 1 */
 270         f25519_mul__distinct(r, x, x);
 271         f25519_mul__distinct(s, r, x);
 272
 273         /* 1 x 248 */
 274         for (i = 0; i < 248; i++) {
 275                 f25519_mul__distinct(r, s, s);
 276                 f25519_mul__distinct(s, r, x);
 277         }
 278
 279         /* 0 */
 280         f25519_mul__distinct(r, s, s);
 281
 282         /* 1 */
 283         f25519_mul__distinct(s, r, r);
 284         f25519_mul__distinct(r, s, x);
 285 }
 286
 287 void f25519_sqrt(uint8_t *r, const uint8_t *a)
 288 {
 289         uint8_t v[F25519_SIZE];
 290         uint8_t i[F25519_SIZE];
 291         uint8_t x[F25519_SIZE];
 292         uint8_t y[F25519_SIZE];
 293
 294         /* v = (2a)^((p-5)/8) [x = 2a] */
 295         f25519_mul_c(x, a, 2);
 296         exp2523(v, x, y);
 297
 298         /* i = 2av^2 - 1 */
 299         f25519_mul__distinct(y, v, v);
 300         f25519_mul__distinct(i, x, y);
 301         f25519_load(y, 1);
 302         f25519_sub(i, i, y);
 303
 304         /* r = avi */
 305         f25519_mul__distinct(x, v, a);
 306         f25519_mul__distinct(r, x, i);
 307 }