_ Git - BoarSSL/blob - Crypto/MutableECPointPrime.cs

   1 /*
   2  * Copyright (c) 2017 Thomas Pornin <pornin@bolet.org>
   3  *
   4  * Permission is hereby granted, free of charge, to any person obtaining
   5  * a copy of this software and associated documentation files (the
   6  * "Software"), to deal in the Software without restriction, including
   7  * without limitation the rights to use, copy, modify, merge, publish,
   8  * distribute, sublicense, and/or sell copies of the Software, and to
   9  * permit persons to whom the Software is furnished to do so, subject to
  10  * the following conditions:
  11  *
  12  * The above copyright notice and this permission notice shall be
  13  * included in all copies or substantial portions of the Software.
  14  *
  15  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  16  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  17  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  18  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  19  * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  20  * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  21  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  22  * SOFTWARE.
  23  */
  24
  25 using System;
  26
  27 namespace Crypto {
  28
  29 /*
  30  * An implementation of MutableECPoint for curves in a prime field.
  31  * Jacobian coordinates are used.
  32  */
  33
  34 internal class MutableECPointPrime : MutableECPoint {
  35
  36         /*
  37          * Internal representation:
  38          *    x = X / (Z^2)
  39          *    y = Y / (Z^3)
  40          * For the point at infinity, Z == 0.
  41          *
  42          * When 'affine' is true (0xFFFFFFFF), then mx and my are in
  43          * normal representation, and mz is either 0 (point at infinity)
  44          * or 1 (otherwise). When 'affine' is false (0x00000000), then
  45          * mx, my and mz are in Montgomery representation.
  46          *
  47          * Note that in affine coordinates, the point at infinity admits
  48          * several equivalent representations. In non-affine
  49          * coordinates, all points have several equivalent
  50          * representations.
  51          */
  52
  53         ECCurvePrime curve;
  54         ModInt mx, my, mz;
  55         uint affine;
  56         ModInt mt1, mt2, mt3, mt4, mt5;
  57         ModInt ms1, ms2, ms3, ms4, ms5, ms6;
  58
  59         /*
  60          * Create a new instance. It is initialized to the point at
  61          * infinity.
  62          */
  63         internal MutableECPointPrime(ECCurvePrime curve)
  64         {
  65                 this.curve = curve;
  66                 mx = curve.mp.Dup();
  67                 my = curve.mp.Dup();
  68                 mz = curve.mp.Dup();
  69                 mt1 = curve.mp.Dup();
  70                 mt2 = curve.mp.Dup();
  71                 mt3 = curve.mp.Dup();
  72                 mt4 = curve.mp.Dup();
  73                 mt5 = curve.mp.Dup();
  74                 ms1 = curve.mp.Dup();
  75                 ms2 = curve.mp.Dup();
  76                 ms3 = curve.mp.Dup();
  77                 ms4 = curve.mp.Dup();
  78                 ms5 = curve.mp.Dup();
  79                 ms6 = curve.mp.Dup();
  80                 affine = 0xFFFFFFFF;
  81         }
  82
  83         internal override ECCurve Curve {
  84                 get {
  85                         return curve;
  86                 }
  87         }
  88
  89         internal override uint IsInfinityCT {
  90                 get {
  91                         return mz.IsZeroCT;
  92                 }
  93         }
  94
  95         internal override void Normalize()
  96         {
  97                 ToAffine();
  98         }
  99
 100         internal override byte[] Encode(bool compressed)
 101         {
 102                 ToAffine();
 103                 if (IsInfinity) {
 104                         return new byte[1];
 105                 }
 106                 if (compressed) {
 107                         byte[] enc = new byte[curve.EncodedLengthCompressed];
 108                         enc[0] = (byte)(0x02 + my.GetLSB());
 109                         mx.Encode(enc, 1, enc.Length - 1);
 110                         return enc;
 111                 } else {
 112                         byte[] enc = new byte[curve.EncodedLength];
 113                         int flen = (enc.Length - 1) >> 1;
 114                         enc[0] = 0x04;
 115                         mx.Encode(enc, 1, flen);
 116                         my.Encode(enc, 1 + flen, flen);
 117                         return enc;
 118                 }
 119         }
 120
 121         internal override uint Encode(byte[] dst, bool compressed)
 122         {
 123                 ToAffine();
 124                 if (compressed) {
 125                         int len = curve.EncodedLengthCompressed;
 126                         if (dst.Length != len) {
 127                                 throw new CryptoException(
 128                                         "invalid output length");
 129                         }
 130                         dst[0] = (byte)(0x02 + my.GetLSB());
 131                         mx.Encode(dst, 1, len - 1);
 132                 } else {
 133                         int len = curve.EncodedLength;
 134                         if (dst.Length != len) {
 135                                 throw new CryptoException(
 136                                         "invalid output length");
 137                         }
 138                         int flen = (len - 1) >> 1;
 139                         dst[0] = 0x04;
 140                         mx.Encode(dst, 1, flen);
 141                         my.Encode(dst, 1 + flen, flen);
 142                 }
 143                 return ~IsInfinityCT;
 144         }
 145
 146         internal override uint DecodeCT(byte[] enc)
 147         {
 148                 /*
 149                  * Format (specified in IEEE P1363, annex E):
 150                  *
 151                  *   0x00             point at infinity
 152                  *   0x02+b <X>       compressed, b = lsb of Y
 153                  *   0x04 <X> <Y>     uncompressed
 154                  *   0x06+b <X> <Y>   uncompressed, b = lsb of Y
 155                  *
 156                  * Coordinates X and Y are in unsigned big-endian
 157                  * notation with exactly the length of the modulus.
 158                  *
 159                  * We want constant-time decoding, up to the encoded
 160                  * length. This means that the four following situations
 161                  * can be differentiated:
 162                  * -- Point is zero (length = 1)
 163                  * -- Point is compressed (length = 1 + flen)
 164                  * -- Point is uncompressed or hybrid (length = 1 + 2*flen)
 165                  * -- Length is neither 1, 1+flen or 1+2*flen.
 166                  */
 167
 168                 int flen = curve.flen;
 169                 uint good = 0xFFFFFFFF;
 170                 if (enc.Length == 1) {
 171                         /*
 172                          * 1-byte encoding is point at infinity; the
 173                          * byte shall have value 0.
 174                          */
 175                         int z = enc[0];
 176                         good &= ~(uint)((z | -z) >> 31);
 177                         SetZero();
 178                 } else if (enc.Length == 1 + flen) {
 179                         /*
 180                          * Compressed encoding. Leading byte is 0x02 or
 181                          * 0x03.
 182                          */
 183                         int z = (enc[0] & 0xFE) - 0x02;
 184                         good &= ~(uint)((z | -z) >> 31);
 185                         uint lsbValue = (uint)(enc[0] & 1);
 186                         good &= mx.Decode(enc, 1, flen);
 187                         RebuildY2();
 188                         if (curve.pMod4 == 3) {
 189                                 good &= my.SqrtBlum();
 190                         } else {
 191                                 /*
 192                                  * Square roots modulo a non-Blum prime
 193                                  * are a bit more complex. We do not
 194                                  * support them yet (TODO).
 195                                  */
 196                                 good = 0x00000000;
 197                         }
 198
 199                         /*
 200                          * Adjust Y depending on LSB.
 201                          */
 202                         mt1.Set(my);
 203                         mt1.Negate();
 204                         uint dn = (uint)-(int)(my.GetLSB() ^ lsbValue);
 205                         my.CondCopy(mt1, dn);
 206
 207                         /*
 208                          * A corner case: LSB adjustment works only if
 209                          * Y != 0. If Y is 0 and requested LSB is 1,
 210                          * then the decoding fails. Note that this case
 211                          * cannot happen with usual prime curves, because
 212                          * they have a prime order, implying that there is
 213                          * no valid point such that Y = 0 (that would be
 214                          * a point of order 2).
 215                          */
 216                         good &= ~(uint)-(int)(my.GetLSB() ^ lsbValue);
 217
 218                         mz.Set(1);
 219                 } else if (enc.Length == 1 + (flen << 1)) {
 220                         /*
 221                          * Uncompressed or hybrid. Leading byte is either
 222                          * 0x04, 0x06 or 0x07. We verify that the X and
 223                          * Y coordinates fulfill the curve equation.
 224                          */
 225                         int fb = enc[0];
 226                         int z = (fb & 0xFC) - 0x04;
 227                         good &= ~(uint)((z | -z) >> 31);
 228                         z = fb - 0x05;
 229                         good &= (uint)((z | -z) >> 31);
 230                         good &= mx.Decode(enc, 1, flen);
 231                         RebuildY2();
 232                         mt1.Set(my);
 233                         mt1.FromMonty();
 234                         good &= my.Decode(enc, 1 + flen, flen);
 235                         mt2.Set(my);
 236                         mt2.MontySquare();
 237                         good &= mt1.EqCT(mt2);
 238
 239                         /*
 240                          * We must check the LSB for hybrid encoding.
 241                          * The check fails if the encoding is marked as
 242                          * hybrid AND the LSB does not match.
 243                          */
 244                         int lm = (fb >> 1) & ((int)my.GetLSB() ^ fb) & 1;
 245                         good &= ~(uint)-lm;
 246
 247                         mz.Set(1);
 248                 } else {
 249                         good = 0x00000000;
 250                 }
 251
 252                 /*
 253                  * If decoding failed, then we force the value to 0.
 254                  * Otherwise, we got a value. Either way, this uses
 255                  * affine coordinates.
 256                  */
 257                 mx.CondCopy(curve.mp, ~good);
 258                 my.CondCopy(curve.mp, ~good);
 259                 mz.CondCopy(curve.mp, ~good);
 260                 affine = 0xFFFFFFFF;
 261                 return good;
 262         }
 263
 264         internal override byte[] X {
 265                 get {
 266                         ToAffine();
 267                         return mx.Encode();
 268                 }
 269         }
 270
 271         internal override byte[] Y {
 272                 get {
 273                         ToAffine();
 274                         return my.Encode();
 275                 }
 276         }
 277
 278         internal override MutableECPoint Dup()
 279         {
 280                 MutableECPointPrime Q = new MutableECPointPrime(curve);
 281                 Q.Set(this);
 282                 return Q;
 283         }
 284
 285         internal void Set(byte[] X, byte[] Y, bool check)
 286         {
 287                 mx.Decode(X);
 288                 my.Decode(Y);
 289                 mz.Set(1);
 290                 affine = 0xFFFFFFFF;
 291                 if (check) {
 292                         CheckEquation();
 293                 }
 294         }
 295
 296         internal void Set(ModInt X, ModInt Y, bool check)
 297         {
 298                 mx.Set(X);
 299                 my.Set(Y);
 300                 mz.Set(1);
 301                 affine = 0xFFFFFFFF;
 302                 if (check) {
 303                         CheckEquation();
 304                 }
 305         }
 306
 307         void CheckEquation()
 308         {
 309                 curve.RebuildY2(mx, mt1, mt2);
 310                 mt2.Set(my);
 311                 mt2.ToMonty();
 312                 mt2.MontyMul(my);
 313                 if (!mt1.Eq(mt2)) {
 314                         throw new CryptoException(
 315                                 "Point is not on the curve");
 316                 }
 317         }
 318
 319         internal override void SetZero()
 320         {
 321                 mx.Set(0);
 322                 my.Set(0);
 323                 mz.Set(0);
 324                 affine = 0xFFFFFFFF;
 325         }
 326
 327         internal override void Set(MutableECPoint Q)
 328         {
 329                 MutableECPointPrime R = SameCurve(Q);
 330                 mx.Set(R.mx);
 331                 my.Set(R.my);
 332                 mz.Set(R.mz);
 333                 affine = R.affine;
 334         }
 335
 336         internal override void Set(MutableECPoint Q, uint ctl)
 337         {
 338                 MutableECPointPrime R = SameCurve(Q);
 339                 mx.CondCopy(R.mx, ctl);
 340                 my.CondCopy(R.my, ctl);
 341                 mz.CondCopy(R.mz, ctl);
 342                 affine ^= ctl & (affine ^ R.affine);
 343         }
 344
 345         internal override void SetMux(uint ctl,
 346                 MutableECPoint P1, MutableECPoint P2)
 347         {
 348                 SetMuxInner(ctl, SameCurve(P1), SameCurve(P2));
 349         }
 350
 351         void SetMuxInner(uint ctl,
 352                 MutableECPointPrime P1, MutableECPointPrime P2)
 353         {
 354                 mx.CopyMux(ctl, P1.mx, P2.mx);
 355                 my.CopyMux(ctl, P1.my, P2.my);
 356                 mz.CopyMux(ctl, P1.mz, P2.mz);
 357                 affine = P2.affine ^ (ctl & (P1.affine ^ P2.affine));
 358         }
 359
 360         internal override void DoubleCT()
 361         {
 362                 ToJacobian();
 363
 364                 /*
 365                  * Formulas are:
 366                  *   S = 4*X*Y^2
 367                  *   M = 3*X^2 + a*Z^4
 368                  *   X' = M^2 - 2*S
 369                  *   Y' = M*(S - X') - 8*Y^4
 370                  *   Z' = 2*Y*Z
 371                  *
 372                  * These formulas also happen to work properly (with our
 373                  * chosen representation) when the source point has
 374                  * order 2 (Y = 0 implies Z' = 0) and when the source
 375                  * point is already the point at infinity (Z = 0 implies
 376                  * Z' = 0).
 377                  *
 378                  * When a = -3, the value of M can be computed with the
 379                  * more efficient formula:
 380                  *   M = 3*(X+Z^2)*(X-Z^2)
 381                  */
 382
 383                 /*
 384                  * Compute M in t1.
 385                  */
 386                 if (curve.aIsM3) {
 387                         /*
 388                          * Set t1 = Z^2.
 389                          */
 390                         mt1.Set(mz);
 391                         mt1.MontySquare();
 392
 393                         /*
 394                          * Set t2 = X-Z^2 and then t1 = X+Z^2.
 395                          */
 396                         mt2.Set(mx);
 397                         mt2.Sub(mt1);
 398                         mt1.Add(mx);
 399
 400                         /*
 401                          * Set t1 = 3*(X+Z^2)*(X-Z^2).
 402                          */
 403                         mt1.MontyMul(mt2);
 404                         mt2.Set(mt1);
 405                         mt1.Add(mt2);
 406                         mt1.Add(mt2);
 407                 } else {
 408                         /*
 409                          * Set t1 = 3*X^2.
 410                          */
 411                         mt1.Set(mx);
 412                         mt1.MontySquare();
 413                         mt2.Set(mt1);
 414                         mt1.Add(mt2);
 415                         mt1.Add(mt2);
 416
 417                         /*
 418                          * Set t2 = a*Z^4.
 419                          */
 420                         mt2.Set(mz);
 421                         mt2.MontySquare();
 422                         mt2.MontySquare();
 423                         mt2.MontyMul(curve.ma);
 424
 425                         /*
 426                          * Set t1 = 3*X^2 + a*Z^4.
 427                          */
 428                         mt1.Add(mt2);
 429                 }
 430
 431                 /*
 432                  * Compute S = 4*X*Y^2 in t2. We also save 2*Y^2 in mt3.
 433                  */
 434                 mt2.Set(my);
 435                 mt2.MontySquare();
 436                 mt2.Add(mt2);
 437                 mt3.Set(mt2);
 438                 mt2.Add(mt2);
 439                 mt2.MontyMul(mx);
 440
 441                 /*
 442                  * Compute X' = M^2 - 2*S.
 443                  */
 444                 mx.Set(mt1);
 445                 mx.MontySquare();
 446                 mx.Sub(mt2);
 447                 mx.Sub(mt2);
 448
 449                 /*
 450                  * Compute Z' = 2*Y*Z.
 451                  */
 452                 mz.MontyMul(my);
 453                 mz.Add(mz);
 454
 455                 /*
 456                  * Compute Y' = M*(S - X') - 8*Y^4. We already have
 457                  * 4*Y^2 in t3.
 458                  */
 459                 mt2.Sub(mx);
 460                 mt2.MontyMul(mt1);
 461                 mt3.MontySquare();
 462                 mt3.Add(mt3);
 463                 my.Set(mt2);
 464                 my.Sub(mt3);
 465         }
 466
 467         internal override uint AddCT(MutableECPoint Q)
 468         {
 469                 MutableECPointPrime P2 = SameCurve(Q);
 470
 471                 if (P2.affine != 0) {
 472                         ms4.Set(P2.mx);
 473                         ms5.Set(P2.my);
 474                         ms6.Set(P2.mz);
 475                         ms4.ToMonty();
 476                         ms5.ToMonty();
 477                         ms6.SetMonty(~ms6.IsZeroCT);
 478                         return AddCTInner(ms4, ms5, ms6, true);
 479                 } else {
 480                         return AddCTInner(P2.mx, P2.my, P2.mz, false);
 481                 }
 482         }
 483
 484         /*
 485          * Inner function for addition. The Jacobian coordinates for
 486          * the operand are provided in Montogomery representation. If
 487          * p2affine is true, then it is guaranteed that p2z is 1
 488          * (converted to Montogomery).
 489          */
 490         uint AddCTInner(ModInt p2x, ModInt p2y, ModInt p2z, bool p2affine)
 491         {
 492                 /*
 493                  * In this comment, the two operands are called P1 and
 494                  * P2. P1 is this instance; P2 is the operand. Coordinates
 495                  * of P1 are (X1,Y1,Z1). Coordinates of P2 are (X2,Y2,Z2).
 496                  *
 497                  * Formulas:
 498                  *   U1 = X1 * Z2^2
 499                  *   U2 = X2 * Z1^2
 500                  *   S1 = Y1 * Z2^3
 501                  *   S2 = Y2 * Z1^3
 502                  *   H = U2 - U1
 503                  *   R = S2 - S1
 504                  *   X3 = R^2 - H^3 - 2*U1*H^2
 505                  *   Y3 = R*(U1*H^2 - X3) - S1*H^3
 506                  *   Z3 = H*Z1*Z2
 507                  *
 508                  * If both P1 and P2 are 0, then the formulas yield 0,
 509                  * which is fine. If one of P1 and P2 is 0 (but not both),
 510                  * then we get 0 as result, which is wrong and must be
 511                  * fixed at the end.
 512                  *
 513                  * If U1 == U2 and S1 == S2 then this means that either
 514                  * P1 or P2 is 0 (or both), or P1 == P2. In the latter
 515                  * case, the formulas are wrong and we must report
 516                  * an error.
 517                  *
 518                  * If U1 == U2 and S1 != S2 then P1 + P2 = 0. We get H = 0,
 519                  * which implies that we obtain the point at infinity,
 520                  * which is fine.
 521                  */
 522
 523                 uint P1IsZero = mz.IsZeroCT;
 524                 uint P2IsZero = p2z.IsZeroCT;
 525
 526                 ToJacobian();
 527
 528                 /*
 529                  * Save this value, in case the operand turns out to
 530                  * be the point at infinity.
 531                  */
 532                 ms1.Set(mx);
 533                 ms2.Set(my);
 534                 ms3.Set(mz);
 535
 536                 /*
 537                  * Compute U1 = X1*Z2^2 in t1, and S1 = Y1*Z2^3 in t3.
 538                  */
 539                 if (p2affine) {
 540                         mt1.Set(mx);
 541                         mt3.Set(my);
 542                 } else {
 543                         mt3.Set(p2z);
 544                         mt3.MontySquare();
 545                         mt1.Set(mx);
 546                         mt1.MontyMul(mt3);
 547                         mt3.MontyMul(p2z);
 548                         mt3.MontyMul(my);
 549                 }
 550                 //PrintMR(" u1 = x1*z2^2", mt1);
 551                 //PrintMR(" s1 = y1*z2^3", mt3);
 552
 553                 /*
 554                  * Compute U2 = X2*Z1^2 in t2, and S2 = Y2*Z1^3 in t4.
 555                  */
 556                 mt4.Set(mz);
 557                 mt4.MontySquare();
 558                 mt2.Set(p2x);
 559                 mt2.MontyMul(mt4);
 560                 mt4.MontyMul(mz);
 561                 mt4.MontyMul(p2y);
 562                 //PrintMR(" u2 = x2*z1^2", mt2);
 563                 //PrintMR(" s2 = y2*z1^3", mt4);
 564
 565                 /*
 566                  * Compute H = U2 - U1 in t2, and R = S2 - S1 in t4.
 567                  */
 568                 mt2.Sub(mt1);
 569                 mt4.Sub(mt3);
 570                 //PrintMR(" h = u2-u1", mt2);
 571                 //PrintMR(" r = s2-s1", mt4);
 572
 573                 /*
 574                  * If both H and R are 0, then we may have a problem
 575                  * (either P1 == P2, or P1 == 0, or P2 == 0).
 576                  */
 577                 uint formProb = mt2.IsZeroCT & mt4.IsZeroCT;
 578
 579                 /*
 580                  * Compute U1*H^2 in t1 and H^3 in t5.
 581                  */
 582                 mt5.Set(mt2);
 583                 mt5.MontySquare();
 584                 mt1.MontyMul(mt5);
 585                 mt5.MontyMul(mt2);
 586                 //PrintMR(" u1*h^2", mt1);
 587                 //PrintMR(" h^3", mt5);
 588
 589                 /*
 590                  * Compute X3 = R^2 - H^3 - 2*U1*H^2.
 591                  */
 592                 mx.Set(mt4);
 593                 mx.MontySquare();
 594                 mx.Sub(mt5);
 595                 mx.Sub(mt1);
 596                 mx.Sub(mt1);
 597                 //PrintMR(" x3 = r^2-h^3-2*u1*h^2", mx);
 598
 599                 /*
 600                  * Compute Y3 = R*(U1*H^2 - X3) - S1*H^3.
 601                  */
 602                 mt1.Sub(mx);
 603                 mt1.MontyMul(mt4);
 604                 mt5.MontyMul(mt3);
 605                 mt1.Sub(mt5);
 606                 my.Set(mt1);
 607                 //PrintMR(" y3 = r*(u1*h^2-x3)-s1*h^3", my);
 608
 609                 /*
 610                  * Compute Z3 = H*Z1*Z2.
 611                  */
 612                 mz.MontyMul(mt2);
 613                 if (!p2affine) {
 614                         mz.MontyMul(p2z);
 615                 }
 616                 //PrintMR(" z3 = h*z1*z2", mz);
 617
 618                 /*
 619                  * Fixup: handle the cases where P1 = 0 or P2 = 0.
 620                  */
 621                 mx.CondCopy(ms1, P2IsZero);
 622                 my.CondCopy(ms2, P2IsZero);
 623                 mz.CondCopy(ms3, P2IsZero);
 624                 mx.CondCopy(p2x, P1IsZero);
 625                 my.CondCopy(p2y, P1IsZero);
 626                 mz.CondCopy(p2z, P1IsZero);
 627
 628                 /*
 629                  * Report failure when P1 == P2, except when one of
 630                  * the points was zero (or both) because that case
 631                  * was properly handled.
 632                  */
 633                 return (~formProb) | P1IsZero | P2IsZero;
 634         }
 635
 636         internal override void NegCT()
 637         {
 638                 my.Negate();
 639         }
 640
 641         internal override uint MulSpecCT(byte[] n)
 642         {
 643                 uint good = 0xFFFFFFFF;
 644
 645                 /*
 646                  * Create and populate window.
 647                  *
 648                  * If this instance is 0, then we only add 0 to 0 and
 649                  * double 0, for which DoubleCT() and AddCT() work
 650                  * properly.
 651                  *
 652                  * If this instance (P) is non-zero, then x*P for all
 653                  * x in the 1..16 range shall be non-zero and distinct,
 654                  * since the subgroup order is prime and at least 17.
 655                  * Thus, we never add two equal points together in the
 656                  * window construction.
 657                  *
 658                  * We MUST ensure that all points are in the same
 659                  * coordinate convention (affine or Jacobian) to ensure
 660                  * constant-time execution. TODO: measure to see which
 661                  * is best: all affine or all Jacobian. All affine implies
 662                  * 14 or 15 extra divisions, but saves a few hundreds of
 663                  * multiplications.
 664                  */
 665                 MutableECPointPrime[] w = new MutableECPointPrime[16];
 666                 w[0] = new MutableECPointPrime(curve);
 667                 w[0].ToJacobian();
 668                 w[1] = new MutableECPointPrime(curve);
 669                 w[1].Set(this);
 670                 w[1].ToJacobian();
 671                 for (int i = 2; (i + 1) < w.Length; i += 2) {
 672                         w[i] = new MutableECPointPrime(curve);
 673                         w[i].Set(w[i >> 1]);
 674                         w[i].DoubleCT();
 675                         w[i + 1] = new MutableECPointPrime(curve);
 676                         w[i + 1].Set(w[i]);
 677                         good &= w[i + 1].AddCT(this);
 678                 }
 679
 680                 /* obsolete
 681                 for (int i = 0; i < w.Length; i ++) {
 682                         w[i].ToAffine();
 683                         w[i].Print("Win " + i);
 684                         w[i].ToJacobian();
 685                 }
 686                 Console.WriteLine("good = {0}", (int)good);
 687                 */
 688
 689                 /*
 690                  * Set this value to 0. We also set it already to
 691                  * Jacobian coordinates, since it will be done that
 692                  * way anyway. This instance will serve as accumulator.
 693                  */
 694                 mx.Set(0);
 695                 my.Set(0);
 696                 mz.Set(0);
 697                 affine = 0x00000000;
 698
 699                 /*
 700                  * We process the multiplier by 4-bit nibbles, starting
 701                  * with the most-significant one (the high nibble of the
 702                  * first byte, since we use big-endian notation).
 703                  *
 704                  * For each nibble, we perform a constant-time lookup
 705                  * in the window, to obtain the point to add to the
 706                  * current value of the accumulator. Thanks to the
 707                  * conditions on the operands (prime subgroup order and
 708                  * so on), all the additions below must work.
 709                  */
 710                 MutableECPointPrime t = new MutableECPointPrime(curve);
 711                 for (int i = (n.Length << 1) - 1; i >= 0; i --) {
 712                         int b = n[n.Length - 1 - (i >> 1)];
 713                         int j = (b >> ((i & 1) << 2)) & 0x0F;
 714                         for (int k = 0; k < 16; k ++) {
 715                                 t.Set(w[k], ~(uint)(((j - k) | (k - j)) >> 31));
 716                         }
 717                         good &= AddCT(t);
 718                         if (i > 0) {
 719                                 DoubleCT();
 720                                 DoubleCT();
 721                                 DoubleCT();
 722                                 DoubleCT();
 723                         }
 724                 }
 725
 726                 return good;
 727         }
 728
 729         internal override uint EqCT(MutableECPoint Q)
 730         {
 731                 MutableECPointPrime R = SameCurve(Q);
 732                 if (affine != 0) {
 733                         if (R.affine != 0) {
 734                                 return mx.EqCT(R.mx)
 735                                         & my.EqCT(R.my)
 736                                         & mz.EqCT(R.mz);
 737                         } else {
 738                                 return EqCTMixed(R, this);
 739                         }
 740                 } else if (R.affine != 0) {
 741                         return EqCTMixed(this, R);
 742                 }
 743
 744                 /*
 745                  * Both points are in Jacobian coordinates.
 746                  * If Z1 and Z2 are non-zero, then equality is
 747                  * achieved if and only if both following equations
 748                  * are true:
 749                  *     X1*(Z2^2) = X2*(Z1^2)
 750                  *     Y1*(Z2^3) = Y2*(Z1^3)
 751                  * If Z1 or Z2 is zero, then equality is achieved
 752                  * if and only if both are zero.
 753                  */
 754                 mt1.Set(mz);
 755                 mt1.MontySquare();
 756                 mt2.Set(R.mz);
 757                 mt2.MontySquare();
 758                 mt3.Set(mx);
 759                 mt3.MontyMul(mt2);
 760                 mt4.Set(R.mx);
 761                 mt4.MontyMul(mt1);
 762                 uint r = mt3.EqCT(mt4);
 763                 mt1.MontyMul(mz);
 764                 mt2.MontyMul(R.mz);
 765                 mt3.Set(my);
 766                 mt3.MontyMul(mt2);
 767                 mt4.Set(R.my);
 768                 mt4.MontyMul(mt1);
 769                 r &= mt3.EqCT(mt4);
 770                 uint z1z = mz.IsZeroCT;
 771                 uint z2z = R.mz.IsZeroCT;
 772                 return (r & ~(z1z | z2z)) ^ (z1z & z2z);
 773         }
 774
 775         /*
 776          * Mixed comparison: P1 is in Jacobian coordinates, P2 is in
 777          * affine coordinates.
 778          */
 779         uint EqCTMixed(MutableECPointPrime P1, MutableECPointPrime P2)
 780         {
 781                 /*
 782                  * If either P1 or P2 is infinity, then they are equal
 783                  * if and only if they both are infinity.
 784                  *
 785                  * If neither is infinity, then we must check the following:
 786                  *    X1 = X2*(Z1^2)
 787                  *    Y1 = Y2*(Z1^3)
 788                  * Beware that X1, Y1 and Z1 are in Montgomery representation,
 789                  * while X2 and Y2 are not.
 790                  */
 791                 mt1.Set(P1.mz);
 792                 mt1.MontySquare();
 793                 mt2.Set(P2.mx);
 794                 mt2.MontyMul(mt1);
 795                 mt3.Set(P1.mx);
 796                 mt3.FromMonty();
 797                 uint r = mt2.EqCT(mt3);
 798                 mt1.MontyMul(P1.mz);
 799                 mt1.MontyMul(P2.my);
 800                 mt2.Set(P1.my);
 801                 mt2.FromMonty();
 802                 r &= mt1.EqCT(mt2);
 803                 uint z1z = P1.mz.IsZeroCT;
 804                 uint z2z = P2.mz.IsZeroCT;
 805                 return (r & ~(z1z | z2z)) ^ (z1z & z2z);
 806         }
 807
 808         MutableECPointPrime SameCurve(MutableECPoint Q)
 809         {
 810                 MutableECPointPrime R = Q as MutableECPointPrime;
 811                 if (R == null || !curve.Equals(R.curve)) {
 812                         throw new CryptoException("Mixed curves");
 813                 }
 814                 return R;
 815         }
 816
 817         /*
 818          * Convert to Jabobian coordinates (if not already done).
 819          */
 820         void ToJacobian()
 821         {
 822                 if (affine == 0) {
 823                         return;
 824                 }
 825
 826                 /*
 827                  * Since Z = 0 or 1 in affine coordinates, we can
 828                  * use SetMonty().
 829                  */
 830                 mx.ToMonty();
 831                 my.ToMonty();
 832                 mz.SetMonty(~mz.IsZeroCT);
 833                 affine = 0x00000000;
 834         }
 835
 836         /*
 837          * Convert to affine coordinates (if not already done).
 838          */
 839         void ToAffine()
 840         {
 841                 if (affine != 0) {
 842                         return;
 843                 }
 844
 845                 /*
 846                  * Divisions are expensive, so we want to make only one,
 847                  * not two. This involves some games with Montgomery
 848                  * representation.
 849                  *
 850                  * A number a in Montgomery representation means that
 851                  * the value we have is equal to aR. Montgomery
 852                  * multiplication of a by b yields ab/R (so, if we
 853                  * apply it to aR and bR, we get abR).
 854                  */
 855
 856                 /* Save Z*R in mt1. */
 857                 mt1.Set(mz);
 858
 859                 /* Compute Z^3 in mz. */
 860                 mz.MontySquare();
 861                 mz.MontyMul(mt1);
 862                 mz.FromMonty();
 863
 864                 /* Compute t2 = 1/Z^3. */
 865                 mt2.Set(mz);
 866                 mt2.Invert();
 867                 uint cc = ~mt2.IsZeroCT;
 868
 869                 /* Compute y. */
 870                 my.MontyMul(mt2);
 871
 872                 /* Compute t2 = 1/Z^2. */
 873                 mt2.MontyMul(mt1);
 874
 875                 /* Compute x. */
 876                 mx.MontyMul(mt2);
 877
 878                 /*
 879                  * If the point is infinity (division by Z^2 failed),
 880                  * then set all coordinates to 0. Otherwise, set mz
 881                  * to exactly 1.
 882                  */
 883                 mx.CondCopy(curve.mp, ~cc);
 884                 my.CondCopy(curve.mp, ~cc);
 885                 mz.Set((int)cc & 1);
 886
 887                 affine = 0xFFFFFFFF;
 888         }
 889
 890         /*
 891          * Compute Y^2 into my, using the value in mx as X. Both values
 892          * are in normal (non-Montgomery) representation.
 893          */
 894         void RebuildY2()
 895         {
 896                 my.Set(mx);
 897                 my.ToMonty();
 898                 mt1.Set(my);
 899                 my.MontySquare();
 900                 my.MontyMul(mx);
 901                 mt1.MontyMul(curve.ma);
 902                 my.Add(mt1);
 903                 my.Add(curve.mb);
 904         }
 905 }
 906
 907 }