tx-signer/src/main/java/com/google/crypto/tink/subtle/Ed25519.java

   1 // Copyright 2017 Google Inc.\r
   2 //\r
   3 // Licensed under the Apache License, Version 2.0 (the "License");\r
   4 // you may not use this file except in compliance with the License.\r
   5 // You may obtain a copy of the License at\r
   6 //\r
   7 //      http://www.apache.org/licenses/LICENSE-2.0\r
   8 //\r
   9 // Unless required by applicable law or agreed to in writing, software\r
  10 // distributed under the License is distributed on an "AS IS" BASIS,\r
  11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\r
  12 // See the License for the specific language governing permissions and\r
  13 // limitations under the License.\r
  14 //\r
  15 ////////////////////////////////////////////////////////////////////////////////\r
  16 \r
  17 package com.google.crypto.tink.subtle;\r
  18 \r
  19 import static com.google.crypto.tink.subtle.Ed25519Constants.B2;\r
  20 import static com.google.crypto.tink.subtle.Ed25519Constants.B_TABLE;\r
  21 import static com.google.crypto.tink.subtle.Ed25519Constants.D;\r
  22 import static com.google.crypto.tink.subtle.Ed25519Constants.D2;\r
  23 import static com.google.crypto.tink.subtle.Ed25519Constants.SQRTM1;\r
  24 import static com.google.crypto.tink.subtle.Field25519.FIELD_LEN;\r
  25 import static com.google.crypto.tink.subtle.Field25519.LIMB_CNT;\r
  26 \r
  27 import java.security.GeneralSecurityException;\r
  28 import java.security.MessageDigest;\r
  29 import java.util.Arrays;\r
  30 \r
  31 /**\r
  32  * This implementation is based on the ed25519/ref10 implementation in NaCl.\r
  33  *\r
  34  * <p>It implements this twisted Edwards curve:\r
  35  *\r
  36  * <pre>\r
  37  * -x^2 + y^2 = 1 + (-121665 / 121666 mod 2^255-19)*x^2*y^2\r
  38  * </pre>\r
  39  *\r
  40  * @see <a href="https://eprint.iacr.org/2008/013.pdf">Bernstein D.J., Birkner P., Joye M., Lange\r
  41  * T., Peters C. (2008) Twisted Edwards Curves</a>\r
  42  * @see <a href="https://eprint.iacr.org/2008/522.pdf">Hisil H., Wong K.KH., Carter G., Dawson E.\r
  43  * (2008) Twisted Edwards Curves Revisited</a>\r
  44  */\r
  45 public final class Ed25519 {\r
  46 \r
  47     public static final int SECRET_KEY_LEN = FIELD_LEN;\r
  48     public static final int PUBLIC_KEY_LEN = FIELD_LEN;\r
  49     public static final int SIGNATURE_LEN = FIELD_LEN * 2;\r
  50 \r
  51     // (x = 0, y = 1) point\r
  52     private static final CachedXYT CACHED_NEUTRAL = new CachedXYT(\r
  53             new long[]{1, 0, 0, 0, 0, 0, 0, 0, 0, 0},\r
  54             new long[]{1, 0, 0, 0, 0, 0, 0, 0, 0, 0},\r
  55             new long[]{0, 0, 0, 0, 0, 0, 0, 0, 0, 0});\r
  56     private static final PartialXYZT NEUTRAL = new PartialXYZT(\r
  57             new XYZ(new long[]{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},\r
  58                     new long[]{1, 0, 0, 0, 0, 0, 0, 0, 0, 0},\r
  59                     new long[]{1, 0, 0, 0, 0, 0, 0, 0, 0, 0}),\r
  60             new long[]{1, 0, 0, 0, 0, 0, 0, 0, 0, 0});\r
  61 \r
  62     /**\r
  63      * Projective point representation (X:Y:Z) satisfying x = X/Z, y = Y/Z\r
  64      * <p>\r
  65      * Note that this is referred as ge_p2 in ref10 impl.\r
  66      * Also note that x = X, y = Y and z = Z below following Java coding style.\r
  67      * <p>\r
  68      * See\r
  69      * Koyama K., Tsuruoka Y. (1993) Speeding up Elliptic Cryptosystems by Using a Signed Binary\r
  70      * Window Method.\r
  71      * <p>\r
  72      * https://hyperelliptic.org/EFD/g1p/auto-twisted-projective.html\r
  73      */\r
  74     private static class XYZ {\r
  75 \r
  76         final long[] x;\r
  77         final long[] y;\r
  78         final long[] z;\r
  79 \r
  80         XYZ() {\r
  81             this(new long[LIMB_CNT], new long[LIMB_CNT], new long[LIMB_CNT]);\r
  82         }\r
  83 \r
  84         XYZ(long[] x, long[] y, long[] z) {\r
  85             this.x = x;\r
  86             this.y = y;\r
  87             this.z = z;\r
  88         }\r
  89 \r
  90         XYZ(XYZ xyz) {\r
  91             x = Arrays.copyOf(xyz.x, LIMB_CNT);\r
  92             y = Arrays.copyOf(xyz.y, LIMB_CNT);\r
  93             z = Arrays.copyOf(xyz.z, LIMB_CNT);\r
  94         }\r
  95 \r
  96         XYZ(PartialXYZT partialXYZT) {\r
  97             this();\r
  98             fromPartialXYZT(this, partialXYZT);\r
  99         }\r
 100 \r
 101         /**\r
 102          * ge_p1p1_to_p2.c\r
 103          */\r
 104         static XYZ fromPartialXYZT(XYZ out, PartialXYZT in) {\r
 105             Field25519.mult(out.x, in.xyz.x, in.t);\r
 106             Field25519.mult(out.y, in.xyz.y, in.xyz.z);\r
 107             Field25519.mult(out.z, in.xyz.z, in.t);\r
 108             return out;\r
 109         }\r
 110 \r
 111         /**\r
 112          * Encodes this point to bytes.\r
 113          */\r
 114         byte[] toBytes() {\r
 115             long[] recip = new long[LIMB_CNT];\r
 116             long[] x = new long[LIMB_CNT];\r
 117             long[] y = new long[LIMB_CNT];\r
 118             Field25519.inverse(recip, z);\r
 119             Field25519.mult(x, this.x, recip);\r
 120             Field25519.mult(y, this.y, recip);\r
 121             byte[] s = Field25519.contract(y);\r
 122             s[31] = (byte) (s[31] ^ (getLsb(x) << 7));\r
 123             return s;\r
 124         }\r
 125 \r
 126         /**\r
 127          * Checks that the point is on curve\r
 128          */\r
 129         boolean isOnCurve() {\r
 130             long[] x2 = new long[LIMB_CNT];\r
 131             Field25519.square(x2, x);\r
 132             long[] y2 = new long[LIMB_CNT];\r
 133             Field25519.square(y2, y);\r
 134             long[] z2 = new long[LIMB_CNT];\r
 135             Field25519.square(z2, z);\r
 136             long[] z4 = new long[LIMB_CNT];\r
 137             Field25519.square(z4, z2);\r
 138             long[] lhs = new long[LIMB_CNT];\r
 139             // lhs = y^2 - x^2\r
 140             Field25519.sub(lhs, y2, x2);\r
 141             // lhs = z^2 * (y2 - x2)\r
 142             Field25519.mult(lhs, lhs, z2);\r
 143             long[] rhs = new long[LIMB_CNT];\r
 144             // rhs = x^2 * y^2\r
 145             Field25519.mult(rhs, x2, y2);\r
 146             // rhs = D * x^2 * y^2\r
 147             Field25519.mult(rhs, rhs, D);\r
 148             // rhs = z^4 + D * x^2 * y^2\r
 149             Field25519.sum(rhs, z4);\r
 150             // z^2 (y^2 - x^2) == z^4 + D * x^2 * y^2\r
 151             return Bytes.equal(Field25519.contract(lhs), Field25519.contract(rhs));\r
 152         }\r
 153     }\r
 154 \r
 155     /**\r
 156      * Represents extended projective point representation (X:Y:Z:T) satisfying x = X/Z, y = Y/Z,\r
 157      * XY = ZT\r
 158      * <p>\r
 159      * Note that this is referred as ge_p3 in ref10 impl.\r
 160      * Also note that t = T below following Java coding style.\r
 161      * <p>\r
 162      * See\r
 163      * Hisil H., Wong K.KH., Carter G., Dawson E. (2008) Twisted Edwards Curves Revisited.\r
 164      * <p>\r
 165      * https://hyperelliptic.org/EFD/g1p/auto-twisted-extended.html\r
 166      */\r
 167     private static class XYZT {\r
 168 \r
 169         final XYZ xyz;\r
 170         final long[] t;\r
 171 \r
 172         XYZT() {\r
 173             this(new XYZ(), new long[LIMB_CNT]);\r
 174         }\r
 175 \r
 176         XYZT(XYZ xyz, long[] t) {\r
 177             this.xyz = xyz;\r
 178             this.t = t;\r
 179         }\r
 180 \r
 181         XYZT(PartialXYZT partialXYZT) {\r
 182             this();\r
 183             fromPartialXYZT(this, partialXYZT);\r
 184         }\r
 185 \r
 186         /**\r
 187          * ge_p1p1_to_p2.c\r
 188          */\r
 189         private static XYZT fromPartialXYZT(XYZT out, PartialXYZT in) {\r
 190             Field25519.mult(out.xyz.x, in.xyz.x, in.t);\r
 191             Field25519.mult(out.xyz.y, in.xyz.y, in.xyz.z);\r
 192             Field25519.mult(out.xyz.z, in.xyz.z, in.t);\r
 193             Field25519.mult(out.t, in.xyz.x, in.xyz.y);\r
 194             return out;\r
 195         }\r
 196 \r
 197         /**\r
 198          * Decodes {@code s} into an extented projective point.\r
 199          * See Section 5.1.3 Decoding in https://tools.ietf.org/html/rfc8032#section-5.1.3\r
 200          */\r
 201         private static XYZT fromBytesNegateVarTime(byte[] s) throws GeneralSecurityException {\r
 202             long[] x = new long[LIMB_CNT];\r
 203             long[] y = Field25519.expand(s);\r
 204             long[] z = new long[LIMB_CNT];\r
 205             z[0] = 1;\r
 206             long[] t = new long[LIMB_CNT];\r
 207             long[] u = new long[LIMB_CNT];\r
 208             long[] v = new long[LIMB_CNT];\r
 209             long[] vxx = new long[LIMB_CNT];\r
 210             long[] check = new long[LIMB_CNT];\r
 211             Field25519.square(u, y);\r
 212             Field25519.mult(v, u, D);\r
 213             Field25519.sub(u, u, z); // u = y^2 - 1\r
 214             Field25519.sum(v, v, z); // v = dy^2 + 1\r
 215 \r
 216             long[] v3 = new long[LIMB_CNT];\r
 217             Field25519.square(v3, v);\r
 218             Field25519.mult(v3, v3, v); // v3 = v^3\r
 219             Field25519.square(x, v3);\r
 220             Field25519.mult(x, x, v);\r
 221             Field25519.mult(x, x, u); // x = uv^7\r
 222 \r
 223             pow2252m3(x, x); // x = (uv^7)^((q-5)/8)\r
 224             Field25519.mult(x, x, v3);\r
 225             Field25519.mult(x, x, u); // x = uv^3(uv^7)^((q-5)/8)\r
 226 \r
 227             Field25519.square(vxx, x);\r
 228             Field25519.mult(vxx, vxx, v);\r
 229             Field25519.sub(check, vxx, u); // vx^2-u\r
 230             if (isNonZeroVarTime(check)) {\r
 231                 Field25519.sum(check, vxx, u); // vx^2+u\r
 232                 if (isNonZeroVarTime(check)) {\r
 233                     throw new GeneralSecurityException("Cannot convert given bytes to extended projective "\r
 234                             + "coordinates. No square root exists for modulo 2^255-19");\r
 235                 }\r
 236                 Field25519.mult(x, x, SQRTM1);\r
 237             }\r
 238 \r
 239             if (!isNonZeroVarTime(x) && (s[31] & 0xff) >> 7 != 0) {\r
 240                 throw new GeneralSecurityException("Cannot convert given bytes to extended projective "\r
 241                         + "coordinates. Computed x is zero and encoded x's least significant bit is not zero");\r
 242             }\r
 243             if (getLsb(x) == ((s[31] & 0xff) >> 7)) {\r
 244                 neg(x, x);\r
 245             }\r
 246 \r
 247             Field25519.mult(t, x, y);\r
 248             return new XYZT(new XYZ(x, y, z), t);\r
 249         }\r
 250     }\r
 251 \r
 252     /**\r
 253      * Partial projective point representation ((X:Z),(Y:T)) satisfying x=X/Z, y=Y/T\r
 254      * <p>\r
 255      * Note that this is referred as complete form in the original ref10 impl (ge_p1p1).\r
 256      * Also note that t = T below following Java coding style.\r
 257      * <p>\r
 258      * Although this has the same types as XYZT, it is redefined to have its own type so that it is\r
 259      * readable and 1:1 corresponds to ref10 impl.\r
 260      * <p>\r
 261      * Can be converted to XYZT as follows:\r
 262      * X1 = X * T = x * Z * T = x * Z1\r
 263      * Y1 = Y * Z = y * T * Z = y * Z1\r
 264      * Z1 = Z * T = Z * T\r
 265      * T1 = X * Y = x * Z * y * T = x * y * Z1 = X1Y1 / Z1\r
 266      */\r
 267     private static class PartialXYZT {\r
 268 \r
 269         final XYZ xyz;\r
 270         final long[] t;\r
 271 \r
 272         PartialXYZT() {\r
 273             this(new XYZ(), new long[LIMB_CNT]);\r
 274         }\r
 275 \r
 276         PartialXYZT(XYZ xyz, long[] t) {\r
 277             this.xyz = xyz;\r
 278             this.t = t;\r
 279         }\r
 280 \r
 281         PartialXYZT(PartialXYZT other) {\r
 282             xyz = new XYZ(other.xyz);\r
 283             t = Arrays.copyOf(other.t, LIMB_CNT);\r
 284         }\r
 285     }\r
 286 \r
 287     /**\r
 288      * Corresponds to the caching mentioned in the last paragraph of Section 3.1 of\r
 289      * Hisil H., Wong K.KH., Carter G., Dawson E. (2008) Twisted Edwards Curves Revisited.\r
 290      * with Z = 1.\r
 291      */\r
 292     static class CachedXYT {\r
 293 \r
 294         final long[] yPlusX;\r
 295         final long[] yMinusX;\r
 296         final long[] t2d;\r
 297 \r
 298         CachedXYT() {\r
 299             this(new long[LIMB_CNT], new long[LIMB_CNT], new long[LIMB_CNT]);\r
 300         }\r
 301 \r
 302         /**\r
 303          * Creates a cached XYZT with Z = 1\r
 304          *\r
 305          * @param yPlusX  y + x\r
 306          * @param yMinusX y - x\r
 307          * @param t2d     2d * xy\r
 308          */\r
 309         CachedXYT(long[] yPlusX, long[] yMinusX, long[] t2d) {\r
 310             this.yPlusX = yPlusX;\r
 311             this.yMinusX = yMinusX;\r
 312             this.t2d = t2d;\r
 313         }\r
 314 \r
 315         CachedXYT(CachedXYT other) {\r
 316             yPlusX = Arrays.copyOf(other.yPlusX, LIMB_CNT);\r
 317             yMinusX = Arrays.copyOf(other.yMinusX, LIMB_CNT);\r
 318             t2d = Arrays.copyOf(other.t2d, LIMB_CNT);\r
 319         }\r
 320 \r
 321         // z is one implicitly, so this just copies {@code in} to {@code output}.\r
 322         void multByZ(long[] output, long[] in) {\r
 323             System.arraycopy(in, 0, output, 0, LIMB_CNT);\r
 324         }\r
 325 \r
 326         /**\r
 327          * If icopy is 1, copies {@code other} into this point. Time invariant wrt to icopy value.\r
 328          */\r
 329         void copyConditional(CachedXYT other, int icopy) {\r
 330             Curve25519.copyConditional(yPlusX, other.yPlusX, icopy);\r
 331             Curve25519.copyConditional(yMinusX, other.yMinusX, icopy);\r
 332             Curve25519.copyConditional(t2d, other.t2d, icopy);\r
 333         }\r
 334     }\r
 335 \r
 336     private static class CachedXYZT extends CachedXYT {\r
 337 \r
 338         private final long[] z;\r
 339 \r
 340         CachedXYZT() {\r
 341             this(new long[LIMB_CNT], new long[LIMB_CNT], new long[LIMB_CNT], new long[LIMB_CNT]);\r
 342         }\r
 343 \r
 344         /**\r
 345          * ge_p3_to_cached.c\r
 346          */\r
 347         CachedXYZT(XYZT xyzt) {\r
 348             this();\r
 349             Field25519.sum(yPlusX, xyzt.xyz.y, xyzt.xyz.x);\r
 350             Field25519.sub(yMinusX, xyzt.xyz.y, xyzt.xyz.x);\r
 351             System.arraycopy(xyzt.xyz.z, 0, z, 0, LIMB_CNT);\r
 352             Field25519.mult(t2d, xyzt.t, D2);\r
 353         }\r
 354 \r
 355         /**\r
 356          * Creates a cached XYZT\r
 357          *\r
 358          * @param yPlusX  Y + X\r
 359          * @param yMinusX Y - X\r
 360          * @param z       Z\r
 361          * @param t2d     2d * (XY/Z)\r
 362          */\r
 363         CachedXYZT(long[] yPlusX, long[] yMinusX, long[] z, long[] t2d) {\r
 364             super(yPlusX, yMinusX, t2d);\r
 365             this.z = z;\r
 366         }\r
 367 \r
 368         @Override\r
 369         public void multByZ(long[] output, long[] in) {\r
 370             Field25519.mult(output, in, z);\r
 371         }\r
 372     }\r
 373 \r
 374     /**\r
 375      * Addition defined in Section 3.1 of\r
 376      * Hisil H., Wong K.KH., Carter G., Dawson E. (2008) Twisted Edwards Curves Revisited.\r
 377      * <p>\r
 378      * Please note that this is a partial of the operation listed there leaving out the final\r
 379      * conversion from PartialXYZT to XYZT.\r
 380      *\r
 381      * @param extended extended projective point input\r
 382      * @param cached   cached projective point input\r
 383      */\r
 384     private static void add(PartialXYZT partialXYZT, XYZT extended, CachedXYT cached) {\r
 385         long[] t = new long[LIMB_CNT];\r
 386 \r
 387         // Y1 + X1\r
 388         Field25519.sum(partialXYZT.xyz.x, extended.xyz.y, extended.xyz.x);\r
 389 \r
 390         // Y1 - X1\r
 391         Field25519.sub(partialXYZT.xyz.y, extended.xyz.y, extended.xyz.x);\r
 392 \r
 393         // A = (Y1 - X1) * (Y2 - X2)\r
 394         Field25519.mult(partialXYZT.xyz.y, partialXYZT.xyz.y, cached.yMinusX);\r
 395 \r
 396         // B = (Y1 + X1) * (Y2 + X2)\r
 397         Field25519.mult(partialXYZT.xyz.z, partialXYZT.xyz.x, cached.yPlusX);\r
 398 \r
 399         // C = T1 * 2d * T2 = 2d * T1 * T2 (2d is written as k in the paper)\r
 400         Field25519.mult(partialXYZT.t, extended.t, cached.t2d);\r
 401 \r
 402         // Z1 * Z2\r
 403         cached.multByZ(partialXYZT.xyz.x, extended.xyz.z);\r
 404 \r
 405         // D = 2 * Z1 * Z2\r
 406         Field25519.sum(t, partialXYZT.xyz.x, partialXYZT.xyz.x);\r
 407 \r
 408         // X3 = B - A\r
 409         Field25519.sub(partialXYZT.xyz.x, partialXYZT.xyz.z, partialXYZT.xyz.y);\r
 410 \r
 411         // Y3 = B + A\r
 412         Field25519.sum(partialXYZT.xyz.y, partialXYZT.xyz.z, partialXYZT.xyz.y);\r
 413 \r
 414         // Z3 = D + C\r
 415         Field25519.sum(partialXYZT.xyz.z, t, partialXYZT.t);\r
 416 \r
 417         // T3 = D - C\r
 418         Field25519.sub(partialXYZT.t, t, partialXYZT.t);\r
 419     }\r
 420 \r
 421     /**\r
 422      * Based on the addition defined in Section 3.1 of\r
 423      * Hisil H., Wong K.KH., Carter G., Dawson E. (2008) Twisted Edwards Curves Revisited.\r
 424      * <p>\r
 425      * Please note that this is a partial of the operation listed there leaving out the final\r
 426      * conversion from PartialXYZT to XYZT.\r
 427      *\r
 428      * @param extended extended projective point input\r
 429      * @param cached   cached projective point input\r
 430      */\r
 431     private static void sub(PartialXYZT partialXYZT, XYZT extended, CachedXYT cached) {\r
 432         long[] t = new long[LIMB_CNT];\r
 433 \r
 434         // Y1 + X1\r
 435         Field25519.sum(partialXYZT.xyz.x, extended.xyz.y, extended.xyz.x);\r
 436 \r
 437         // Y1 - X1\r
 438         Field25519.sub(partialXYZT.xyz.y, extended.xyz.y, extended.xyz.x);\r
 439 \r
 440         // A = (Y1 - X1) * (Y2 + X2)\r
 441         Field25519.mult(partialXYZT.xyz.y, partialXYZT.xyz.y, cached.yPlusX);\r
 442 \r
 443         // B = (Y1 + X1) * (Y2 - X2)\r
 444         Field25519.mult(partialXYZT.xyz.z, partialXYZT.xyz.x, cached.yMinusX);\r
 445 \r
 446         // C = T1 * 2d * T2 = 2d * T1 * T2 (2d is written as k in the paper)\r
 447         Field25519.mult(partialXYZT.t, extended.t, cached.t2d);\r
 448 \r
 449         // Z1 * Z2\r
 450         cached.multByZ(partialXYZT.xyz.x, extended.xyz.z);\r
 451 \r
 452         // D = 2 * Z1 * Z2\r
 453         Field25519.sum(t, partialXYZT.xyz.x, partialXYZT.xyz.x);\r
 454 \r
 455         // X3 = B - A\r
 456         Field25519.sub(partialXYZT.xyz.x, partialXYZT.xyz.z, partialXYZT.xyz.y);\r
 457 \r
 458         // Y3 = B + A\r
 459         Field25519.sum(partialXYZT.xyz.y, partialXYZT.xyz.z, partialXYZT.xyz.y);\r
 460 \r
 461         // Z3 = D - C\r
 462         Field25519.sub(partialXYZT.xyz.z, t, partialXYZT.t);\r
 463 \r
 464         // T3 = D + C\r
 465         Field25519.sum(partialXYZT.t, t, partialXYZT.t);\r
 466     }\r
 467 \r
 468     /**\r
 469      * Doubles {@code p} and puts the result into this PartialXYZT.\r
 470      * <p>\r
 471      * This is based on the addition defined in formula 7 in Section 3.3 of\r
 472      * Hisil H., Wong K.KH., Carter G., Dawson E. (2008) Twisted Edwards Curves Revisited.\r
 473      * <p>\r
 474      * Please note that this is a partial of the operation listed there leaving out the final\r
 475      * conversion from PartialXYZT to XYZT and also this fixes a typo in calculation of Y3 and T3 in\r
 476      * the paper, H should be replaced with A+B.\r
 477      */\r
 478     private static void doubleXYZ(PartialXYZT partialXYZT, XYZ p) {\r
 479         long[] t0 = new long[LIMB_CNT];\r
 480 \r
 481         // XX = X1^2\r
 482         Field25519.square(partialXYZT.xyz.x, p.x);\r
 483 \r
 484         // YY = Y1^2\r
 485         Field25519.square(partialXYZT.xyz.z, p.y);\r
 486 \r
 487         // B' = Z1^2\r
 488         Field25519.square(partialXYZT.t, p.z);\r
 489 \r
 490         // B = 2 * B'\r
 491         Field25519.sum(partialXYZT.t, partialXYZT.t, partialXYZT.t);\r
 492 \r
 493         // A = X1 + Y1\r
 494         Field25519.sum(partialXYZT.xyz.y, p.x, p.y);\r
 495 \r
 496         // AA = A^2\r
 497         Field25519.square(t0, partialXYZT.xyz.y);\r
 498 \r
 499         // Y3 = YY + XX\r
 500         Field25519.sum(partialXYZT.xyz.y, partialXYZT.xyz.z, partialXYZT.xyz.x);\r
 501 \r
 502         // Z3 = YY - XX\r
 503         Field25519.sub(partialXYZT.xyz.z, partialXYZT.xyz.z, partialXYZT.xyz.x);\r
 504 \r
 505         // X3 = AA - Y3\r
 506         Field25519.sub(partialXYZT.xyz.x, t0, partialXYZT.xyz.y);\r
 507 \r
 508         // T3 = B - Z3\r
 509         Field25519.sub(partialXYZT.t, partialXYZT.t, partialXYZT.xyz.z);\r
 510     }\r
 511 \r
 512     /**\r
 513      * Doubles {@code p} and puts the result into this PartialXYZT.\r
 514      */\r
 515     private static void doubleXYZT(PartialXYZT partialXYZT, XYZT p) {\r
 516         doubleXYZ(partialXYZT, p.xyz);\r
 517     }\r
 518 \r
 519     /**\r
 520      * Compares two byte values in constant time.\r
 521      * <p>\r
 522      * Please note that this doesn't reuse {@link Curve25519#eq} method since the below inputs are\r
 523      * byte values.\r
 524      */\r
 525     private static int eq(int a, int b) {\r
 526         int r = ~(a ^ b) & 0xff;\r
 527         r &= r << 4;\r
 528         r &= r << 2;\r
 529         r &= r << 1;\r
 530         return (r >> 7) & 1;\r
 531     }\r
 532 \r
 533     /**\r
 534      * This is a constant time operation where point b*B*256^pos is stored in {@code t}.\r
 535      * When b is 0, t remains the same (i.e., neutral point).\r
 536      * <p>\r
 537      * Although B_TABLE[32][8] (B_TABLE[i][j] = (j+1)*B*256^i) has j values in [0, 7], the select\r
 538      * method negates the corresponding point if b is negative (which is straight forward in elliptic\r
 539      * curves by just negating y coordinate). Therefore we can get multiples of B with the half of\r
 540      * memory requirements.\r
 541      *\r
 542      * @param t   neutral element (i.e., point 0), also serves as output.\r
 543      * @param pos in B[pos][j] = (j+1)*B*256^pos\r
 544      * @param b   value in [-8, 8] range.\r
 545      */\r
 546     private static void select(CachedXYT t, int pos, byte b) {\r
 547         int bnegative = (b & 0xff) >> 7;\r
 548         int babs = b - (((-bnegative) & b) << 1);\r
 549 \r
 550         t.copyConditional(B_TABLE[pos][0], eq(babs, 1));\r
 551         t.copyConditional(B_TABLE[pos][1], eq(babs, 2));\r
 552         t.copyConditional(B_TABLE[pos][2], eq(babs, 3));\r
 553         t.copyConditional(B_TABLE[pos][3], eq(babs, 4));\r
 554         t.copyConditional(B_TABLE[pos][4], eq(babs, 5));\r
 555         t.copyConditional(B_TABLE[pos][5], eq(babs, 6));\r
 556         t.copyConditional(B_TABLE[pos][6], eq(babs, 7));\r
 557         t.copyConditional(B_TABLE[pos][7], eq(babs, 8));\r
 558 \r
 559         long[] yPlusX = Arrays.copyOf(t.yMinusX, LIMB_CNT);\r
 560         long[] yMinusX = Arrays.copyOf(t.yPlusX, LIMB_CNT);\r
 561         long[] t2d = Arrays.copyOf(t.t2d, LIMB_CNT);\r
 562         neg(t2d, t2d);\r
 563         CachedXYT minust = new CachedXYT(yPlusX, yMinusX, t2d);\r
 564         t.copyConditional(minust, bnegative);\r
 565     }\r
 566 \r
 567     /**\r
 568      * Computes {@code a}*B\r
 569      * where a = a[0]+256*a[1]+...+256^31 a[31] and\r
 570      * B is the Ed25519 base point (x,4/5) with x positive.\r
 571      * <p>\r
 572      * Preconditions:\r
 573      * a[31] <= 127\r
 574      *\r
 575      * @throws IllegalStateException iff there is arithmetic error.\r
 576      */\r
 577     @SuppressWarnings("NarrowingCompoundAssignment")\r
 578     private static XYZ scalarMultWithBase(byte[] a) {\r
 579         byte[] e = new byte[2 * FIELD_LEN];\r
 580         for (int i = 0; i < FIELD_LEN; i++) {\r
 581             e[2 * i + 0] = (byte) (((a[i] & 0xff) >> 0) & 0xf);\r
 582             e[2 * i + 1] = (byte) (((a[i] & 0xff) >> 4) & 0xf);\r
 583         }\r
 584         // each e[i] is between 0 and 15\r
 585         // e[63] is between 0 and 7\r
 586 \r
 587         // Rewrite e in a way that each e[i] is in [-8, 8].\r
 588         // This can be done since a[63] is in [0, 7], the carry-over onto the most significant byte\r
 589         // a[63] can be at most 1.\r
 590         int carry = 0;\r
 591         for (int i = 0; i < e.length - 1; i++) {\r
 592             e[i] += carry;\r
 593             carry = e[i] + 8;\r
 594             carry >>= 4;\r
 595             e[i] -= carry << 4;\r
 596         }\r
 597         e[e.length - 1] += carry;\r
 598 \r
 599         PartialXYZT ret = new PartialXYZT(NEUTRAL);\r
 600         XYZT xyzt = new XYZT();\r
 601         // Although B_TABLE's i can be at most 31 (stores only 32 4bit multiples of B) and we have 64\r
 602         // 4bit values in e array, the below for loop adds cached values by iterating e by two in odd\r
 603         // indices. After the result, we can double the result point 4 times to shift the multiplication\r
 604         // scalar by 4 bits.\r
 605         for (int i = 1; i < e.length; i += 2) {\r
 606             CachedXYT t = new CachedXYT(CACHED_NEUTRAL);\r
 607             select(t, i / 2, e[i]);\r
 608             add(ret, XYZT.fromPartialXYZT(xyzt, ret), t);\r
 609         }\r
 610 \r
 611         // Doubles the result 4 times to shift the multiplication scalar 4 bits to get the actual result\r
 612         // for the odd indices in e.\r
 613         XYZ xyz = new XYZ();\r
 614         doubleXYZ(ret, XYZ.fromPartialXYZT(xyz, ret));\r
 615         doubleXYZ(ret, XYZ.fromPartialXYZT(xyz, ret));\r
 616         doubleXYZ(ret, XYZ.fromPartialXYZT(xyz, ret));\r
 617         doubleXYZ(ret, XYZ.fromPartialXYZT(xyz, ret));\r
 618 \r
 619         // Add multiples of B for even indices of e.\r
 620         for (int i = 0; i < e.length; i += 2) {\r
 621             CachedXYT t = new CachedXYT(CACHED_NEUTRAL);\r
 622             select(t, i / 2, e[i]);\r
 623             add(ret, XYZT.fromPartialXYZT(xyzt, ret), t);\r
 624         }\r
 625 \r
 626         // This check is to protect against flaws, i.e. if there is a computation error through a\r
 627         // faulty CPU or if the implementation contains a bug.\r
 628         XYZ result = new XYZ(ret);\r
 629         if (!result.isOnCurve()) {\r
 630             throw new IllegalStateException("arithmetic error in scalar multiplication");\r
 631         }\r
 632         return result;\r
 633     }\r
 634 \r
 635     /**\r
 636      * Computes {@code a}*B\r
 637      * where a = a[0]+256*a[1]+...+256^31 a[31] and\r
 638      * B is the Ed25519 base point (x,4/5) with x positive.\r
 639      * <p>\r
 640      * Preconditions:\r
 641      * a[31] <= 127\r
 642      */\r
 643     public static byte[] scalarMultWithBaseToBytes(byte[] a) {\r
 644         return scalarMultWithBase(a).toBytes();\r
 645     }\r
 646 \r
 647     @SuppressWarnings("NarrowingCompoundAssignment")\r
 648     private static byte[] slide(byte[] a) {\r
 649         byte[] r = new byte[256];\r
 650         // Writes each bit in a[0..31] into r[0..255]:\r
 651         // a = a[0]+256*a[1]+...+256^31*a[31] is equal to\r
 652         // r = r[0]+2*r[1]+...+2^255*r[255]\r
 653         for (int i = 0; i < 256; i++) {\r
 654             r[i] = (byte) (1 & ((a[i >> 3] & 0xff) >> (i & 7)));\r
 655         }\r
 656 \r
 657         // Transforms r[i] as odd values in [-15, 15]\r
 658         for (int i = 0; i < 256; i++) {\r
 659             if (r[i] != 0) {\r
 660                 for (int b = 1; b <= 6 && i + b < 256; b++) {\r
 661                     if (r[i + b] != 0) {\r
 662                         if (r[i] + (r[i + b] << b) <= 15) {\r
 663                             r[i] += r[i + b] << b;\r
 664                             r[i + b] = 0;\r
 665                         } else if (r[i] - (r[i + b] << b) >= -15) {\r
 666                             r[i] -= r[i + b] << b;\r
 667                             for (int k = i + b; k < 256; k++) {\r
 668                                 if (r[k] == 0) {\r
 669                                     r[k] = 1;\r
 670                                     break;\r
 671                                 }\r
 672                                 r[k] = 0;\r
 673                             }\r
 674                         } else {\r
 675                             break;\r
 676                         }\r
 677                     }\r
 678                 }\r
 679             }\r
 680         }\r
 681         return r;\r
 682     }\r
 683 \r
 684     /**\r
 685      * Computes {@code a}*{@code pointA}+{@code b}*B\r
 686      * where a = a[0]+256*a[1]+...+256^31*a[31].\r
 687      * and b = b[0]+256*b[1]+...+256^31*b[31].\r
 688      * B is the Ed25519 base point (x,4/5) with x positive.\r
 689      * <p>\r
 690      * Note that execution time varies based on the input since this will only be used in verification\r
 691      * of signatures.\r
 692      */\r
 693     private static XYZ doubleScalarMultVarTime(byte[] a, XYZT pointA, byte[] b) {\r
 694         // pointA, 3*pointA, 5*pointA, 7*pointA, 9*pointA, 11*pointA, 13*pointA, 15*pointA\r
 695         CachedXYZT[] pointAArray = new CachedXYZT[8];\r
 696         pointAArray[0] = new CachedXYZT(pointA);\r
 697         PartialXYZT t = new PartialXYZT();\r
 698         doubleXYZT(t, pointA);\r
 699         XYZT doubleA = new XYZT(t);\r
 700         for (int i = 1; i < pointAArray.length; i++) {\r
 701             add(t, doubleA, pointAArray[i - 1]);\r
 702             pointAArray[i] = new CachedXYZT(new XYZT(t));\r
 703         }\r
 704 \r
 705         byte[] aSlide = slide(a);\r
 706         byte[] bSlide = slide(b);\r
 707         t = new PartialXYZT(NEUTRAL);\r
 708         XYZT u = new XYZT();\r
 709         int i = 255;\r
 710         for (; i >= 0; i--) {\r
 711             if (aSlide[i] != 0 || bSlide[i] != 0) {\r
 712                 break;\r
 713             }\r
 714         }\r
 715         for (; i >= 0; i--) {\r
 716             doubleXYZ(t, new XYZ(t));\r
 717             if (aSlide[i] > 0) {\r
 718                 add(t, XYZT.fromPartialXYZT(u, t), pointAArray[aSlide[i] / 2]);\r
 719             } else if (aSlide[i] < 0) {\r
 720                 sub(t, XYZT.fromPartialXYZT(u, t), pointAArray[-aSlide[i] / 2]);\r
 721             }\r
 722             if (bSlide[i] > 0) {\r
 723                 add(t, XYZT.fromPartialXYZT(u, t), B2[bSlide[i] / 2]);\r
 724             } else if (bSlide[i] < 0) {\r
 725                 sub(t, XYZT.fromPartialXYZT(u, t), B2[-bSlide[i] / 2]);\r
 726             }\r
 727         }\r
 728 \r
 729         return new XYZ(t);\r
 730     }\r
 731 \r
 732     /**\r
 733      * Returns true if {@code in} is nonzero.\r
 734      * <p>\r
 735      * Note that execution time might depend on the input {@code in}.\r
 736      */\r
 737     private static boolean isNonZeroVarTime(long[] in) {\r
 738         long[] inCopy = new long[in.length + 1];\r
 739         System.arraycopy(in, 0, inCopy, 0, in.length);\r
 740         Field25519.reduceCoefficients(inCopy);\r
 741         byte[] bytes = Field25519.contract(inCopy);\r
 742         for (byte b : bytes) {\r
 743             if (b != 0) {\r
 744                 return true;\r
 745             }\r
 746         }\r
 747         return false;\r
 748     }\r
 749 \r
 750     /**\r
 751      * Returns the least significant bit of {@code in}.\r
 752      */\r
 753     private static int getLsb(long[] in) {\r
 754         return Field25519.contract(in)[0] & 1;\r
 755     }\r
 756 \r
 757     /**\r
 758      * Negates all values in {@code in} and store it in {@code out}.\r
 759      */\r
 760     private static void neg(long[] out, long[] in) {\r
 761         for (int i = 0; i < in.length; i++) {\r
 762             out[i] = -in[i];\r
 763         }\r
 764     }\r
 765 \r
 766     /**\r
 767      * Computes {@code in}^(2^252-3) mod 2^255-19 and puts the result in {@code out}.\r
 768      */\r
 769     private static void pow2252m3(long[] out, long[] in) {\r
 770         long[] t0 = new long[LIMB_CNT];\r
 771         long[] t1 = new long[LIMB_CNT];\r
 772         long[] t2 = new long[LIMB_CNT];\r
 773 \r
 774         // z2 = z1^2^1\r
 775         Field25519.square(t0, in);\r
 776 \r
 777         // z8 = z2^2^2\r
 778         Field25519.square(t1, t0);\r
 779         for (int i = 1; i < 2; i++) {\r
 780             Field25519.square(t1, t1);\r
 781         }\r
 782 \r
 783         // z9 = z1*z8\r
 784         Field25519.mult(t1, in, t1);\r
 785 \r
 786         // z11 = z2*z9\r
 787         Field25519.mult(t0, t0, t1);\r
 788 \r
 789         // z22 = z11^2^1\r
 790         Field25519.square(t0, t0);\r
 791 \r
 792         // z_5_0 = z9*z22\r
 793         Field25519.mult(t0, t1, t0);\r
 794 \r
 795         // z_10_5 = z_5_0^2^5\r
 796         Field25519.square(t1, t0);\r
 797         for (int i = 1; i < 5; i++) {\r
 798             Field25519.square(t1, t1);\r
 799         }\r
 800 \r
 801         // z_10_0 = z_10_5*z_5_0\r
 802         Field25519.mult(t0, t1, t0);\r
 803 \r
 804         // z_20_10 = z_10_0^2^10\r
 805         Field25519.square(t1, t0);\r
 806         for (int i = 1; i < 10; i++) {\r
 807             Field25519.square(t1, t1);\r
 808         }\r
 809 \r
 810         // z_20_0 = z_20_10*z_10_0\r
 811         Field25519.mult(t1, t1, t0);\r
 812 \r
 813         // z_40_20 = z_20_0^2^20\r
 814         Field25519.square(t2, t1);\r
 815         for (int i = 1; i < 20; i++) {\r
 816             Field25519.square(t2, t2);\r
 817         }\r
 818 \r
 819         // z_40_0 = z_40_20*z_20_0\r
 820         Field25519.mult(t1, t2, t1);\r
 821 \r
 822         // z_50_10 = z_40_0^2^10\r
 823         Field25519.square(t1, t1);\r
 824         for (int i = 1; i < 10; i++) {\r
 825             Field25519.square(t1, t1);\r
 826         }\r
 827 \r
 828         // z_50_0 = z_50_10*z_10_0\r
 829         Field25519.mult(t0, t1, t0);\r
 830 \r
 831         // z_100_50 = z_50_0^2^50\r
 832         Field25519.square(t1, t0);\r
 833         for (int i = 1; i < 50; i++) {\r
 834             Field25519.square(t1, t1);\r
 835         }\r
 836 \r
 837         // z_100_0 = z_100_50*z_50_0\r
 838         Field25519.mult(t1, t1, t0);\r
 839 \r
 840         // z_200_100 = z_100_0^2^100\r
 841         Field25519.square(t2, t1);\r
 842         for (int i = 1; i < 100; i++) {\r
 843             Field25519.square(t2, t2);\r
 844         }\r
 845 \r
 846         // z_200_0 = z_200_100*z_100_0\r
 847         Field25519.mult(t1, t2, t1);\r
 848 \r
 849         // z_250_50 = z_200_0^2^50\r
 850         Field25519.square(t1, t1);\r
 851         for (int i = 1; i < 50; i++) {\r
 852             Field25519.square(t1, t1);\r
 853         }\r
 854 \r
 855         // z_250_0 = z_250_50*z_50_0\r
 856         Field25519.mult(t0, t1, t0);\r
 857 \r
 858         // z_252_2 = z_250_0^2^2\r
 859         Field25519.square(t0, t0);\r
 860         for (int i = 1; i < 2; i++) {\r
 861             Field25519.square(t0, t0);\r
 862         }\r
 863 \r
 864         // z_252_3 = z_252_2*z1\r
 865         Field25519.mult(out, t0, in);\r
 866     }\r
 867 \r
 868     /**\r
 869      * Returns 3 bytes of {@code in} starting from {@code idx} in Little-Endian format.\r
 870      */\r
 871     private static long load3(byte[] in, int idx) {\r
 872         long result;\r
 873         result = (long) in[idx] & 0xff;\r
 874         result |= (long) (in[idx + 1] & 0xff) << 8;\r
 875         result |= (long) (in[idx + 2] & 0xff) << 16;\r
 876         return result;\r
 877     }\r
 878 \r
 879     /**\r
 880      * Returns 4 bytes of {@code in} starting from {@code idx} in Little-Endian format.\r
 881      */\r
 882     private static long load4(byte[] in, int idx) {\r
 883         long result = load3(in, idx);\r
 884         result |= (long) (in[idx + 3] & 0xff) << 24;\r
 885         return result;\r
 886     }\r
 887 \r
 888     /**\r
 889      * Input:\r
 890      * s[0]+256*s[1]+...+256^63*s[63] = s\r
 891      * <p>\r
 892      * Output:\r
 893      * s[0]+256*s[1]+...+256^31*s[31] = s mod l\r
 894      * where l = 2^252 + 27742317777372353535851937790883648493.\r
 895      * Overwrites s in place.\r
 896      */\r
 897     public static void reduce(byte[] s) {\r
 898         // Observation:\r
 899         // 2^252 mod l is equivalent to -27742317777372353535851937790883648493 mod l\r
 900         // Let m = -27742317777372353535851937790883648493\r
 901         // Thus a*2^252+b mod l is equivalent to a*m+b mod l\r
 902         //\r
 903         // First s is divided into chunks of 21 bits as follows:\r
 904         // s0+2^21*s1+2^42*s3+...+2^462*s23 = s[0]+256*s[1]+...+256^63*s[63]\r
 905         long s0 = 2097151 & load3(s, 0);\r
 906         long s1 = 2097151 & (load4(s, 2) >> 5);\r
 907         long s2 = 2097151 & (load3(s, 5) >> 2);\r
 908         long s3 = 2097151 & (load4(s, 7) >> 7);\r
 909         long s4 = 2097151 & (load4(s, 10) >> 4);\r
 910         long s5 = 2097151 & (load3(s, 13) >> 1);\r
 911         long s6 = 2097151 & (load4(s, 15) >> 6);\r
 912         long s7 = 2097151 & (load3(s, 18) >> 3);\r
 913         long s8 = 2097151 & load3(s, 21);\r
 914         long s9 = 2097151 & (load4(s, 23) >> 5);\r
 915         long s10 = 2097151 & (load3(s, 26) >> 2);\r
 916         long s11 = 2097151 & (load4(s, 28) >> 7);\r
 917         long s12 = 2097151 & (load4(s, 31) >> 4);\r
 918         long s13 = 2097151 & (load3(s, 34) >> 1);\r
 919         long s14 = 2097151 & (load4(s, 36) >> 6);\r
 920         long s15 = 2097151 & (load3(s, 39) >> 3);\r
 921         long s16 = 2097151 & load3(s, 42);\r
 922         long s17 = 2097151 & (load4(s, 44) >> 5);\r
 923         long s18 = 2097151 & (load3(s, 47) >> 2);\r
 924         long s19 = 2097151 & (load4(s, 49) >> 7);\r
 925         long s20 = 2097151 & (load4(s, 52) >> 4);\r
 926         long s21 = 2097151 & (load3(s, 55) >> 1);\r
 927         long s22 = 2097151 & (load4(s, 57) >> 6);\r
 928         long s23 = (load4(s, 60) >> 3);\r
 929         long carry0;\r
 930         long carry1;\r
 931         long carry2;\r
 932         long carry3;\r
 933         long carry4;\r
 934         long carry5;\r
 935         long carry6;\r
 936         long carry7;\r
 937         long carry8;\r
 938         long carry9;\r
 939         long carry10;\r
 940         long carry11;\r
 941         long carry12;\r
 942         long carry13;\r
 943         long carry14;\r
 944         long carry15;\r
 945         long carry16;\r
 946 \r
 947         // s23*2^462 = s23*2^210*2^252 is equivalent to s23*2^210*m in mod l\r
 948         // As m is a 125 bit number, the result needs to scattered to 6 limbs (125/21 ceil is 6)\r
 949         // starting from s11 (s11*2^210)\r
 950         // m = [666643, 470296, 654183, -997805, 136657, -683901] in 21-bit limbs\r
 951         s11 += s23 * 666643;\r
 952         s12 += s23 * 470296;\r
 953         s13 += s23 * 654183;\r
 954         s14 -= s23 * 997805;\r
 955         s15 += s23 * 136657;\r
 956         s16 -= s23 * 683901;\r
 957         // s23 = 0;\r
 958 \r
 959         s10 += s22 * 666643;\r
 960         s11 += s22 * 470296;\r
 961         s12 += s22 * 654183;\r
 962         s13 -= s22 * 997805;\r
 963         s14 += s22 * 136657;\r
 964         s15 -= s22 * 683901;\r
 965         // s22 = 0;\r
 966 \r
 967         s9 += s21 * 666643;\r
 968         s10 += s21 * 470296;\r
 969         s11 += s21 * 654183;\r
 970         s12 -= s21 * 997805;\r
 971         s13 += s21 * 136657;\r
 972         s14 -= s21 * 683901;\r
 973         // s21 = 0;\r
 974 \r
 975         s8 += s20 * 666643;\r
 976         s9 += s20 * 470296;\r
 977         s10 += s20 * 654183;\r
 978         s11 -= s20 * 997805;\r
 979         s12 += s20 * 136657;\r
 980         s13 -= s20 * 683901;\r
 981         // s20 = 0;\r
 982 \r
 983         s7 += s19 * 666643;\r
 984         s8 += s19 * 470296;\r
 985         s9 += s19 * 654183;\r
 986         s10 -= s19 * 997805;\r
 987         s11 += s19 * 136657;\r
 988         s12 -= s19 * 683901;\r
 989         // s19 = 0;\r
 990 \r
 991         s6 += s18 * 666643;\r
 992         s7 += s18 * 470296;\r
 993         s8 += s18 * 654183;\r
 994         s9 -= s18 * 997805;\r
 995         s10 += s18 * 136657;\r
 996         s11 -= s18 * 683901;\r
 997         // s18 = 0;\r
 998 \r
 999         // Reduce the bit length of limbs from s6 to s15 to 21-bits.\r
1000         carry6 = (s6 + (1 << 20)) >> 21;\r
1001         s7 += carry6;\r
1002         s6 -= carry6 << 21;\r
1003         carry8 = (s8 + (1 << 20)) >> 21;\r
1004         s9 += carry8;\r
1005         s8 -= carry8 << 21;\r
1006         carry10 = (s10 + (1 << 20)) >> 21;\r
1007         s11 += carry10;\r
1008         s10 -= carry10 << 21;\r
1009         carry12 = (s12 + (1 << 20)) >> 21;\r
1010         s13 += carry12;\r
1011         s12 -= carry12 << 21;\r
1012         carry14 = (s14 + (1 << 20)) >> 21;\r
1013         s15 += carry14;\r
1014         s14 -= carry14 << 21;\r
1015         carry16 = (s16 + (1 << 20)) >> 21;\r
1016         s17 += carry16;\r
1017         s16 -= carry16 << 21;\r
1018 \r
1019         carry7 = (s7 + (1 << 20)) >> 21;\r
1020         s8 += carry7;\r
1021         s7 -= carry7 << 21;\r
1022         carry9 = (s9 + (1 << 20)) >> 21;\r
1023         s10 += carry9;\r
1024         s9 -= carry9 << 21;\r
1025         carry11 = (s11 + (1 << 20)) >> 21;\r
1026         s12 += carry11;\r
1027         s11 -= carry11 << 21;\r
1028         carry13 = (s13 + (1 << 20)) >> 21;\r
1029         s14 += carry13;\r
1030         s13 -= carry13 << 21;\r
1031         carry15 = (s15 + (1 << 20)) >> 21;\r
1032         s16 += carry15;\r
1033         s15 -= carry15 << 21;\r
1034 \r
1035         // Resume reduction where we left off.\r
1036         s5 += s17 * 666643;\r
1037         s6 += s17 * 470296;\r
1038         s7 += s17 * 654183;\r
1039         s8 -= s17 * 997805;\r
1040         s9 += s17 * 136657;\r
1041         s10 -= s17 * 683901;\r
1042         // s17 = 0;\r
1043 \r
1044         s4 += s16 * 666643;\r
1045         s5 += s16 * 470296;\r
1046         s6 += s16 * 654183;\r
1047         s7 -= s16 * 997805;\r
1048         s8 += s16 * 136657;\r
1049         s9 -= s16 * 683901;\r
1050         // s16 = 0;\r
1051 \r
1052         s3 += s15 * 666643;\r
1053         s4 += s15 * 470296;\r
1054         s5 += s15 * 654183;\r
1055         s6 -= s15 * 997805;\r
1056         s7 += s15 * 136657;\r
1057         s8 -= s15 * 683901;\r
1058         // s15 = 0;\r
1059 \r
1060         s2 += s14 * 666643;\r
1061         s3 += s14 * 470296;\r
1062         s4 += s14 * 654183;\r
1063         s5 -= s14 * 997805;\r
1064         s6 += s14 * 136657;\r
1065         s7 -= s14 * 683901;\r
1066         // s14 = 0;\r
1067 \r
1068         s1 += s13 * 666643;\r
1069         s2 += s13 * 470296;\r
1070         s3 += s13 * 654183;\r
1071         s4 -= s13 * 997805;\r
1072         s5 += s13 * 136657;\r
1073         s6 -= s13 * 683901;\r
1074         // s13 = 0;\r
1075 \r
1076         s0 += s12 * 666643;\r
1077         s1 += s12 * 470296;\r
1078         s2 += s12 * 654183;\r
1079         s3 -= s12 * 997805;\r
1080         s4 += s12 * 136657;\r
1081         s5 -= s12 * 683901;\r
1082         s12 = 0;\r
1083 \r
1084         // Reduce the range of limbs from s0 to s11 to 21-bits.\r
1085         carry0 = (s0 + (1 << 20)) >> 21;\r
1086         s1 += carry0;\r
1087         s0 -= carry0 << 21;\r
1088         carry2 = (s2 + (1 << 20)) >> 21;\r
1089         s3 += carry2;\r
1090         s2 -= carry2 << 21;\r
1091         carry4 = (s4 + (1 << 20)) >> 21;\r
1092         s5 += carry4;\r
1093         s4 -= carry4 << 21;\r
1094         carry6 = (s6 + (1 << 20)) >> 21;\r
1095         s7 += carry6;\r
1096         s6 -= carry6 << 21;\r
1097         carry8 = (s8 + (1 << 20)) >> 21;\r
1098         s9 += carry8;\r
1099         s8 -= carry8 << 21;\r
1100         carry10 = (s10 + (1 << 20)) >> 21;\r
1101         s11 += carry10;\r
1102         s10 -= carry10 << 21;\r
1103 \r
1104         carry1 = (s1 + (1 << 20)) >> 21;\r
1105         s2 += carry1;\r
1106         s1 -= carry1 << 21;\r
1107         carry3 = (s3 + (1 << 20)) >> 21;\r
1108         s4 += carry3;\r
1109         s3 -= carry3 << 21;\r
1110         carry5 = (s5 + (1 << 20)) >> 21;\r
1111         s6 += carry5;\r
1112         s5 -= carry5 << 21;\r
1113         carry7 = (s7 + (1 << 20)) >> 21;\r
1114         s8 += carry7;\r
1115         s7 -= carry7 << 21;\r
1116         carry9 = (s9 + (1 << 20)) >> 21;\r
1117         s10 += carry9;\r
1118         s9 -= carry9 << 21;\r
1119         carry11 = (s11 + (1 << 20)) >> 21;\r
1120         s12 += carry11;\r
1121         s11 -= carry11 << 21;\r
1122 \r
1123         s0 += s12 * 666643;\r
1124         s1 += s12 * 470296;\r
1125         s2 += s12 * 654183;\r
1126         s3 -= s12 * 997805;\r
1127         s4 += s12 * 136657;\r
1128         s5 -= s12 * 683901;\r
1129         s12 = 0;\r
1130 \r
1131         // Carry chain reduction to propagate excess bits from s0 to s5 to the most significant limbs.\r
1132         carry0 = s0 >> 21;\r
1133         s1 += carry0;\r
1134         s0 -= carry0 << 21;\r
1135         carry1 = s1 >> 21;\r
1136         s2 += carry1;\r
1137         s1 -= carry1 << 21;\r
1138         carry2 = s2 >> 21;\r
1139         s3 += carry2;\r
1140         s2 -= carry2 << 21;\r
1141         carry3 = s3 >> 21;\r
1142         s4 += carry3;\r
1143         s3 -= carry3 << 21;\r
1144         carry4 = s4 >> 21;\r
1145         s5 += carry4;\r
1146         s4 -= carry4 << 21;\r
1147         carry5 = s5 >> 21;\r
1148         s6 += carry5;\r
1149         s5 -= carry5 << 21;\r
1150         carry6 = s6 >> 21;\r
1151         s7 += carry6;\r
1152         s6 -= carry6 << 21;\r
1153         carry7 = s7 >> 21;\r
1154         s8 += carry7;\r
1155         s7 -= carry7 << 21;\r
1156         carry8 = s8 >> 21;\r
1157         s9 += carry8;\r
1158         s8 -= carry8 << 21;\r
1159         carry9 = s9 >> 21;\r
1160         s10 += carry9;\r
1161         s9 -= carry9 << 21;\r
1162         carry10 = s10 >> 21;\r
1163         s11 += carry10;\r
1164         s10 -= carry10 << 21;\r
1165         carry11 = s11 >> 21;\r
1166         s12 += carry11;\r
1167         s11 -= carry11 << 21;\r
1168 \r
1169         // Do one last reduction as s12 might be 1.\r
1170         s0 += s12 * 666643;\r
1171         s1 += s12 * 470296;\r
1172         s2 += s12 * 654183;\r
1173         s3 -= s12 * 997805;\r
1174         s4 += s12 * 136657;\r
1175         s5 -= s12 * 683901;\r
1176         // s12 = 0;\r
1177 \r
1178         carry0 = s0 >> 21;\r
1179         s1 += carry0;\r
1180         s0 -= carry0 << 21;\r
1181         carry1 = s1 >> 21;\r
1182         s2 += carry1;\r
1183         s1 -= carry1 << 21;\r
1184         carry2 = s2 >> 21;\r
1185         s3 += carry2;\r
1186         s2 -= carry2 << 21;\r
1187         carry3 = s3 >> 21;\r
1188         s4 += carry3;\r
1189         s3 -= carry3 << 21;\r
1190         carry4 = s4 >> 21;\r
1191         s5 += carry4;\r
1192         s4 -= carry4 << 21;\r
1193         carry5 = s5 >> 21;\r
1194         s6 += carry5;\r
1195         s5 -= carry5 << 21;\r
1196         carry6 = s6 >> 21;\r
1197         s7 += carry6;\r
1198         s6 -= carry6 << 21;\r
1199         carry7 = s7 >> 21;\r
1200         s8 += carry7;\r
1201         s7 -= carry7 << 21;\r
1202         carry8 = s8 >> 21;\r
1203         s9 += carry8;\r
1204         s8 -= carry8 << 21;\r
1205         carry9 = s9 >> 21;\r
1206         s10 += carry9;\r
1207         s9 -= carry9 << 21;\r
1208         carry10 = s10 >> 21;\r
1209         s11 += carry10;\r
1210         s10 -= carry10 << 21;\r
1211 \r
1212         // Serialize the result into the s.\r
1213         s[0] = (byte) s0;\r
1214         s[1] = (byte) (s0 >> 8);\r
1215         s[2] = (byte) ((s0 >> 16) | (s1 << 5));\r
1216         s[3] = (byte) (s1 >> 3);\r
1217         s[4] = (byte) (s1 >> 11);\r
1218         s[5] = (byte) ((s1 >> 19) | (s2 << 2));\r
1219         s[6] = (byte) (s2 >> 6);\r
1220         s[7] = (byte) ((s2 >> 14) | (s3 << 7));\r
1221         s[8] = (byte) (s3 >> 1);\r
1222         s[9] = (byte) (s3 >> 9);\r
1223         s[10] = (byte) ((s3 >> 17) | (s4 << 4));\r
1224         s[11] = (byte) (s4 >> 4);\r
1225         s[12] = (byte) (s4 >> 12);\r
1226         s[13] = (byte) ((s4 >> 20) | (s5 << 1));\r
1227         s[14] = (byte) (s5 >> 7);\r
1228         s[15] = (byte) ((s5 >> 15) | (s6 << 6));\r
1229         s[16] = (byte) (s6 >> 2);\r
1230         s[17] = (byte) (s6 >> 10);\r
1231         s[18] = (byte) ((s6 >> 18) | (s7 << 3));\r
1232         s[19] = (byte) (s7 >> 5);\r
1233         s[20] = (byte) (s7 >> 13);\r
1234         s[21] = (byte) s8;\r
1235         s[22] = (byte) (s8 >> 8);\r
1236         s[23] = (byte) ((s8 >> 16) | (s9 << 5));\r
1237         s[24] = (byte) (s9 >> 3);\r
1238         s[25] = (byte) (s9 >> 11);\r
1239         s[26] = (byte) ((s9 >> 19) | (s10 << 2));\r
1240         s[27] = (byte) (s10 >> 6);\r
1241         s[28] = (byte) ((s10 >> 14) | (s11 << 7));\r
1242         s[29] = (byte) (s11 >> 1);\r
1243         s[30] = (byte) (s11 >> 9);\r
1244         s[31] = (byte) (s11 >> 17);\r
1245     }\r
1246 \r
1247     /**\r
1248      * Input:\r
1249      * a[0]+256*a[1]+...+256^31*a[31] = a\r
1250      * b[0]+256*b[1]+...+256^31*b[31] = b\r
1251      * c[0]+256*c[1]+...+256^31*c[31] = c\r
1252      * <p>\r
1253      * Output:\r
1254      * s[0]+256*s[1]+...+256^31*s[31] = (ab+c) mod l\r
1255      * where l = 2^252 + 27742317777372353535851937790883648493.\r
1256      */\r
1257     public static void mulAdd(byte[] s, byte[] a, byte[] b, byte[] c) {\r
1258         // This is very similar to Ed25519.reduce, the difference in here is that it computes ab+c\r
1259         // See Ed25519.reduce for related comments.\r
1260         long a0 = 2097151 & load3(a, 0);\r
1261         long a1 = 2097151 & (load4(a, 2) >> 5);\r
1262         long a2 = 2097151 & (load3(a, 5) >> 2);\r
1263         long a3 = 2097151 & (load4(a, 7) >> 7);\r
1264         long a4 = 2097151 & (load4(a, 10) >> 4);\r
1265         long a5 = 2097151 & (load3(a, 13) >> 1);\r
1266         long a6 = 2097151 & (load4(a, 15) >> 6);\r
1267         long a7 = 2097151 & (load3(a, 18) >> 3);\r
1268         long a8 = 2097151 & load3(a, 21);\r
1269         long a9 = 2097151 & (load4(a, 23) >> 5);\r
1270         long a10 = 2097151 & (load3(a, 26) >> 2);\r
1271         long a11 = (load4(a, 28) >> 7);\r
1272         long b0 = 2097151 & load3(b, 0);\r
1273         long b1 = 2097151 & (load4(b, 2) >> 5);\r
1274         long b2 = 2097151 & (load3(b, 5) >> 2);\r
1275         long b3 = 2097151 & (load4(b, 7) >> 7);\r
1276         long b4 = 2097151 & (load4(b, 10) >> 4);\r
1277         long b5 = 2097151 & (load3(b, 13) >> 1);\r
1278         long b6 = 2097151 & (load4(b, 15) >> 6);\r
1279         long b7 = 2097151 & (load3(b, 18) >> 3);\r
1280         long b8 = 2097151 & load3(b, 21);\r
1281         long b9 = 2097151 & (load4(b, 23) >> 5);\r
1282         long b10 = 2097151 & (load3(b, 26) >> 2);\r
1283         long b11 = (load4(b, 28) >> 7);\r
1284         long c0 = 2097151 & load3(c, 0);\r
1285         long c1 = 2097151 & (load4(c, 2) >> 5);\r
1286         long c2 = 2097151 & (load3(c, 5) >> 2);\r
1287         long c3 = 2097151 & (load4(c, 7) >> 7);\r
1288         long c4 = 2097151 & (load4(c, 10) >> 4);\r
1289         long c5 = 2097151 & (load3(c, 13) >> 1);\r
1290         long c6 = 2097151 & (load4(c, 15) >> 6);\r
1291         long c7 = 2097151 & (load3(c, 18) >> 3);\r
1292         long c8 = 2097151 & load3(c, 21);\r
1293         long c9 = 2097151 & (load4(c, 23) >> 5);\r
1294         long c10 = 2097151 & (load3(c, 26) >> 2);\r
1295         long c11 = (load4(c, 28) >> 7);\r
1296         long s0;\r
1297         long s1;\r
1298         long s2;\r
1299         long s3;\r
1300         long s4;\r
1301         long s5;\r
1302         long s6;\r
1303         long s7;\r
1304         long s8;\r
1305         long s9;\r
1306         long s10;\r
1307         long s11;\r
1308         long s12;\r
1309         long s13;\r
1310         long s14;\r
1311         long s15;\r
1312         long s16;\r
1313         long s17;\r
1314         long s18;\r
1315         long s19;\r
1316         long s20;\r
1317         long s21;\r
1318         long s22;\r
1319         long s23;\r
1320         long carry0;\r
1321         long carry1;\r
1322         long carry2;\r
1323         long carry3;\r
1324         long carry4;\r
1325         long carry5;\r
1326         long carry6;\r
1327         long carry7;\r
1328         long carry8;\r
1329         long carry9;\r
1330         long carry10;\r
1331         long carry11;\r
1332         long carry12;\r
1333         long carry13;\r
1334         long carry14;\r
1335         long carry15;\r
1336         long carry16;\r
1337         long carry17;\r
1338         long carry18;\r
1339         long carry19;\r
1340         long carry20;\r
1341         long carry21;\r
1342         long carry22;\r
1343 \r
1344         s0 = c0 + a0 * b0;\r
1345         s1 = c1 + a0 * b1 + a1 * b0;\r
1346         s2 = c2 + a0 * b2 + a1 * b1 + a2 * b0;\r
1347         s3 = c3 + a0 * b3 + a1 * b2 + a2 * b1 + a3 * b0;\r
1348         s4 = c4 + a0 * b4 + a1 * b3 + a2 * b2 + a3 * b1 + a4 * b0;\r
1349         s5 = c5 + a0 * b5 + a1 * b4 + a2 * b3 + a3 * b2 + a4 * b1 + a5 * b0;\r
1350         s6 = c6 + a0 * b6 + a1 * b5 + a2 * b4 + a3 * b3 + a4 * b2 + a5 * b1 + a6 * b0;\r
1351         s7 = c7 + a0 * b7 + a1 * b6 + a2 * b5 + a3 * b4 + a4 * b3 + a5 * b2 + a6 * b1 + a7 * b0;\r
1352         s8 = c8 + a0 * b8 + a1 * b7 + a2 * b6 + a3 * b5 + a4 * b4 + a5 * b3 + a6 * b2 + a7 * b1\r
1353                 + a8 * b0;\r
1354         s9 = c9 + a0 * b9 + a1 * b8 + a2 * b7 + a3 * b6 + a4 * b5 + a5 * b4 + a6 * b3 + a7 * b2\r
1355                 + a8 * b1 + a9 * b0;\r
1356         s10 = c10 + a0 * b10 + a1 * b9 + a2 * b8 + a3 * b7 + a4 * b6 + a5 * b5 + a6 * b4 + a7 * b3\r
1357                 + a8 * b2 + a9 * b1 + a10 * b0;\r
1358         s11 = c11 + a0 * b11 + a1 * b10 + a2 * b9 + a3 * b8 + a4 * b7 + a5 * b6 + a6 * b5 + a7 * b4\r
1359                 + a8 * b3 + a9 * b2 + a10 * b1 + a11 * b0;\r
1360         s12 = a1 * b11 + a2 * b10 + a3 * b9 + a4 * b8 + a5 * b7 + a6 * b6 + a7 * b5 + a8 * b4 + a9 * b3\r
1361                 + a10 * b2 + a11 * b1;\r
1362         s13 = a2 * b11 + a3 * b10 + a4 * b9 + a5 * b8 + a6 * b7 + a7 * b6 + a8 * b5 + a9 * b4 + a10 * b3\r
1363                 + a11 * b2;\r
1364         s14 = a3 * b11 + a4 * b10 + a5 * b9 + a6 * b8 + a7 * b7 + a8 * b6 + a9 * b5 + a10 * b4\r
1365                 + a11 * b3;\r
1366         s15 = a4 * b11 + a5 * b10 + a6 * b9 + a7 * b8 + a8 * b7 + a9 * b6 + a10 * b5 + a11 * b4;\r
1367         s16 = a5 * b11 + a6 * b10 + a7 * b9 + a8 * b8 + a9 * b7 + a10 * b6 + a11 * b5;\r
1368         s17 = a6 * b11 + a7 * b10 + a8 * b9 + a9 * b8 + a10 * b7 + a11 * b6;\r
1369         s18 = a7 * b11 + a8 * b10 + a9 * b9 + a10 * b8 + a11 * b7;\r
1370         s19 = a8 * b11 + a9 * b10 + a10 * b9 + a11 * b8;\r
1371         s20 = a9 * b11 + a10 * b10 + a11 * b9;\r
1372         s21 = a10 * b11 + a11 * b10;\r
1373         s22 = a11 * b11;\r
1374         s23 = 0;\r
1375 \r
1376         carry0 = (s0 + (1 << 20)) >> 21;\r
1377         s1 += carry0;\r
1378         s0 -= carry0 << 21;\r
1379         carry2 = (s2 + (1 << 20)) >> 21;\r
1380         s3 += carry2;\r
1381         s2 -= carry2 << 21;\r
1382         carry4 = (s4 + (1 << 20)) >> 21;\r
1383         s5 += carry4;\r
1384         s4 -= carry4 << 21;\r
1385         carry6 = (s6 + (1 << 20)) >> 21;\r
1386         s7 += carry6;\r
1387         s6 -= carry6 << 21;\r
1388         carry8 = (s8 + (1 << 20)) >> 21;\r
1389         s9 += carry8;\r
1390         s8 -= carry8 << 21;\r
1391         carry10 = (s10 + (1 << 20)) >> 21;\r
1392         s11 += carry10;\r
1393         s10 -= carry10 << 21;\r
1394         carry12 = (s12 + (1 << 20)) >> 21;\r
1395         s13 += carry12;\r
1396         s12 -= carry12 << 21;\r
1397         carry14 = (s14 + (1 << 20)) >> 21;\r
1398         s15 += carry14;\r
1399         s14 -= carry14 << 21;\r
1400         carry16 = (s16 + (1 << 20)) >> 21;\r
1401         s17 += carry16;\r
1402         s16 -= carry16 << 21;\r
1403         carry18 = (s18 + (1 << 20)) >> 21;\r
1404         s19 += carry18;\r
1405         s18 -= carry18 << 21;\r
1406         carry20 = (s20 + (1 << 20)) >> 21;\r
1407         s21 += carry20;\r
1408         s20 -= carry20 << 21;\r
1409         carry22 = (s22 + (1 << 20)) >> 21;\r
1410         s23 += carry22;\r
1411         s22 -= carry22 << 21;\r
1412 \r
1413         carry1 = (s1 + (1 << 20)) >> 21;\r
1414         s2 += carry1;\r
1415         s1 -= carry1 << 21;\r
1416         carry3 = (s3 + (1 << 20)) >> 21;\r
1417         s4 += carry3;\r
1418         s3 -= carry3 << 21;\r
1419         carry5 = (s5 + (1 << 20)) >> 21;\r
1420         s6 += carry5;\r
1421         s5 -= carry5 << 21;\r
1422         carry7 = (s7 + (1 << 20)) >> 21;\r
1423         s8 += carry7;\r
1424         s7 -= carry7 << 21;\r
1425         carry9 = (s9 + (1 << 20)) >> 21;\r
1426         s10 += carry9;\r
1427         s9 -= carry9 << 21;\r
1428         carry11 = (s11 + (1 << 20)) >> 21;\r
1429         s12 += carry11;\r
1430         s11 -= carry11 << 21;\r
1431         carry13 = (s13 + (1 << 20)) >> 21;\r
1432         s14 += carry13;\r
1433         s13 -= carry13 << 21;\r
1434         carry15 = (s15 + (1 << 20)) >> 21;\r
1435         s16 += carry15;\r
1436         s15 -= carry15 << 21;\r
1437         carry17 = (s17 + (1 << 20)) >> 21;\r
1438         s18 += carry17;\r
1439         s17 -= carry17 << 21;\r
1440         carry19 = (s19 + (1 << 20)) >> 21;\r
1441         s20 += carry19;\r
1442         s19 -= carry19 << 21;\r
1443         carry21 = (s21 + (1 << 20)) >> 21;\r
1444         s22 += carry21;\r
1445         s21 -= carry21 << 21;\r
1446 \r
1447         s11 += s23 * 666643;\r
1448         s12 += s23 * 470296;\r
1449         s13 += s23 * 654183;\r
1450         s14 -= s23 * 997805;\r
1451         s15 += s23 * 136657;\r
1452         s16 -= s23 * 683901;\r
1453         // s23 = 0;\r
1454 \r
1455         s10 += s22 * 666643;\r
1456         s11 += s22 * 470296;\r
1457         s12 += s22 * 654183;\r
1458         s13 -= s22 * 997805;\r
1459         s14 += s22 * 136657;\r
1460         s15 -= s22 * 683901;\r
1461         // s22 = 0;\r
1462 \r
1463         s9 += s21 * 666643;\r
1464         s10 += s21 * 470296;\r
1465         s11 += s21 * 654183;\r
1466         s12 -= s21 * 997805;\r
1467         s13 += s21 * 136657;\r
1468         s14 -= s21 * 683901;\r
1469         // s21 = 0;\r
1470 \r
1471         s8 += s20 * 666643;\r
1472         s9 += s20 * 470296;\r
1473         s10 += s20 * 654183;\r
1474         s11 -= s20 * 997805;\r
1475         s12 += s20 * 136657;\r
1476         s13 -= s20 * 683901;\r
1477         // s20 = 0;\r
1478 \r
1479         s7 += s19 * 666643;\r
1480         s8 += s19 * 470296;\r
1481         s9 += s19 * 654183;\r
1482         s10 -= s19 * 997805;\r
1483         s11 += s19 * 136657;\r
1484         s12 -= s19 * 683901;\r
1485         // s19 = 0;\r
1486 \r
1487         s6 += s18 * 666643;\r
1488         s7 += s18 * 470296;\r
1489         s8 += s18 * 654183;\r
1490         s9 -= s18 * 997805;\r
1491         s10 += s18 * 136657;\r
1492         s11 -= s18 * 683901;\r
1493         // s18 = 0;\r
1494 \r
1495         carry6 = (s6 + (1 << 20)) >> 21;\r
1496         s7 += carry6;\r
1497         s6 -= carry6 << 21;\r
1498         carry8 = (s8 + (1 << 20)) >> 21;\r
1499         s9 += carry8;\r
1500         s8 -= carry8 << 21;\r
1501         carry10 = (s10 + (1 << 20)) >> 21;\r
1502         s11 += carry10;\r
1503         s10 -= carry10 << 21;\r
1504         carry12 = (s12 + (1 << 20)) >> 21;\r
1505         s13 += carry12;\r
1506         s12 -= carry12 << 21;\r
1507         carry14 = (s14 + (1 << 20)) >> 21;\r
1508         s15 += carry14;\r
1509         s14 -= carry14 << 21;\r
1510         carry16 = (s16 + (1 << 20)) >> 21;\r
1511         s17 += carry16;\r
1512         s16 -= carry16 << 21;\r
1513 \r
1514         carry7 = (s7 + (1 << 20)) >> 21;\r
1515         s8 += carry7;\r
1516         s7 -= carry7 << 21;\r
1517         carry9 = (s9 + (1 << 20)) >> 21;\r
1518         s10 += carry9;\r
1519         s9 -= carry9 << 21;\r
1520         carry11 = (s11 + (1 << 20)) >> 21;\r
1521         s12 += carry11;\r
1522         s11 -= carry11 << 21;\r
1523         carry13 = (s13 + (1 << 20)) >> 21;\r
1524         s14 += carry13;\r
1525         s13 -= carry13 << 21;\r
1526         carry15 = (s15 + (1 << 20)) >> 21;\r
1527         s16 += carry15;\r
1528         s15 -= carry15 << 21;\r
1529 \r
1530         s5 += s17 * 666643;\r
1531         s6 += s17 * 470296;\r
1532         s7 += s17 * 654183;\r
1533         s8 -= s17 * 997805;\r
1534         s9 += s17 * 136657;\r
1535         s10 -= s17 * 683901;\r
1536         // s17 = 0;\r
1537 \r
1538         s4 += s16 * 666643;\r
1539         s5 += s16 * 470296;\r
1540         s6 += s16 * 654183;\r
1541         s7 -= s16 * 997805;\r
1542         s8 += s16 * 136657;\r
1543         s9 -= s16 * 683901;\r
1544         // s16 = 0;\r
1545 \r
1546         s3 += s15 * 666643;\r
1547         s4 += s15 * 470296;\r
1548         s5 += s15 * 654183;\r
1549         s6 -= s15 * 997805;\r
1550         s7 += s15 * 136657;\r
1551         s8 -= s15 * 683901;\r
1552         // s15 = 0;\r
1553 \r
1554         s2 += s14 * 666643;\r
1555         s3 += s14 * 470296;\r
1556         s4 += s14 * 654183;\r
1557         s5 -= s14 * 997805;\r
1558         s6 += s14 * 136657;\r
1559         s7 -= s14 * 683901;\r
1560         // s14 = 0;\r
1561 \r
1562         s1 += s13 * 666643;\r
1563         s2 += s13 * 470296;\r
1564         s3 += s13 * 654183;\r
1565         s4 -= s13 * 997805;\r
1566         s5 += s13 * 136657;\r
1567         s6 -= s13 * 683901;\r
1568         // s13 = 0;\r
1569 \r
1570         s0 += s12 * 666643;\r
1571         s1 += s12 * 470296;\r
1572         s2 += s12 * 654183;\r
1573         s3 -= s12 * 997805;\r
1574         s4 += s12 * 136657;\r
1575         s5 -= s12 * 683901;\r
1576         s12 = 0;\r
1577 \r
1578         carry0 = (s0 + (1 << 20)) >> 21;\r
1579         s1 += carry0;\r
1580         s0 -= carry0 << 21;\r
1581         carry2 = (s2 + (1 << 20)) >> 21;\r
1582         s3 += carry2;\r
1583         s2 -= carry2 << 21;\r
1584         carry4 = (s4 + (1 << 20)) >> 21;\r
1585         s5 += carry4;\r
1586         s4 -= carry4 << 21;\r
1587         carry6 = (s6 + (1 << 20)) >> 21;\r
1588         s7 += carry6;\r
1589         s6 -= carry6 << 21;\r
1590         carry8 = (s8 + (1 << 20)) >> 21;\r
1591         s9 += carry8;\r
1592         s8 -= carry8 << 21;\r
1593         carry10 = (s10 + (1 << 20)) >> 21;\r
1594         s11 += carry10;\r
1595         s10 -= carry10 << 21;\r
1596 \r
1597         carry1 = (s1 + (1 << 20)) >> 21;\r
1598         s2 += carry1;\r
1599         s1 -= carry1 << 21;\r
1600         carry3 = (s3 + (1 << 20)) >> 21;\r
1601         s4 += carry3;\r
1602         s3 -= carry3 << 21;\r
1603         carry5 = (s5 + (1 << 20)) >> 21;\r
1604         s6 += carry5;\r
1605         s5 -= carry5 << 21;\r
1606         carry7 = (s7 + (1 << 20)) >> 21;\r
1607         s8 += carry7;\r
1608         s7 -= carry7 << 21;\r
1609         carry9 = (s9 + (1 << 20)) >> 21;\r
1610         s10 += carry9;\r
1611         s9 -= carry9 << 21;\r
1612         carry11 = (s11 + (1 << 20)) >> 21;\r
1613         s12 += carry11;\r
1614         s11 -= carry11 << 21;\r
1615 \r
1616         s0 += s12 * 666643;\r
1617         s1 += s12 * 470296;\r
1618         s2 += s12 * 654183;\r
1619         s3 -= s12 * 997805;\r
1620         s4 += s12 * 136657;\r
1621         s5 -= s12 * 683901;\r
1622         s12 = 0;\r
1623 \r
1624         carry0 = s0 >> 21;\r
1625         s1 += carry0;\r
1626         s0 -= carry0 << 21;\r
1627         carry1 = s1 >> 21;\r
1628         s2 += carry1;\r
1629         s1 -= carry1 << 21;\r
1630         carry2 = s2 >> 21;\r
1631         s3 += carry2;\r
1632         s2 -= carry2 << 21;\r
1633         carry3 = s3 >> 21;\r
1634         s4 += carry3;\r
1635         s3 -= carry3 << 21;\r
1636         carry4 = s4 >> 21;\r
1637         s5 += carry4;\r
1638         s4 -= carry4 << 21;\r
1639         carry5 = s5 >> 21;\r
1640         s6 += carry5;\r
1641         s5 -= carry5 << 21;\r
1642         carry6 = s6 >> 21;\r
1643         s7 += carry6;\r
1644         s6 -= carry6 << 21;\r
1645         carry7 = s7 >> 21;\r
1646         s8 += carry7;\r
1647         s7 -= carry7 << 21;\r
1648         carry8 = s8 >> 21;\r
1649         s9 += carry8;\r
1650         s8 -= carry8 << 21;\r
1651         carry9 = s9 >> 21;\r
1652         s10 += carry9;\r
1653         s9 -= carry9 << 21;\r
1654         carry10 = s10 >> 21;\r
1655         s11 += carry10;\r
1656         s10 -= carry10 << 21;\r
1657         carry11 = s11 >> 21;\r
1658         s12 += carry11;\r
1659         s11 -= carry11 << 21;\r
1660 \r
1661         s0 += s12 * 666643;\r
1662         s1 += s12 * 470296;\r
1663         s2 += s12 * 654183;\r
1664         s3 -= s12 * 997805;\r
1665         s4 += s12 * 136657;\r
1666         s5 -= s12 * 683901;\r
1667         // s12 = 0;\r
1668 \r
1669         carry0 = s0 >> 21;\r
1670         s1 += carry0;\r
1671         s0 -= carry0 << 21;\r
1672         carry1 = s1 >> 21;\r
1673         s2 += carry1;\r
1674         s1 -= carry1 << 21;\r
1675         carry2 = s2 >> 21;\r
1676         s3 += carry2;\r
1677         s2 -= carry2 << 21;\r
1678         carry3 = s3 >> 21;\r
1679         s4 += carry3;\r
1680         s3 -= carry3 << 21;\r
1681         carry4 = s4 >> 21;\r
1682         s5 += carry4;\r
1683         s4 -= carry4 << 21;\r
1684         carry5 = s5 >> 21;\r
1685         s6 += carry5;\r
1686         s5 -= carry5 << 21;\r
1687         carry6 = s6 >> 21;\r
1688         s7 += carry6;\r
1689         s6 -= carry6 << 21;\r
1690         carry7 = s7 >> 21;\r
1691         s8 += carry7;\r
1692         s7 -= carry7 << 21;\r
1693         carry8 = s8 >> 21;\r
1694         s9 += carry8;\r
1695         s8 -= carry8 << 21;\r
1696         carry9 = s9 >> 21;\r
1697         s10 += carry9;\r
1698         s9 -= carry9 << 21;\r
1699         carry10 = s10 >> 21;\r
1700         s11 += carry10;\r
1701         s10 -= carry10 << 21;\r
1702 \r
1703         s[0] = (byte) s0;\r
1704         s[1] = (byte) (s0 >> 8);\r
1705         s[2] = (byte) ((s0 >> 16) | (s1 << 5));\r
1706         s[3] = (byte) (s1 >> 3);\r
1707         s[4] = (byte) (s1 >> 11);\r
1708         s[5] = (byte) ((s1 >> 19) | (s2 << 2));\r
1709         s[6] = (byte) (s2 >> 6);\r
1710         s[7] = (byte) ((s2 >> 14) | (s3 << 7));\r
1711         s[8] = (byte) (s3 >> 1);\r
1712         s[9] = (byte) (s3 >> 9);\r
1713         s[10] = (byte) ((s3 >> 17) | (s4 << 4));\r
1714         s[11] = (byte) (s4 >> 4);\r
1715         s[12] = (byte) (s4 >> 12);\r
1716         s[13] = (byte) ((s4 >> 20) | (s5 << 1));\r
1717         s[14] = (byte) (s5 >> 7);\r
1718         s[15] = (byte) ((s5 >> 15) | (s6 << 6));\r
1719         s[16] = (byte) (s6 >> 2);\r
1720         s[17] = (byte) (s6 >> 10);\r
1721         s[18] = (byte) ((s6 >> 18) | (s7 << 3));\r
1722         s[19] = (byte) (s7 >> 5);\r
1723         s[20] = (byte) (s7 >> 13);\r
1724         s[21] = (byte) s8;\r
1725         s[22] = (byte) (s8 >> 8);\r
1726         s[23] = (byte) ((s8 >> 16) | (s9 << 5));\r
1727         s[24] = (byte) (s9 >> 3);\r
1728         s[25] = (byte) (s9 >> 11);\r
1729         s[26] = (byte) ((s9 >> 19) | (s10 << 2));\r
1730         s[27] = (byte) (s10 >> 6);\r
1731         s[28] = (byte) ((s10 >> 14) | (s11 << 7));\r
1732         s[29] = (byte) (s11 >> 1);\r
1733         s[30] = (byte) (s11 >> 9);\r
1734         s[31] = (byte) (s11 >> 17);\r
1735     }\r
1736 \r
1737     static byte[] getHashedScalar(final byte[] privateKey)\r
1738             throws GeneralSecurityException {\r
1739         MessageDigest digest = EngineFactory.MESSAGE_DIGEST.getInstance("SHA-512");\r
1740         digest.update(privateKey, 0, FIELD_LEN);\r
1741         byte[] h = digest.digest();\r
1742         // https://tools.ietf.org/html/rfc8032#section-5.1.2.\r
1743         // Clear the lowest three bits of the first octet.\r
1744         h[0] = (byte) (h[0] & 248);\r
1745         // Clear the highest bit of the last octet.\r
1746         h[31] = (byte) (h[31] & 127);\r
1747         // Set the second highest bit if the last octet.\r
1748         h[31] = (byte) (h[31] | 64);\r
1749         return h;\r
1750     }\r
1751 \r
1752     /**\r
1753      * Returns the EdDSA signature for the {@code message} based on the {@code hashedPrivateKey}.\r
1754      *\r
1755      * @param message          to sign\r
1756      * @param hashedPrivateKey {@link Ed25519#getHashedScalar(byte[])} of the private key\r
1757      * @return signature for the {@code message}.\r
1758      * @throws GeneralSecurityException if there is no SHA-512 algorithm defined in\r
1759      *                                  {@link EngineFactory}.MESSAGE_DIGEST.\r
1760      */\r
1761     public static byte[] sign(final byte[] message, final byte[] publicKey, final byte[] hashedPrivateKey)\r
1762             throws GeneralSecurityException {\r
1763         // Copying the message to make it thread-safe. Otherwise, if the caller modifies the message\r
1764         // between the first and the second hash then it might leak the private key.\r
1765         byte[] messageCopy = Arrays.copyOfRange(message, 0, message.length);\r
1766         MessageDigest digest = EngineFactory.MESSAGE_DIGEST.getInstance("SHA-512");\r
1767         digest.update(hashedPrivateKey, FIELD_LEN, FIELD_LEN);\r
1768         digest.update(messageCopy);\r
1769         byte[] r = digest.digest();\r
1770         reduce(r);\r
1771 \r
1772         byte[] rB = Arrays.copyOfRange(scalarMultWithBase(r).toBytes(), 0, FIELD_LEN);\r
1773         digest.reset();\r
1774         digest.update(rB);\r
1775         digest.update(publicKey);\r
1776         digest.update(messageCopy);\r
1777         byte[] hram = digest.digest();\r
1778         reduce(hram);\r
1779         byte[] s = new byte[FIELD_LEN];\r
1780         mulAdd(s, hram, hashedPrivateKey, r);\r
1781         return Bytes.concat(rB, s);\r
1782     }\r
1783 \r
1784 \r
1785     // The order of the generator as unsigned bytes in little endian order.\r
1786     // (2^252 + 0x14def9dea2f79cd65812631a5cf5d3ed, cf. RFC 7748)\r
1787     static final byte[] GROUP_ORDER = new byte[]{\r
1788             (byte) 0xed, (byte) 0xd3, (byte) 0xf5, (byte) 0x5c,\r
1789             (byte) 0x1a, (byte) 0x63, (byte) 0x12, (byte) 0x58,\r
1790             (byte) 0xd6, (byte) 0x9c, (byte) 0xf7, (byte) 0xa2,\r
1791             (byte) 0xde, (byte) 0xf9, (byte) 0xde, (byte) 0x14,\r
1792             (byte) 0x00, (byte) 0x00, (byte) 0x00, (byte) 0x00,\r
1793             (byte) 0x00, (byte) 0x00, (byte) 0x00, (byte) 0x00,\r
1794             (byte) 0x00, (byte) 0x00, (byte) 0x00, (byte) 0x00,\r
1795             (byte) 0x00, (byte) 0x00, (byte) 0x00, (byte) 0x10};\r
1796 \r
1797     // Checks whether s represents an integer smaller than the order of the group.\r
1798     // This is needed to ensure that EdDSA signatures are non-malleable, as failing to check\r
1799     // the range of S allows to modify signatures (cf. RFC 8032, Section 5.2.7 and Section 8.4.)\r
1800     // @param s an integer in little-endian order.\r
1801     private static boolean isSmallerThanGroupOrder(byte[] s) {\r
1802         for (int j = FIELD_LEN - 1; j >= 0; j--) {\r
1803             // compare unsigned bytes\r
1804             int a = s[j] & 0xff;\r
1805             int b = GROUP_ORDER[j] & 0xff;\r
1806             if (a != b) {\r
1807                 return a < b;\r
1808             }\r
1809         }\r
1810         return false;\r
1811     }\r
1812 \r
1813     /**\r
1814      * Returns true if the EdDSA {@code signature} with {@code message}, can be verified with\r
1815      * {@code publicKey}.\r
1816      *\r
1817      * @throws GeneralSecurityException if there is no SHA-512 algorithm defined in\r
1818      *                                  {@link EngineFactory}.MESSAGE_DIGEST.\r
1819      */\r
1820     static boolean verify(final byte[] message, final byte[] signature,\r
1821                           final byte[] publicKey) throws GeneralSecurityException {\r
1822         if (signature.length != SIGNATURE_LEN) {\r
1823             return false;\r
1824         }\r
1825         byte[] s = Arrays.copyOfRange(signature, FIELD_LEN, SIGNATURE_LEN);\r
1826         if (!isSmallerThanGroupOrder(s)) {\r
1827             return false;\r
1828         }\r
1829         MessageDigest digest = EngineFactory.MESSAGE_DIGEST.getInstance("SHA-512");\r
1830         digest.update(signature, 0, FIELD_LEN);\r
1831         digest.update(publicKey);\r
1832         digest.update(message);\r
1833         byte[] h = digest.digest();\r
1834         reduce(h);\r
1835 \r
1836         XYZT negPublicKey = XYZT.fromBytesNegateVarTime(publicKey);\r
1837         XYZ xyz = doubleScalarMultVarTime(h, negPublicKey, s);\r
1838         byte[] expectedR = xyz.toBytes();\r
1839         for (int i = 0; i < FIELD_LEN; i++) {\r
1840             if (expectedR[i] != signature[i]) {\r
1841                 return false;\r
1842             }\r
1843         }\r
1844         return true;\r
1845     }\r
1846 }\r