decaf.tmpl.c 46.6 KB
Newer Older
1
/** @brief Decaf high-level functions. */
Michael Hamburg's avatar
Michael Hamburg committed
2

3 4 5
#define _XOPEN_SOURCE 600 /* for posix_memalign */
#include "word.h"
#include "field.h"
Michael Hamburg's avatar
Michael Hamburg committed
6

7
#include <decaf.h>
8
#include <decaf/ed$(gf_bits).h>
9

10
/* Template stuff */
Michael Hamburg's avatar
Michael Hamburg committed
11 12
#define API_NS(_id) $(c_ns)_##_id
#define SCALAR_BITS $(C_NS)_SCALAR_BITS
13
#define SCALAR_SER_BYTES $(C_NS)_SCALAR_BYTES
Michael Hamburg's avatar
Michael Hamburg committed
14 15 16 17 18 19 20
#define SCALAR_LIMBS $(C_NS)_SCALAR_LIMBS
#define scalar_t API_NS(scalar_t)
#define point_t API_NS(point_t)
#define precomputed_s API_NS(precomputed_s)
#define IMAGINE_TWIST $(imagine_twist)
#define COFACTOR $(cofactor)

21
/* Comb config: number of combs, n, t, s. */
22 23 24 25 26 27 28
#define COMBS_N $(combs.n)
#define COMBS_T $(combs.t)
#define COMBS_S $(combs.s)
#define DECAF_WINDOW_BITS $(window_bits)
#define DECAF_WNAF_FIXED_TABLE_BITS $(wnaf.fixed)
#define DECAF_WNAF_VAR_TABLE_BITS $(wnaf.var)

Michael Hamburg's avatar
Michael Hamburg committed
29 30
#define EDDSA_USE_SIGMA_ISOGENY $(eddsa_sigma_iso)

Michael Hamburg's avatar
Michael Hamburg committed
31
static const int EDWARDS_D = $(d);
32
static const scalar_t point_scalarmul_adjustment = {{{
33 34 35 36
    $(ser((2**(scalar_bits-1+window_bits - ((scalar_bits-1)%window_bits)) - 1) % q,64,"SC_LIMB"))
}}}, precomputed_scalarmul_adjustment = {{{
    $(ser((2**(combs.n*combs.t*combs.s) - 1) % q,64,"SC_LIMB"))
}}};
Michael Hamburg's avatar
Michael Hamburg committed
37

38
const uint8_t decaf_x$(gf_shortname)_base_point[DECAF_X$(gf_shortname)_PUBLIC_BYTES] = { $(ser(mont_base,8)) };
39

40 41
#define RISTRETTO_FACTOR $(C_NS)_RISTRETTO_FACTOR
const gf RISTRETTO_FACTOR = {{{
42
    $(ser(msqrt(d-1 if imagine_twist else -d,modulus,lo_bit_clear=True),gf_lit_limb_bits))
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
}}};

#if IMAGINE_TWIST
#define TWISTED_D (-(EDWARDS_D))
#else
#define TWISTED_D ((EDWARDS_D)-1)
#endif

#if TWISTED_D < 0
#define EFF_D (-(TWISTED_D))
#define NEG_D 1
#else
#define EFF_D TWISTED_D
#define NEG_D 0
#endif

59 60
/* End of template stuff */

61
/* Sanity */
62
#if (COFACTOR == 8) && !IMAGINE_TWIST && !UNSAFE_CURVE_HAS_POINTS_AT_INFINITY
63
/* FUTURE MAGIC: Curve41417 doesn't have these properties. */
64 65 66 67 68
#error "Currently require IMAGINE_TWIST (and thus p=5 mod 8) for cofactor 8"
        /* OK, but why?
         * Two reasons: #1: There are bugs when COFACTOR == && IMAGINE_TWIST
         # #2: 
         */
69
#endif
70 71

#if IMAGINE_TWIST && (P_MOD_8 != 5)
72
    #error "Cannot use IMAGINE_TWIST except for p == 5 mod 8"
73 74
#endif

75
#if (COFACTOR != 8) && (COFACTOR != 4)
76
    #error "COFACTOR must be 4 or 8"
77 78 79
#endif
 
#if IMAGINE_TWIST
80
    extern const gf SQRT_MINUS_ONE;
81 82
#endif

83
#define WBITS DECAF_WORD_BITS /* NB this may be different from ARCH_WORD_BITS */
84

85
extern const point_t API_NS(point_base);
Michael Hamburg's avatar
Michael Hamburg committed
86

87 88
/* Projective Niels coordinates */
typedef struct { gf a, b, c; } niels_s, niels_t[1];
89
typedef struct { niels_t n; gf z; } VECTOR_ALIGNED pniels_s, pniels_t[1];
Michael Hamburg's avatar
Michael Hamburg committed
90

91
/* Precomputed base */
92
struct precomputed_s { niels_t table [COMBS_N<<(COMBS_T-1)]; };
Michael Hamburg's avatar
Michael Hamburg committed
93

94 95 96
extern const gf API_NS(precomputed_base_as_fe)[];
const precomputed_s *API_NS(precomputed_base) =
    (const precomputed_s *) &API_NS(precomputed_base_as_fe);
Michael Hamburg's avatar
Michael Hamburg committed
97

Michael Hamburg's avatar
Michael Hamburg committed
98 99
const size_t API_NS(sizeof_precomputed_s) = sizeof(precomputed_s);
const size_t API_NS(alignof_precomputed_s) = sizeof(big_register_t);
100 101 102

/** Inverse. */
static void
103
gf_invert(gf y, const gf x, int assert_nonzero) {
104 105
    gf t1, t2;
    gf_sqr(t1, x); // o^2
106
    mask_t ret = gf_isr(t2, t1); // +-1/sqrt(o^2) = +-1/o
107 108
    (void)ret;
    if (assert_nonzero) assert(ret);
109 110
    gf_sqr(t1, t2);
    gf_mul(t2, t1, x); // not direct to y in case of alias.
111
    gf_copy(y, t2);
112 113
}

114
/** identity = (0,1) */
115
const point_t API_NS(point_identity) = {{{{{0}}},{{{1}}},{{{1}}},{{{0}}}}};
Michael Hamburg's avatar
Michael Hamburg committed
116

117
/* Predeclare because not static: called by elligator */
118
void API_NS(deisogenize) (
119
    gf_s *__restrict__ s,
120 121
    gf_s *__restrict__ inv_el_sum,
    gf_s *__restrict__ inv_el_m1,
122
    const point_t p,
123
    mask_t toggle_s,
124
    mask_t toggle_altx,
125 126 127
    mask_t toggle_rotation
);

128
void API_NS(deisogenize) (
129
    gf_s *__restrict__ s,
130 131
    gf_s *__restrict__ inv_el_sum,
    gf_s *__restrict__ inv_el_m1,
132
    const point_t p,
133
    mask_t toggle_s,
134
    mask_t toggle_altx,
135
    mask_t toggle_rotation
136
) {
137
#if COFACTOR == 4 && !IMAGINE_TWIST
138
    (void)toggle_rotation; /* Only applies to cofactor 8 */
139 140
    gf t1;
    gf_s *t2 = s, *t3=inv_el_sum, *t4=inv_el_m1;
141
    
142 143 144 145 146 147 148 149
    gf_add(t1,p->x,p->t);
    gf_sub(t2,p->x,p->t);
    gf_mul(t3,t1,t2); /* t3 = num */
    gf_sqr(t2,p->x);
    gf_mul(t1,t2,t3);
    gf_mulw(t2,t1,-1-TWISTED_D); /* -x^2 * (a-d) * num */
    gf_isr(t1,t2);    /* t1 = isr */
    gf_mul(t2,t1,t3); /* t2 = ratio */
150
    gf_mul(t4,t2,RISTRETTO_FACTOR);
151
    mask_t negx = gf_lobit(t4) ^ toggle_altx;
152 153 154 155
    gf_cond_neg(t2, negx);
    gf_mul(t3,t2,p->z);
    gf_sub(t3,t3,p->t);
    gf_mul(t2,t3,p->x);
156 157 158 159 160 161 162 163 164
    gf_mulw(t4,t2,-1-TWISTED_D);
    gf_mul(s,t4,t1);
    mask_t lobs = gf_lobit(s);
    gf_cond_neg(s,lobs);
    gf_copy(inv_el_m1,p->x);
    gf_cond_neg(inv_el_m1,~lobs^negx^toggle_s);
    gf_add(inv_el_m1,inv_el_m1,p->t);
    
#elif COFACTOR == 8 && IMAGINE_TWIST
165
    /* More complicated because of rotation */
166
    gf t1,t2,t3,t4,t5;
167 168 169 170 171 172 173 174 175
    gf_add(t1,p->z,p->y);
    gf_sub(t2,p->z,p->y);
    gf_mul(t3,t1,t2);      /* t3 = num */
    gf_mul(t2,p->x,p->y);  /* t2 = den */
    gf_sqr(t1,t2);
    gf_mul(t4,t1,t3);
    gf_mulw(t1,t4,-1-TWISTED_D);
    gf_isr(t4,t1);         /* isqrt(num*(a-d)*den^2) */
    gf_mul(t1,t2,t4);
176
    gf_mul(t2,t1,RISTRETTO_FACTOR); /* t2 = "iden" in ristretto.sage */
177 178
    gf_mul(t1,t3,t4);                 /* t1 = "inum" in ristretto.sage */

179
    /* Calculate altxy = iden*inum*i*t^2*(d-a) */
180
    gf_mul(t3,t1,t2);
181
    gf_mul_i(t4,t3);
182 183 184 185
    gf_mul(t3,t4,p->t);
    gf_mul(t4,t3,p->t);
    gf_mulw(t3,t4,TWISTED_D+1);      /* iden*inum*i*t^2*(d-1) */
    mask_t rotate = toggle_rotation ^ gf_lobit(t3);
186
    
187
    /* Rotate if altxy is negative */
188
    gf_cond_swap(t1,t2,rotate);
189
    gf_mul_i(t4,p->x);
190 191
    gf_cond_sel(t4,p->y,t4,rotate);  /* t4 = "fac" = ix if rotate, else y */
    
192
    gf_mul_i(t5,RISTRETTO_FACTOR); /* t5 = imi */
193 194 195 196 197
    gf_mul(t3,t5,t2);                /* iden * imi */
    gf_mul(t2,t5,t1);
    gf_mul(t5,t2,p->t);              /* "altx" = iden*imi*t */
    mask_t negx = gf_lobit(t5) ^ toggle_altx;
    
198
    gf_cond_neg(t1,negx^rotate);
199 200
    gf_mul(t2,t1,p->z);
    gf_add(t2,t2,ONE);
201 202 203
    gf_mul(inv_el_sum,t2,t4);
    gf_mul(s,inv_el_sum,t3);
    
204 205 206 207 208 209
    mask_t negs = gf_lobit(s);
    gf_cond_neg(s,negs);
    
    mask_t negz = ~negs ^ toggle_s ^ negx;
    gf_copy(inv_el_m1,p->z);
    gf_cond_neg(inv_el_m1,negz);
210
    gf_sub(inv_el_m1,inv_el_m1,t4);
211
#else
212
#error "Cofactor must be 4 (with no IMAGINE_TWIST) or 8 (with IMAGINE_TWIST)"
213 214 215 216
#endif
}

void API_NS(point_encode)( unsigned char ser[SER_BYTES], const point_t p ) {
217
    gf s,ie1,ie2;
218
    API_NS(deisogenize)(s,ie1,ie2,p,0,0,0);
219
    gf_serialize(ser,s,1);
220
}
221 222 223 224 225 226

decaf_error_t API_NS(point_decode) (
    point_t p,
    const unsigned char ser[SER_BYTES],
    decaf_bool_t allow_identity
) {
227 228 229 230 231 232 233
    gf s, s2, num, tmp;
    gf_s *tmp2=s2, *ynum=p->z, *isr=p->x, *den=p->t;
    
    mask_t succ = gf_deserialize(s, ser, 1);
    succ &= bool_to_mask(allow_identity) | ~gf_eq(s, ZERO);
    succ &= ~gf_lobit(s);
    
234
    gf_sqr(s2,s);                  /* s^2 = -as^2 */
235
#if IMAGINE_TWIST
236
    gf_sub(s2,ZERO,s2);            /* -as^2 */
237
#endif
238 239 240 241 242 243 244 245 246 247 248 249 250
    gf_sub(den,ONE,s2);            /* 1+as^2 */
    gf_add(ynum,ONE,s2);           /* 1-as^2 */
    gf_mulw(num,s2,-4*TWISTED_D);
    gf_sqr(tmp,den);               /* tmp = den^2 */
    gf_add(num,tmp,num);           /* num = den^2 - 4*d*s^2 */
    gf_mul(tmp2,num,tmp);          /* tmp2 = num*den^2 */
    succ &= gf_isr(isr,tmp2);      /* isr = 1/sqrt(num*den^2) */
    gf_mul(tmp,isr,den);           /* isr*den */
    gf_mul(p->y,tmp,ynum);         /* isr*den*(1-as^2) */
    gf_mul(tmp2,tmp,s);            /* s*isr*den */
    gf_add(tmp2,tmp2,tmp2);        /* 2*s*isr*den */
    gf_mul(tmp,tmp2,isr);          /* 2*s*isr^2*den */
    gf_mul(p->x,tmp,num);          /* 2*s*isr^2*den*num */
251
    gf_mul(tmp,tmp2,RISTRETTO_FACTOR); /* 2*s*isr*den*magic */
252
    gf_cond_neg(p->x,gf_lobit(tmp)); /* flip x */
253
    
254 255 256 257
#if COFACTOR==8
    /* Additionally check y != 0 and x*y*isomagic nonegative */
    succ &= ~gf_eq(p->y,ZERO);
    gf_mul(tmp,p->x,p->y);
258
    gf_mul(tmp2,tmp,RISTRETTO_FACTOR);
259
    succ &= ~gf_lobit(tmp2);
260 261 262
#endif

#if IMAGINE_TWIST
263
    gf_copy(tmp,p->x);
264
    gf_mul_i(p->x,tmp);
265
#endif
266 267 268 269

    /* Fill in z and t */
    gf_copy(p->z,ONE);
    gf_mul(p->t,p->x,p->y);
270 271
    
    assert(API_NS(point_valid)(p) | ~succ);
272
    return decaf_succeed_if(mask_to_bool(succ));
273 274 275 276 277 278
}

void API_NS(point_sub) (
    point_t p,
    const point_t q,
    const point_t r
279 280
) {
    gf a, b, c, d;
281 282 283
    gf_sub_nr ( b, q->y, q->x ); /* 3+e */
    gf_sub_nr ( d, r->y, r->x ); /* 3+e */
    gf_add_nr ( c, r->y, r->x ); /* 2+e */
284
    gf_mul ( a, c, b );
285
    gf_add_nr ( b, q->y, q->x ); /* 2+e */
286 287
    gf_mul ( p->y, d, b );
    gf_mul ( b, r->t, q->t );
288
    gf_mulw ( p->x, b, 2*EFF_D );
289 290
    gf_add_nr ( b, a, p->y );    /* 2+e */
    gf_sub_nr ( c, p->y, a );    /* 3+e */
291
    gf_mul ( a, q->z, r->z );
292 293
    gf_add_nr ( a, a, a );       /* 2+e */
    if (GF_HEADROOM <= 3) gf_weak_reduce(a); /* or 1+e */
294
#if NEG_D
295 296
    gf_sub_nr ( p->y, a, p->x ); /* 4+e or 3+e */
    gf_add_nr ( a, a, p->x );    /* 3+e or 2+e */
297
#else
298 299
    gf_add_nr ( p->y, a, p->x ); /* 3+e or 2+e */
    gf_sub_nr ( a, a, p->x );    /* 4+e or 3+e */
300
#endif
301 302 303 304
    gf_mul ( p->z, a, p->y );
    gf_mul ( p->x, p->y, c );
    gf_mul ( p->y, a, b );
    gf_mul ( p->t, b, c );
305
}
306
    
307 308 309 310
void API_NS(point_add) (
    point_t p,
    const point_t q,
    const point_t r
311
) {
312
    gf a, b, c, d;
313 314 315
    gf_sub_nr ( b, q->y, q->x ); /* 3+e */
    gf_sub_nr ( c, r->y, r->x ); /* 3+e */
    gf_add_nr ( d, r->y, r->x ); /* 2+e */
316
    gf_mul ( a, c, b );
317
    gf_add_nr ( b, q->y, q->x ); /* 2+e */
318 319
    gf_mul ( p->y, d, b );
    gf_mul ( b, r->t, q->t );
320
    gf_mulw ( p->x, b, 2*EFF_D );
321 322
    gf_add_nr ( b, a, p->y );    /* 2+e */
    gf_sub_nr ( c, p->y, a );    /* 3+e */
323
    gf_mul ( a, q->z, r->z );
324 325
    gf_add_nr ( a, a, a );       /* 2+e */
    if (GF_HEADROOM <= 3) gf_weak_reduce(a); /* or 1+e */
326
#if NEG_D
327 328
    gf_add_nr ( p->y, a, p->x ); /* 3+e or 2+e */
    gf_sub_nr ( a, a, p->x );    /* 4+e or 3+e */
329
#else
330 331
    gf_sub_nr ( p->y, a, p->x ); /* 4+e or 3+e */
    gf_add_nr ( a, a, p->x );    /* 3+e or 2+e */
332 333 334 335 336
#endif
    gf_mul ( p->z, a, p->y );
    gf_mul ( p->x, p->y, c );
    gf_mul ( p->y, a, b );
    gf_mul ( p->t, b, c );
Michael Hamburg's avatar
Michael Hamburg committed
337
}
338

339
static DECAF_NOINLINE void
340 341 342
point_double_internal (
    point_t p,
    const point_t q,
343
    int before_double
344 345 346 347
) {
    gf a, b, c, d;
    gf_sqr ( c, q->x );
    gf_sqr ( a, q->y );
348 349
    gf_add_nr ( d, c, a );             /* 2+e */
    gf_add_nr ( p->t, q->y, q->x );    /* 2+e */
350
    gf_sqr ( b, p->t );
351 352
    gf_subx_nr ( b, b, d, 3 );         /* 4+e */
    gf_sub_nr ( p->t, a, c );          /* 3+e */
353
    gf_sqr ( p->x, q->z );
354 355 356
    gf_add_nr ( p->z, p->x, p->x );    /* 2+e */
    gf_subx_nr ( a, p->z, p->t, 4 );   /* 6+e */
    if (GF_HEADROOM == 5) gf_weak_reduce(a); /* or 1+e */
357 358 359 360
    gf_mul ( p->x, a, b );
    gf_mul ( p->z, p->t, a );
    gf_mul ( p->y, p->t, d );
    if (!before_double) gf_mul ( p->t, b, d );
361 362
}

363 364
void API_NS(point_double)(point_t p, const point_t q) {
    point_double_internal(p,q,0);
365
}
366

367 368 369
void API_NS(point_negate) (
   point_t nega,
   const point_t a
Mike Hamburg's avatar
Mike Hamburg committed
370 371
) {
    gf_sub(nega->x, ZERO, a->x);
372 373
    gf_copy(nega->y, a->y);
    gf_copy(nega->z, a->z);
Mike Hamburg's avatar
Mike Hamburg committed
374 375 376
    gf_sub(nega->t, ZERO, a->t);
}

377
/* Operations on [p]niels */
378
static DECAF_INLINE void
379 380
cond_neg_niels (
    niels_t n,
381
    mask_t neg
Michael Hamburg's avatar
Michael Hamburg committed
382
) {
383 384
    gf_cond_swap(n->a, n->b, neg);
    gf_cond_neg(n->c, neg);
385 386
}

387
static DECAF_NOINLINE void pt_to_pniels (
388 389 390 391 392
    pniels_t b,
    const point_t a
) {
    gf_sub ( b->n->a, a->y, a->x );
    gf_add ( b->n->b, a->x, a->y );
393
    gf_mulw ( b->n->c, a->t, 2*TWISTED_D );
394 395 396
    gf_add ( b->z, a->z, a->z );
}

397
static DECAF_NOINLINE void pniels_to_pt (
398 399 400 401 402 403 404 405 406 407 408 409
    point_t e,
    const pniels_t d
) {
    gf eu;
    gf_add ( eu, d->n->b, d->n->a );
    gf_sub ( e->y, d->n->b, d->n->a );
    gf_mul ( e->t, e->y, eu);
    gf_mul ( e->x, d->z, e->y );
    gf_mul ( e->y, d->z, eu );
    gf_sqr ( e->z, d->z );
}

410
static DECAF_NOINLINE void
411 412 413 414 415 416 417
niels_to_pt (
    point_t e,
    const niels_t n
) {
    gf_add ( e->y, n->b, n->a );
    gf_sub ( e->x, n->b, n->a );
    gf_mul ( e->t, e->y, e->x );
418
    gf_copy ( e->z, ONE );
419 420
}

421
static DECAF_NOINLINE void
422 423 424
add_niels_to_pt (
    point_t d,
    const niels_t e,
425
    int before_double
426 427
) {
    gf a, b, c;
428
    gf_sub_nr ( b, d->y, d->x ); /* 3+e */
429
    gf_mul ( a, e->a, b );
430
    gf_add_nr ( b, d->x, d->y ); /* 2+e */
431 432
    gf_mul ( d->y, e->b, b );
    gf_mul ( d->x, e->c, d->t );
433 434 435 436
    gf_add_nr ( c, a, d->y );    /* 2+e */
    gf_sub_nr ( b, d->y, a );    /* 3+e */
    gf_sub_nr ( d->y, d->z, d->x ); /* 3+e */
    gf_add_nr ( a, d->x, d->z ); /* 2+e */
437 438 439 440 441 442
    gf_mul ( d->z, a, d->y );
    gf_mul ( d->x, d->y, b );
    gf_mul ( d->y, a, c );
    if (!before_double) gf_mul ( d->t, b, c );
}

443
static DECAF_NOINLINE void
444 445 446
sub_niels_from_pt (
    point_t d,
    const niels_t e,
447
    int before_double
448 449
) {
    gf a, b, c;
450
    gf_sub_nr ( b, d->y, d->x ); /* 3+e */
451
    gf_mul ( a, e->b, b );
452
    gf_add_nr ( b, d->x, d->y ); /* 2+e */
453 454
    gf_mul ( d->y, e->a, b );
    gf_mul ( d->x, e->c, d->t );
455 456 457 458
    gf_add_nr ( c, a, d->y );    /* 2+e */
    gf_sub_nr ( b, d->y, a );    /* 3+e */
    gf_add_nr ( d->y, d->z, d->x ); /* 2+e */
    gf_sub_nr ( a, d->z, d->x ); /* 3+e */
459 460 461 462 463 464 465 466 467 468
    gf_mul ( d->z, a, d->y );
    gf_mul ( d->x, d->y, b );
    gf_mul ( d->y, a, c );
    if (!before_double) gf_mul ( d->t, b, c );
}

static void
add_pniels_to_pt (
    point_t p,
    const pniels_t pn,
469
    int before_double
470 471 472
) {
    gf L0;
    gf_mul ( L0, p->z, pn->z );
473
    gf_copy ( p->z, L0 );
474 475 476 477 478 479 480
    add_niels_to_pt( p, pn->n, before_double );
}

static void
sub_pniels_from_pt (
    point_t p,
    const pniels_t pn,
481
    int before_double
482 483 484
) {
    gf L0;
    gf_mul ( L0, p->z, pn->z );
485
    gf_copy ( p->z, L0 );
486 487 488
    sub_niels_from_pt( p, pn->n, before_double );
}

489
static DECAF_NOINLINE void
490 491 492 493 494 495 496
prepare_fixed_window(
    pniels_t *multiples,
    const point_t b,
    int ntable
) {
    point_t tmp;
    pniels_t pn;
497
    int i;
498 499 500 501 502 503 504 505
    
    point_double_internal(tmp, b, 0);
    pt_to_pniels(pn, tmp);
    pt_to_pniels(multiples[0], b);
    API_NS(point_copy)(tmp, b);
    for (i=1; i<ntable; i++) {
        add_pniels_to_pt(tmp, pn, 0);
        pt_to_pniels(multiples[i], tmp);
506
    }
507 508 509 510 511 512 513 514 515
    
    decaf_bzero(pn,sizeof(pn));
    decaf_bzero(tmp,sizeof(tmp));
}

void API_NS(point_scalarmul) (
    point_t a,
    const point_t b,
    const scalar_t scalar
516
) {
517 518 519 520 521 522
    const int WINDOW = DECAF_WINDOW_BITS,
        WINDOW_MASK = (1<<WINDOW)-1,
        WINDOW_T_MASK = WINDOW_MASK >> 1,
        NTABLE = 1<<(WINDOW-1);
        
    scalar_t scalar1x;
523
    API_NS(scalar_add)(scalar1x, scalar, point_scalarmul_adjustment);
524
    API_NS(scalar_halve)(scalar1x,scalar1x);
525 526 527 528 529 530 531 532 533 534 535 536
    
    /* Set up a precomputed table with odd multiples of b. */
    pniels_t pn, multiples[NTABLE];
    point_t tmp;
    prepare_fixed_window(multiples, b, NTABLE);

    /* Initialize. */
    int i,j,first=1;
    i = SCALAR_BITS - ((SCALAR_BITS-1) % WINDOW) - 1;

    for (; i>=0; i-=WINDOW) {
        /* Fetch another block of bits */
537
        word_t bits = scalar1x->limb[i/WBITS] >> (i%WBITS);
538 539 540 541
        if (i%WBITS >= WBITS-WINDOW && i/WBITS<SCALAR_LIMBS-1) {
            bits ^= scalar1x->limb[i/WBITS+1] << (WBITS - (i%WBITS));
        }
        bits &= WINDOW_MASK;
542
        mask_t inv = (bits>>(WINDOW-1))-1;
543 544 545
        bits ^= inv;
    
        /* Add in from table.  Compute t only on last iteration. */
546
        constant_time_lookup(pn, multiples, sizeof(pn), NTABLE, bits & WINDOW_T_MASK);
547 548 549 550 551 552 553 554 555 556 557 558 559 560
        cond_neg_niels(pn->n, inv);
        if (first) {
            pniels_to_pt(tmp, pn);
            first = 0;
        } else {
           /* Using Hisil et al's lookahead method instead of extensible here
            * for no particular reason.  Double WINDOW times, but only compute t on
            * the last one.
            */
            for (j=0; j<WINDOW-1; j++)
                point_double_internal(tmp, tmp, -1);
            point_double_internal(tmp, tmp, 0);
            add_pniels_to_pt(tmp, pn, i ? -1 : 0);
        }
561
    }
562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584
    
    /* Write out the answer */
    API_NS(point_copy)(a,tmp);
    
    decaf_bzero(scalar1x,sizeof(scalar1x));
    decaf_bzero(pn,sizeof(pn));
    decaf_bzero(multiples,sizeof(multiples));
    decaf_bzero(tmp,sizeof(tmp));
}

void API_NS(point_double_scalarmul) (
    point_t a,
    const point_t b,
    const scalar_t scalarb,
    const point_t c,
    const scalar_t scalarc
) {
    const int WINDOW = DECAF_WINDOW_BITS,
        WINDOW_MASK = (1<<WINDOW)-1,
        WINDOW_T_MASK = WINDOW_MASK >> 1,
        NTABLE = 1<<(WINDOW-1);
        
    scalar_t scalar1x, scalar2x;
585
    API_NS(scalar_add)(scalar1x, scalarb, point_scalarmul_adjustment);
586
    API_NS(scalar_halve)(scalar1x,scalar1x);
587
    API_NS(scalar_add)(scalar2x, scalarc, point_scalarmul_adjustment);
588
    API_NS(scalar_halve)(scalar2x,scalar2x);
589 590 591 592 593 594 595 596 597 598 599 600 601
    
    /* Set up a precomputed table with odd multiples of b. */
    pniels_t pn, multiples1[NTABLE], multiples2[NTABLE];
    point_t tmp;
    prepare_fixed_window(multiples1, b, NTABLE);
    prepare_fixed_window(multiples2, c, NTABLE);

    /* Initialize. */
    int i,j,first=1;
    i = SCALAR_BITS - ((SCALAR_BITS-1) % WINDOW) - 1;

    for (; i>=0; i-=WINDOW) {
        /* Fetch another block of bits */
602
        word_t bits1 = scalar1x->limb[i/WBITS] >> (i%WBITS),
603 604 605 606 607 608 609
                     bits2 = scalar2x->limb[i/WBITS] >> (i%WBITS);
        if (i%WBITS >= WBITS-WINDOW && i/WBITS<SCALAR_LIMBS-1) {
            bits1 ^= scalar1x->limb[i/WBITS+1] << (WBITS - (i%WBITS));
            bits2 ^= scalar2x->limb[i/WBITS+1] << (WBITS - (i%WBITS));
        }
        bits1 &= WINDOW_MASK;
        bits2 &= WINDOW_MASK;
610 611
        mask_t inv1 = (bits1>>(WINDOW-1))-1;
        mask_t inv2 = (bits2>>(WINDOW-1))-1;
612 613 614 615
        bits1 ^= inv1;
        bits2 ^= inv2;
    
        /* Add in from table.  Compute t only on last iteration. */
616
        constant_time_lookup(pn, multiples1, sizeof(pn), NTABLE, bits1 & WINDOW_T_MASK);
617 618 619 620 621 622 623 624 625 626 627 628 629 630
        cond_neg_niels(pn->n, inv1);
        if (first) {
            pniels_to_pt(tmp, pn);
            first = 0;
        } else {
           /* Using Hisil et al's lookahead method instead of extensible here
            * for no particular reason.  Double WINDOW times, but only compute t on
            * the last one.
            */
            for (j=0; j<WINDOW-1; j++)
                point_double_internal(tmp, tmp, -1);
            point_double_internal(tmp, tmp, 0);
            add_pniels_to_pt(tmp, pn, 0);
        }
631
        constant_time_lookup(pn, multiples2, sizeof(pn), NTABLE, bits2 & WINDOW_T_MASK);
632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653
        cond_neg_niels(pn->n, inv2);
        add_pniels_to_pt(tmp, pn, i?-1:0);
    }
    
    /* Write out the answer */
    API_NS(point_copy)(a,tmp);
    

    decaf_bzero(scalar1x,sizeof(scalar1x));
    decaf_bzero(scalar2x,sizeof(scalar2x));
    decaf_bzero(pn,sizeof(pn));
    decaf_bzero(multiples1,sizeof(multiples1));
    decaf_bzero(multiples2,sizeof(multiples2));
    decaf_bzero(tmp,sizeof(tmp));
}

void API_NS(point_dual_scalarmul) (
    point_t a1,
    point_t a2,
    const point_t b,
    const scalar_t scalar1,
    const scalar_t scalar2
654
) {
655 656 657 658 659 660
    const int WINDOW = DECAF_WINDOW_BITS,
        WINDOW_MASK = (1<<WINDOW)-1,
        WINDOW_T_MASK = WINDOW_MASK >> 1,
        NTABLE = 1<<(WINDOW-1);
        
    scalar_t scalar1x, scalar2x;
661
    API_NS(scalar_add)(scalar1x, scalar1, point_scalarmul_adjustment);
662
    API_NS(scalar_halve)(scalar1x,scalar1x);
663
    API_NS(scalar_add)(scalar2x, scalar2, point_scalarmul_adjustment);
664
    API_NS(scalar_halve)(scalar2x,scalar2x);
665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687
    
    /* Set up a precomputed table with odd multiples of b. */
    point_t multiples1[NTABLE], multiples2[NTABLE], working, tmp;
    pniels_t pn;
    
    API_NS(point_copy)(working, b);

    /* Initialize. */
    int i,j;
    
    for (i=0; i<NTABLE; i++) {
        API_NS(point_copy)(multiples1[i], API_NS(point_identity));
        API_NS(point_copy)(multiples2[i], API_NS(point_identity));
    }

    for (i=0; i<SCALAR_BITS; i+=WINDOW) {   
        if (i) {
            for (j=0; j<WINDOW-1; j++)
                point_double_internal(working, working, -1);
            point_double_internal(working, working, 0);
        }
        
        /* Fetch another block of bits */
688 689
        word_t bits1 = scalar1x->limb[i/WBITS] >> (i%WBITS),
               bits2 = scalar2x->limb[i/WBITS] >> (i%WBITS);
690 691 692 693 694 695
        if (i%WBITS >= WBITS-WINDOW && i/WBITS<SCALAR_LIMBS-1) {
            bits1 ^= scalar1x->limb[i/WBITS+1] << (WBITS - (i%WBITS));
            bits2 ^= scalar2x->limb[i/WBITS+1] << (WBITS - (i%WBITS));
        }
        bits1 &= WINDOW_MASK;
        bits2 &= WINDOW_MASK;
696 697
        mask_t inv1 = (bits1>>(WINDOW-1))-1;
        mask_t inv2 = (bits2>>(WINDOW-1))-1;
698 699 700 701 702
        bits1 ^= inv1;
        bits2 ^= inv2;
        
        pt_to_pniels(pn, working);

703
        constant_time_lookup(tmp, multiples1, sizeof(tmp), NTABLE, bits1 & WINDOW_T_MASK);
704 705 706 707 708 709
        cond_neg_niels(pn->n, inv1);
        /* add_pniels_to_pt(multiples1[bits1 & WINDOW_T_MASK], pn, 0); */
        add_pniels_to_pt(tmp, pn, 0);
        constant_time_insert(multiples1, tmp, sizeof(tmp), NTABLE, bits1 & WINDOW_T_MASK);
        
        
710
        constant_time_lookup(tmp, multiples2, sizeof(tmp), NTABLE, bits2 & WINDOW_T_MASK);
711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745
        cond_neg_niels(pn->n, inv1^inv2);
        /* add_pniels_to_pt(multiples2[bits2 & WINDOW_T_MASK], pn, 0); */
        add_pniels_to_pt(tmp, pn, 0);
        constant_time_insert(multiples2, tmp, sizeof(tmp), NTABLE, bits2 & WINDOW_T_MASK);
    }
    
    if (NTABLE > 1) {
        API_NS(point_copy)(working, multiples1[NTABLE-1]);
        API_NS(point_copy)(tmp    , multiples2[NTABLE-1]);
    
        for (i=NTABLE-1; i>1; i--) {
            API_NS(point_add)(multiples1[i-1], multiples1[i-1], multiples1[i]);
            API_NS(point_add)(multiples2[i-1], multiples2[i-1], multiples2[i]);
            API_NS(point_add)(working, working, multiples1[i-1]);
            API_NS(point_add)(tmp,     tmp,     multiples2[i-1]);
        }
    
        API_NS(point_add)(multiples1[0], multiples1[0], multiples1[1]);
        API_NS(point_add)(multiples2[0], multiples2[0], multiples2[1]);
        point_double_internal(working, working, 0);
        point_double_internal(tmp,         tmp, 0);
        API_NS(point_add)(a1, working, multiples1[0]);
        API_NS(point_add)(a2, tmp,     multiples2[0]);
    } else {
        API_NS(point_copy)(a1, multiples1[0]);
        API_NS(point_copy)(a2, multiples2[0]);
    }

    decaf_bzero(scalar1x,sizeof(scalar1x));
    decaf_bzero(scalar2x,sizeof(scalar2x));
    decaf_bzero(pn,sizeof(pn));
    decaf_bzero(multiples1,sizeof(multiples1));
    decaf_bzero(multiples2,sizeof(multiples2));
    decaf_bzero(tmp,sizeof(tmp));
    decaf_bzero(working,sizeof(working));
746 747
}

748
decaf_bool_t API_NS(point_eq) ( const point_t p, const point_t q ) {
749
    /* equality mod 2-torsion compares x/y */
Michael Hamburg's avatar
Michael Hamburg committed
750 751 752
    gf a, b;
    gf_mul ( a, p->y, q->x );
    gf_mul ( b, q->y, p->x );
753
    mask_t succ = gf_eq(a,b);
754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770
    
    #if (COFACTOR == 8) && IMAGINE_TWIST
        gf_mul ( a, p->y, q->y );
        gf_mul ( b, q->x, p->x );
        #if !(IMAGINE_TWIST)
            gf_sub ( a, ZERO, a );
        #else
           /* Interesting note: the 4tor would normally be rotation.
            * But because of the *i twist, it's actually
            * (x,y) <-> (iy,ix)
            */
    
           /* No code, just a comment. */
        #endif
        succ |= gf_eq(a,b);
    #endif
    
771
    return mask_to_bool(succ);
772 773
}

774 775
decaf_bool_t API_NS(point_valid) (
    const point_t p
Michael Hamburg's avatar
Michael Hamburg committed
776 777 778 779
) {
    gf a,b,c;
    gf_mul(a,p->x,p->y);
    gf_mul(b,p->z,p->t);
780
    mask_t out = gf_eq(a,b);
Michael Hamburg's avatar
Michael Hamburg committed
781 782 783 784
    gf_sqr(a,p->x);
    gf_sqr(b,p->y);
    gf_sub(a,b,a);
    gf_sqr(b,p->t);
785
    gf_mulw(c,b,TWISTED_D);
Michael Hamburg's avatar
Michael Hamburg committed
786
    gf_sqr(b,p->z);
787
    gf_add(b,b,c);
Michael Hamburg's avatar
Michael Hamburg committed
788
    out &= gf_eq(a,b);
789
    out &= ~gf_eq(p->z,ZERO);
790
    return mask_to_bool(out);
Michael Hamburg's avatar
Michael Hamburg committed
791
}
792

793 794 795 796
void API_NS(point_debugging_torque) (
    point_t q,
    const point_t p
) {
797
#if COFACTOR == 8 && IMAGINE_TWIST
798 799 800
    gf tmp;
    gf_mul(tmp,p->x,SQRT_MINUS_ONE);
    gf_mul(q->x,p->y,SQRT_MINUS_ONE);
801 802
    gf_copy(q->y,tmp);
    gf_copy(q->z,p->z);
803 804 805 806
    gf_sub(q->t,ZERO,p->t);
#else
    gf_sub(q->x,ZERO,p->x);
    gf_sub(q->y,ZERO,p->y);
807 808
    gf_copy(q->z,p->z);
    gf_copy(q->t,p->t);
809 810 811 812 813 814 815 816 817
#endif
}

void API_NS(point_debugging_pscale) (
    point_t q,
    const point_t p,
    const uint8_t factor[SER_BYTES]
) {
    gf gfac,tmp;
818 819
    /* NB this means you'll never pscale by negative numbers for p521 */
    ignore_result(gf_deserialize(gfac,factor,0));
820
    gf_cond_sel(gfac,gfac,ONE,gf_eq(gfac,ZERO));
821
    gf_mul(tmp,p->x,gfac);
822
    gf_copy(q->x,tmp);
823
    gf_mul(tmp,p->y,gfac);
824
    gf_copy(q->y,tmp);
825
    gf_mul(tmp,p->z,gfac);
826
    gf_copy(q->z,tmp);
827
    gf_mul(tmp,p->t,gfac);
828
    gf_copy(q->t,tmp);
829 830 831 832 833 834
}

static void gf_batch_invert (
    gf *__restrict__ out,
    const gf *in,
    unsigned int n
835
) {
836 837 838
    gf t1;
    assert(n>1);
  
839
    gf_copy(out[1], in[0]);
840 841 842 843 844 845
    int i;
    for (i=1; i<(int) (n-1); i++) {
        gf_mul(out[i+1], out[i], in[i]);
    }
    gf_mul(out[0], out[n-1], in[n-1]);

846
    gf_invert(out[0], out[0], 1);
847 848 849

    for (i=n-1; i>0; i--) {
        gf_mul(t1, out[i], out[0]);
850
        gf_copy(out[i], t1);
851
        gf_mul(t1, out[0], in[i]);
852
        gf_copy(out[0], t1);
853
    }
854 855
}

856 857 858 859 860 861 862 863 864 865 866 867 868
static void batch_normalize_niels (
    niels_t *table,
    const gf *zs,
    gf *__restrict__ zis,
    int n
) {
    int i;
    gf product;
    gf_batch_invert(zis, zs, n);

    for (i=0; i<n; i++) {
        gf_mul(product, table[i]->a, zis[i]);
        gf_strong_reduce(product);
869
        gf_copy(table[i]->a, product);
870 871 872
        
        gf_mul(product, table[i]->b, zis[i]);
        gf_strong_reduce(product);
873
        gf_copy(table[i]->b, product);
874 875 876
        
        gf_mul(product, table[i]->c, zis[i]);
        gf_strong_reduce(product);
877
        gf_copy(table[i]->c, product);
878 879 880 881 882 883 884 885 886
    }
    
    decaf_bzero(product,sizeof(product));
}

void API_NS(precompute) (
    precomputed_s *table,
    const point_t base
) { 
887
    const unsigned int n = COMBS_N, t = COMBS_T, s = COMBS_S;
888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921
    assert(n*t*s >= SCALAR_BITS);
  
    point_t working, start, doubles[t-1];
    API_NS(point_copy)(working, base);
    pniels_t pn_tmp;
  
    gf zs[n<<(t-1)], zis[n<<(t-1)];
  
    unsigned int i,j,k;
    
    /* Compute n tables */
    for (i=0; i<n; i++) {

        /* Doubling phase */
        for (j=0; j<t; j++) {
            if (j) API_NS(point_add)(start, start, working);
            else API_NS(point_copy)(start, working);

            if (j==t-1 && i==n-1) break;

            point_double_internal(working, working,0);
            if (j<t-1) API_NS(point_copy)(doubles[j], working);

            for (k=0; k<s-1; k++)
                point_double_internal(working, working, k<s-2);
        }

        /* Gray-code phase */
        for (j=0;; j++) {
            int gray = j ^ (j>>1);
            int idx = (((i+1)<<(t-1))-1) ^ gray;

            pt_to_pniels(pn_tmp, start);
            memcpy(table->table[idx], pn_tmp->n, sizeof(pn_tmp->n));
922
            gf_copy(zs[idx], pn_tmp->z);
923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947
			
            if (j >= (1u<<(t-1)) - 1) break;
            int delta = (j+1) ^ ((j+1)>>1) ^ gray;

            for (k=0; delta>1; k++)
                delta >>=1;
            
            if (gray & (1<<k)) {
                API_NS(point_add)(start, start, doubles[k]);
            } else {
                API_NS(point_sub)(start, start, doubles[k]);
            }
        }
    }
    
    batch_normalize_niels(table->table,(const gf *)zs,zis,n<<(t-1));
    
    decaf_bzero(zs,sizeof(zs));
    decaf_bzero(zis,sizeof(zis));
    decaf_bzero(pn_tmp,sizeof(pn_tmp));
    decaf_bzero(working,sizeof(working));
    decaf_bzero(start,sizeof(start));
    decaf_bzero(doubles,sizeof(doubles));
}

948
static DECAF_INLINE void
949
constant_time_lookup_niels (
950 951 952 953 954
    niels_s *__restrict__ ni,
    const niels_t *table,
    int nelts,
    int idx
) {
955
    constant_time_lookup(ni, table, sizeof(niels_s), nelts, idx);
956 957 958 959 960 961 962 963 964
}

void API_NS(precomputed_scalarmul) (
    point_t out,
    const precomputed_s *table,
    const scalar_t scalar
) {
    int i;
    unsigned j,k;
965
    const unsigned int n = COMBS_N, t = COMBS_T, s = COMBS_S;
966 967
    
    scalar_t scalar1x;
968
    API_NS(scalar_add)(scalar1x, scalar, precomputed_scalarmul_adjustment);
969
    API_NS(scalar_halve)(scalar1x,scalar1x);
970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985
    
    niels_t ni;
    
    for (i=s-1; i>=0; i--) {
        if (i != (int)s-1) point_double_internal(out,out,0);
        
        for (j=0; j<n; j++) {
            int tab = 0;
         
            for (k=0; k<t; k++) {
                unsigned int bit = i + s*(k + j*t);
                if (bit < SCALAR_BITS) {
                    tab |= (scalar1x->limb[bit/WBITS] >> (bit%WBITS) & 1) << k;
                }
            }
            
986
            mask_t invert = (tab>>(t-1))-1;
987 988 989
            tab ^= invert;
            tab &= (1<<(t-1)) - 1;

990
            constant_time_lookup_niels(ni, &table->table[j<<(t-1)], 1<<(t-1), tab);
991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010

            cond_neg_niels(ni, invert);
            if ((i!=(int)s-1)||j) {
                add_niels_to_pt(out, ni, j==n-1 && i);
            } else {
                niels_to_pt(out, ni);
            }
        }
    }
    
    decaf_bzero(ni,sizeof(ni));
    decaf_bzero(scalar1x,sizeof(scalar1x));
}

void API_NS(point_cond_sel) (
    point_t out,
    const point_t a,
    const point_t b,
    decaf_bool_t pick_b
) {
1011
    constant_time_select(out,a,b,sizeof(point_t),bool_to_mask(pick_b),0);
1012 1013 1014 1015 1016 1017 1018
}

/* FUTURE: restore Curve25519 Montgomery ladder? */
decaf_error_t API_NS(direct_scalarmul) (
    uint8_t scaled[SER_BYTES],
    const uint8_t base[SER_BYTES],
    const scalar_t scalar,
1019 1020 1021
    decaf_bool_t allow_identity,
    decaf_bool_t short_circuit
) {
1022
    point_t basep;
1023 1024
    decaf_error_t succ = API_NS(point_decode)(basep, base, allow_identity);
    if (short_circuit && succ != DECAF_SUCCESS) return succ;
1025 1026 1027 1028
    API_NS(point_cond_sel)(basep, API_NS(point_base), basep, succ);
    API_NS(point_scalarmul)(basep, basep, scalar);
    API_NS(point_encode)(scaled, basep);
    API_NS(point_destroy)(basep);
1029
    return succ;
1030 1031
}

1032
void API_NS(point_mul_by_cofactor_and_encode_like_eddsa) (
1033
    uint8_t enc[DECAF_EDDSA_$(gf_shortname)_PUBLIC_BYTES],
1034 1035 1036 1037 1038
    const point_t p
) {
    
    /* The point is now on the twisted curve.  Move it to untwisted. */
    gf x, y, z, t;
1039 1040 1041 1042 1043 1044
    point_t q;
#if COFACTOR == 8
    API_NS(point_double)(q,p);
#else
    API_NS(point_copy)(q,p);
#endif
Michael Hamburg's avatar
Michael Hamburg committed
1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060
    
#if EDDSA_USE_SIGMA_ISOGENY
    {
        /* Use 4-isogeny like ed25519:
         *   2*x*y*sqrt(d/a-1)/(ax^2 + y^2 - 2)
         *   (y^2 - ax^2)/(y^2 + ax^2)
         * with a = -1, d = -EDWARDS_D:
         *   -2xysqrt(EDWARDS_D-1)/(2z^2-y^2+x^2)
         *   (y^2+x^2)/(y^2-x^2)
         */
        gf u;
        gf_sqr ( x, q->x ); // x^2
        gf_sqr ( t, q->y ); // y^2
        gf_add( u, x, t ); // x^2 + y^2
        gf_add( z, q->y, q->x );
        gf_sqr ( y, z);
1061
        gf_sub ( y, u, y ); // -2xy
Michael Hamburg's avatar
Michael Hamburg committed
1062 1063 1064 1065 1066 1067 1068 1069
        gf_sub ( z, t, x ); // y^2 - x^2
        gf_sqr ( x, q->z );
        gf_add ( t, x, x);
        gf_sub ( t, t, z);  // 2z^2 - y^2 + x^2
        gf_mul ( x, y, z ); // 2xy(y^2-x^2)
        gf_mul ( y, u, t ); // (x^2+y^2)(2z^2-y^2+x^2)
        gf_mul ( u, z, t );
        gf_copy( z, u );
1070
        gf_mul ( u, x, RISTRETTO_FACTOR );
1071
#if IMAGINE_TWIST
1072
        gf_mul_i( x, u );
1073 1074
#else
#error "... probably wrong"
Michael Hamburg's avatar
Michael Hamburg committed
1075
        gf_copy( x, u );
1076
#endif
Michael Hamburg's avatar
Michael Hamburg committed
1077 1078 1079
        decaf_bzero(u,sizeof(u));
    }
#elif IMAGINE_TWIST
1080
    {
1081
        API_NS(point_double)(q,q);
1082
        API_NS(point_double)(q,q);
1083
        gf_mul_i(x, q->x);
1084 1085 1086 1087 1088 1089 1090
        gf_copy(y, q->y);
        gf_copy(z, q->z);
    }
#else
    {
        /* 4-isogeny: 2xy/(y^+x^2), (y^2-x^2)/(2z^2-y^2+x^2) */
        gf u;
1091 1092
        gf_sqr ( x, q->x );
        gf_sqr ( t, q->y );
1093
        gf_add( u, x, t );
1094
        gf_add( z, q->y, q->x );
1095
        gf_sqr ( y, z);
1096
        gf_sub ( y, y, u );
1097
        gf_sub ( z, t, x );
1098
        gf_sqr ( x, q->z );
1099 1100 1101 1102 1103
        gf_add ( t, x, x); 
        gf_sub ( t, t, z);
        gf_mul ( x, t, y );
        gf_mul ( y, z, u );
        gf_mul ( z, u, t );
Michael Hamburg's avatar
Michael Hamburg committed
1104
        decaf_bzero(u,sizeof(u));
1105 1106 1107
    }
#endif
    /* Affinize */
1108
    gf_invert(z,z,1);
1109 1110 1111 1112
    gf_mul(t,x,z);
    gf_mul(x,y,z);
    
    /* Encode */
1113
    enc[DECAF_EDDSA_$(gf_shortname)_PRIVATE_BYTES-1] = 0;
1114
    gf_serialize(enc, x, 1);
1115
    enc[DECAF_EDDSA_$(gf_shortname)_PRIVATE_BYTES-1] |= 0x80 & gf_lobit(t);
1116 1117 1118 1119 1120

    decaf_bzero(x,sizeof(x));
    decaf_bzero(y,sizeof(y));
    decaf_bzero(z,sizeof(z));
    decaf_bzero(t,sizeof(t));
1121
    API_NS(point_destroy)(q);
1122 1123
}

1124

1125
decaf_error_t API_NS(point_decode_like_eddsa_and_ignore_cofactor) (
1126
    point_t p,
1127
    const uint8_t enc[DECAF_EDDSA_$(gf_shortname)_PUBLIC_BYTES]
1128
) {
1129
    uint8_t enc2[DECAF_EDDSA_$(gf_shortname)_PUBLIC_BYTES];
1130 1131
    memcpy(enc2,enc,sizeof(enc2));

1132 1133
    mask_t low = ~word_is_zero(enc2[DECAF_EDDSA_$(gf_shortname)_PRIVATE_BYTES-1] & 0x80);
    enc2[DECAF_EDDSA_$(gf_shortname)_PRIVATE_BYTES-1] &= ~0x80;
1134
    
1135
    mask_t succ = gf_deserialize(p->y, enc2, 1);
1136
#if $(gf_bits % 8) == 0
1137
    succ &= word_is_zero(enc2[DECAF_EDDSA_$(gf_shortname)_PRIVATE_BYTES-1]);
1138
#endif
Michael Hamburg's avatar
Michael Hamburg committed
1139 1140 1141

    gf_sqr(p->x,p->y);
    gf_sub(p->z,ONE,p->x); /* num = 1-y^2 */
1142 1143 1144 1145 1146 1147 1148
    #if EDDSA_USE_SIGMA_ISOGENY
        gf_mulw(p->t,p->z,EDWARDS_D); /* d-dy^2 */
        gf_mulw(p->x,p->z,EDWARDS_D-1); /* num = (1-y^2)(d-1) */
        gf_copy(p->z,p->x);
    #else
        gf_mulw(p->t,p->x,EDWARDS_D); /* dy^2 */
    #endif
Michael Hamburg's avatar
Michael Hamburg committed
1149 1150
    gf_sub(p->t,ONE,p->t); /* denom = 1-dy^2 or 1-d + dy^2 */
    
1151
    gf_mul(p->x,p->z,p->t);
Michael Hamburg's avatar
Michael Hamburg committed
1152 1153 1154
    succ &= gf_isr(p->t,p->x); /* 1/sqrt(num * denom) */
    
    gf_mul(p->x,p->t,p->z); /* sqrt(num / denom) */
1155
    gf_cond_neg(p->x,gf_lobit(p->x)^low);
1156
    gf_copy(p->z,ONE);
Michael Hamburg's avatar
Michael Hamburg committed
1157
  
1158
    #if EDDSA_USE_SIGMA_ISOGENY
Michael Hamburg's avatar
Michael Hamburg committed
1159 1160
    {
       /* Use 4-isogeny like ed25519:
1161
        *   2*x*y/sqrt(1-d/a)/(ax^2 + y^2 - 2)
Michael Hamburg's avatar
Michael Hamburg committed
1162
        *   (y^2 - ax^2)/(y^2 + ax^2)
1163 1164 1165
        * (MAGIC: above formula may be off by a factor of -a
        * or something somewhere; check it for other a)
        *
Michael Hamburg's avatar
Michael Hamburg committed
1166
        * with a = -1, d = -EDWARDS_D:
1167
        *   -2xy/sqrt(1-EDWARDS_D)/(2z^2-y^2+x^2)
Michael Hamburg's avatar
Michael Hamburg committed
1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179
        *   (y^2+x^2)/(y^2-x^2)
        */
        gf a, b, c, d;
        gf_sqr ( c, p->x );
        gf_sqr ( a, p->y );
        gf_add ( d, c, a ); // x^2 + y^2
        gf_add ( p->t, p->y, p->x );
        gf_sqr ( b, p->t );
        gf_sub ( b, b, d ); // 2xy
        gf_sub ( p->t, a, c ); // y^2 - x^2
        gf_sqr ( p->x, p->z );
        gf_add ( p->z, p->x, p->x );
1180
        gf_sub ( c, p->z, p->t ); // 2z^2 - y^2 + x^2
1181 1182
        gf_div_i ( a, c );
        gf_mul ( c, a, RISTRETTO_FACTOR );
Michael Hamburg's avatar
Michael Hamburg committed
1183 1184 1185 1186 1187 1188 1189 1190 1191
        gf_mul ( p->x, b, p->t); // (2xy)(y^2-x^2)
        gf_mul ( p->z, p->t, c ); // (y^2-x^2)sd(2z^2 - y^2 + x^2)
        gf_mul ( p->y, d, c ); // (y^2+x^2)sd(2z^2 - y^2 + x^2)
        gf_mul ( p->t, d, b );
        decaf_bzero(a,sizeof(a));
        decaf_bzero(b,sizeof(b));
        decaf_bzero(c,sizeof(c));
        decaf_bzero(d,sizeof(d));
    } 
1192
    #elif IMAGINE_TWIST
1193 1194 1195
    {
        gf_mul(p->t,p->x,SQRT_MINUS_ONE);
        gf_copy(p->x,p->t);
1196
        gf_mul(p->t,p->x,p->y);
1197
    }
1198
    #else
1199
    {
1200
        /* 4-isogeny 2xy/(y^2-ax^2), (y^2+ax^2)/(2-y^2-ax^2) */
1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220
        gf a, b, c, d;
        gf_sqr ( c, p->x );
        gf_sqr ( a, p->y );
        gf_add ( d, c, a );
        gf_add ( p->t, p->y, p->x );
        gf_sqr ( b, p->t );
        gf_sub ( b, b, d );
        gf_sub ( p->t, a, c );
        gf_sqr ( p->x, p->z );
        gf_add ( p->z, p->x, p->x );
        gf_sub ( a, p->z, d );
        gf_mul ( p->x, a, b );
        gf_mul ( p->z, p->t, a );
        gf_mul ( p->y, p->t, d );
        gf_mul ( p->t, b, d );
        decaf_bzero(a,sizeof(a));
        decaf_bzero(b,sizeof(b));
        decaf_bzero(c,sizeof(c));
        decaf_bzero(d,sizeof(d));
    }
1221
    #endif
1222 1223 1224
    
    decaf_bzero(enc2,sizeof(enc2));
    assert(API_NS(point_valid)(p) || ~succ);
1225
    
1226
    return decaf_succeed_if(mask_to_bool(succ));
1227 1228
}

1229
decaf_error_t decaf_x$(gf_shortname) (
1230 1231 1232 1233 1234
    uint8_t out[X_PUBLIC_BYTES],
    const uint8_t base[X_PUBLIC_BYTES],
    const uint8_t scalar[X_PRIVATE_BYTES]
) {
    gf x1, x2, z2, x3, z3, t1, t2;
1235
    ignore_result(gf_deserialize(x1,base,1));
1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254
    gf_copy(x2,ONE);
    gf_copy(z2,ZERO);
    gf_copy(x3,x1);
    gf_copy(z3,ONE);
    
    int t;
    mask_t swap = 0;
    
    for (t = X_PRIVATE_BITS-1; t>=0; t--) {
        uint8_t sb = scalar[t/8];
        
        /* Scalar conditioning */
        if (t/8==0) sb &= -(uint8_t)COFACTOR;
        else if (t == X_PRIVATE_BITS-1) sb = -1;
        
        mask_t k_t = (sb>>(t%8)) & 1;
        k_t = -k_t; /* set to all 0s or all 1s */
        
        swap ^= k_t;
1255 1256
        gf_cond_swap(x2,x3,swap);
        gf_cond_swap(z2,z3,swap);
1257 1258
        swap = k_t;
        
1259 1260 1261
        gf_add_nr(t1,x2,z2); /* A = x2 + z2 */        /* 2+e */
        gf_sub_nr(t2,x2,z2); /* B = x2 - z2 */        /* 3+e */
        gf_sub_nr(z2,x3,z3); /* D = x3 - z3 */        /* 3+e */
1262
        gf_mul(x2,t1,z2);    /* DA */
1263
        gf_add_nr(z2,z3,x3); /* C = x3 + z3 */        /* 2+e */
1264
        gf_mul(x3,t2,z2);    /* CB */
1265
        gf_sub_nr(z3,x2,x3); /* DA-CB */              /* 3+e */
1266 1267
        gf_sqr(z2,z3);       /* (DA-CB)^2 */
        gf_mul(z3,x1,z2);    /* z3 = x1(DA-CB)^2 */
1268
        gf_add_nr(z2,x2,x3); /* (DA+CB) */            /* 2+e */
1269
        gf_sqr(x3,z2);       /* x3 = (DA+CB)^2 */
1270
        
1271 1272 1273
        gf_sqr(z2,t1);       /* AA = A^2 */
        gf_sqr(t1,t2);       /* BB = B^2 */
        gf_mul(x2,z2,t1);    /* x2 = AA*BB */
1274
        gf_sub_nr(t2,z2,t1); /* E = AA-BB */          /* 3+e */
1275