Field elements are written as an array of signed, 64-bit limbs (an array of longs), least + * significant first. The value of the field element is: + * + *
+ * x[0] + 2^26·x[1] + 2^51·x[2] + 2^77·x[3] + 2^102·x[4] + 2^128·x[5] + 2^153·x[6] + 2^179·x[7] + + * 2^204·x[8] + 2^230·x[9], + *+ * + *
i.e. the limbs are 26, 25, 26, 25, ... bits wide. + */ +@Alpha +final class Field25519 { + /** + * During Field25519 computation, the mixed radix representation may be in different forms: + *
+ * TODO(quannguyen): + *
+ * On entry: in1, in2 are in reduced-size form. + */ + static void sum(long[] output, long[] in1, long[] in2) { + for (int i = 0; i < LIMB_CNT; i++) { + output[i] = in1[i] + in2[i]; + } + } + + /** + * Sums two numbers: output += in + *
+ * On entry: in is in reduced-size form. + */ + static void sum(long[] output, long[] in) { + sum(output, output, in); + } + + /** + * Find the difference of two numbers: output = in1 - in2 + * (note the order of the arguments!). + *
+ * On entry: in1, in2 are in reduced-size form. + */ + static void sub(long[] output, long[] in1, long[] in2) { + for (int i = 0; i < LIMB_CNT; i++) { + output[i] = in1[i] - in2[i]; + } + } + + /** + * Find the difference of two numbers: output = in - output + * (note the order of the arguments!). + *
+ * On entry: in, output are in reduced-size form. + */ + static void sub(long[] output, long[] in) { + sub(output, in, output); + } + + /** + * Multiply a number by a scalar: output = in * scalar + */ + static void scalarProduct(long[] output, long[] in, long scalar) { + for (int i = 0; i < LIMB_CNT; i++) { + output[i] = in[i] * scalar; + } + } + + /** + * Multiply two numbers: out = in2 * in + *
+ * output must be distinct to both inputs. The inputs are reduced coefficient form, + * the output is not. + *
+ * out[x] <= 14 * the largest product of the input limbs. + */ + static void product(long[] out, long[] in2, long[] in) { + out[0] = in2[0] * in[0]; + out[1] = in2[0] * in[1] + + in2[1] * in[0]; + out[2] = 2 * in2[1] * in[1] + + in2[0] * in[2] + + in2[2] * in[0]; + out[3] = in2[1] * in[2] + + in2[2] * in[1] + + in2[0] * in[3] + + in2[3] * in[0]; + out[4] = in2[2] * in[2] + + 2 * (in2[1] * in[3] + in2[3] * in[1]) + + in2[0] * in[4] + + in2[4] * in[0]; + out[5] = in2[2] * in[3] + + in2[3] * in[2] + + in2[1] * in[4] + + in2[4] * in[1] + + in2[0] * in[5] + + in2[5] * in[0]; + out[6] = 2 * (in2[3] * in[3] + in2[1] * in[5] + in2[5] * in[1]) + + in2[2] * in[4] + + in2[4] * in[2] + + in2[0] * in[6] + + in2[6] * in[0]; + out[7] = in2[3] * in[4] + + in2[4] * in[3] + + in2[2] * in[5] + + in2[5] * in[2] + + in2[1] * in[6] + + in2[6] * in[1] + + in2[0] * in[7] + + in2[7] * in[0]; + out[8] = in2[4] * in[4] + + 2 * (in2[3] * in[5] + in2[5] * in[3] + in2[1] * in[7] + in2[7] * in[1]) + + in2[2] * in[6] + + in2[6] * in[2] + + in2[0] * in[8] + + in2[8] * in[0]; + out[9] = in2[4] * in[5] + + in2[5] * in[4] + + in2[3] * in[6] + + in2[6] * in[3] + + in2[2] * in[7] + + in2[7] * in[2] + + in2[1] * in[8] + + in2[8] * in[1] + + in2[0] * in[9] + + in2[9] * in[0]; + out[10] = + 2 * (in2[5] * in[5] + in2[3] * in[7] + in2[7] * in[3] + in2[1] * in[9] + in2[9] * in[1]) + + in2[4] * in[6] + + in2[6] * in[4] + + in2[2] * in[8] + + in2[8] * in[2]; + out[11] = in2[5] * in[6] + + in2[6] * in[5] + + in2[4] * in[7] + + in2[7] * in[4] + + in2[3] * in[8] + + in2[8] * in[3] + + in2[2] * in[9] + + in2[9] * in[2]; + out[12] = in2[6] * in[6] + + 2 * (in2[5] * in[7] + in2[7] * in[5] + in2[3] * in[9] + in2[9] * in[3]) + + in2[4] * in[8] + + in2[8] * in[4]; + out[13] = in2[6] * in[7] + + in2[7] * in[6] + + in2[5] * in[8] + + in2[8] * in[5] + + in2[4] * in[9] + + in2[9] * in[4]; + out[14] = 2 * (in2[7] * in[7] + in2[5] * in[9] + in2[9] * in[5]) + + in2[6] * in[8] + + in2[8] * in[6]; + out[15] = in2[7] * in[8] + + in2[8] * in[7] + + in2[6] * in[9] + + in2[9] * in[6]; + out[16] = in2[8] * in[8] + + 2 * (in2[7] * in[9] + in2[9] * in[7]); + out[17] = in2[8] * in[9] + + in2[9] * in[8]; + out[18] = 2 * in2[9] * in[9]; + } + + /** + * Reduce a long form to a reduced-size form by taking the input mod 2^255 - 19. + *
+ * On entry: |output[i]| < 14*2^54 + * On exit: |output[0..8]| < 280*2^54 + */ + static void reduceSizeByModularReduction(long[] output) { + // The coefficients x[10], x[11],..., x[18] are eliminated by reduction modulo 2^255 - 19. + // For example, the coefficient x[18] is multiplied by 19 and added to the coefficient x[8]. + // + // Each of these shifts and adds ends up multiplying the value by 19. + // + // For output[0..8], the absolute entry value is < 14*2^54 and we add, at most, 19*14*2^54 thus, + // on exit, |output[0..8]| < 280*2^54. + output[8] += output[18] << 4; + output[8] += output[18] << 1; + output[8] += output[18]; + output[7] += output[17] << 4; + output[7] += output[17] << 1; + output[7] += output[17]; + output[6] += output[16] << 4; + output[6] += output[16] << 1; + output[6] += output[16]; + output[5] += output[15] << 4; + output[5] += output[15] << 1; + output[5] += output[15]; + output[4] += output[14] << 4; + output[4] += output[14] << 1; + output[4] += output[14]; + output[3] += output[13] << 4; + output[3] += output[13] << 1; + output[3] += output[13]; + output[2] += output[12] << 4; + output[2] += output[12] << 1; + output[2] += output[12]; + output[1] += output[11] << 4; + output[1] += output[11] << 1; + output[1] += output[11]; + output[0] += output[10] << 4; + output[0] += output[10] << 1; + output[0] += output[10]; + } + + /** + * Reduce all coefficients of the short form input so that |x| < 2^26. + *
+ * On entry: |output[i]| < 280*2^54 + */ + static void reduceCoefficients(long[] output) { + output[10] = 0; + + for (int i = 0; i < LIMB_CNT; i += 2) { + long over = output[i] / TWO_TO_26; + // The entry condition (that |output[i]| < 280*2^54) means that over is, at most, 280*2^28 in + // the first iteration of this loop. This is added to the next limb and we can approximate the + // resulting bound of that limb by 281*2^54. + output[i] -= over << 26; + output[i + 1] += over; + + // For the first iteration, |output[i+1]| < 281*2^54, thus |over| < 281*2^29. When this is + // added to the next limb, the resulting bound can be approximated as 281*2^54. + // + // For subsequent iterations of the loop, 281*2^54 remains a conservative bound and no + // overflow occurs. + over = output[i + 1] / TWO_TO_25; + output[i + 1] -= over << 25; + output[i + 2] += over; + } + // Now |output[10]| < 281*2^29 and all other coefficients are reduced. + output[0] += output[10] << 4; + output[0] += output[10] << 1; + output[0] += output[10]; + + output[10] = 0; + // Now output[1..9] are reduced, and |output[0]| < 2^26 + 19*281*2^29 so |over| will be no more + // than 2^16. + long over = output[0] / TWO_TO_26; + output[0] -= over << 26; + output[1] += over; + // Now output[0,2..9] are reduced, and |output[1]| < 2^25 + 2^16 < 2^26. The bound on + // |output[1]| is sufficient to meet our needs. + } + + /** + * A helpful wrapper around {@ref Field25519#product}: output = in * in2. + *
+ * On entry: |in[i]| < 2^27 and |in2[i]| < 2^27. + *
+ * The output is reduced degree (indeed, one need only provide storage for 10 limbs) and + * |output[i]| < 2^26. + */ + static void mult(long[] output, long[] in, long[] in2) { + long[] t = new long[19]; + product(t, in, in2); + // |t[i]| < 14*2^54 + reduceSizeByModularReduction(t); + reduceCoefficients(t); + // |t[i]| < 2^26 + System.arraycopy(t, 0, output, 0, LIMB_CNT); + } + + /** + * Square a number: out = in**2 + *
+ * output must be distinct from the input. The inputs are reduced coefficient form, the output is + * not. + *
+ * out[x] <= 14 * the largest product of the input limbs. + */ + private static void squareInner(long[] out, long[] in) { + out[0] = in[0] * in[0]; + out[1] = 2 * in[0] * in[1]; + out[2] = 2 * (in[1] * in[1] + in[0] * in[2]); + out[3] = 2 * (in[1] * in[2] + in[0] * in[3]); + out[4] = in[2] * in[2] + + 4 * in[1] * in[3] + + 2 * in[0] * in[4]; + out[5] = 2 * (in[2] * in[3] + in[1] * in[4] + in[0] * in[5]); + out[6] = 2 * (in[3] * in[3] + in[2] * in[4] + in[0] * in[6] + 2 * in[1] * in[5]); + out[7] = 2 * (in[3] * in[4] + in[2] * in[5] + in[1] * in[6] + in[0] * in[7]); + out[8] = in[4] * in[4] + + 2 * (in[2] * in[6] + in[0] * in[8] + 2 * (in[1] * in[7] + in[3] * in[5])); + out[9] = 2 * (in[4] * in[5] + in[3] * in[6] + in[2] * in[7] + in[1] * in[8] + in[0] * in[9]); + out[10] = 2 * (in[5] * in[5] + + in[4] * in[6] + + in[2] * in[8] + + 2 * (in[3] * in[7] + in[1] * in[9])); + out[11] = 2 * (in[5] * in[6] + in[4] * in[7] + in[3] * in[8] + in[2] * in[9]); + out[12] = in[6] * in[6] + + 2 * (in[4] * in[8] + 2 * (in[5] * in[7] + in[3] * in[9])); + out[13] = 2 * (in[6] * in[7] + in[5] * in[8] + in[4] * in[9]); + out[14] = 2 * (in[7] * in[7] + in[6] * in[8] + 2 * in[5] * in[9]); + out[15] = 2 * (in[7] * in[8] + in[6] * in[9]); + out[16] = in[8] * in[8] + 4 * in[7] * in[9]; + out[17] = 2 * in[8] * in[9]; + out[18] = 2 * in[9] * in[9]; + } + + /** + * Returns in^2. + *
+ * On entry: The |in| argument is in reduced coefficients form and |in[i]| < 2^27. + *
+ * On exit: The |output| argument is in reduced coefficients form (indeed, one need only provide + * storage for 10 limbs) and |out[i]| < 2^26. + */ + static void square(long[] output, long[] in) { + long[] t = new long[19]; + squareInner(t, in); + // |t[i]| < 14*2^54 because the largest product of two limbs will be < 2^(27+27) and SquareInner + // adds together, at most, 14 of those products. + reduceSizeByModularReduction(t); + reduceCoefficients(t); + // |t[i]| < 2^26 + System.arraycopy(t, 0, output, 0, LIMB_CNT); + } + + /** + * Takes a little-endian, 32-byte number and expands it into mixed radix form. + */ + static long[] expand(byte[] input) { + long[] output = new long[LIMB_CNT]; + for (int i = 0; i < LIMB_CNT; i++) { + output[i] = ((((long) (input[EXPAND_START[i]] & 0xff)) + | ((long) (input[EXPAND_START[i] + 1] & 0xff)) << 8 + | ((long) (input[EXPAND_START[i] + 2] & 0xff)) << 16 + | ((long) (input[EXPAND_START[i] + 3] & 0xff)) << 24) >> EXPAND_SHIFT[i]) & MASK[i & 1]; + } + return output; + } + + /** + * Takes a fully reduced mixed radix form number and contract it into a little-endian, 32-byte + * array. + *
+ * On entry: |input_limbs[i]| < 2^26 + */ + @SuppressWarnings("NarrowingCompoundAssignment") + static byte[] contract(long[] inputLimbs) { + long[] input = Arrays.copyOf(inputLimbs, LIMB_CNT); + for (int j = 0; j < 2; j++) { + for (int i = 0; i < 9; i++) { + // This calculation is a time-invariant way to make input[i] non-negative by borrowing + // from the next-larger limb. + int carry = -(int) ((input[i] & (input[i] >> 31)) >> SHIFT[i & 1]); + input[i] = input[i] + (carry << SHIFT[i & 1]); + input[i + 1] -= carry; + } + + // There's no greater limb for input[9] to borrow from, but we can multiply by 19 and borrow + // from input[0], which is valid mod 2^255-19. + { + int carry = -(int) ((input[9] & (input[9] >> 31)) >> 25); + input[9] += (carry << 25); + input[0] -= (carry * 19); + } + + // After the first iteration, input[1..9] are non-negative and fit within 25 or 26 bits, + // depending on position. However, input[0] may be negative. + } + + // The first borrow-propagation pass above ended with every limb except (possibly) input[0] + // non-negative. + // + // If input[0] was negative after the first pass, then it was because of a carry from input[9]. + // On entry, input[9] < 2^26 so the carry was, at most, one, since (2**26-1) >> 25 = 1. Thus + // input[0] >= -19. + // + // In the second pass, each limb is decreased by at most one. Thus the second borrow-propagation + // pass could only have wrapped around to decrease input[0] again if the first pass left + // input[0] negative *and* input[1] through input[9] were all zero. In that case, input[1] is + // now 2^25 - 1, and this last borrow-propagation step will leave input[1] non-negative. + { + int carry = -(int) ((input[0] & (input[0] >> 31)) >> 26); + input[0] += (carry << 26); + input[1] -= carry; + } + + // All input[i] are now non-negative. However, there might be values between 2^25 and 2^26 in a + // limb which is, nominally, 25 bits wide. + for (int j = 0; j < 2; j++) { + for (int i = 0; i < 9; i++) { + int carry = (int) (input[i] >> SHIFT[i & 1]); + input[i] &= MASK[i & 1]; + input[i + 1] += carry; + } + } + + { + int carry = (int) (input[9] >> 25); + input[9] &= 0x1ffffff; + input[0] += 19 * carry; + } + + // If the first carry-chain pass, just above, ended up with a carry from input[9], and that + // caused input[0] to be out-of-bounds, then input[0] was < 2^26 + 2*19, because the carry was, + // at most, two. + // + // If the second pass carried from input[9] again then input[0] is < 2*19 and the input[9] -> + // input[0] carry didn't push input[0] out of bounds. + + // It still remains the case that input might be between 2^255-19 and 2^255. In this case, + // input[1..9] must take their maximum value and input[0] must be >= (2^255-19) & 0x3ffffff, + // which is 0x3ffffed. + int mask = gte((int) input[0], 0x3ffffed); + for (int i = 1; i < LIMB_CNT; i++) { + mask &= eq((int) input[i], MASK[i & 1]); + } + + // mask is either 0xffffffff (if input >= 2^255-19) and zero otherwise. Thus this conditionally + // subtracts 2^255-19. + input[0] -= mask & 0x3ffffed; + input[1] -= mask & 0x1ffffff; + for (int i = 2; i < LIMB_CNT; i += 2) { + input[i] -= mask & 0x3ffffff; + input[i + 1] -= mask & 0x1ffffff; + } + + for (int i = 0; i < LIMB_CNT; i++) { + input[i] <<= EXPAND_SHIFT[i]; + } + byte[] output = new byte[FIELD_LEN]; + for (int i = 0; i < LIMB_CNT; i++) { + output[EXPAND_START[i]] |= input[i] & 0xff; + output[EXPAND_START[i] + 1] |= (input[i] >> 8) & 0xff; + output[EXPAND_START[i] + 2] |= (input[i] >> 16) & 0xff; + output[EXPAND_START[i] + 3] |= (input[i] >> 24) & 0xff; + } + return output; + } + + /** + * Computes inverse of z = z(2^255 - 21) + *
+ * Shamelessly copied from agl's code which was shamelessly copied from djb's code. Only the + * comment format and the variable namings are different from those. + */ + static void inverse(long[] out, long[] z) { + long[] z2 = new long[Field25519.LIMB_CNT]; + long[] z9 = new long[Field25519.LIMB_CNT]; + long[] z11 = new long[Field25519.LIMB_CNT]; + long[] z2To5Minus1 = new long[Field25519.LIMB_CNT]; + long[] z2To10Minus1 = new long[Field25519.LIMB_CNT]; + long[] z2To20Minus1 = new long[Field25519.LIMB_CNT]; + long[] z2To50Minus1 = new long[Field25519.LIMB_CNT]; + long[] z2To100Minus1 = new long[Field25519.LIMB_CNT]; + long[] t0 = new long[Field25519.LIMB_CNT]; + long[] t1 = new long[Field25519.LIMB_CNT]; + + square(z2, z); // 2 + square(t1, z2); // 4 + square(t0, t1); // 8 + mult(z9, t0, z); // 9 + mult(z11, z9, z2); // 11 + square(t0, z11); // 22 + mult(z2To5Minus1, t0, z9); // 2^5 - 2^0 = 31 + + square(t0, z2To5Minus1); // 2^6 - 2^1 + square(t1, t0); // 2^7 - 2^2 + square(t0, t1); // 2^8 - 2^3 + square(t1, t0); // 2^9 - 2^4 + square(t0, t1); // 2^10 - 2^5 + mult(z2To10Minus1, t0, z2To5Minus1); // 2^10 - 2^0 + + square(t0, z2To10Minus1); // 2^11 - 2^1 + square(t1, t0); // 2^12 - 2^2 + for (int i = 2; i < 10; i += 2) { // 2^20 - 2^10 + square(t0, t1); + square(t1, t0); + } + mult(z2To20Minus1, t1, z2To10Minus1); // 2^20 - 2^0 + + square(t0, z2To20Minus1); // 2^21 - 2^1 + square(t1, t0); // 2^22 - 2^2 + for (int i = 2; i < 20; i += 2) { // 2^40 - 2^20 + square(t0, t1); + square(t1, t0); + } + mult(t0, t1, z2To20Minus1); // 2^40 - 2^0 + + square(t1, t0); // 2^41 - 2^1 + square(t0, t1); // 2^42 - 2^2 + for (int i = 2; i < 10; i += 2) { // 2^50 - 2^10 + square(t1, t0); + square(t0, t1); + } + mult(z2To50Minus1, t0, z2To10Minus1); // 2^50 - 2^0 + + square(t0, z2To50Minus1); // 2^51 - 2^1 + square(t1, t0); // 2^52 - 2^2 + for (int i = 2; i < 50; i += 2) { // 2^100 - 2^50 + square(t0, t1); + square(t1, t0); + } + mult(z2To100Minus1, t1, z2To50Minus1); // 2^100 - 2^0 + + square(t1, z2To100Minus1); // 2^101 - 2^1 + square(t0, t1); // 2^102 - 2^2 + for (int i = 2; i < 100; i += 2) { // 2^200 - 2^100 + square(t1, t0); + square(t0, t1); + } + mult(t1, t0, z2To100Minus1); // 2^200 - 2^0 + + square(t0, t1); // 2^201 - 2^1 + square(t1, t0); // 2^202 - 2^2 + for (int i = 2; i < 50; i += 2) { // 2^250 - 2^50 + square(t0, t1); + square(t1, t0); + } + mult(t0, t1, z2To50Minus1); // 2^250 - 2^0 + + square(t1, t0); // 2^251 - 2^1 + square(t0, t1); // 2^252 - 2^2 + square(t1, t0); // 2^253 - 2^3 + square(t0, t1); // 2^254 - 2^4 + square(t1, t0); // 2^255 - 2^5 + mult(out, t1, z11); // 2^255 - 21 + } + + + /** + * Returns 0xffffffff iff a == b and zero otherwise. + */ + private static int eq(int a, int b) { + a = ~(a ^ b); + a &= a << 16; + a &= a << 8; + a &= a << 4; + a &= a << 2; + a &= a << 1; + return a >> 31; + } + + /** + * returns 0xffffffff if a >= b and zero otherwise, where a and b are both non-negative. + */ + private static int gte(int a, int b) { + a -= b; + // a >= 0 iff a >= b. + return ~(a >> 31); + } +}