mozilla/security/nss/lib/freebl/mpi/mpi.c - Issue 14249009: Change the NSS and NSPR source tree to the new directory structure to be

Side by Side Diff: mozilla/security/nss/lib/freebl/mpi/mpi.c

Issue 14249009: Change the NSS and NSPR source tree to the new directory structure to be (Closed) Base URL: svn://svn.chromium.org/chrome/trunk/deps/third_party/nss/

Patch Set: Created 7 years, 8 months ago

Use n/p to move between diff chunks; N/P to move between comments. Draft comments are only viewable by you.

Jump to:

View unified diff | Download patch | Annotate | Revision Log

OLD	NEW
	(Empty)
1 /*

2 * mpi.c

3 *

4 * Arbitrary precision integer arithmetic library

5 *

6 * This Source Code Form is subject to the terms of the Mozilla Public

7 * License, v. 2.0. If a copy of the MPL was not distributed with this

8 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

9 /* $Id: mpi.c,v 1.51 2012/04/25 14:49:50 gerv%gerv.net Exp $ */

10

11 #include "mpi-priv.h"

12 #if defined(OSF1)

13 #include <c_asm.h>

14 #endif

15

16 #if defined(__arm__) && \

17 ((defined(__thumb__) && !defined(__thumb2__)) \|\| defined(__ARM_ARCH_3__))

18 /* 16-bit thumb or ARM v3 doesn't work inlined assember version */

19 #undef MP_ASSEMBLY_MULTIPLY

20 #undef MP_ASSEMBLY_SQUARE

21 #endif

22

23 #if MP_LOGTAB

24 /*

25 A table of the logs of 2 for various bases (the 0 and 1 entries of

26 this table are meaningless and should not be referenced).

27

28 This table is used to compute output lengths for the mp_toradix()

29 function. Since a number n in radix r takes up about log_r(n)

30 digits, we estimate the output size by taking the least integer

31 greater than log_r(n), where:

32

33 log_r(n) = log_2(n) * log_r(2)

34

35 This table, therefore, is a table of log_r(2) for 2 <= r <= 36,

36 which are the output bases supported.

37 */

38 #include "logtab.h"

39 #endif

40

41 /* {{{ Constant strings */

42

43 /* Constant strings returned by mp_strerror() */

44 static const char *mp_err_string[] = {

45 "unknown result code", /* say what? */

46 "boolean true", /* MP_OKAY, MP_YES */

47 "boolean false", /* MP_NO */

48 "out of memory", /* MP_MEM */

49 "argument out of range", /* MP_RANGE */

50 "invalid input parameter", /* MP_BADARG */

51 "result is undefined" /* MP_UNDEF */

52 };

53

54 /* Value to digit maps for radix conversion */

55

56 /* s_dmap_1 - standard digits and letters */

57 static const char *s_dmap_1 =

58 "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";

59

60 /* }}} */

61

62 unsigned long mp_allocs;

63 unsigned long mp_frees;

64 unsigned long mp_copies;

65

66 /* {{{ Default precision manipulation */

67

68 /* Default precision for newly created mp_int's */

69 static mp_size s_mp_defprec = MP_DEFPREC;

70

71 mp_size mp_get_prec(void)

72 {

73 return s_mp_defprec;

74

75 } /* end mp_get_prec() */

76

77 void mp_set_prec(mp_size prec)

78 {

79 if(prec == 0)

80 s_mp_defprec = MP_DEFPREC;

81 else

82 s_mp_defprec = prec;

83

84 } /* end mp_set_prec() */

85

86 /* }}} */

87

88 /------------------------------------------------------------------------/

89 /* {{{ mp_init(mp) */

90

91 /*

92 mp_init(mp)

93

94 Initialize a new zero-valued mp_int. Returns MP_OKAY if successful,

95 MP_MEM if memory could not be allocated for the structure.

96 */

97

98 mp_err mp_init(mp_int *mp)

99 {

100 return mp_init_size(mp, s_mp_defprec);

101

102 } /* end mp_init() */

103

104 /* }}} */

105

106 /* {{{ mp_init_size(mp, prec) */

107

108 /*

109 mp_init_size(mp, prec)

110

111 Initialize a new zero-valued mp_int with at least the given

112 precision; returns MP_OKAY if successful, or MP_MEM if memory could

113 not be allocated for the structure.

114 */

115

116 mp_err mp_init_size(mp_int *mp, mp_size prec)

117 {

118 ARGCHK(mp != NULL && prec > 0, MP_BADARG);

119

120 prec = MP_ROUNDUP(prec, s_mp_defprec);

121 if((DIGITS(mp) = s_mp_alloc(prec, sizeof(mp_digit))) == NULL)

122 return MP_MEM;

123

124 SIGN(mp) = ZPOS;

125 USED(mp) = 1;

126 ALLOC(mp) = prec;

127

128 return MP_OKAY;

129

130 } /* end mp_init_size() */

131

132 /* }}} */

133

134 /* {{{ mp_init_copy(mp, from) */

135

136 /*

137 mp_init_copy(mp, from)

138

139 Initialize mp as an exact copy of from. Returns MP_OKAY if

140 successful, MP_MEM if memory could not be allocated for the new

141 structure.

142 */

143

144 mp_err mp_init_copy(mp_int mp, const mp_int from)

145 {

146 ARGCHK(mp != NULL && from != NULL, MP_BADARG);

147

148 if(mp == from)

149 return MP_OKAY;

150

151 if((DIGITS(mp) = s_mp_alloc(ALLOC(from), sizeof(mp_digit))) == NULL)

152 return MP_MEM;

153

154 s_mp_copy(DIGITS(from), DIGITS(mp), USED(from));

155 USED(mp) = USED(from);

156 ALLOC(mp) = ALLOC(from);

157 SIGN(mp) = SIGN(from);

158

159 return MP_OKAY;

160

161 } /* end mp_init_copy() */

162

163 /* }}} */

164

165 /* {{{ mp_copy(from, to) */

166

167 /*

168 mp_copy(from, to)

169

170 Copies the mp_int 'from' to the mp_int 'to'. It is presumed that

171 'to' has already been initialized (if not, use mp_init_copy()

172 instead). If 'from' and 'to' are identical, nothing happens.

173 */

174

175 mp_err mp_copy(const mp_int from, mp_int to)

176 {

177 ARGCHK(from != NULL && to != NULL, MP_BADARG);

178

179 if(from == to)

180 return MP_OKAY;

181

182 { /* copy */

183 mp_digit *tmp;

184

185 /*

186 If the allocated buffer in 'to' already has enough space to hold

187 all the used digits of 'from', we'll re-use it to avoid hitting

188 the memory allocater more than necessary; otherwise, we'd have

189 to grow anyway, so we just allocate a hunk and make the copy as

190 usual

191 */

192 if(ALLOC(to) >= USED(from)) {

193 s_mp_setz(DIGITS(to) + USED(from), ALLOC(to) - USED(from));

194 s_mp_copy(DIGITS(from), DIGITS(to), USED(from));

195

196 } else {

197 if((tmp = s_mp_alloc(ALLOC(from), sizeof(mp_digit))) == NULL)

198 return MP_MEM;

199

200 s_mp_copy(DIGITS(from), tmp, USED(from));

201

202 if(DIGITS(to) != NULL) {

203 #if MP_CRYPTO

204 s_mp_setz(DIGITS(to), ALLOC(to));

205 #endif

206 s_mp_free(DIGITS(to));

207 }

208

209 DIGITS(to) = tmp;

210 ALLOC(to) = ALLOC(from);

211 }

212

213 /* Copy the precision and sign from the original */

214 USED(to) = USED(from);

215 SIGN(to) = SIGN(from);

216 } /* end copy */

217

218 return MP_OKAY;

219

220 } /* end mp_copy() */

221

222 /* }}} */

223

224 /* {{{ mp_exch(mp1, mp2) */

225

226 /*

227 mp_exch(mp1, mp2)

228

229 Exchange mp1 and mp2 without allocating any intermediate memory

230 (well, unless you count the stack space needed for this call and the

231 locals it creates...). This cannot fail.

232 */

233

234 void mp_exch(mp_int mp1, mp_int mp2)

235 {

236 #if MP_ARGCHK == 2

237 assert(mp1 != NULL && mp2 != NULL);

238 #else

239 if(mp1 == NULL \|\| mp2 == NULL)

240 return;

241 #endif

242

243 s_mp_exch(mp1, mp2);

244

245 } /* end mp_exch() */

246

247 /* }}} */

248

249 /* {{{ mp_clear(mp) */

250

251 /*

252 mp_clear(mp)

253

254 Release the storage used by an mp_int, and void its fields so that

255 if someone calls mp_clear() again for the same int later, we won't

256 get tollchocked.

257 */

258

259 void mp_clear(mp_int *mp)

260 {

261 if(mp == NULL)

262 return;

263

264 if(DIGITS(mp) != NULL) {

265 #if MP_CRYPTO

266 s_mp_setz(DIGITS(mp), ALLOC(mp));

267 #endif

268 s_mp_free(DIGITS(mp));

269 DIGITS(mp) = NULL;

270 }

271

272 USED(mp) = 0;

273 ALLOC(mp) = 0;

274

275 } /* end mp_clear() */

276

277 /* }}} */

278

279 /* {{{ mp_zero(mp) */

280

281 /*

282 mp_zero(mp)

283

284 Set mp to zero. Does not change the allocated size of the structure,

285 and therefore cannot fail (except on a bad argument, which we ignore)

286 */

287 void mp_zero(mp_int *mp)

288 {

289 if(mp == NULL)

290 return;

291

292 s_mp_setz(DIGITS(mp), ALLOC(mp));

293 USED(mp) = 1;

294 SIGN(mp) = ZPOS;

295

296 } /* end mp_zero() */

297

298 /* }}} */

299

300 /* {{{ mp_set(mp, d) */

301

302 void mp_set(mp_int *mp, mp_digit d)

303 {

304 if(mp == NULL)

305 return;

306

307 mp_zero(mp);

308 DIGIT(mp, 0) = d;

309

310 } /* end mp_set() */

311

312 /* }}} */

313

314 /* {{{ mp_set_int(mp, z) */

315

316 mp_err mp_set_int(mp_int *mp, long z)

317 {

318 int ix;

319 unsigned long v = labs(z);

320 mp_err res;

321

322 ARGCHK(mp != NULL, MP_BADARG);

323

324 mp_zero(mp);

325 if(z == 0)

326 return MP_OKAY; /* shortcut for zero */

327

328 if (sizeof v <= sizeof(mp_digit)) {

329 DIGIT(mp,0) = v;

330 } else {

331 for (ix = sizeof(long) - 1; ix >= 0; ix--) {

332 if ((res = s_mp_mul_d(mp, (UCHAR_MAX + 1))) != MP_OKAY)

333 return res;

334

335 res = s_mp_add_d(mp, (mp_digit)((v >> (ix * CHAR_BIT)) & UCHAR_MAX));

336 if (res != MP_OKAY)

337 return res;

338 }

339 }

340 if(z < 0)

341 SIGN(mp) = NEG;

342

343 return MP_OKAY;

344

345 } /* end mp_set_int() */

346

347 /* }}} */

348

349 /* {{{ mp_set_ulong(mp, z) */

350

351 mp_err mp_set_ulong(mp_int *mp, unsigned long z)

352 {

353 int ix;

354 mp_err res;

355

356 ARGCHK(mp != NULL, MP_BADARG);

357

358 mp_zero(mp);

359 if(z == 0)

360 return MP_OKAY; /* shortcut for zero */

361

362 if (sizeof z <= sizeof(mp_digit)) {

363 DIGIT(mp,0) = z;

364 } else {

365 for (ix = sizeof(long) - 1; ix >= 0; ix--) {

366 if ((res = s_mp_mul_d(mp, (UCHAR_MAX + 1))) != MP_OKAY)

367 return res;

368

369 res = s_mp_add_d(mp, (mp_digit)((z >> (ix * CHAR_BIT)) & UCHAR_MAX));

370 if (res != MP_OKAY)

371 return res;

372 }

373 }

374 return MP_OKAY;

375 } /* end mp_set_ulong() */

376

377 /* }}} */

378

379 /------------------------------------------------------------------------/

380 /* {{{ Digit arithmetic */

381

382 /* {{{ mp_add_d(a, d, b) */

383

384 /*

385 mp_add_d(a, d, b)

386

387 Compute the sum b = a + d, for a single digit d. Respects the sign of

388 its primary addend (single digits are unsigned anyway).

389 */

390

391 mp_err mp_add_d(const mp_int a, mp_digit d, mp_int b)

392 {

393 mp_int tmp;

394 mp_err res;

395

396 ARGCHK(a != NULL && b != NULL, MP_BADARG);

397

398 if((res = mp_init_copy(&tmp, a)) != MP_OKAY)

399 return res;

400

401 if(SIGN(&tmp) == ZPOS) {

402 if((res = s_mp_add_d(&tmp, d)) != MP_OKAY)

403 goto CLEANUP;

404 } else if(s_mp_cmp_d(&tmp, d) >= 0) {

405 if((res = s_mp_sub_d(&tmp, d)) != MP_OKAY)

406 goto CLEANUP;

407 } else {

408 mp_neg(&tmp, &tmp);

409

410 DIGIT(&tmp, 0) = d - DIGIT(&tmp, 0);

411 }

412

413 if(s_mp_cmp_d(&tmp, 0) == 0)

414 SIGN(&tmp) = ZPOS;

415

416 s_mp_exch(&tmp, b);

417

418 CLEANUP:

419 mp_clear(&tmp);

420 return res;

421

422 } /* end mp_add_d() */

423

424 /* }}} */

425

426 /* {{{ mp_sub_d(a, d, b) */

427

428 /*

429 mp_sub_d(a, d, b)

430

431 Compute the difference b = a - d, for a single digit d. Respects the

432 sign of its subtrahend (single digits are unsigned anyway).

433 */

434

435 mp_err mp_sub_d(const mp_int a, mp_digit d, mp_int b)

436 {

437 mp_int tmp;

438 mp_err res;

439

440 ARGCHK(a != NULL && b != NULL, MP_BADARG);

441

442 if((res = mp_init_copy(&tmp, a)) != MP_OKAY)

443 return res;

444

445 if(SIGN(&tmp) == NEG) {

446 if((res = s_mp_add_d(&tmp, d)) != MP_OKAY)

447 goto CLEANUP;

448 } else if(s_mp_cmp_d(&tmp, d) >= 0) {

449 if((res = s_mp_sub_d(&tmp, d)) != MP_OKAY)

450 goto CLEANUP;

451 } else {

452 mp_neg(&tmp, &tmp);

453

454 DIGIT(&tmp, 0) = d - DIGIT(&tmp, 0);

455 SIGN(&tmp) = NEG;

456 }

457

458 if(s_mp_cmp_d(&tmp, 0) == 0)

459 SIGN(&tmp) = ZPOS;

460

461 s_mp_exch(&tmp, b);

462

463 CLEANUP:

464 mp_clear(&tmp);

465 return res;

466

467 } /* end mp_sub_d() */

468

469 /* }}} */

470

471 /* {{{ mp_mul_d(a, d, b) */

472

473 /*

474 mp_mul_d(a, d, b)

475

476 Compute the product b = a * d, for a single digit d. Respects the sign

477 of its multiplicand (single digits are unsigned anyway)

478 */

479

480 mp_err mp_mul_d(const mp_int a, mp_digit d, mp_int b)

481 {

482 mp_err res;

483

484 ARGCHK(a != NULL && b != NULL, MP_BADARG);

485

486 if(d == 0) {

487 mp_zero(b);

488 return MP_OKAY;

489 }

490

491 if((res = mp_copy(a, b)) != MP_OKAY)

492 return res;

493

494 res = s_mp_mul_d(b, d);

495

496 return res;

497

498 } /* end mp_mul_d() */

499

500 /* }}} */

501

502 /* {{{ mp_mul_2(a, c) */

503

504 mp_err mp_mul_2(const mp_int a, mp_int c)

505 {

506 mp_err res;

507

508 ARGCHK(a != NULL && c != NULL, MP_BADARG);

509

510 if((res = mp_copy(a, c)) != MP_OKAY)

511 return res;

512

513 return s_mp_mul_2(c);

514

515 } /* end mp_mul_2() */

516

517 /* }}} */

518

519 /* {{{ mp_div_d(a, d, q, r) */

520

521 /*

522 mp_div_d(a, d, q, r)

523

524 Compute the quotient q = a / d and remainder r = a mod d, for a

525 single digit d. Respects the sign of its divisor (single digits are

526 unsigned anyway).

527 */

528

529 mp_err mp_div_d(const mp_int a, mp_digit d, mp_int q, mp_digit *r)

530 {

531 mp_err res;

532 mp_int qp;

533 mp_digit rem;

534 int pow;

535

536 ARGCHK(a != NULL, MP_BADARG);

537

538 if(d == 0)

539 return MP_RANGE;

540

541 /* Shortcut for powers of two ... */

542 if((pow = s_mp_ispow2d(d)) >= 0) {

543 mp_digit mask;

544

545 mask = ((mp_digit)1 << pow) - 1;

546 rem = DIGIT(a, 0) & mask;

547

548 if(q) {

549 mp_copy(a, q);

550 s_mp_div_2d(q, pow);

551 }

552

553 if(r)

554 *r = rem;

555

556 return MP_OKAY;

557 }

558

559 if((res = mp_init_copy(&qp, a)) != MP_OKAY)

560 return res;

561

562 res = s_mp_div_d(&qp, d, &rem);

563

564 if(s_mp_cmp_d(&qp, 0) == 0)

565 SIGN(q) = ZPOS;

566

567 if(r)

568 *r = rem;

569

570 if(q)

571 s_mp_exch(&qp, q);

572

573 mp_clear(&qp);

574 return res;

575

576 } /* end mp_div_d() */

577

578 /* }}} */

579

580 /* {{{ mp_div_2(a, c) */

581

582 /*

583 mp_div_2(a, c)

584

585 Compute c = a / 2, disregarding the remainder.

586 */

587

588 mp_err mp_div_2(const mp_int a, mp_int c)

589 {

590 mp_err res;

591

592 ARGCHK(a != NULL && c != NULL, MP_BADARG);

593

594 if((res = mp_copy(a, c)) != MP_OKAY)

595 return res;

596

597 s_mp_div_2(c);

598

599 return MP_OKAY;

600

601 } /* end mp_div_2() */

602

603 /* }}} */

604

605 /* {{{ mp_expt_d(a, d, b) */

606

607 mp_err mp_expt_d(const mp_int a, mp_digit d, mp_int c)

608 {

609 mp_int s, x;

610 mp_err res;

611

612 ARGCHK(a != NULL && c != NULL, MP_BADARG);

613

614 if((res = mp_init(&s)) != MP_OKAY)

615 return res;

616 if((res = mp_init_copy(&x, a)) != MP_OKAY)

617 goto X;

618

619 DIGIT(&s, 0) = 1;

620

621 while(d != 0) {

622 if(d & 1) {

623 if((res = s_mp_mul(&s, &x)) != MP_OKAY)

624 goto CLEANUP;

625 }

626

627 d /= 2;

628

629 if((res = s_mp_sqr(&x)) != MP_OKAY)

630 goto CLEANUP;

631 }

632

633 s_mp_exch(&s, c);

634

635 CLEANUP:

636 mp_clear(&x);

637 X:

638 mp_clear(&s);

639

640 return res;

641

642 } /* end mp_expt_d() */

643

644 /* }}} */

645

646 /* }}} */

647

648 /------------------------------------------------------------------------/

649 /* {{{ Full arithmetic */

650

651 /* {{{ mp_abs(a, b) */

652

653 /*

654 mp_abs(a, b)

655

656 Compute b = \|a\|. 'a' and 'b' may be identical.

657 */

658

659 mp_err mp_abs(const mp_int a, mp_int b)

660 {

661 mp_err res;

662

663 ARGCHK(a != NULL && b != NULL, MP_BADARG);

664

665 if((res = mp_copy(a, b)) != MP_OKAY)

666 return res;

667

668 SIGN(b) = ZPOS;

669

670 return MP_OKAY;

671

672 } /* end mp_abs() */

673

674 /* }}} */

675

676 /* {{{ mp_neg(a, b) */

677

678 /*

679 mp_neg(a, b)

680

681 Compute b = -a. 'a' and 'b' may be identical.

682 */

683

684 mp_err mp_neg(const mp_int a, mp_int b)

685 {

686 mp_err res;

687

688 ARGCHK(a != NULL && b != NULL, MP_BADARG);

689

690 if((res = mp_copy(a, b)) != MP_OKAY)

691 return res;

692

693 if(s_mp_cmp_d(b, 0) == MP_EQ)

694 SIGN(b) = ZPOS;

695 else

696 SIGN(b) = (SIGN(b) == NEG) ? ZPOS : NEG;

697

698 return MP_OKAY;

699

700 } /* end mp_neg() */

701

702 /* }}} */

703

704 /* {{{ mp_add(a, b, c) */

705

706 /*

707 mp_add(a, b, c)

708

709 Compute c = a + b. All parameters may be identical.

710 */

711

712 mp_err mp_add(const mp_int a, const mp_int b, mp_int *c)

713 {

714 mp_err res;

715

716 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

717

718 if(SIGN(a) == SIGN(b)) { /* same sign: add values, keep sign */

719 MP_CHECKOK( s_mp_add_3arg(a, b, c) );

720 } else if(s_mp_cmp(a, b) >= 0) { /* different sign: \|a\| >= \|b\| */

721 MP_CHECKOK( s_mp_sub_3arg(a, b, c) );

722 } else { /* different sign: \|a\| < \|b\| */

723 MP_CHECKOK( s_mp_sub_3arg(b, a, c) );

724 }

725

726 if (s_mp_cmp_d(c, 0) == MP_EQ)

727 SIGN(c) = ZPOS;

728

729 CLEANUP:

730 return res;

731

732 } /* end mp_add() */

733

734 /* }}} */

735

736 /* {{{ mp_sub(a, b, c) */

737

738 /*

739 mp_sub(a, b, c)

740

741 Compute c = a - b. All parameters may be identical.

742 */

743

744 mp_err mp_sub(const mp_int a, const mp_int b, mp_int *c)

745 {

746 mp_err res;

747 int magDiff;

748

749 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

750

751 if (a == b) {

752 mp_zero(c);

753 return MP_OKAY;

754 }

755

756 if (MP_SIGN(a) != MP_SIGN(b)) {

757 MP_CHECKOK( s_mp_add_3arg(a, b, c) );

758 } else if (!(magDiff = s_mp_cmp(a, b))) {

759 mp_zero(c);

760 res = MP_OKAY;

761 } else if (magDiff > 0) {

762 MP_CHECKOK( s_mp_sub_3arg(a, b, c) );

763 } else {

764 MP_CHECKOK( s_mp_sub_3arg(b, a, c) );

765 MP_SIGN(c) = !MP_SIGN(a);

766 }

767

768 if (s_mp_cmp_d(c, 0) == MP_EQ)

769 MP_SIGN(c) = MP_ZPOS;

770

771 CLEANUP:

772 return res;

773

774 } /* end mp_sub() */

775

776 /* }}} */

777

778 /* {{{ mp_mul(a, b, c) */

779

780 /*

781 mp_mul(a, b, c)

782

783 Compute c = a * b. All parameters may be identical.

784 */

785 mp_err mp_mul(const mp_int a, const mp_int b, mp_int * c)

786 {

787 mp_digit *pb;

788 mp_int tmp;

789 mp_err res;

790 mp_size ib;

791 mp_size useda, usedb;

792

793 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

794

795 if (a == c) {

796 if ((res = mp_init_copy(&tmp, a)) != MP_OKAY)

797 return res;

798 if (a == b)

799 b = &tmp;

800 a = &tmp;

801 } else if (b == c) {

802 if ((res = mp_init_copy(&tmp, b)) != MP_OKAY)

803 return res;

804 b = &tmp;

805 } else {

806 MP_DIGITS(&tmp) = 0;

807 }

808

809 if (MP_USED(a) < MP_USED(b)) {

810 const mp_int xch = b; / switch a and b, to do fewer outer loops */

811 b = a;

812 a = xch;

813 }

814

815 MP_USED(c) = 1; MP_DIGIT(c, 0) = 0;

816 if((res = s_mp_pad(c, USED(a) + USED(b))) != MP_OKAY)

817 goto CLEANUP;

818

819 #ifdef NSS_USE_COMBA

820 if ((MP_USED(a) == MP_USED(b)) && IS_POWER_OF_2(MP_USED(b))) {

821 if (MP_USED(a) == 4) {

822 s_mp_mul_comba_4(a, b, c);

823 goto CLEANUP;

824 }

825 if (MP_USED(a) == 8) {

826 s_mp_mul_comba_8(a, b, c);

827 goto CLEANUP;

828 }

829 if (MP_USED(a) == 16) {

830 s_mp_mul_comba_16(a, b, c);

831 goto CLEANUP;

832 }

833 if (MP_USED(a) == 32) {

834 s_mp_mul_comba_32(a, b, c);

835 goto CLEANUP;

836 }

837 }

838 #endif

839

840 pb = MP_DIGITS(b);

841 s_mpv_mul_d(MP_DIGITS(a), MP_USED(a), *pb++, MP_DIGITS(c));

842

843 /* Outer loop: Digits of b */

844 useda = MP_USED(a);

845 usedb = MP_USED(b);

846 for (ib = 1; ib < usedb; ib++) {

847 mp_digit b_i = *pb++;

848

849 /* Inner product: Digits of a */

850 if (b_i)

851 s_mpv_mul_d_add(MP_DIGITS(a), useda, b_i, MP_DIGITS(c) + ib);

852 else

853 MP_DIGIT(c, ib + useda) = b_i;

854 }

855

856 s_mp_clamp(c);

857

858 if(SIGN(a) == SIGN(b) \|\| s_mp_cmp_d(c, 0) == MP_EQ)

859 SIGN(c) = ZPOS;

860 else

861 SIGN(c) = NEG;

862

863 CLEANUP:

864 mp_clear(&tmp);

865 return res;

866 } /* end mp_mul() */

867

868 /* }}} */

869

870 /* {{{ mp_sqr(a, sqr) */

871

872 #if MP_SQUARE

873 /*

874 Computes the square of a. This can be done more

875 efficiently than a general multiplication, because many of the

876 computation steps are redundant when squaring. The inner product

877 step is a bit more complicated, but we save a fair number of

878 iterations of the multiplication loop.

879 */

880

881 /* sqr = a^2; Caller provides both a and tmp; */

882 mp_err mp_sqr(const mp_int a, mp_int sqr)

883 {

884 mp_digit *pa;

885 mp_digit d;

886 mp_err res;

887 mp_size ix;

888 mp_int tmp;

889 int count;

890

891 ARGCHK(a != NULL && sqr != NULL, MP_BADARG);

892

893 if (a == sqr) {

894 if((res = mp_init_copy(&tmp, a)) != MP_OKAY)

895 return res;

896 a = &tmp;

897 } else {

898 DIGITS(&tmp) = 0;

899 res = MP_OKAY;

900 }

901

902 ix = 2 * MP_USED(a);

903 if (ix > MP_ALLOC(sqr)) {

904 MP_USED(sqr) = 1;

905 MP_CHECKOK( s_mp_grow(sqr, ix) );

906 }

907 MP_USED(sqr) = ix;

908 MP_DIGIT(sqr, 0) = 0;

909

910 #ifdef NSS_USE_COMBA

911 if (IS_POWER_OF_2(MP_USED(a))) {

912 if (MP_USED(a) == 4) {

913 s_mp_sqr_comba_4(a, sqr);

914 goto CLEANUP;

915 }

916 if (MP_USED(a) == 8) {

917 s_mp_sqr_comba_8(a, sqr);

918 goto CLEANUP;

919 }

920 if (MP_USED(a) == 16) {

921 s_mp_sqr_comba_16(a, sqr);

922 goto CLEANUP;

923 }

924 if (MP_USED(a) == 32) {

925 s_mp_sqr_comba_32(a, sqr);

926 goto CLEANUP;

927 }

928 }

929 #endif

930

931 pa = MP_DIGITS(a);

932 count = MP_USED(a) - 1;

933 if (count > 0) {

934 d = *pa++;

935 s_mpv_mul_d(pa, count, d, MP_DIGITS(sqr) + 1);

936 for (ix = 3; --count > 0; ix += 2) {

937 d = *pa++;

938 s_mpv_mul_d_add(pa, count, d, MP_DIGITS(sqr) + ix);

939 } /* for(ix ...) */

940 MP_DIGIT(sqr, MP_USED(sqr)-1) = 0; /* above loop stopped short of this. */

941

942 /* now sqr = 2 /

943 s_mp_mul_2(sqr);

944 } else {

945 MP_DIGIT(sqr, 1) = 0;

946 }

947

948 /* now add the squares of the digits of a to sqr. */

949 s_mpv_sqr_add_prop(MP_DIGITS(a), MP_USED(a), MP_DIGITS(sqr));

950

951 SIGN(sqr) = ZPOS;

952 s_mp_clamp(sqr);

953

954 CLEANUP:

955 mp_clear(&tmp);

956 return res;

957

958 } /* end mp_sqr() */

959 #endif

960

961 /* }}} */

962

963 /* {{{ mp_div(a, b, q, r) */

964

965 /*

966 mp_div(a, b, q, r)

967

968 Compute q = a / b and r = a mod b. Input parameters may be re-used

969 as output parameters. If q or r is NULL, that portion of the

970 computation will be discarded (although it will still be computed)

971 */

972 mp_err mp_div(const mp_int a, const mp_int b, mp_int q, mp_int r)

973 {

974 mp_err res;

975 mp_int pQ, pR;

976 mp_int qtmp, rtmp, btmp;

977 int cmp;

978 mp_sign signA;

979 mp_sign signB;

980

981 ARGCHK(a != NULL && b != NULL, MP_BADARG);

982

983 signA = MP_SIGN(a);

984 signB = MP_SIGN(b);

985

986 if(mp_cmp_z(b) == MP_EQ)

987 return MP_RANGE;

988

989 DIGITS(&qtmp) = 0;

990 DIGITS(&rtmp) = 0;

991 DIGITS(&btmp) = 0;

992

993 /* Set up some temporaries... */

994 if (!r \|\| r == a \|\| r == b) {

995 MP_CHECKOK( mp_init_copy(&rtmp, a) );

996 pR = &rtmp;

997 } else {

998 MP_CHECKOK( mp_copy(a, r) );

999 pR = r;

1000 }

1001

1002 if (!q \|\| q == a \|\| q == b) {

1003 MP_CHECKOK( mp_init_size(&qtmp, MP_USED(a)) );

1004 pQ = &qtmp;

1005 } else {

1006 MP_CHECKOK( s_mp_pad(q, MP_USED(a)) );

1007 pQ = q;

1008 mp_zero(pQ);

1009 }

1010

1011 /*

1012 If \|a\| <= \|b\|, we can compute the solution without division;

1013 otherwise, we actually do the work required.

1014 */

1015 if ((cmp = s_mp_cmp(a, b)) <= 0) {

1016 if (cmp) {

1017 /* r was set to a above. */

1018 mp_zero(pQ);

1019 } else {

1020 mp_set(pQ, 1);

1021 mp_zero(pR);

1022 }

1023 } else {

1024 MP_CHECKOK( mp_init_copy(&btmp, b) );

1025 MP_CHECKOK( s_mp_div(pR, &btmp, pQ) );

1026 }

1027

1028 /* Compute the signs for the output */

1029 MP_SIGN(pR) = signA; /* Sr = Sa */

1030 /* Sq = ZPOS if Sa == Sb / / Sq = NEG if Sa != Sb */

1031 MP_SIGN(pQ) = (signA == signB) ? ZPOS : NEG;

1032

1033 if(s_mp_cmp_d(pQ, 0) == MP_EQ)

1034 SIGN(pQ) = ZPOS;

1035 if(s_mp_cmp_d(pR, 0) == MP_EQ)

1036 SIGN(pR) = ZPOS;

1037

1038 /* Copy output, if it is needed */

1039 if(q && q != pQ)

1040 s_mp_exch(pQ, q);

1041

1042 if(r && r != pR)

1043 s_mp_exch(pR, r);

1044

1045 CLEANUP:

1046 mp_clear(&btmp);

1047 mp_clear(&rtmp);

1048 mp_clear(&qtmp);

1049

1050 return res;

1051

1052 } /* end mp_div() */

1053

1054 /* }}} */

1055

1056 /* {{{ mp_div_2d(a, d, q, r) */

1057

1058 mp_err mp_div_2d(const mp_int a, mp_digit d, mp_int q, mp_int *r)

1059 {

1060 mp_err res;

1061

1062 ARGCHK(a != NULL, MP_BADARG);

1063

1064 if(q) {

1065 if((res = mp_copy(a, q)) != MP_OKAY)

1066 return res;

1067 }

1068 if(r) {

1069 if((res = mp_copy(a, r)) != MP_OKAY)

1070 return res;

1071 }

1072 if(q) {

1073 s_mp_div_2d(q, d);

1074 }

1075 if(r) {

1076 s_mp_mod_2d(r, d);

1077 }

1078

1079 return MP_OKAY;

1080

1081 } /* end mp_div_2d() */

1082

1083 /* }}} */

1084

1085 /* {{{ mp_expt(a, b, c) */

1086

1087 /*

1088 mp_expt(a, b, c)

1089

1090 Compute c = a ** b, that is, raise a to the b power. Uses a

1091 standard iterative square-and-multiply technique.

1092 */

1093

1094 mp_err mp_expt(mp_int a, mp_int b, mp_int *c)

1095 {

1096 mp_int s, x;

1097 mp_err res;

1098 mp_digit d;

1099 int dig, bit;

1100

1101 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

1102

1103 if(mp_cmp_z(b) < 0)

1104 return MP_RANGE;

1105

1106 if((res = mp_init(&s)) != MP_OKAY)

1107 return res;

1108

1109 mp_set(&s, 1);

1110

1111 if((res = mp_init_copy(&x, a)) != MP_OKAY)

1112 goto X;

1113

1114 /* Loop over low-order digits in ascending order */

1115 for(dig = 0; dig < (USED(b) - 1); dig++) {

1116 d = DIGIT(b, dig);

1117

1118 /* Loop over bits of each non-maximal digit */

1119 for(bit = 0; bit < DIGIT_BIT; bit++) {

1120 if(d & 1) {

1121 if((res = s_mp_mul(&s, &x)) != MP_OKAY)

1122 goto CLEANUP;

1123 }

1124

1125 d >>= 1;

1126

1127 if((res = s_mp_sqr(&x)) != MP_OKAY)

1128 goto CLEANUP;

1129 }

1130 }

1131

1132 /* Consider now the last digit... */

1133 d = DIGIT(b, dig);

1134

1135 while(d) {

1136 if(d & 1) {

1137 if((res = s_mp_mul(&s, &x)) != MP_OKAY)

1138 goto CLEANUP;

1139 }

1140

1141 d >>= 1;

1142

1143 if((res = s_mp_sqr(&x)) != MP_OKAY)

1144 goto CLEANUP;

1145 }

1146

1147 if(mp_iseven(b))

1148 SIGN(&s) = SIGN(a);

1149

1150 res = mp_copy(&s, c);

1151

1152 CLEANUP:

1153 mp_clear(&x);

1154 X:

1155 mp_clear(&s);

1156

1157 return res;

1158

1159 } /* end mp_expt() */

1160

1161 /* }}} */

1162

1163 /* {{{ mp_2expt(a, k) */

1164

1165 /* Compute a = 2^k */

1166

1167 mp_err mp_2expt(mp_int *a, mp_digit k)

1168 {

1169 ARGCHK(a != NULL, MP_BADARG);

1170

1171 return s_mp_2expt(a, k);

1172

1173 } /* end mp_2expt() */

1174

1175 /* }}} */

1176

1177 /* {{{ mp_mod(a, m, c) */

1178

1179 /*

1180 mp_mod(a, m, c)

1181

1182 Compute c = a (mod m). Result will always be 0 <= c < m.

1183 */

1184

1185 mp_err mp_mod(const mp_int a, const mp_int m, mp_int *c)

1186 {

1187 mp_err res;

1188 int mag;

1189

1190 ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);

1191

1192 if(SIGN(m) == NEG)

1193 return MP_RANGE;

1194

1195 /*

1196 If \|a\| > m, we need to divide to get the remainder and take the

1197 absolute value.

1198

1199 If \|a\| < m, we don't need to do any division, just copy and adjust

1200 the sign (if a is negative).

1201

1202 If \|a\| == m, we can simply set the result to zero.

1203

1204 This order is intended to minimize the average path length of the

1205 comparison chain on common workloads -- the most frequent cases are

1206 that \|a\| != m, so we do those first.

1207 */

1208 if((mag = s_mp_cmp(a, m)) > 0) {

1209 if((res = mp_div(a, m, NULL, c)) != MP_OKAY)

1210 return res;

1211

1212 if(SIGN(c) == NEG) {

1213 if((res = mp_add(c, m, c)) != MP_OKAY)

1214 return res;

1215 }

1216

1217 } else if(mag < 0) {

1218 if((res = mp_copy(a, c)) != MP_OKAY)

1219 return res;

1220

1221 if(mp_cmp_z(a) < 0) {

1222 if((res = mp_add(c, m, c)) != MP_OKAY)

1223 return res;

1224

1225 }

1226

1227 } else {

1228 mp_zero(c);

1229

1230 }

1231

1232 return MP_OKAY;

1233

1234 } /* end mp_mod() */

1235

1236 /* }}} */

1237

1238 /* {{{ mp_mod_d(a, d, c) */

1239

1240 /*

1241 mp_mod_d(a, d, c)

1242

1243 Compute c = a (mod d). Result will always be 0 <= c < d

1244 */

1245 mp_err mp_mod_d(const mp_int a, mp_digit d, mp_digit c)

1246 {

1247 mp_err res;

1248 mp_digit rem;

1249

1250 ARGCHK(a != NULL && c != NULL, MP_BADARG);

1251

1252 if(s_mp_cmp_d(a, d) > 0) {

1253 if((res = mp_div_d(a, d, NULL, &rem)) != MP_OKAY)

1254 return res;

1255

1256 } else {

1257 if(SIGN(a) == NEG)

1258 rem = d - DIGIT(a, 0);

1259 else

1260 rem = DIGIT(a, 0);

1261 }

1262

1263 if(c)

1264 *c = rem;

1265

1266 return MP_OKAY;

1267

1268 } /* end mp_mod_d() */

1269

1270 /* }}} */

1271

1272 /* {{{ mp_sqrt(a, b) */

1273

1274 /*

1275 mp_sqrt(a, b)

1276

1277 Compute the integer square root of a, and store the result in b.

1278 Uses an integer-arithmetic version of Newton's iterative linear

1279 approximation technique to determine this value; the result has the

1280 following two properties:

1281

1282 b^2 <= a

1283 (b+1)^2 >= a

1284

1285 It is a range error to pass a negative value.

1286 */

1287 mp_err mp_sqrt(const mp_int a, mp_int b)

1288 {

1289 mp_int x, t;

1290 mp_err res;

1291 mp_size used;

1292

1293 ARGCHK(a != NULL && b != NULL, MP_BADARG);

1294

1295 /* Cannot take square root of a negative value */

1296 if(SIGN(a) == NEG)

1297 return MP_RANGE;

1298

1299 /* Special cases for zero and one, trivial */

1300 if(mp_cmp_d(a, 1) <= 0)

1301 return mp_copy(a, b);

1302

1303 /* Initialize the temporaries we'll use below */

1304 if((res = mp_init_size(&t, USED(a))) != MP_OKAY)

1305 return res;

1306

1307 /* Compute an initial guess for the iteration as a itself */

1308 if((res = mp_init_copy(&x, a)) != MP_OKAY)

1309 goto X;

1310

1311 used = MP_USED(&x);

1312 if (used > 1) {

1313 s_mp_rshd(&x, used / 2);

1314 }

1315

1316 for(;;) {

1317 /* t = (x * x) - a */

1318 mp_copy(&x, &t); /* can't fail, t is big enough for original x */

1319 if((res = mp_sqr(&t, &t)) != MP_OKAY \|\|

1320 (res = mp_sub(&t, a, &t)) != MP_OKAY)

1321 goto CLEANUP;

1322

1323 /* t = t / 2x */

1324 s_mp_mul_2(&x);

1325 if((res = mp_div(&t, &x, &t, NULL)) != MP_OKAY)

1326 goto CLEANUP;

1327 s_mp_div_2(&x);

1328

1329 /* Terminate the loop, if the quotient is zero */

1330 if(mp_cmp_z(&t) == MP_EQ)

1331 break;

1332

1333 /* x = x - t */

1334 if((res = mp_sub(&x, &t, &x)) != MP_OKAY)

1335 goto CLEANUP;

1336

1337 }

1338

1339 /* Copy result to output parameter */

1340 mp_sub_d(&x, 1, &x);

1341 s_mp_exch(&x, b);

1342

1343 CLEANUP:

1344 mp_clear(&x);

1345 X:

1346 mp_clear(&t);

1347

1348 return res;

1349

1350 } /* end mp_sqrt() */

1351

1352 /* }}} */

1353

1354 /* }}} */

1355

1356 /------------------------------------------------------------------------/

1357 /* {{{ Modular arithmetic */

1358

1359 #if MP_MODARITH

1360 /* {{{ mp_addmod(a, b, m, c) */

1361

1362 /*

1363 mp_addmod(a, b, m, c)

1364

1365 Compute c = (a + b) mod m

1366 */

1367

1368 mp_err mp_addmod(const mp_int a, const mp_int b, const mp_int m, mp_int c)

1369 {

1370 mp_err res;

1371

1372 ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);

1373

1374 if((res = mp_add(a, b, c)) != MP_OKAY)

1375 return res;

1376 if((res = mp_mod(c, m, c)) != MP_OKAY)

1377 return res;

1378

1379 return MP_OKAY;

1380

1381 }

1382

1383 /* }}} */

1384

1385 /* {{{ mp_submod(a, b, m, c) */

1386

1387 /*

1388 mp_submod(a, b, m, c)

1389

1390 Compute c = (a - b) mod m

1391 */

1392

1393 mp_err mp_submod(const mp_int a, const mp_int b, const mp_int m, mp_int c)

1394 {

1395 mp_err res;

1396

1397 ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);

1398

1399 if((res = mp_sub(a, b, c)) != MP_OKAY)

1400 return res;

1401 if((res = mp_mod(c, m, c)) != MP_OKAY)

1402 return res;

1403

1404 return MP_OKAY;

1405

1406 }

1407

1408 /* }}} */

1409

1410 /* {{{ mp_mulmod(a, b, m, c) */

1411

1412 /*

1413 mp_mulmod(a, b, m, c)

1414

1415 Compute c = (a * b) mod m

1416 */

1417

1418 mp_err mp_mulmod(const mp_int a, const mp_int b, const mp_int m, mp_int c)

1419 {

1420 mp_err res;

1421

1422 ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);

1423

1424 if((res = mp_mul(a, b, c)) != MP_OKAY)

1425 return res;

1426 if((res = mp_mod(c, m, c)) != MP_OKAY)

1427 return res;

1428

1429 return MP_OKAY;

1430

1431 }

1432

1433 /* }}} */

1434

1435 /* {{{ mp_sqrmod(a, m, c) */

1436

1437 #if MP_SQUARE

1438 mp_err mp_sqrmod(const mp_int a, const mp_int m, mp_int *c)

1439 {

1440 mp_err res;

1441

1442 ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);

1443

1444 if((res = mp_sqr(a, c)) != MP_OKAY)

1445 return res;

1446 if((res = mp_mod(c, m, c)) != MP_OKAY)

1447 return res;

1448

1449 return MP_OKAY;

1450

1451 } /* end mp_sqrmod() */

1452 #endif

1453

1454 /* }}} */

1455

1456 /* {{{ s_mp_exptmod(a, b, m, c) */

1457

1458 /*

1459 s_mp_exptmod(a, b, m, c)

1460

1461 Compute c = (a ** b) mod m. Uses a standard square-and-multiply

1462 method with modular reductions at each step. (This is basically the

1463 same code as mp_expt(), except for the addition of the reductions)

1464

1465 The modular reductions are done using Barrett's algorithm (see

1466 s_mp_reduce() below for details)

1467 */

1468

1469 mp_err s_mp_exptmod(const mp_int a, const mp_int b, const mp_int m, mp_int c )

1470 {

1471 mp_int s, x, mu;

1472 mp_err res;

1473 mp_digit d;

1474 int dig, bit;

1475

1476 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

1477

1478 if(mp_cmp_z(b) < 0 \|\| mp_cmp_z(m) <= 0)

1479 return MP_RANGE;

1480

1481 if((res = mp_init(&s)) != MP_OKAY)

1482 return res;

1483 if((res = mp_init_copy(&x, a)) != MP_OKAY \|\|

1484 (res = mp_mod(&x, m, &x)) != MP_OKAY)

1485 goto X;

1486 if((res = mp_init(&mu)) != MP_OKAY)

1487 goto MU;

1488

1489 mp_set(&s, 1);

1490

1491 /* mu = b^2k / m */

1492 s_mp_add_d(&mu, 1);

1493 s_mp_lshd(&mu, 2 * USED(m));

1494 if((res = mp_div(&mu, m, &mu, NULL)) != MP_OKAY)

1495 goto CLEANUP;

1496

1497 /* Loop over digits of b in ascending order, except highest order */

1498 for(dig = 0; dig < (USED(b) - 1); dig++) {

1499 d = DIGIT(b, dig);

1500

1501 /* Loop over the bits of the lower-order digits */

1502 for(bit = 0; bit < DIGIT_BIT; bit++) {

1503 if(d & 1) {

1504 if((res = s_mp_mul(&s, &x)) != MP_OKAY)

1505 goto CLEANUP;

1506 if((res = s_mp_reduce(&s, m, &mu)) != MP_OKAY)

1507 goto CLEANUP;

1508 }

1509

1510 d >>= 1;

1511

1512 if((res = s_mp_sqr(&x)) != MP_OKAY)

1513 goto CLEANUP;

1514 if((res = s_mp_reduce(&x, m, &mu)) != MP_OKAY)

1515 goto CLEANUP;

1516 }

1517 }

1518

1519 /* Now do the last digit... */

1520 d = DIGIT(b, dig);

1521

1522 while(d) {

1523 if(d & 1) {

1524 if((res = s_mp_mul(&s, &x)) != MP_OKAY)

1525 goto CLEANUP;

1526 if((res = s_mp_reduce(&s, m, &mu)) != MP_OKAY)

1527 goto CLEANUP;

1528 }

1529

1530 d >>= 1;

1531

1532 if((res = s_mp_sqr(&x)) != MP_OKAY)

1533 goto CLEANUP;

1534 if((res = s_mp_reduce(&x, m, &mu)) != MP_OKAY)

1535 goto CLEANUP;

1536 }

1537

1538 s_mp_exch(&s, c);

1539

1540 CLEANUP:

1541 mp_clear(&mu);

1542 MU:

1543 mp_clear(&x);

1544 X:

1545 mp_clear(&s);

1546

1547 return res;

1548

1549 } /* end s_mp_exptmod() */

1550

1551 /* }}} */

1552

1553 /* {{{ mp_exptmod_d(a, d, m, c) */

1554

1555 mp_err mp_exptmod_d(const mp_int a, mp_digit d, const mp_int m, mp_int *c)

1556 {

1557 mp_int s, x;

1558 mp_err res;

1559

1560 ARGCHK(a != NULL && c != NULL, MP_BADARG);

1561

1562 if((res = mp_init(&s)) != MP_OKAY)

1563 return res;

1564 if((res = mp_init_copy(&x, a)) != MP_OKAY)

1565 goto X;

1566

1567 mp_set(&s, 1);

1568

1569 while(d != 0) {

1570 if(d & 1) {

1571 if((res = s_mp_mul(&s, &x)) != MP_OKAY \|\|

1572 (res = mp_mod(&s, m, &s)) != MP_OKAY)

1573 goto CLEANUP;

1574 }

1575

1576 d /= 2;

1577

1578 if((res = s_mp_sqr(&x)) != MP_OKAY \|\|

1579 (res = mp_mod(&x, m, &x)) != MP_OKAY)

1580 goto CLEANUP;

1581 }

1582

1583 s_mp_exch(&s, c);

1584

1585 CLEANUP:

1586 mp_clear(&x);

1587 X:

1588 mp_clear(&s);

1589

1590 return res;

1591

1592 } /* end mp_exptmod_d() */

1593

1594 /* }}} */

1595 #endif /* if MP_MODARITH */

1596

1597 /* }}} */

1598

1599 /------------------------------------------------------------------------/

1600 /* {{{ Comparison functions */

1601

1602 /* {{{ mp_cmp_z(a) */

1603

1604 /*

1605 mp_cmp_z(a)

1606

1607 Compare a <=> 0. Returns <0 if a<0, 0 if a=0, >0 if a>0.

1608 */

1609

1610 int mp_cmp_z(const mp_int *a)

1611 {

1612 if(SIGN(a) == NEG)

1613 return MP_LT;

1614 else if(USED(a) == 1 && DIGIT(a, 0) == 0)

1615 return MP_EQ;

1616 else

1617 return MP_GT;

1618

1619 } /* end mp_cmp_z() */

1620

1621 /* }}} */

1622

1623 /* {{{ mp_cmp_d(a, d) */

1624

1625 /*

1626 mp_cmp_d(a, d)

1627

1628 Compare a <=> d. Returns <0 if a<d, 0 if a=d, >0 if a>d

1629 */

1630

1631 int mp_cmp_d(const mp_int *a, mp_digit d)

1632 {

1633 ARGCHK(a != NULL, MP_EQ);

1634

1635 if(SIGN(a) == NEG)

1636 return MP_LT;

1637

1638 return s_mp_cmp_d(a, d);

1639

1640 } /* end mp_cmp_d() */

1641

1642 /* }}} */

1643

1644 /* {{{ mp_cmp(a, b) */

1645

1646 int mp_cmp(const mp_int a, const mp_int b)

1647 {

1648 ARGCHK(a != NULL && b != NULL, MP_EQ);

1649

1650 if(SIGN(a) == SIGN(b)) {

1651 int mag;

1652

1653 if((mag = s_mp_cmp(a, b)) == MP_EQ)

1654 return MP_EQ;

1655

1656 if(SIGN(a) == ZPOS)

1657 return mag;

1658 else

1659 return -mag;

1660

1661 } else if(SIGN(a) == ZPOS) {

1662 return MP_GT;

1663 } else {

1664 return MP_LT;

1665 }

1666

1667 } /* end mp_cmp() */

1668

1669 /* }}} */

1670

1671 /* {{{ mp_cmp_mag(a, b) */

1672

1673 /*

1674 mp_cmp_mag(a, b)

1675

1676 Compares \|a\| <=> \|b\|, and returns an appropriate comparison result

1677 */

1678

1679 int mp_cmp_mag(mp_int a, mp_int b)

1680 {

1681 ARGCHK(a != NULL && b != NULL, MP_EQ);

1682

1683 return s_mp_cmp(a, b);

1684

1685 } /* end mp_cmp_mag() */

1686

1687 /* }}} */

1688

1689 /* {{{ mp_cmp_int(a, z) */

1690

1691 /*

1692 This just converts z to an mp_int, and uses the existing comparison

1693 routines. This is sort of inefficient, but it's not clear to me how

1694 frequently this wil get used anyway. For small positive constants,

1695 you can always use mp_cmp_d(), and for zero, there is mp_cmp_z().

1696 */

1697 int mp_cmp_int(const mp_int *a, long z)

1698 {

1699 mp_int tmp;

1700 int out;

1701

1702 ARGCHK(a != NULL, MP_EQ);

1703

1704 mp_init(&tmp); mp_set_int(&tmp, z);

1705 out = mp_cmp(a, &tmp);

1706 mp_clear(&tmp);

1707

1708 return out;

1709

1710 } /* end mp_cmp_int() */

1711

1712 /* }}} */

1713

1714 /* {{{ mp_isodd(a) */

1715

1716 /*

1717 mp_isodd(a)

1718

1719 Returns a true (non-zero) value if a is odd, false (zero) otherwise.

1720 */

1721 int mp_isodd(const mp_int *a)

1722 {

1723 ARGCHK(a != NULL, 0);

1724

1725 return (int)(DIGIT(a, 0) & 1);

1726

1727 } /* end mp_isodd() */

1728

1729 /* }}} */

1730

1731 /* {{{ mp_iseven(a) */

1732

1733 int mp_iseven(const mp_int *a)

1734 {

1735 return !mp_isodd(a);

1736

1737 } /* end mp_iseven() */

1738

1739 /* }}} */

1740

1741 /* }}} */

1742

1743 /------------------------------------------------------------------------/

1744 /* {{{ Number theoretic functions */

1745

1746 #if MP_NUMTH

1747 /* {{{ mp_gcd(a, b, c) */

1748

1749 /*

1750 Like the old mp_gcd() function, except computes the GCD using the

1751 binary algorithm due to Josef Stein in 1961 (via Knuth).

1752 */

1753 mp_err mp_gcd(mp_int a, mp_int b, mp_int *c)

1754 {

1755 mp_err res;

1756 mp_int u, v, t;

1757 mp_size k = 0;

1758

1759 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

1760

1761 if(mp_cmp_z(a) == MP_EQ && mp_cmp_z(b) == MP_EQ)

1762 return MP_RANGE;

1763 if(mp_cmp_z(a) == MP_EQ) {

1764 return mp_copy(b, c);

1765 } else if(mp_cmp_z(b) == MP_EQ) {

1766 return mp_copy(a, c);

1767 }

1768

1769 if((res = mp_init(&t)) != MP_OKAY)

1770 return res;

1771 if((res = mp_init_copy(&u, a)) != MP_OKAY)

1772 goto U;

1773 if((res = mp_init_copy(&v, b)) != MP_OKAY)

1774 goto V;

1775

1776 SIGN(&u) = ZPOS;

1777 SIGN(&v) = ZPOS;

1778

1779 /* Divide out common factors of 2 until at least 1 of a, b is even */

1780 while(mp_iseven(&u) && mp_iseven(&v)) {

1781 s_mp_div_2(&u);

1782 s_mp_div_2(&v);

1783 ++k;

1784 }

1785

1786 /* Initialize t */

1787 if(mp_isodd(&u)) {

1788 if((res = mp_copy(&v, &t)) != MP_OKAY)

1789 goto CLEANUP;

1790

1791 /* t = -v */

1792 if(SIGN(&v) == ZPOS)

1793 SIGN(&t) = NEG;

1794 else

1795 SIGN(&t) = ZPOS;

1796

1797 } else {

1798 if((res = mp_copy(&u, &t)) != MP_OKAY)

1799 goto CLEANUP;

1800

1801 }

1802

1803 for(;;) {

1804 while(mp_iseven(&t)) {

1805 s_mp_div_2(&t);

1806 }

1807

1808 if(mp_cmp_z(&t) == MP_GT) {

1809 if((res = mp_copy(&t, &u)) != MP_OKAY)

1810 goto CLEANUP;

1811

1812 } else {

1813 if((res = mp_copy(&t, &v)) != MP_OKAY)

1814 goto CLEANUP;

1815

1816 /* v = -t */

1817 if(SIGN(&t) == ZPOS)

1818 SIGN(&v) = NEG;

1819 else

1820 SIGN(&v) = ZPOS;

1821 }

1822

1823 if((res = mp_sub(&u, &v, &t)) != MP_OKAY)

1824 goto CLEANUP;

1825

1826 if(s_mp_cmp_d(&t, 0) == MP_EQ)

1827 break;

1828 }

1829

1830 s_mp_2expt(&v, k); /* v = 2^k */

1831 res = mp_mul(&u, &v, c); /* c = u * v */

1832

1833 CLEANUP:

1834 mp_clear(&v);

1835 V:

1836 mp_clear(&u);

1837 U:

1838 mp_clear(&t);

1839

1840 return res;

1841

1842 } /* end mp_gcd() */

1843

1844 /* }}} */

1845

1846 /* {{{ mp_lcm(a, b, c) */

1847

1848 /* We compute the least common multiple using the rule:

1849

1850 ab = [a, b](a, b)

1851

1852 ... by computing the product, and dividing out the gcd.

1853 */

1854

1855 mp_err mp_lcm(mp_int a, mp_int b, mp_int *c)

1856 {

1857 mp_int gcd, prod;

1858 mp_err res;

1859

1860 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

1861

1862 /* Set up temporaries */

1863 if((res = mp_init(&gcd)) != MP_OKAY)

1864 return res;

1865 if((res = mp_init(&prod)) != MP_OKAY)

1866 goto GCD;

1867

1868 if((res = mp_mul(a, b, &prod)) != MP_OKAY)

1869 goto CLEANUP;

1870 if((res = mp_gcd(a, b, &gcd)) != MP_OKAY)

1871 goto CLEANUP;

1872

1873 res = mp_div(&prod, &gcd, c, NULL);

1874

1875 CLEANUP:

1876 mp_clear(&prod);

1877 GCD:

1878 mp_clear(&gcd);

1879

1880 return res;

1881

1882 } /* end mp_lcm() */

1883

1884 /* }}} */

1885

1886 /* {{{ mp_xgcd(a, b, g, x, y) */

1887

1888 /*

1889 mp_xgcd(a, b, g, x, y)

1890

1891 Compute g = (a, b) and values x and y satisfying Bezout's identity

1892 (that is, ax + by = g). This uses the binary extended GCD algorithm

1893 based on the Stein algorithm used for mp_gcd()

1894 See algorithm 14.61 in Handbook of Applied Cryptogrpahy.

1895 */

1896

1897 mp_err mp_xgcd(const mp_int a, const mp_int b, mp_int g, mp_int x, mp_int *y )

1898 {

1899 mp_int gx, xc, yc, u, v, A, B, C, D;

1900 mp_int *clean[9];

1901 mp_err res;

1902 int last = -1;

1903

1904 if(mp_cmp_z(b) == 0)

1905 return MP_RANGE;

1906

1907 /* Initialize all these variables we need */

1908 MP_CHECKOK( mp_init(&u) );

1909 clean[++last] = &u;

1910 MP_CHECKOK( mp_init(&v) );

1911 clean[++last] = &v;

1912 MP_CHECKOK( mp_init(&gx) );

1913 clean[++last] = &gx;

1914 MP_CHECKOK( mp_init(&A) );

1915 clean[++last] = &A;

1916 MP_CHECKOK( mp_init(&B) );

1917 clean[++last] = &B;

1918 MP_CHECKOK( mp_init(&C) );

1919 clean[++last] = &C;

1920 MP_CHECKOK( mp_init(&D) );

1921 clean[++last] = &D;

1922 MP_CHECKOK( mp_init_copy(&xc, a) );

1923 clean[++last] = &xc;

1924 mp_abs(&xc, &xc);

1925 MP_CHECKOK( mp_init_copy(&yc, b) );

1926 clean[++last] = &yc;

1927 mp_abs(&yc, &yc);

1928

1929 mp_set(&gx, 1);

1930

1931 /* Divide by two until at least one of them is odd */

1932 while(mp_iseven(&xc) && mp_iseven(&yc)) {

1933 mp_size nx = mp_trailing_zeros(&xc);

1934 mp_size ny = mp_trailing_zeros(&yc);

1935 mp_size n = MP_MIN(nx, ny);

1936 s_mp_div_2d(&xc,n);

1937 s_mp_div_2d(&yc,n);

1938 MP_CHECKOK( s_mp_mul_2d(&gx,n) );

1939 }

1940

1941 mp_copy(&xc, &u);

1942 mp_copy(&yc, &v);

1943 mp_set(&A, 1); mp_set(&D, 1);

1944

1945 /* Loop through binary GCD algorithm */

1946 do {

1947 while(mp_iseven(&u)) {

1948 s_mp_div_2(&u);

1949

1950 if(mp_iseven(&A) && mp_iseven(&B)) {

1951 s_mp_div_2(&A); s_mp_div_2(&B);

1952 } else {

1953 MP_CHECKOK( mp_add(&A, &yc, &A) );

1954 s_mp_div_2(&A);

1955 MP_CHECKOK( mp_sub(&B, &xc, &B) );

1956 s_mp_div_2(&B);

1957 }

1958 }

1959

1960 while(mp_iseven(&v)) {

1961 s_mp_div_2(&v);

1962

1963 if(mp_iseven(&C) && mp_iseven(&D)) {

1964 s_mp_div_2(&C); s_mp_div_2(&D);

1965 } else {

1966 MP_CHECKOK( mp_add(&C, &yc, &C) );

1967 s_mp_div_2(&C);

1968 MP_CHECKOK( mp_sub(&D, &xc, &D) );

1969 s_mp_div_2(&D);

1970 }

1971 }

1972

1973 if(mp_cmp(&u, &v) >= 0) {

1974 MP_CHECKOK( mp_sub(&u, &v, &u) );

1975 MP_CHECKOK( mp_sub(&A, &C, &A) );

1976 MP_CHECKOK( mp_sub(&B, &D, &B) );

1977 } else {

1978 MP_CHECKOK( mp_sub(&v, &u, &v) );

1979 MP_CHECKOK( mp_sub(&C, &A, &C) );

1980 MP_CHECKOK( mp_sub(&D, &B, &D) );

1981 }

1982 } while (mp_cmp_z(&u) != 0);

1983

1984 /* copy results to output */

1985 if(x)

1986 MP_CHECKOK( mp_copy(&C, x) );

1987

1988 if(y)

1989 MP_CHECKOK( mp_copy(&D, y) );

1990

1991 if(g)

1992 MP_CHECKOK( mp_mul(&gx, &v, g) );

1993

1994 CLEANUP:

1995 while(last >= 0)

1996 mp_clear(clean[last--]);

1997

1998 return res;

1999

2000 } /* end mp_xgcd() */

2001

2002 /* }}} */

2003

2004 mp_size mp_trailing_zeros(const mp_int *mp)

2005 {

2006 mp_digit d;

2007 mp_size n = 0;

2008 int ix;

2009

2010 if (!mp \|\| !MP_DIGITS(mp) \|\| !mp_cmp_z(mp))

2011 return n;

2012

2013 for (ix = 0; !(d = MP_DIGIT(mp,ix)) && (ix < MP_USED(mp)); ++ix)

2014 n += MP_DIGIT_BIT;

2015 if (!d)

2016 return 0; /* shouldn't happen, but ... */

2017 #if !defined(MP_USE_UINT_DIGIT)

2018 if (!(d & 0xffffffffU)) {

2019 d >>= 32;

2020 n += 32;

2021 }

2022 #endif

2023 if (!(d & 0xffffU)) {

2024 d >>= 16;

2025 n += 16;

2026 }

2027 if (!(d & 0xffU)) {

2028 d >>= 8;

2029 n += 8;

2030 }

2031 if (!(d & 0xfU)) {

2032 d >>= 4;

2033 n += 4;

2034 }

2035 if (!(d & 0x3U)) {

2036 d >>= 2;

2037 n += 2;

2038 }

2039 if (!(d & 0x1U)) {

2040 d >>= 1;

2041 n += 1;

2042 }

2043 #if MP_ARGCHK == 2

2044 assert(0 != (d & 1));

2045 #endif

2046 return n;

2047 }

2048

2049 /* Given a and prime p, computes c and k such that ac == 2*k (mod p).

2050 ** Returns k (positive) or error (negative).

2051 ** This technique from the paper "Fast Modular Reciprocals" (unpublished)

2052 ** by Richard Schroeppel (a.k.a. Captain Nemo).

2053 */

2054 mp_err s_mp_almost_inverse(const mp_int a, const mp_int p, mp_int *c)

2055 {

2056 mp_err res;

2057 mp_err k = 0;

2058 mp_int d, f, g;

2059

2060 ARGCHK(a && p && c, MP_BADARG);

2061

2062 MP_DIGITS(&d) = 0;

2063 MP_DIGITS(&f) = 0;

2064 MP_DIGITS(&g) = 0;

2065 MP_CHECKOK( mp_init(&d) );

2066 MP_CHECKOK( mp_init_copy(&f, a) ); /* f = a */

2067 MP_CHECKOK( mp_init_copy(&g, p) ); /* g = p */

2068

2069 mp_set(c, 1);

2070 mp_zero(&d);

2071

2072 if (mp_cmp_z(&f) == 0) {

2073 res = MP_UNDEF;

2074 } else

2075 for (;;) {

2076 int diff_sign;

2077 while (mp_iseven(&f)) {

2078 mp_size n = mp_trailing_zeros(&f);

2079 if (!n) {

2080 res = MP_UNDEF;

2081 goto CLEANUP;

2082 }

2083 s_mp_div_2d(&f, n);

2084 MP_CHECKOK( s_mp_mul_2d(&d, n) );

2085 k += n;

2086 }

2087 if (mp_cmp_d(&f, 1) == MP_EQ) { /* f == 1 */

2088 res = k;

2089 break;

2090 }

2091 diff_sign = mp_cmp(&f, &g);

2092 if (diff_sign < 0) { /* f < g */

2093 s_mp_exch(&f, &g);

2094 s_mp_exch(c, &d);

2095 } else if (diff_sign == 0) { /* f == g */

2096 res = MP_UNDEF; /* a and p are not relatively prime */

2097 break;

2098 }

2099 if ((MP_DIGIT(&f,0) % 4) == (MP_DIGIT(&g,0) % 4)) {

2100 MP_CHECKOK( mp_sub(&f, &g, &f) ); /* f = f - g */

2101 MP_CHECKOK( mp_sub(c, &d, c) ); /* c = c - d */

2102 } else {

2103 MP_CHECKOK( mp_add(&f, &g, &f) ); /* f = f + g */

2104 MP_CHECKOK( mp_add(c, &d, c) ); /* c = c + d */

2105 }

2106 }

2107 if (res >= 0) {

2108 while (MP_SIGN(c) != MP_ZPOS) {

2109 MP_CHECKOK( mp_add(c, p, c) );

2110 }

2111 res = k;

2112 }

2113

2114 CLEANUP:

2115 mp_clear(&d);

2116 mp_clear(&f);

2117 mp_clear(&g);

2118 return res;

2119 }

2120

2121 /* Compute T = (P ** -1) mod MP_RADIX. Also works for 16-bit mp_digits.

2122 ** This technique from the paper "Fast Modular Reciprocals" (unpublished)

2123 ** by Richard Schroeppel (a.k.a. Captain Nemo).

2124 */

2125 mp_digit s_mp_invmod_radix(mp_digit P)

2126 {

2127 mp_digit T = P;

2128 T = 2 - (P T);

2129 T = 2 - (P T);

2130 T = 2 - (P T);

2131 T = 2 - (P T);

2132 #if !defined(MP_USE_UINT_DIGIT)

2133 T = 2 - (P T);

2134 T = 2 - (P T);

2135 #endif

2136 return T;

2137 }

2138

2139 /* Given c, k, and prime p, where ac == 2*k (mod p),

2140 Compute x = (a -1) mod p. This is similar to Montgomery reduction.

2141 ** This technique from the paper "Fast Modular Reciprocals" (unpublished)

2142 ** by Richard Schroeppel (a.k.a. Captain Nemo).

2143 */

2144 mp_err s_mp_fixup_reciprocal(const mp_int c, const mp_int p, int k, mp_int *x )

2145 {

2146 int k_orig = k;

2147 mp_digit r;

2148 mp_size ix;

2149 mp_err res;

2150

2151 if (mp_cmp_z(c) < 0) { /* c < 0 */

2152 MP_CHECKOK( mp_add(c, p, x) ); /* x = c + p */

2153 } else {

2154 MP_CHECKOK( mp_copy(c, x) ); /* x = c */

2155 }

2156

2157 /* make sure x is large enough */

2158 ix = MP_HOWMANY(k, MP_DIGIT_BIT) + MP_USED(p) + 1;

2159 ix = MP_MAX(ix, MP_USED(x));

2160 MP_CHECKOK( s_mp_pad(x, ix) );

2161

2162 r = 0 - s_mp_invmod_radix(MP_DIGIT(p,0));

2163

2164 for (ix = 0; k > 0; ix++) {

2165 int j = MP_MIN(k, MP_DIGIT_BIT);

2166 mp_digit v = r * MP_DIGIT(x, ix);

2167 if (j < MP_DIGIT_BIT) {

2168 v &= ((mp_digit)1 << j) - 1; /* v = v mod (2 ** j) */

2169 }

2170 s_mp_mul_d_add_offset(p, v, x, ix); /* x += p * v * (RADIX ** ix) */

2171 k -= j;

2172 }

2173 s_mp_clamp(x);

2174 s_mp_div_2d(x, k_orig);

2175 res = MP_OKAY;

2176

2177 CLEANUP:

2178 return res;

2179 }

2180

2181 /* compute mod inverse using Schroeppel's method, only if m is odd */

2182 mp_err s_mp_invmod_odd_m(const mp_int a, const mp_int m, mp_int *c)

2183 {

2184 int k;

2185 mp_err res;

2186 mp_int x;

2187

2188 ARGCHK(a && m && c, MP_BADARG);

2189

2190 if(mp_cmp_z(a) == 0 \|\| mp_cmp_z(m) == 0)

2191 return MP_RANGE;

2192 if (mp_iseven(m))

2193 return MP_UNDEF;

2194

2195 MP_DIGITS(&x) = 0;

2196

2197 if (a == c) {

2198 if ((res = mp_init_copy(&x, a)) != MP_OKAY)

2199 return res;

2200 if (a == m)

2201 m = &x;

2202 a = &x;

2203 } else if (m == c) {

2204 if ((res = mp_init_copy(&x, m)) != MP_OKAY)

2205 return res;

2206 m = &x;

2207 } else {

2208 MP_DIGITS(&x) = 0;

2209 }

2210

2211 MP_CHECKOK( s_mp_almost_inverse(a, m, c) );

2212 k = res;

2213 MP_CHECKOK( s_mp_fixup_reciprocal(c, m, k, c) );

2214 CLEANUP:

2215 mp_clear(&x);

2216 return res;

2217 }

2218

2219 /* Known good algorithm for computing modular inverse. But slow. */

2220 mp_err mp_invmod_xgcd(const mp_int a, const mp_int m, mp_int *c)

2221 {

2222 mp_int g, x;

2223 mp_err res;

2224

2225 ARGCHK(a && m && c, MP_BADARG);

2226

2227 if(mp_cmp_z(a) == 0 \|\| mp_cmp_z(m) == 0)

2228 return MP_RANGE;

2229

2230 MP_DIGITS(&g) = 0;

2231 MP_DIGITS(&x) = 0;

2232 MP_CHECKOK( mp_init(&x) );

2233 MP_CHECKOK( mp_init(&g) );

2234

2235 MP_CHECKOK( mp_xgcd(a, m, &g, &x, NULL) );

2236

2237 if (mp_cmp_d(&g, 1) != MP_EQ) {

2238 res = MP_UNDEF;

2239 goto CLEANUP;

2240 }

2241

2242 res = mp_mod(&x, m, c);

2243 SIGN(c) = SIGN(a);

2244

2245 CLEANUP:

2246 mp_clear(&x);

2247 mp_clear(&g);

2248

2249 return res;

2250 }

2251

2252 /* modular inverse where modulus is 2*k. /

2253 /* c = a-1 mod 2k */

2254 mp_err s_mp_invmod_2d(const mp_int a, mp_size k, mp_int c)

2255 {

2256 mp_err res;

2257 mp_size ix = k + 4;

2258 mp_int t0, t1, val, tmp, two2k;

2259

2260 static const mp_digit d2 = 2;

2261 static const mp_int two = { MP_ZPOS, 1, 1, (mp_digit *)&d2 };

2262

2263 if (mp_iseven(a))

2264 return MP_UNDEF;

2265 if (k <= MP_DIGIT_BIT) {

2266 mp_digit i = s_mp_invmod_radix(MP_DIGIT(a,0));

2267 if (k < MP_DIGIT_BIT)

2268 i &= ((mp_digit)1 << k) - (mp_digit)1;

2269 mp_set(c, i);

2270 return MP_OKAY;

2271 }

2272 MP_DIGITS(&t0) = 0;

2273 MP_DIGITS(&t1) = 0;

2274 MP_DIGITS(&val) = 0;

2275 MP_DIGITS(&tmp) = 0;

2276 MP_DIGITS(&two2k) = 0;

2277 MP_CHECKOK( mp_init_copy(&val, a) );

2278 s_mp_mod_2d(&val, k);

2279 MP_CHECKOK( mp_init_copy(&t0, &val) );

2280 MP_CHECKOK( mp_init_copy(&t1, &t0) );

2281 MP_CHECKOK( mp_init(&tmp) );

2282 MP_CHECKOK( mp_init(&two2k) );

2283 MP_CHECKOK( s_mp_2expt(&two2k, k) );

2284 do {

2285 MP_CHECKOK( mp_mul(&val, &t1, &tmp) );

2286 MP_CHECKOK( mp_sub(&two, &tmp, &tmp) );

2287 MP_CHECKOK( mp_mul(&t1, &tmp, &t1) );

2288 s_mp_mod_2d(&t1, k);

2289 while (MP_SIGN(&t1) != MP_ZPOS) {

2290 MP_CHECKOK( mp_add(&t1, &two2k, &t1) );

2291 }

2292 if (mp_cmp(&t1, &t0) == MP_EQ)

2293 break;

2294 MP_CHECKOK( mp_copy(&t1, &t0) );

2295 } while (--ix > 0);

2296 if (!ix) {

2297 res = MP_UNDEF;

2298 } else {

2299 mp_exch(c, &t1);

2300 }

2301

2302 CLEANUP:

2303 mp_clear(&t0);

2304 mp_clear(&t1);

2305 mp_clear(&val);

2306 mp_clear(&tmp);

2307 mp_clear(&two2k);

2308 return res;

2309 }

2310

2311 mp_err s_mp_invmod_even_m(const mp_int a, const mp_int m, mp_int *c)

2312 {

2313 mp_err res;

2314 mp_size k;

2315 mp_int oddFactor, evenFactor; /* factors of the modulus */

2316 mp_int oddPart, evenPart; /* parts to combine via CRT. */

2317 mp_int C2, tmp1, tmp2;

2318

2319 /static const mp_digit d1 = 1; /

2320 /static const mp_int one = { MP_ZPOS, 1, 1, (mp_digit )&d1 }; */

2321

2322 if ((res = s_mp_ispow2(m)) >= 0) {

2323 k = res;

2324 return s_mp_invmod_2d(a, k, c);

2325 }

2326 MP_DIGITS(&oddFactor) = 0;

2327 MP_DIGITS(&evenFactor) = 0;

2328 MP_DIGITS(&oddPart) = 0;

2329 MP_DIGITS(&evenPart) = 0;

2330 MP_DIGITS(&C2) = 0;

2331 MP_DIGITS(&tmp1) = 0;

2332 MP_DIGITS(&tmp2) = 0;

2333

2334 MP_CHECKOK( mp_init_copy(&oddFactor, m) ); /* oddFactor = m */

2335 MP_CHECKOK( mp_init(&evenFactor) );

2336 MP_CHECKOK( mp_init(&oddPart) );

2337 MP_CHECKOK( mp_init(&evenPart) );

2338 MP_CHECKOK( mp_init(&C2) );

2339 MP_CHECKOK( mp_init(&tmp1) );

2340 MP_CHECKOK( mp_init(&tmp2) );

2341

2342 k = mp_trailing_zeros(m);

2343 s_mp_div_2d(&oddFactor, k);

2344 MP_CHECKOK( s_mp_2expt(&evenFactor, k) );

2345

2346 /* compute a*-1 mod oddFactor. /

2347 MP_CHECKOK( s_mp_invmod_odd_m(a, &oddFactor, &oddPart) );

2348 /* compute a-1 mod evenFactor, where evenFactor == 2k. */

2349 MP_CHECKOK( s_mp_invmod_2d( a, k, &evenPart) );

2350

2351 /* Use Chinese Remainer theorem to compute a*-1 mod m. /

2352 /* let m1 = oddFactor, v1 = oddPart,

2353 * let m2 = evenFactor, v2 = evenPart.

2354 */

2355

2356 /* Compute C2 = m1*-1 mod m2. /

2357 MP_CHECKOK( s_mp_invmod_2d(&oddFactor, k, &C2) );

2358

2359 /* compute u = (v2 - v1)C2 mod m2 /

2360 MP_CHECKOK( mp_sub(&evenPart, &oddPart, &tmp1) );

2361 MP_CHECKOK( mp_mul(&tmp1, &C2, &tmp2) );

2362 s_mp_mod_2d(&tmp2, k);

2363 while (MP_SIGN(&tmp2) != MP_ZPOS) {

2364 MP_CHECKOK( mp_add(&tmp2, &evenFactor, &tmp2) );

2365 }

2366

2367 /* compute answer = v1 + um1 /

2368 MP_CHECKOK( mp_mul(&tmp2, &oddFactor, c) );

2369 MP_CHECKOK( mp_add(&oddPart, c, c) );

2370 /* not sure this is necessary, but it's low cost if not. */

2371 MP_CHECKOK( mp_mod(c, m, c) );

2372

2373 CLEANUP:

2374 mp_clear(&oddFactor);

2375 mp_clear(&evenFactor);

2376 mp_clear(&oddPart);

2377 mp_clear(&evenPart);

2378 mp_clear(&C2);

2379 mp_clear(&tmp1);

2380 mp_clear(&tmp2);

2381 return res;

2382 }

2383

2384

2385 /* {{{ mp_invmod(a, m, c) */

2386

2387 /*

2388 mp_invmod(a, m, c)

2389

2390 Compute c = a^-1 (mod m), if there is an inverse for a (mod m).

2391 This is equivalent to the question of whether (a, m) = 1. If not,

2392 MP_UNDEF is returned, and there is no inverse.

2393 */

2394

2395 mp_err mp_invmod(const mp_int a, const mp_int m, mp_int *c)

2396 {

2397

2398 ARGCHK(a && m && c, MP_BADARG);

2399

2400 if(mp_cmp_z(a) == 0 \|\| mp_cmp_z(m) == 0)

2401 return MP_RANGE;

2402

2403 if (mp_isodd(m)) {

2404 return s_mp_invmod_odd_m(a, m, c);

2405 }

2406 if (mp_iseven(a))

2407 return MP_UNDEF; /* not invertable */

2408

2409 return s_mp_invmod_even_m(a, m, c);

2410

2411 } /* end mp_invmod() */

2412

2413 /* }}} */

2414 #endif /* if MP_NUMTH */

2415

2416 /* }}} */

2417

2418 /------------------------------------------------------------------------/

2419 /* {{{ mp_print(mp, ofp) */

2420

2421 #if MP_IOFUNC

2422 /*

2423 mp_print(mp, ofp)

2424

2425 Print a textual representation of the given mp_int on the output

2426 stream 'ofp'. Output is generated using the internal radix.

2427 */

2428

2429 void mp_print(mp_int mp, FILE ofp)

2430 {

2431 int ix;

2432

2433 if(mp == NULL \|\| ofp == NULL)

2434 return;

2435

2436 fputc((SIGN(mp) == NEG) ? '-' : '+', ofp);

2437

2438 for(ix = USED(mp) - 1; ix >= 0; ix--) {

2439 fprintf(ofp, DIGIT_FMT, DIGIT(mp, ix));

2440 }

2441

2442 } /* end mp_print() */

2443

2444 #endif /* if MP_IOFUNC */

2445

2446 /* }}} */

2447

2448 /------------------------------------------------------------------------/

2449 /* {{{ More I/O Functions */

2450

2451 /* {{{ mp_read_raw(mp, str, len) */

2452

2453 /*

2454 mp_read_raw(mp, str, len)

2455

2456 Read in a raw value (base 256) into the given mp_int

2457 */

2458

2459 mp_err mp_read_raw(mp_int mp, char str, int len)

2460 {

2461 int ix;

2462 mp_err res;

2463 unsigned char ustr = (unsigned char )str;

2464

2465 ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);

2466

2467 mp_zero(mp);

2468

2469 /* Get sign from first byte */

2470 if(ustr[0])

2471 SIGN(mp) = NEG;

2472 else

2473 SIGN(mp) = ZPOS;

2474

2475 /* Read the rest of the digits */

2476 for(ix = 1; ix < len; ix++) {

2477 if((res = mp_mul_d(mp, 256, mp)) != MP_OKAY)

2478 return res;

2479 if((res = mp_add_d(mp, ustr[ix], mp)) != MP_OKAY)

2480 return res;

2481 }

2482

2483 return MP_OKAY;

2484

2485 } /* end mp_read_raw() */

2486

2487 /* }}} */

2488

2489 /* {{{ mp_raw_size(mp) */

2490

2491 int mp_raw_size(mp_int *mp)

2492 {

2493 ARGCHK(mp != NULL, 0);

2494

2495 return (USED(mp) * sizeof(mp_digit)) + 1;

2496

2497 } /* end mp_raw_size() */

2498

2499 /* }}} */

2500

2501 /* {{{ mp_toraw(mp, str) */

2502

2503 mp_err mp_toraw(mp_int mp, char str)

2504 {

2505 int ix, jx, pos = 1;

2506

2507 ARGCHK(mp != NULL && str != NULL, MP_BADARG);

2508

2509 str[0] = (char)SIGN(mp);

2510

2511 /* Iterate over each digit... */

2512 for(ix = USED(mp) - 1; ix >= 0; ix--) {

2513 mp_digit d = DIGIT(mp, ix);

2514

2515 /* Unpack digit bytes, high order first */

2516 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {

2517 str[pos++] = (char)(d >> (jx * CHAR_BIT));

2518 }

2519 }

2520

2521 return MP_OKAY;

2522

2523 } /* end mp_toraw() */

2524

2525 /* }}} */

2526

2527 /* {{{ mp_read_radix(mp, str, radix) */

2528

2529 /*

2530 mp_read_radix(mp, str, radix)

2531

2532 Read an integer from the given string, and set mp to the resulting

2533 value. The input is presumed to be in base 10. Leading non-digit

2534 characters are ignored, and the function reads until a non-digit

2535 character or the end of the string.

2536 */

2537

2538 mp_err mp_read_radix(mp_int mp, const char str, int radix)

2539 {

2540 int ix = 0, val = 0;

2541 mp_err res;

2542 mp_sign sig = ZPOS;

2543

2544 ARGCHK(mp != NULL && str != NULL && radix >= 2 && radix <= MAX_RADIX,

2545 MP_BADARG);

2546

2547 mp_zero(mp);

2548

2549 /* Skip leading non-digit characters until a digit or '-' or '+' */

2550 while(str[ix] &&

2551 (s_mp_tovalue(str[ix], radix) < 0) &&

2552 str[ix] != '-' &&

2553 str[ix] != '+') {

2554 ++ix;

2555 }

2556

2557 if(str[ix] == '-') {

2558 sig = NEG;

2559 ++ix;

2560 } else if(str[ix] == '+') {

2561 sig = ZPOS; /* this is the default anyway... */

2562 ++ix;

2563 }

2564

2565 while((val = s_mp_tovalue(str[ix], radix)) >= 0) {

2566 if((res = s_mp_mul_d(mp, radix)) != MP_OKAY)

2567 return res;

2568 if((res = s_mp_add_d(mp, val)) != MP_OKAY)

2569 return res;

2570 ++ix;

2571 }

2572

2573 if(s_mp_cmp_d(mp, 0) == MP_EQ)

2574 SIGN(mp) = ZPOS;

2575 else

2576 SIGN(mp) = sig;

2577

2578 return MP_OKAY;

2579

2580 } /* end mp_read_radix() */

2581

2582 mp_err mp_read_variable_radix(mp_int a, const char str, int default_radix)

2583 {

2584 int radix = default_radix;

2585 int cx;

2586 mp_sign sig = ZPOS;

2587 mp_err res;

2588

2589 /* Skip leading non-digit characters until a digit or '-' or '+' */

2590 while ((cx = *str) != 0 &&

2591 (s_mp_tovalue(cx, radix) < 0) &&

2592 cx != '-' &&

2593 cx != '+') {

2594 ++str;

2595 }

2596

2597 if (cx == '-') {

2598 sig = NEG;

2599 ++str;

2600 } else if (cx == '+') {

2601 sig = ZPOS; /* this is the default anyway... */

2602 ++str;

2603 }

2604

2605 if (str[0] == '0') {

2606 if ((str[1] \| 0x20) == 'x') {

2607 radix = 16;

2608 str += 2;

2609 } else {

2610 radix = 8;

2611 str++;

2612 }

2613 }

2614 res = mp_read_radix(a, str, radix);

2615 if (res == MP_OKAY) {

2616 MP_SIGN(a) = (s_mp_cmp_d(a, 0) == MP_EQ) ? ZPOS : sig;

2617 }

2618 return res;

2619 }

2620

2621 /* }}} */

2622

2623 /* {{{ mp_radix_size(mp, radix) */

2624

2625 int mp_radix_size(mp_int *mp, int radix)

2626 {

2627 int bits;

2628

2629 if(!mp \|\| radix < 2 \|\| radix > MAX_RADIX)

2630 return 0;

2631

2632 bits = USED(mp) * DIGIT_BIT - 1;

2633

2634 return s_mp_outlen(bits, radix);

2635

2636 } /* end mp_radix_size() */

2637

2638 /* }}} */

2639

2640 /* {{{ mp_toradix(mp, str, radix) */

2641

2642 mp_err mp_toradix(mp_int mp, char str, int radix)

2643 {

2644 int ix, pos = 0;

2645

2646 ARGCHK(mp != NULL && str != NULL, MP_BADARG);

2647 ARGCHK(radix > 1 && radix <= MAX_RADIX, MP_RANGE);

2648

2649 if(mp_cmp_z(mp) == MP_EQ) {

2650 str[0] = '0';

2651 str[1] = '\0';

2652 } else {

2653 mp_err res;

2654 mp_int tmp;

2655 mp_sign sgn;

2656 mp_digit rem, rdx = (mp_digit)radix;

2657 char ch;

2658

2659 if((res = mp_init_copy(&tmp, mp)) != MP_OKAY)

2660 return res;

2661

2662 /* Save sign for later, and take absolute value */

2663 sgn = SIGN(&tmp); SIGN(&tmp) = ZPOS;

2664

2665 /* Generate output digits in reverse order */

2666 while(mp_cmp_z(&tmp) != 0) {

2667 if((res = mp_div_d(&tmp, rdx, &tmp, &rem)) != MP_OKAY) {

2668 mp_clear(&tmp);

2669 return res;

2670 }

2671

2672 /* Generate digits, use capital letters */

2673 ch = s_mp_todigit(rem, radix, 0);

2674

2675 str[pos++] = ch;

2676 }

2677

2678 /* Add - sign if original value was negative */

2679 if(sgn == NEG)

2680 str[pos++] = '-';

2681

2682 /* Add trailing NUL to end the string */

2683 str[pos--] = '\0';

2684

2685 /* Reverse the digits and sign indicator */

2686 ix = 0;

2687 while(ix < pos) {

2688 char tmp = str[ix];

2689

2690 str[ix] = str[pos];

2691 str[pos] = tmp;

2692 ++ix;

2693 --pos;

2694 }

2695

2696 mp_clear(&tmp);

2697 }

2698

2699 return MP_OKAY;

2700

2701 } /* end mp_toradix() */

2702

2703 /* }}} */

2704

2705 /* {{{ mp_tovalue(ch, r) */

2706

2707 int mp_tovalue(char ch, int r)

2708 {

2709 return s_mp_tovalue(ch, r);

2710

2711 } /* end mp_tovalue() */

2712

2713 /* }}} */

2714

2715 /* }}} */

2716

2717 /* {{{ mp_strerror(ec) */

2718

2719 /*

2720 mp_strerror(ec)

2721

2722 Return a string describing the meaning of error code 'ec'. The

2723 string returned is allocated in static memory, so the caller should

2724 not attempt to modify or free the memory associated with this

2725 string.

2726 */

2727 const char *mp_strerror(mp_err ec)

2728 {

2729 int aec = (ec < 0) ? -ec : ec;

2730

2731 /* Code values are negative, so the senses of these comparisons

2732 are accurate */

2733 if(ec < MP_LAST_CODE \|\| ec > MP_OKAY) {

2734 return mp_err_string[0]; /* unknown error code */

2735 } else {

2736 return mp_err_string[aec + 1];

2737 }

2738

2739 } /* end mp_strerror() */

2740

2741 /* }}} */

2742

2743 /========================================================================/

2744 /------------------------------------------------------------------------/

2745 /* Static function definitions (internal use only) */

2746

2747 /* {{{ Memory management */

2748

2749 /* {{{ s_mp_grow(mp, min) */

2750

2751 /* Make sure there are at least 'min' digits allocated to mp */

2752 mp_err s_mp_grow(mp_int *mp, mp_size min)

2753 {

2754 if(min > ALLOC(mp)) {

2755 mp_digit *tmp;

2756

2757 /* Set min to next nearest default precision block size */

2758 min = MP_ROUNDUP(min, s_mp_defprec);

2759

2760 if((tmp = s_mp_alloc(min, sizeof(mp_digit))) == NULL)

2761 return MP_MEM;

2762

2763 s_mp_copy(DIGITS(mp), tmp, USED(mp));

2764

2765 #if MP_CRYPTO

2766 s_mp_setz(DIGITS(mp), ALLOC(mp));

2767 #endif

2768 s_mp_free(DIGITS(mp));

2769 DIGITS(mp) = tmp;

2770 ALLOC(mp) = min;

2771 }

2772

2773 return MP_OKAY;

2774

2775 } /* end s_mp_grow() */

2776

2777 /* }}} */

2778

2779 /* {{{ s_mp_pad(mp, min) */

2780

2781 /* Make sure the used size of mp is at least 'min', growing if needed */

2782 mp_err s_mp_pad(mp_int *mp, mp_size min)

2783 {

2784 if(min > USED(mp)) {

2785 mp_err res;

2786

2787 /* Make sure there is room to increase precision */

2788 if (min > ALLOC(mp)) {

2789 if ((res = s_mp_grow(mp, min)) != MP_OKAY)

2790 return res;

2791 } else {

2792 s_mp_setz(DIGITS(mp) + USED(mp), min - USED(mp));

2793 }

2794

2795 /* Increase precision; should already be 0-filled */

2796 USED(mp) = min;

2797 }

2798

2799 return MP_OKAY;

2800

2801 } /* end s_mp_pad() */

2802

2803 /* }}} */

2804

2805 /* {{{ s_mp_setz(dp, count) */

2806

2807 #if MP_MACRO == 0

2808 /* Set 'count' digits pointed to by dp to be zeroes */

2809 void s_mp_setz(mp_digit *dp, mp_size count)

2810 {

2811 #if MP_MEMSET == 0

2812 int ix;

2813

2814 for(ix = 0; ix < count; ix++)

2815 dp[ix] = 0;

2816 #else

2817 memset(dp, 0, count * sizeof(mp_digit));

2818 #endif

2819

2820 } /* end s_mp_setz() */

2821 #endif

2822

2823 /* }}} */

2824

2825 /* {{{ s_mp_copy(sp, dp, count) */

2826

2827 #if MP_MACRO == 0

2828 /* Copy 'count' digits from sp to dp */

2829 void s_mp_copy(const mp_digit sp, mp_digit dp, mp_size count)

2830 {

2831 #if MP_MEMCPY == 0

2832 int ix;

2833

2834 for(ix = 0; ix < count; ix++)

2835 dp[ix] = sp[ix];

2836 #else

2837 memcpy(dp, sp, count * sizeof(mp_digit));

2838 #endif

2839 ++mp_copies;

2840

2841 } /* end s_mp_copy() */

2842 #endif

2843

2844 /* }}} */

2845

2846 /* {{{ s_mp_alloc(nb, ni) */

2847

2848 #if MP_MACRO == 0

2849 /* Allocate ni records of nb bytes each, and return a pointer to that */

2850 void *s_mp_alloc(size_t nb, size_t ni)

2851 {

2852 ++mp_allocs;

2853 return calloc(nb, ni);

2854

2855 } /* end s_mp_alloc() */

2856 #endif

2857

2858 /* }}} */

2859

2860 /* {{{ s_mp_free(ptr) */

2861

2862 #if MP_MACRO == 0

2863 /* Free the memory pointed to by ptr */

2864 void s_mp_free(void *ptr)

2865 {

2866 if(ptr) {

2867 ++mp_frees;

2868 free(ptr);

2869 }

2870 } /* end s_mp_free() */

2871 #endif

2872

2873 /* }}} */

2874

2875 /* {{{ s_mp_clamp(mp) */

2876

2877 #if MP_MACRO == 0

2878 /* Remove leading zeroes from the given value */

2879 void s_mp_clamp(mp_int *mp)

2880 {

2881 mp_size used = MP_USED(mp);

2882 while (used > 1 && DIGIT(mp, used - 1) == 0)

2883 --used;

2884 MP_USED(mp) = used;

2885 } /* end s_mp_clamp() */

2886 #endif

2887

2888 /* }}} */

2889

2890 /* {{{ s_mp_exch(a, b) */

2891

2892 /* Exchange the data for a and b; (b, a) = (a, b) */

2893 void s_mp_exch(mp_int a, mp_int b)

2894 {

2895 mp_int tmp;

2896

2897 tmp = *a;

2898 a = b;

2899 *b = tmp;

2900

2901 } /* end s_mp_exch() */

2902

2903 /* }}} */

2904

2905 /* }}} */

2906

2907 /* {{{ Arithmetic helpers */

2908

2909 /* {{{ s_mp_lshd(mp, p) */

2910

2911 /*

2912 Shift mp leftward by p digits, growing if needed, and zero-filling

2913 the in-shifted digits at the right end. This is a convenient

2914 alternative to multiplication by powers of the radix

2915 */

2916

2917 mp_err s_mp_lshd(mp_int *mp, mp_size p)

2918 {

2919 mp_err res;

2920 mp_size pos;

2921 int ix;

2922

2923 if(p == 0)

2924 return MP_OKAY;

2925

2926 if (MP_USED(mp) == 1 && MP_DIGIT(mp, 0) == 0)

2927 return MP_OKAY;

2928

2929 if((res = s_mp_pad(mp, USED(mp) + p)) != MP_OKAY)

2930 return res;

2931

2932 pos = USED(mp) - 1;

2933

2934 /* Shift all the significant figures over as needed */

2935 for(ix = pos - p; ix >= 0; ix--)

2936 DIGIT(mp, ix + p) = DIGIT(mp, ix);

2937

2938 /* Fill the bottom digits with zeroes */

2939 for(ix = 0; ix < p; ix++)

2940 DIGIT(mp, ix) = 0;

2941

2942 return MP_OKAY;

2943

2944 } /* end s_mp_lshd() */

2945

2946 /* }}} */

2947

2948 /* {{{ s_mp_mul_2d(mp, d) */

2949

2950 /*

2951 Multiply the integer by 2^d, where d is a number of bits. This

2952 amounts to a bitwise shift of the value.

2953 */

2954 mp_err s_mp_mul_2d(mp_int *mp, mp_digit d)

2955 {

2956 mp_err res;

2957 mp_digit dshift, bshift;

2958 mp_digit mask;

2959

2960 ARGCHK(mp != NULL, MP_BADARG);

2961

2962 dshift = d / MP_DIGIT_BIT;

2963 bshift = d % MP_DIGIT_BIT;

2964 /* bits to be shifted out of the top word */

2965 mask = ((mp_digit)~0 << (MP_DIGIT_BIT - bshift));

2966 mask &= MP_DIGIT(mp, MP_USED(mp) - 1);

2967

2968 if (MP_OKAY != (res = s_mp_pad(mp, MP_USED(mp) + dshift + (mask != 0) )))

2969 return res;

2970

2971 if (dshift && MP_OKAY != (res = s_mp_lshd(mp, dshift)))

2972 return res;

2973

2974 if (bshift) {

2975 mp_digit *pa = MP_DIGITS(mp);

2976 mp_digit *alim = pa + MP_USED(mp);

2977 mp_digit prev = 0;

2978

2979 for (pa += dshift; pa < alim; ) {

2980 mp_digit x = *pa;

2981 *pa++ = (x << bshift) \| prev;

2982 prev = x >> (DIGIT_BIT - bshift);

2983 }

2984 }

2985

2986 s_mp_clamp(mp);

2987 return MP_OKAY;

2988 } /* end s_mp_mul_2d() */

2989

2990 /* {{{ s_mp_rshd(mp, p) */

2991

2992 /*

2993 Shift mp rightward by p digits. Maintains the invariant that

2994 digits above the precision are all zero. Digits shifted off the

2995 end are lost. Cannot fail.

2996 */

2997

2998 void s_mp_rshd(mp_int *mp, mp_size p)

2999 {

3000 mp_size ix;

3001 mp_digit src, dst;

3002

3003 if(p == 0)

3004 return;

3005

3006 /* Shortcut when all digits are to be shifted off */

3007 if(p >= USED(mp)) {

3008 s_mp_setz(DIGITS(mp), ALLOC(mp));

3009 USED(mp) = 1;

3010 SIGN(mp) = ZPOS;

3011 return;

3012 }

3013

3014 /* Shift all the significant figures over as needed */

3015 dst = MP_DIGITS(mp);

3016 src = dst + p;

3017 for (ix = USED(mp) - p; ix > 0; ix--)

3018 dst++ = src++;

3019

3020 MP_USED(mp) -= p;

3021 /* Fill the top digits with zeroes */

3022 while (p-- > 0)

3023 *dst++ = 0;

3024

3025 #if 0

3026 /* Strip off any leading zeroes */

3027 s_mp_clamp(mp);

3028 #endif

3029

3030 } /* end s_mp_rshd() */

3031

3032 /* }}} */

3033

3034 /* {{{ s_mp_div_2(mp) */

3035

3036 /* Divide by two -- take advantage of radix properties to do it fast */

3037 void s_mp_div_2(mp_int *mp)

3038 {

3039 s_mp_div_2d(mp, 1);

3040

3041 } /* end s_mp_div_2() */

3042

3043 /* }}} */

3044

3045 /* {{{ s_mp_mul_2(mp) */

3046

3047 mp_err s_mp_mul_2(mp_int *mp)

3048 {

3049 mp_digit *pd;

3050 int ix, used;

3051 mp_digit kin = 0;

3052

3053 /* Shift digits leftward by 1 bit */

3054 used = MP_USED(mp);

3055 pd = MP_DIGITS(mp);

3056 for (ix = 0; ix < used; ix++) {

3057 mp_digit d = *pd;

3058 *pd++ = (d << 1) \| kin;

3059 kin = (d >> (DIGIT_BIT - 1));

3060 }

3061

3062 /* Deal with rollover from last digit */

3063 if (kin) {

3064 if (ix >= ALLOC(mp)) {

3065 mp_err res;

3066 if((res = s_mp_grow(mp, ALLOC(mp) + 1)) != MP_OKAY)

3067 return res;

3068 }

3069

3070 DIGIT(mp, ix) = kin;

3071 USED(mp) += 1;

3072 }

3073

3074 return MP_OKAY;

3075

3076 } /* end s_mp_mul_2() */

3077

3078 /* }}} */

3079

3080 /* {{{ s_mp_mod_2d(mp, d) */

3081

3082 /*

3083 Remainder the integer by 2^d, where d is a number of bits. This

3084 amounts to a bitwise AND of the value, and does not require the full

3085 division code

3086 */

3087 void s_mp_mod_2d(mp_int *mp, mp_digit d)

3088 {

3089 mp_size ndig = (d / DIGIT_BIT), nbit = (d % DIGIT_BIT);

3090 mp_size ix;

3091 mp_digit dmask;

3092

3093 if(ndig >= USED(mp))

3094 return;

3095

3096 /* Flush all the bits above 2^d in its digit */

3097 dmask = ((mp_digit)1 << nbit) - 1;

3098 DIGIT(mp, ndig) &= dmask;

3099

3100 /* Flush all digits above the one with 2^d in it */

3101 for(ix = ndig + 1; ix < USED(mp); ix++)

3102 DIGIT(mp, ix) = 0;

3103

3104 s_mp_clamp(mp);

3105

3106 } /* end s_mp_mod_2d() */

3107

3108 /* }}} */

3109

3110 /* {{{ s_mp_div_2d(mp, d) */

3111

3112 /*

3113 Divide the integer by 2^d, where d is a number of bits. This

3114 amounts to a bitwise shift of the value, and does not require the

3115 full division code (used in Barrett reduction, see below)

3116 */

3117 void s_mp_div_2d(mp_int *mp, mp_digit d)

3118 {

3119 int ix;

3120 mp_digit save, next, mask;

3121

3122 s_mp_rshd(mp, d / DIGIT_BIT);

3123 d %= DIGIT_BIT;

3124 if (d) {

3125 mask = ((mp_digit)1 << d) - 1;

3126 save = 0;

3127 for(ix = USED(mp) - 1; ix >= 0; ix--) {

3128 next = DIGIT(mp, ix) & mask;

3129 DIGIT(mp, ix) = (DIGIT(mp, ix) >> d) \| (save << (DIGIT_BIT - d));

3130 save = next;

3131 }

3132 }

3133 s_mp_clamp(mp);

3134

3135 } /* end s_mp_div_2d() */

3136

3137 /* }}} */

3138

3139 /* {{{ s_mp_norm(a, b, d) /

3140

3141 /*

3142 s_mp_norm(a, b, *d)

3143

3144 Normalize a and b for division, where b is the divisor. In order

3145 that we might make good guesses for quotient digits, we want the

3146 leading digit of b to be at least half the radix, which we

3147 accomplish by multiplying a and b by a power of 2. The exponent

3148 (shift count) is placed in *pd, so that the remainder can be shifted

3149 back at the end of the division process.

3150 */

3151

3152 mp_err s_mp_norm(mp_int a, mp_int b, mp_digit *pd)

3153 {

3154 mp_digit d;

3155 mp_digit mask;

3156 mp_digit b_msd;

3157 mp_err res = MP_OKAY;

3158

3159 d = 0;

3160 mask = DIGIT_MAX & ~(DIGIT_MAX >> 1); /* mask is msb of digit */

3161 b_msd = DIGIT(b, USED(b) - 1);

3162 while (!(b_msd & mask)) {

3163 b_msd <<= 1;

3164 ++d;

3165 }

3166

3167 if (d) {

3168 MP_CHECKOK( s_mp_mul_2d(a, d) );

3169 MP_CHECKOK( s_mp_mul_2d(b, d) );

3170 }

3171

3172 *pd = d;

3173 CLEANUP:

3174 return res;

3175

3176 } /* end s_mp_norm() */

3177

3178 /* }}} */

3179

3180 /* }}} */

3181

3182 /* {{{ Primitive digit arithmetic */

3183

3184 /* {{{ s_mp_add_d(mp, d) */

3185

3186 /* Add d to \|mp\| in place */

3187 mp_err s_mp_add_d(mp_int mp, mp_digit d) / unsigned digit addition */

3188 {

3189 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3190 mp_word w, k = 0;

3191 mp_size ix = 1;

3192

3193 w = (mp_word)DIGIT(mp, 0) + d;

3194 DIGIT(mp, 0) = ACCUM(w);

3195 k = CARRYOUT(w);

3196

3197 while(ix < USED(mp) && k) {

3198 w = (mp_word)DIGIT(mp, ix) + k;

3199 DIGIT(mp, ix) = ACCUM(w);

3200 k = CARRYOUT(w);

3201 ++ix;

3202 }

3203

3204 if(k != 0) {

3205 mp_err res;

3206

3207 if((res = s_mp_pad(mp, USED(mp) + 1)) != MP_OKAY)

3208 return res;

3209

3210 DIGIT(mp, ix) = (mp_digit)k;

3211 }

3212

3213 return MP_OKAY;

3214 #else

3215 mp_digit * pmp = MP_DIGITS(mp);

3216 mp_digit sum, mp_i, carry = 0;

3217 mp_err res = MP_OKAY;

3218 int used = (int)MP_USED(mp);

3219

3220 mp_i = *pmp;

3221 *pmp++ = sum = d + mp_i;

3222 carry = (sum < d);

3223 while (carry && --used > 0) {

3224 mp_i = *pmp;

3225 *pmp++ = sum = carry + mp_i;

3226 carry = !sum;

3227 }

3228 if (carry && !used) {

3229 /* mp is growing */

3230 used = MP_USED(mp);

3231 MP_CHECKOK( s_mp_pad(mp, used + 1) );

3232 MP_DIGIT(mp, used) = carry;

3233 }

3234 CLEANUP:

3235 return res;

3236 #endif

3237 } /* end s_mp_add_d() */

3238

3239 /* }}} */

3240

3241 /* {{{ s_mp_sub_d(mp, d) */

3242

3243 /* Subtract d from \|mp\| in place, assumes \|mp\| > d */

3244 mp_err s_mp_sub_d(mp_int mp, mp_digit d) / unsigned digit subtract */

3245 {

3246 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)

3247 mp_word w, b = 0;

3248 mp_size ix = 1;

3249

3250 /* Compute initial subtraction */

3251 w = (RADIX + (mp_word)DIGIT(mp, 0)) - d;

3252 b = CARRYOUT(w) ? 0 : 1;

3253 DIGIT(mp, 0) = ACCUM(w);

3254

3255 /* Propagate borrows leftward */

3256 while(b && ix < USED(mp)) {

3257 w = (RADIX + (mp_word)DIGIT(mp, ix)) - b;

3258 b = CARRYOUT(w) ? 0 : 1;

3259 DIGIT(mp, ix) = ACCUM(w);

3260 ++ix;

3261 }

3262

3263 /* Remove leading zeroes */

3264 s_mp_clamp(mp);

3265

3266 /* If we have a borrow out, it's a violation of the input invariant */

3267 if(b)

3268 return MP_RANGE;

3269 else

3270 return MP_OKAY;

3271 #else

3272 mp_digit *pmp = MP_DIGITS(mp);

3273 mp_digit mp_i, diff, borrow;

3274 mp_size used = MP_USED(mp);

3275

3276 mp_i = *pmp;

3277 *pmp++ = diff = mp_i - d;

3278 borrow = (diff > mp_i);

3279 while (borrow && --used) {

3280 mp_i = *pmp;

3281 *pmp++ = diff = mp_i - borrow;

3282 borrow = (diff > mp_i);

3283 }

3284 s_mp_clamp(mp);

3285 return (borrow && !used) ? MP_RANGE : MP_OKAY;

3286 #endif

3287 } /* end s_mp_sub_d() */

3288

3289 /* }}} */

3290

3291 /* {{{ s_mp_mul_d(a, d) */

3292

3293 /* Compute a = a * d, single digit multiplication */

3294 mp_err s_mp_mul_d(mp_int *a, mp_digit d)

3295 {

3296 mp_err res;

3297 mp_size used;

3298 int pow;

3299

3300 if (!d) {

3301 mp_zero(a);

3302 return MP_OKAY;

3303 }

3304 if (d == 1)

3305 return MP_OKAY;

3306 if (0 <= (pow = s_mp_ispow2d(d))) {

3307 return s_mp_mul_2d(a, (mp_digit)pow);

3308 }

3309

3310 used = MP_USED(a);

3311 MP_CHECKOK( s_mp_pad(a, used + 1) );

3312

3313 s_mpv_mul_d(MP_DIGITS(a), used, d, MP_DIGITS(a));

3314

3315 s_mp_clamp(a);

3316

3317 CLEANUP:

3318 return res;

3319

3320 } /* end s_mp_mul_d() */

3321

3322 /* }}} */

3323

3324 /* {{{ s_mp_div_d(mp, d, r) */

3325

3326 /*

3327 s_mp_div_d(mp, d, r)

3328

3329 Compute the quotient mp = mp / d and remainder r = mp mod d, for a

3330 single digit d. If r is null, the remainder will be discarded.

3331 */

3332

3333 mp_err s_mp_div_d(mp_int mp, mp_digit d, mp_digit r)

3334 {

3335 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)

3336 mp_word w = 0, q;

3337 #else

3338 mp_digit w, q;

3339 #endif

3340 int ix;

3341 mp_err res;

3342 mp_int quot;

3343 mp_int rem;

3344

3345 if(d == 0)

3346 return MP_RANGE;

3347 if (d == 1) {

3348 if (r)

3349 *r = 0;

3350 return MP_OKAY;

3351 }

3352 /* could check for power of 2 here, but mp_div_d does that. */

3353 if (MP_USED(mp) == 1) {

3354 mp_digit n = MP_DIGIT(mp,0);

3355 mp_digit rem;

3356

3357 q = n / d;

3358 rem = n % d;

3359 MP_DIGIT(mp,0) = q;

3360 if (r)

3361 *r = rem;

3362 return MP_OKAY;

3363 }

3364

3365 MP_DIGITS(&rem) = 0;

3366 MP_DIGITS(&quot) = 0;

3367 /* Make room for the quotient */

3368 MP_CHECKOK( mp_init_size(&quot, USED(mp)) );

3369

3370 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)

3371 for(ix = USED(mp) - 1; ix >= 0; ix--) {

3372 w = (w << DIGIT_BIT) \| DIGIT(mp, ix);

3373

3374 if(w >= d) {

3375 q = w / d;

3376 w = w % d;

3377 } else {

3378 q = 0;

3379 }

3380

3381 s_mp_lshd(&quot, 1);

3382 DIGIT(&quot, 0) = (mp_digit)q;

3383 }

3384 #else

3385 {

3386 mp_digit p;

3387 #if !defined(MP_ASSEMBLY_DIV_2DX1D)

3388 mp_digit norm;

3389 #endif

3390

3391 MP_CHECKOK( mp_init_copy(&rem, mp) );

3392

3393 #if !defined(MP_ASSEMBLY_DIV_2DX1D)

3394 MP_DIGIT(&quot, 0) = d;

3395 MP_CHECKOK( s_mp_norm(&rem, &quot, &norm) );

3396 if (norm)

3397 d <<= norm;

3398 MP_DIGIT(&quot, 0) = 0;

3399 #endif

3400

3401 p = 0;

3402 for (ix = USED(&rem) - 1; ix >= 0; ix--) {

3403 w = DIGIT(&rem, ix);

3404

3405 if (p) {

3406 MP_CHECKOK( s_mpv_div_2dx1d(p, w, d, &q, &w) );

3407 } else if (w >= d) {

3408 q = w / d;

3409 w = w % d;

3410 } else {

3411 q = 0;

3412 }

3413

3414 MP_CHECKOK( s_mp_lshd(&quot, 1) );

3415 DIGIT(&quot, 0) = q;

3416 p = w;

3417 }

3418 #if !defined(MP_ASSEMBLY_DIV_2DX1D)

3419 if (norm)

3420 w >>= norm;

3421 #endif

3422 }

3423 #endif

3424

3425 /* Deliver the remainder, if desired */

3426 if(r)

3427 *r = (mp_digit)w;

3428

3429 s_mp_clamp(&quot);

3430 mp_exch(&quot, mp);

3431 CLEANUP:

3432 mp_clear(&quot);

3433 mp_clear(&rem);

3434

3435 return res;

3436 } /* end s_mp_div_d() */

3437

3438 /* }}} */

3439

3440

3441 /* }}} */

3442

3443 /* {{{ Primitive full arithmetic */

3444

3445 /* {{{ s_mp_add(a, b) */

3446

3447 /* Compute a = \|a\| + \|b\| */

3448 mp_err s_mp_add(mp_int a, const mp_int b) /* magnitude addition */

3449 {

3450 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3451 mp_word w = 0;

3452 #else

3453 mp_digit d, sum, carry = 0;

3454 #endif

3455 mp_digit pa, pb;

3456 mp_size ix;

3457 mp_size used;

3458 mp_err res;

3459

3460 /* Make sure a has enough precision for the output value */

3461 if((USED(b) > USED(a)) && (res = s_mp_pad(a, USED(b))) != MP_OKAY)

3462 return res;

3463

3464 /*

3465 Add up all digits up to the precision of b. If b had initially

3466 the same precision as a, or greater, we took care of it by the

3467 padding step above, so there is no problem. If b had initially

3468 less precision, we'll have to make sure the carry out is duly

3469 propagated upward among the higher-order digits of the sum.

3470 */

3471 pa = MP_DIGITS(a);

3472 pb = MP_DIGITS(b);

3473 used = MP_USED(b);

3474 for(ix = 0; ix < used; ix++) {

3475 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3476 w = w + pa + pb++;

3477 *pa++ = ACCUM(w);

3478 w = CARRYOUT(w);

3479 #else

3480 d = *pa;

3481 sum = d + *pb++;

3482 d = (sum < d); /* detect overflow */

3483 *pa++ = sum += carry;

3484 carry = d + (sum < carry); /* detect overflow */

3485 #endif

3486 }

3487

3488 /* If we run out of 'b' digits before we're actually done, make

3489 sure the carries get propagated upward...

3490 */

3491 used = MP_USED(a);

3492 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3493 while (w && ix < used) {

3494 w = w + *pa;

3495 *pa++ = ACCUM(w);

3496 w = CARRYOUT(w);

3497 ++ix;

3498 }

3499 #else

3500 while (carry && ix < used) {

3501 sum = carry + *pa;

3502 *pa++ = sum;

3503 carry = !sum;

3504 ++ix;

3505 }

3506 #endif

3507

3508 /* If there's an overall carry out, increase precision and include

3509 it. We could have done this initially, but why touch the memory

3510 allocator unless we're sure we have to?

3511 */

3512 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3513 if (w) {

3514 if((res = s_mp_pad(a, used + 1)) != MP_OKAY)

3515 return res;

3516

3517 DIGIT(a, ix) = (mp_digit)w;

3518 }

3519 #else

3520 if (carry) {

3521 if((res = s_mp_pad(a, used + 1)) != MP_OKAY)

3522 return res;

3523

3524 DIGIT(a, used) = carry;

3525 }

3526 #endif

3527

3528 return MP_OKAY;

3529 } /* end s_mp_add() */

3530

3531 /* }}} */

3532

3533 /* Compute c = \|a\| + \|b\| / / magnitude addition */

3534 mp_err s_mp_add_3arg(const mp_int a, const mp_int b, mp_int *c)

3535 {

3536 mp_digit pa, pb, *pc;

3537 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3538 mp_word w = 0;

3539 #else

3540 mp_digit sum, carry = 0, d;

3541 #endif

3542 mp_size ix;

3543 mp_size used;

3544 mp_err res;

3545

3546 MP_SIGN(c) = MP_SIGN(a);

3547 if (MP_USED(a) < MP_USED(b)) {

3548 const mp_int *xch = a;

3549 a = b;

3550 b = xch;

3551 }

3552

3553 /* Make sure a has enough precision for the output value */

3554 if (MP_OKAY != (res = s_mp_pad(c, MP_USED(a))))

3555 return res;

3556

3557 /*

3558 Add up all digits up to the precision of b. If b had initially

3559 the same precision as a, or greater, we took care of it by the

3560 exchange step above, so there is no problem. If b had initially

3561 less precision, we'll have to make sure the carry out is duly

3562 propagated upward among the higher-order digits of the sum.

3563 */

3564 pa = MP_DIGITS(a);

3565 pb = MP_DIGITS(b);

3566 pc = MP_DIGITS(c);

3567 used = MP_USED(b);

3568 for (ix = 0; ix < used; ix++) {

3569 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3570 w = w + pa++ + pb++;

3571 *pc++ = ACCUM(w);

3572 w = CARRYOUT(w);

3573 #else

3574 d = *pa++;

3575 sum = d + *pb++;

3576 d = (sum < d); /* detect overflow */

3577 *pc++ = sum += carry;

3578 carry = d + (sum < carry); /* detect overflow */

3579 #endif

3580 }

3581

3582 /* If we run out of 'b' digits before we're actually done, make

3583 sure the carries get propagated upward...

3584 */

3585 for (used = MP_USED(a); ix < used; ++ix) {

3586 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3587 w = w + *pa++;

3588 *pc++ = ACCUM(w);

3589 w = CARRYOUT(w);

3590 #else

3591 pc++ = sum = carry + pa++;

3592 carry = (sum < carry);

3593 #endif

3594 }

3595

3596 /* If there's an overall carry out, increase precision and include

3597 it. We could have done this initially, but why touch the memory

3598 allocator unless we're sure we have to?

3599 */

3600 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3601 if (w) {

3602 if((res = s_mp_pad(c, used + 1)) != MP_OKAY)

3603 return res;

3604

3605 DIGIT(c, used) = (mp_digit)w;

3606 ++used;

3607 }

3608 #else

3609 if (carry) {

3610 if((res = s_mp_pad(c, used + 1)) != MP_OKAY)

3611 return res;

3612

3613 DIGIT(c, used) = carry;

3614 ++used;

3615 }

3616 #endif

3617 MP_USED(c) = used;

3618 return MP_OKAY;

3619 }

3620 /* {{{ s_mp_add_offset(a, b, offset) */

3621

3622 /* Compute a = \|a\| + ( \|b\| * (RADIX ** offset) ) */

3623 mp_err s_mp_add_offset(mp_int a, mp_int b, mp_size offset)

3624 {

3625 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3626 mp_word w, k = 0;

3627 #else

3628 mp_digit d, sum, carry = 0;

3629 #endif

3630 mp_size ib;

3631 mp_size ia;

3632 mp_size lim;

3633 mp_err res;

3634

3635 /* Make sure a has enough precision for the output value */

3636 lim = MP_USED(b) + offset;

3637 if((lim > USED(a)) && (res = s_mp_pad(a, lim)) != MP_OKAY)

3638 return res;

3639

3640 /*

3641 Add up all digits up to the precision of b. If b had initially

3642 the same precision as a, or greater, we took care of it by the

3643 padding step above, so there is no problem. If b had initially

3644 less precision, we'll have to make sure the carry out is duly

3645 propagated upward among the higher-order digits of the sum.

3646 */

3647 lim = USED(b);

3648 for(ib = 0, ia = offset; ib < lim; ib++, ia++) {

3649 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3650 w = (mp_word)DIGIT(a, ia) + DIGIT(b, ib) + k;

3651 DIGIT(a, ia) = ACCUM(w);

3652 k = CARRYOUT(w);

3653 #else

3654 d = MP_DIGIT(a, ia);

3655 sum = d + MP_DIGIT(b, ib);

3656 d = (sum < d);

3657 MP_DIGIT(a,ia) = sum += carry;

3658 carry = d + (sum < carry);

3659 #endif

3660 }

3661

3662 /* If we run out of 'b' digits before we're actually done, make

3663 sure the carries get propagated upward...

3664 */

3665 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3666 for (lim = MP_USED(a); k && (ia < lim); ++ia) {

3667 w = (mp_word)DIGIT(a, ia) + k;

3668 DIGIT(a, ia) = ACCUM(w);

3669 k = CARRYOUT(w);

3670 }

3671 #else

3672 for (lim = MP_USED(a); carry && (ia < lim); ++ia) {

3673 d = MP_DIGIT(a, ia);

3674 MP_DIGIT(a,ia) = sum = d + carry;

3675 carry = (sum < d);

3676 }

3677 #endif

3678

3679 /* If there's an overall carry out, increase precision and include

3680 it. We could have done this initially, but why touch the memory

3681 allocator unless we're sure we have to?

3682 */

3683 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)

3684 if(k) {

3685 if((res = s_mp_pad(a, USED(a) + 1)) != MP_OKAY)

3686 return res;

3687

3688 DIGIT(a, ia) = (mp_digit)k;

3689 }

3690 #else

3691 if (carry) {

3692 if((res = s_mp_pad(a, lim + 1)) != MP_OKAY)

3693 return res;

3694

3695 DIGIT(a, lim) = carry;

3696 }

3697 #endif

3698 s_mp_clamp(a);

3699

3700 return MP_OKAY;

3701

3702 } /* end s_mp_add_offset() */

3703

3704 /* }}} */

3705

3706 /* {{{ s_mp_sub(a, b) */

3707

3708 /* Compute a = \|a\| - \|b\|, assumes \|a\| >= \|b\| */

3709 mp_err s_mp_sub(mp_int a, const mp_int b) /* magnitude subtract */

3710 {

3711 mp_digit pa, pb, *limit;

3712 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)

3713 mp_sword w = 0;

3714 #else

3715 mp_digit d, diff, borrow = 0;

3716 #endif

3717

3718 /*

3719 Subtract and propagate borrow. Up to the precision of b, this

3720 accounts for the digits of b; after that, we just make sure the

3721 carries get to the right place. This saves having to pad b out to

3722 the precision of a just to make the loops work right...

3723 */

3724 pa = MP_DIGITS(a);

3725 pb = MP_DIGITS(b);

3726 limit = pb + MP_USED(b);

3727 while (pb < limit) {

3728 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)

3729 w = w + pa - pb++;

3730 *pa++ = ACCUM(w);

3731 w >>= MP_DIGIT_BIT;

3732 #else

3733 d = *pa;

3734 diff = d - *pb++;

3735 d = (diff > d); /* detect borrow */

3736 if (borrow && --diff == MP_DIGIT_MAX)

3737 ++d;

3738 *pa++ = diff;

3739 borrow = d;

3740 #endif

3741 }

3742 limit = MP_DIGITS(a) + MP_USED(a);

3743 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)

3744 while (w && pa < limit) {

3745 w = w + *pa;

3746 *pa++ = ACCUM(w);

3747 w >>= MP_DIGIT_BIT;

3748 }

3749 #else

3750 while (borrow && pa < limit) {

3751 d = *pa;

3752 *pa++ = diff = d - borrow;

3753 borrow = (diff > d);

3754 }

3755 #endif

3756

3757 /* Clobber any leading zeroes we created */

3758 s_mp_clamp(a);

3759

3760 /*

3761 If there was a borrow out, then \|b\| > \|a\| in violation

3762 of our input invariant. We've already done the work,

3763 but we'll at least complain about it...

3764 */

3765 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)

3766 return w ? MP_RANGE : MP_OKAY;

3767 #else

3768 return borrow ? MP_RANGE : MP_OKAY;

3769 #endif

3770 } /* end s_mp_sub() */

3771

3772 /* }}} */

3773

3774 /* Compute c = \|a\| - \|b\|, assumes \|a\| >= \|b\| / / magnitude subtract */

3775 mp_err s_mp_sub_3arg(const mp_int a, const mp_int b, mp_int *c)

3776 {

3777 mp_digit pa, pb, *pc;

3778 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)

3779 mp_sword w = 0;

3780 #else

3781 mp_digit d, diff, borrow = 0;

3782 #endif

3783 int ix, limit;

3784 mp_err res;

3785

3786 MP_SIGN(c) = MP_SIGN(a);

3787

3788 /* Make sure a has enough precision for the output value */

3789 if (MP_OKAY != (res = s_mp_pad(c, MP_USED(a))))

3790 return res;

3791

3792 /*

3793 Subtract and propagate borrow. Up to the precision of b, this

3794 accounts for the digits of b; after that, we just make sure the

3795 carries get to the right place. This saves having to pad b out to

3796 the precision of a just to make the loops work right...

3797 */

3798 pa = MP_DIGITS(a);

3799 pb = MP_DIGITS(b);

3800 pc = MP_DIGITS(c);

3801 limit = MP_USED(b);

3802 for (ix = 0; ix < limit; ++ix) {

3803 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)

3804 w = w + pa++ - pb++;

3805 *pc++ = ACCUM(w);

3806 w >>= MP_DIGIT_BIT;

3807 #else

3808 d = *pa++;

3809 diff = d - *pb++;

3810 d = (diff > d);

3811 if (borrow && --diff == MP_DIGIT_MAX)

3812 ++d;

3813 *pc++ = diff;

3814 borrow = d;

3815 #endif

3816 }

3817 for (limit = MP_USED(a); ix < limit; ++ix) {

3818 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)

3819 w = w + *pa++;

3820 *pc++ = ACCUM(w);

3821 w >>= MP_DIGIT_BIT;

3822 #else

3823 d = *pa++;

3824 *pc++ = diff = d - borrow;

3825 borrow = (diff > d);

3826 #endif

3827 }

3828

3829 /* Clobber any leading zeroes we created */

3830 MP_USED(c) = ix;

3831 s_mp_clamp(c);

3832

3833 /*

3834 If there was a borrow out, then \|b\| > \|a\| in violation

3835 of our input invariant. We've already done the work,

3836 but we'll at least complain about it...

3837 */

3838 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)

3839 return w ? MP_RANGE : MP_OKAY;

3840 #else

3841 return borrow ? MP_RANGE : MP_OKAY;

3842 #endif

3843 }

3844 /* {{{ s_mp_mul(a, b) */

3845

3846 /* Compute a = \|a\| * \|b\| */

3847 mp_err s_mp_mul(mp_int a, const mp_int b)

3848 {

3849 return mp_mul(a, b, a);

3850 } /* end s_mp_mul() */

3851

3852 /* }}} */

3853

3854 #if defined(MP_USE_UINT_DIGIT) && defined(MP_USE_LONG_LONG_MULTIPLY)

3855 /* This trick works on Sparc V8 CPUs with the Workshop compilers. */

3856 #define MP_MUL_DxD(a, b, Phi, Plo) \

3857 { unsigned long long product = (unsigned long long)a * b; \

3858 Plo = (mp_digit)product; \

3859 Phi = (mp_digit)(product >> MP_DIGIT_BIT); }

3860 #elif defined(OSF1)

3861 #define MP_MUL_DxD(a, b, Phi, Plo) \

3862 { Plo = asm ("mulq %a0, %a1, %v0", a, b);\

3863 Phi = asm ("umulh %a0, %a1, %v0", a, b); }

3864 #else

3865 #define MP_MUL_DxD(a, b, Phi, Plo) \

3866 { mp_digit a0b1, a1b0; \

3867 Plo = (a & MP_HALF_DIGIT_MAX) * (b & MP_HALF_DIGIT_MAX); \

3868 Phi = (a >> MP_HALF_DIGIT_BIT) * (b >> MP_HALF_DIGIT_BIT); \

3869 a0b1 = (a & MP_HALF_DIGIT_MAX) * (b >> MP_HALF_DIGIT_BIT); \

3870 a1b0 = (a >> MP_HALF_DIGIT_BIT) * (b & MP_HALF_DIGIT_MAX); \

3871 a1b0 += a0b1; \

3872 Phi += a1b0 >> MP_HALF_DIGIT_BIT; \

3873 if (a1b0 < a0b1) \

3874 Phi += MP_HALF_RADIX; \

3875 a1b0 <<= MP_HALF_DIGIT_BIT; \

3876 Plo += a1b0; \

3877 if (Plo < a1b0) \

3878 ++Phi; \

3879 }

3880 #endif

3881

3882 #if !defined(MP_ASSEMBLY_MULTIPLY)

3883 /* c = a * b */

3884 void s_mpv_mul_d(const mp_digit a, mp_size a_len, mp_digit b, mp_digit c)

3885 {

3886 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)

3887 mp_digit d = 0;

3888

3889 /* Inner product: Digits of a */

3890 while (a_len--) {

3891 mp_word w = ((mp_word)b * *a++) + d;

3892 *c++ = ACCUM(w);

3893 d = CARRYOUT(w);

3894 }

3895 *c = d;

3896 #else

3897 mp_digit carry = 0;

3898 while (a_len--) {

3899 mp_digit a_i = *a++;

3900 mp_digit a0b0, a1b1;

3901

3902 MP_MUL_DxD(a_i, b, a1b1, a0b0);

3903

3904 a0b0 += carry;

3905 if (a0b0 < carry)

3906 ++a1b1;

3907 *c++ = a0b0;

3908 carry = a1b1;

3909 }

3910 *c = carry;

3911 #endif

3912 }

3913

3914 /* c += a * b */

3915 void s_mpv_mul_d_add(const mp_digit *a, mp_size a_len, mp_digit b,

3916 mp_digit *c)

3917 {

3918 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)

3919 mp_digit d = 0;

3920

3921 /* Inner product: Digits of a */

3922 while (a_len--) {

3923 mp_word w = ((mp_word)b * a++) + c + d;

3924 *c++ = ACCUM(w);

3925 d = CARRYOUT(w);

3926 }

3927 *c = d;

3928 #else

3929 mp_digit carry = 0;

3930 while (a_len--) {

3931 mp_digit a_i = *a++;

3932 mp_digit a0b0, a1b1;

3933

3934 MP_MUL_DxD(a_i, b, a1b1, a0b0);

3935

3936 a0b0 += carry;

3937 if (a0b0 < carry)

3938 ++a1b1;

3939 a0b0 += a_i = *c;

3940 if (a0b0 < a_i)

3941 ++a1b1;

3942 *c++ = a0b0;

3943 carry = a1b1;

3944 }

3945 *c = carry;

3946 #endif

3947 }

3948

3949 /* Presently, this is only used by the Montgomery arithmetic code. */

3950 /* c += a * b */

3951 void s_mpv_mul_d_add_prop(const mp_digit a, mp_size a_len, mp_digit b, mp_digit c)

3952 {

3953 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)

3954 mp_digit d = 0;

3955

3956 /* Inner product: Digits of a */

3957 while (a_len--) {

3958 mp_word w = ((mp_word)b * a++) + c + d;

3959 *c++ = ACCUM(w);

3960 d = CARRYOUT(w);

3961 }

3962

3963 while (d) {

3964 mp_word w = (mp_word)*c + d;

3965 *c++ = ACCUM(w);

3966 d = CARRYOUT(w);

3967 }

3968 #else

3969 mp_digit carry = 0;

3970 while (a_len--) {

3971 mp_digit a_i = *a++;

3972 mp_digit a0b0, a1b1;

3973

3974 MP_MUL_DxD(a_i, b, a1b1, a0b0);

3975

3976 a0b0 += carry;

3977 if (a0b0 < carry)

3978 ++a1b1;

3979

3980 a0b0 += a_i = *c;

3981 if (a0b0 < a_i)

3982 ++a1b1;

3983

3984 *c++ = a0b0;

3985 carry = a1b1;

3986 }

3987 while (carry) {

3988 mp_digit c_i = *c;

3989 carry += c_i;

3990 *c++ = carry;

3991 carry = carry < c_i;

3992 }

3993 #endif

3994 }

3995 #endif

3996

3997 #if defined(MP_USE_UINT_DIGIT) && defined(MP_USE_LONG_LONG_MULTIPLY)

3998 /* This trick works on Sparc V8 CPUs with the Workshop compilers. */

3999 #define MP_SQR_D(a, Phi, Plo) \

4000 { unsigned long long square = (unsigned long long)a * a; \

4001 Plo = (mp_digit)square; \

4002 Phi = (mp_digit)(square >> MP_DIGIT_BIT); }

4003 #elif defined(OSF1)

4004 #define MP_SQR_D(a, Phi, Plo) \

4005 { Plo = asm ("mulq %a0, %a0, %v0", a);\

4006 Phi = asm ("umulh %a0, %a0, %v0", a); }

4007 #else

4008 #define MP_SQR_D(a, Phi, Plo) \

4009 { mp_digit Pmid; \

4010 Plo = (a & MP_HALF_DIGIT_MAX) * (a & MP_HALF_DIGIT_MAX); \

4011 Phi = (a >> MP_HALF_DIGIT_BIT) * (a >> MP_HALF_DIGIT_BIT); \

4012 Pmid = (a & MP_HALF_DIGIT_MAX) * (a >> MP_HALF_DIGIT_BIT); \

4013 Phi += Pmid >> (MP_HALF_DIGIT_BIT - 1); \

4014 Pmid <<= (MP_HALF_DIGIT_BIT + 1); \

4015 Plo += Pmid; \

4016 if (Plo < Pmid) \

4017 ++Phi; \

4018 }

4019 #endif

4020

4021 #if !defined(MP_ASSEMBLY_SQUARE)

4022 /* Add the squares of the digits of a to the digits of b. */

4023 void s_mpv_sqr_add_prop(const mp_digit pa, mp_size a_len, mp_digit ps)

4024 {

4025 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)

4026 mp_word w;

4027 mp_digit d;

4028 mp_size ix;

4029

4030 w = 0;

4031 #define ADD_SQUARE(n) \

4032 d = pa[n]; \

4033 w += (d * (mp_word)d) + ps[2*n]; \

4034 ps[2*n] = ACCUM(w); \

4035 w = (w >> DIGIT_BIT) + ps[2*n+1]; \

4036 ps[2*n+1] = ACCUM(w); \

4037 w = (w >> DIGIT_BIT)

4038

4039 for (ix = a_len; ix >= 4; ix -= 4) {

4040 ADD_SQUARE(0);

4041 ADD_SQUARE(1);

4042 ADD_SQUARE(2);

4043 ADD_SQUARE(3);

4044 pa += 4;

4045 ps += 8;

4046 }

4047 if (ix) {

4048 ps += 2*ix;

4049 pa += ix;

4050 switch (ix) {

4051 case 3: ADD_SQUARE(-3); /* FALLTHRU */

4052 case 2: ADD_SQUARE(-2); /* FALLTHRU */

4053 case 1: ADD_SQUARE(-1); /* FALLTHRU */

4054 case 0: break;

4055 }

4056 }

4057 while (w) {

4058 w += *ps;

4059 *ps++ = ACCUM(w);

4060 w = (w >> DIGIT_BIT);

4061 }

4062 #else

4063 mp_digit carry = 0;

4064 while (a_len--) {

4065 mp_digit a_i = *pa++;

4066 mp_digit a0a0, a1a1;

4067

4068 MP_SQR_D(a_i, a1a1, a0a0);

4069

4070 /* here a1a1 and a0a0 constitute a_i ** 2 */

4071 a0a0 += carry;

4072 if (a0a0 < carry)

4073 ++a1a1;

4074

4075 /* now add to ps */

4076 a0a0 += a_i = *ps;

4077 if (a0a0 < a_i)

4078 ++a1a1;

4079 *ps++ = a0a0;

4080 a1a1 += a_i = *ps;

4081 carry = (a1a1 < a_i);

4082 *ps++ = a1a1;

4083 }

4084 while (carry) {

4085 mp_digit s_i = *ps;

4086 carry += s_i;

4087 *ps++ = carry;

4088 carry = carry < s_i;

4089 }

4090 #endif

4091 }

4092 #endif

4093

4094 #if (defined(MP_NO_MP_WORD) \|\| defined(MP_NO_DIV_WORD)) \

4095 && !defined(MP_ASSEMBLY_DIV_2DX1D)

4096 /*

4097 ** Divide 64-bit (Nhi,Nlo) by 32-bit divisor, which must be normalized

4098 ** so its high bit is 1. This code is from NSPR.

4099 */

4100 mp_err s_mpv_div_2dx1d(mp_digit Nhi, mp_digit Nlo, mp_digit divisor,

4101 mp_digit qp, mp_digit rp)

4102 {

4103 mp_digit d1, d0, q1, q0;

4104 mp_digit r1, r0, m;

4105

4106 d1 = divisor >> MP_HALF_DIGIT_BIT;

4107 d0 = divisor & MP_HALF_DIGIT_MAX;

4108 r1 = Nhi % d1;

4109 q1 = Nhi / d1;

4110 m = q1 * d0;

4111 r1 = (r1 << MP_HALF_DIGIT_BIT) \| (Nlo >> MP_HALF_DIGIT_BIT);

4112 if (r1 < m) {

4113 q1--, r1 += divisor;

4114 if (r1 >= divisor && r1 < m) {

4115 q1--, r1 += divisor;

4116 }

4117 }

4118 r1 -= m;

4119 r0 = r1 % d1;

4120 q0 = r1 / d1;

4121 m = q0 * d0;

4122 r0 = (r0 << MP_HALF_DIGIT_BIT) \| (Nlo & MP_HALF_DIGIT_MAX);

4123 if (r0 < m) {

4124 q0--, r0 += divisor;

4125 if (r0 >= divisor && r0 < m) {

4126 q0--, r0 += divisor;

4127 }

4128 }

4129 if (qp)

4130 *qp = (q1 << MP_HALF_DIGIT_BIT) \| q0;

4131 if (rp)

4132 *rp = r0 - m;

4133 return MP_OKAY;

4134 }

4135 #endif

4136

4137 #if MP_SQUARE

4138 /* {{{ s_mp_sqr(a) */

4139

4140 mp_err s_mp_sqr(mp_int *a)

4141 {

4142 mp_err res;

4143 mp_int tmp;

4144

4145 if((res = mp_init_size(&tmp, 2 * USED(a))) != MP_OKAY)

4146 return res;

4147 res = mp_sqr(a, &tmp);

4148 if (res == MP_OKAY) {

4149 s_mp_exch(&tmp, a);

4150 }

4151 mp_clear(&tmp);

4152 return res;

4153 }

4154

4155 /* }}} */

4156 #endif

4157

4158 /* {{{ s_mp_div(a, b) */

4159

4160 /*

4161 s_mp_div(a, b)

4162

4163 Compute a = a / b and b = a mod b. Assumes b > a.

4164 */

4165

4166 mp_err s_mp_div(mp_int rem, / i: dividend, o: remainder */

4167 mp_int div, / i: divisor */

4168 mp_int quot) / i: 0; o: quotient */

4169 {

4170 mp_int part, t;

4171 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)

4172 mp_word q_msd;

4173 #else

4174 mp_digit q_msd;

4175 #endif

4176 mp_err res;

4177 mp_digit d;

4178 mp_digit div_msd;

4179 int ix;

4180

4181 if(mp_cmp_z(div) == 0)

4182 return MP_RANGE;

4183

4184 DIGITS(&t) = 0;

4185 /* Shortcut if divisor is power of two */

4186 if((ix = s_mp_ispow2(div)) >= 0) {

4187 MP_CHECKOK( mp_copy(rem, quot) );

4188 s_mp_div_2d(quot, (mp_digit)ix);

4189 s_mp_mod_2d(rem, (mp_digit)ix);

4190

4191 return MP_OKAY;

4192 }

4193

4194 MP_SIGN(rem) = ZPOS;

4195 MP_SIGN(div) = ZPOS;

4196

4197 /* A working temporary for division */

4198 MP_CHECKOK( mp_init_size(&t, MP_ALLOC(rem)));

4199

4200 /* Normalize to optimize guessing */

4201 MP_CHECKOK( s_mp_norm(rem, div, &d) );

4202

4203 part = *rem;

4204

4205 /* Perform the division itself...woo! */

4206 MP_USED(quot) = MP_ALLOC(quot);

4207

4208 /* Find a partial substring of rem which is at least div */

4209 /* If we didn't find one, we're finished dividing */

4210 while (MP_USED(rem) > MP_USED(div) \|\| s_mp_cmp(rem, div) >= 0) {

4211 int i;

4212 int unusedRem;

4213

4214 unusedRem = MP_USED(rem) - MP_USED(div);

4215 MP_DIGITS(&part) = MP_DIGITS(rem) + unusedRem;

4216 MP_ALLOC(&part) = MP_ALLOC(rem) - unusedRem;

4217 MP_USED(&part) = MP_USED(div);

4218 if (s_mp_cmp(&part, div) < 0) {

4219 -- unusedRem;

4220 #if MP_ARGCHK == 2

4221 assert(unusedRem >= 0);

4222 #endif

4223 -- MP_DIGITS(&part);

4224 ++ MP_USED(&part);

4225 ++ MP_ALLOC(&part);

4226 }

4227

4228 /* Compute a guess for the next quotient digit */

4229 q_msd = MP_DIGIT(&part, MP_USED(&part) - 1);

4230 div_msd = MP_DIGIT(div, MP_USED(div) - 1);

4231 if (q_msd >= div_msd) {

4232 q_msd = 1;

4233 } else if (MP_USED(&part) > 1) {

4234 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)

4235 q_msd = (q_msd << MP_DIGIT_BIT) \| MP_DIGIT(&part, MP_USED(&part) - 2);

4236 q_msd /= div_msd;

4237 if (q_msd == RADIX)

4238 --q_msd;

4239 #else

4240 mp_digit r;

4241 MP_CHECKOK( s_mpv_div_2dx1d(q_msd, MP_DIGIT(&part, MP_USED(&part) - 2),

4242 div_msd, &q_msd, &r) );

4243 #endif

4244 } else {

4245 q_msd = 0;

4246 }

4247 #if MP_ARGCHK == 2

4248 assert(q_msd > 0); /* This case should never occur any more. */

4249 #endif

4250 if (q_msd <= 0)

4251 break;

4252

4253 /* See what that multiplies out to */

4254 mp_copy(div, &t);

4255 MP_CHECKOK( s_mp_mul_d(&t, (mp_digit)q_msd) );

4256

4257 /*

4258 If it's too big, back it off. We should not have to do this

4259 more than once, or, in rare cases, twice. Knuth describes a

4260 method by which this could be reduced to a maximum of once, but

4261 I didn't implement that here.

4262 * When using s_mpv_div_2dx1d, we may have to do this 3 times.

4263 */

4264 for (i = 4; s_mp_cmp(&t, &part) > 0 && i > 0; --i) {

4265 --q_msd;

4266 s_mp_sub(&t, div); /* t -= div */

4267 }

4268 if (i < 0) {

4269 res = MP_RANGE;

4270 goto CLEANUP;

4271 }

4272

4273 /* At this point, q_msd should be the right next digit */

4274 MP_CHECKOK( s_mp_sub(&part, &t) ); /* part -= t */

4275 s_mp_clamp(rem);

4276

4277 /*

4278 Include the digit in the quotient. We allocated enough memory

4279 for any quotient we could ever possibly get, so we should not

4280 have to check for failures here

4281 */

4282 MP_DIGIT(quot, unusedRem) = (mp_digit)q_msd;

4283 }

4284

4285 /* Denormalize remainder */

4286 if (d) {

4287 s_mp_div_2d(rem, d);

4288 }

4289

4290 s_mp_clamp(quot);

4291

4292 CLEANUP:

4293 mp_clear(&t);

4294

4295 return res;

4296

4297 } /* end s_mp_div() */

4298

4299

4300 /* }}} */

4301

4302 /* {{{ s_mp_2expt(a, k) */

4303

4304 mp_err s_mp_2expt(mp_int *a, mp_digit k)

4305 {

4306 mp_err res;

4307 mp_size dig, bit;

4308

4309 dig = k / DIGIT_BIT;

4310 bit = k % DIGIT_BIT;

4311

4312 mp_zero(a);

4313 if((res = s_mp_pad(a, dig + 1)) != MP_OKAY)

4314 return res;

4315

4316 DIGIT(a, dig) \|= ((mp_digit)1 << bit);

4317

4318 return MP_OKAY;

4319

4320 } /* end s_mp_2expt() */

4321

4322 /* }}} */

4323

4324 /* {{{ s_mp_reduce(x, m, mu) */

4325

4326 /*

4327 Compute Barrett reduction, x (mod m), given a precomputed value for

4328 mu = b^2k / m, where b = RADIX and k = #digits(m). This should be

4329 faster than straight division, when many reductions by the same

4330 value of m are required (such as in modular exponentiation). This

4331 can nearly halve the time required to do modular exponentiation,

4332 as compared to using the full integer divide to reduce.

4333

4334 This algorithm was derived from the _Handbook of Applied

4335 Cryptography_ by Menezes, Oorschot and VanStone, Ch. 14,

4336 pp. 603-604.

4337 */

4338

4339 mp_err s_mp_reduce(mp_int x, const mp_int m, const mp_int *mu)

4340 {

4341 mp_int q;

4342 mp_err res;

4343

4344 if((res = mp_init_copy(&q, x)) != MP_OKAY)

4345 return res;

4346

4347 s_mp_rshd(&q, USED(m) - 1); /* q1 = x / b^(k-1) */

4348 s_mp_mul(&q, mu); /* q2 = q1 * mu */

4349 s_mp_rshd(&q, USED(m) + 1); /* q3 = q2 / b^(k+1) */

4350

4351 /* x = x mod b^(k+1), quick (no division) */

4352 s_mp_mod_2d(x, DIGIT_BIT * (USED(m) + 1));

4353

4354 /* q = q * m mod b^(k+1), quick (no division) */

4355 s_mp_mul(&q, m);

4356 s_mp_mod_2d(&q, DIGIT_BIT * (USED(m) + 1));

4357

4358 /* x = x - q */

4359 if((res = mp_sub(x, &q, x)) != MP_OKAY)

4360 goto CLEANUP;

4361

4362 /* If x < 0, add b^(k+1) to it */

4363 if(mp_cmp_z(x) < 0) {

4364 mp_set(&q, 1);

4365 if((res = s_mp_lshd(&q, USED(m) + 1)) != MP_OKAY)

4366 goto CLEANUP;

4367 if((res = mp_add(x, &q, x)) != MP_OKAY)

4368 goto CLEANUP;

4369 }

4370

4371 /* Back off if it's too big */

4372 while(mp_cmp(x, m) >= 0) {

4373 if((res = s_mp_sub(x, m)) != MP_OKAY)

4374 break;

4375 }

4376

4377 CLEANUP:

4378 mp_clear(&q);

4379

4380 return res;

4381

4382 } /* end s_mp_reduce() */

4383

4384 /* }}} */

4385

4386 /* }}} */

4387

4388 /* {{{ Primitive comparisons */

4389

4390 /* {{{ s_mp_cmp(a, b) */

4391

4392 /* Compare \|a\| <=> \|b\|, return 0 if equal, <0 if a<b, >0 if a>b */

4393 int s_mp_cmp(const mp_int a, const mp_int b)

4394 {

4395 mp_size used_a = MP_USED(a);

4396 {

4397 mp_size used_b = MP_USED(b);

4398

4399 if (used_a > used_b)

4400 goto IS_GT;

4401 if (used_a < used_b)

4402 goto IS_LT;

4403 }

4404 {

4405 mp_digit pa, pb;

4406 mp_digit da = 0, db = 0;

4407

4408 #define CMP_AB(n) if ((da = pa[n]) != (db = pb[n])) goto done

4409

4410 pa = MP_DIGITS(a) + used_a;

4411 pb = MP_DIGITS(b) + used_a;

4412 while (used_a >= 4) {

4413 pa -= 4;

4414 pb -= 4;

4415 used_a -= 4;

4416 CMP_AB(3);

4417 CMP_AB(2);

4418 CMP_AB(1);

4419 CMP_AB(0);

4420 }

4421 while (used_a-- > 0 && ((da = --pa) == (db = --pb)))

4422 /* do nothing */;

4423 done:

4424 if (da > db)

4425 goto IS_GT;

4426 if (da < db)

4427 goto IS_LT;

4428 }

4429 return MP_EQ;

4430 IS_LT:

4431 return MP_LT;

4432 IS_GT:

4433 return MP_GT;

4434 } /* end s_mp_cmp() */

4435

4436 /* }}} */

4437

4438 /* {{{ s_mp_cmp_d(a, d) */

4439

4440 /* Compare \|a\| <=> d, return 0 if equal, <0 if a<d, >0 if a>d */

4441 int s_mp_cmp_d(const mp_int *a, mp_digit d)

4442 {

4443 if(USED(a) > 1)

4444 return MP_GT;

4445

4446 if(DIGIT(a, 0) < d)

4447 return MP_LT;

4448 else if(DIGIT(a, 0) > d)

4449 return MP_GT;

4450 else

4451 return MP_EQ;

4452

4453 } /* end s_mp_cmp_d() */

4454

4455 /* }}} */

4456

4457 /* {{{ s_mp_ispow2(v) */

4458

4459 /*

4460 Returns -1 if the value is not a power of two; otherwise, it returns

4461 k such that v = 2^k, i.e. lg(v).

4462 */

4463 int s_mp_ispow2(const mp_int *v)

4464 {

4465 mp_digit d;

4466 int extra = 0, ix;

4467

4468 ix = MP_USED(v) - 1;

4469 d = MP_DIGIT(v, ix); /* most significant digit of v */

4470

4471 extra = s_mp_ispow2d(d);

4472 if (extra < 0 \|\| ix == 0)

4473 return extra;

4474

4475 while (--ix >= 0) {

4476 if (DIGIT(v, ix) != 0)

4477 return -1; /* not a power of two */

4478 extra += MP_DIGIT_BIT;

4479 }

4480

4481 return extra;

4482

4483 } /* end s_mp_ispow2() */

4484

4485 /* }}} */

4486

4487 /* {{{ s_mp_ispow2d(d) */

4488

4489 int s_mp_ispow2d(mp_digit d)

4490 {

4491 if ((d != 0) && ((d & (d-1)) == 0)) { /* d is a power of 2 */

4492 int pow = 0;

4493 #if defined (MP_USE_UINT_DIGIT)

4494 if (d & 0xffff0000U)

4495 pow += 16;

4496 if (d & 0xff00ff00U)

4497 pow += 8;

4498 if (d & 0xf0f0f0f0U)

4499 pow += 4;

4500 if (d & 0xccccccccU)

4501 pow += 2;

4502 if (d & 0xaaaaaaaaU)

4503 pow += 1;

4504 #elif defined(MP_USE_LONG_LONG_DIGIT)

4505 if (d & 0xffffffff00000000ULL)

4506 pow += 32;

4507 if (d & 0xffff0000ffff0000ULL)

4508 pow += 16;

4509 if (d & 0xff00ff00ff00ff00ULL)

4510 pow += 8;

4511 if (d & 0xf0f0f0f0f0f0f0f0ULL)

4512 pow += 4;

4513 if (d & 0xccccccccccccccccULL)

4514 pow += 2;

4515 if (d & 0xaaaaaaaaaaaaaaaaULL)

4516 pow += 1;

4517 #elif defined(MP_USE_LONG_DIGIT)

4518 if (d & 0xffffffff00000000UL)

4519 pow += 32;

4520 if (d & 0xffff0000ffff0000UL)

4521 pow += 16;

4522 if (d & 0xff00ff00ff00ff00UL)

4523 pow += 8;

4524 if (d & 0xf0f0f0f0f0f0f0f0UL)

4525 pow += 4;

4526 if (d & 0xccccccccccccccccUL)

4527 pow += 2;

4528 if (d & 0xaaaaaaaaaaaaaaaaUL)

4529 pow += 1;

4530 #else

4531 #error "unknown type for mp_digit"

4532 #endif

4533 return pow;

4534 }

4535 return -1;

4536

4537 } /* end s_mp_ispow2d() */

4538

4539 /* }}} */

4540

4541 /* }}} */

4542

4543 /* {{{ Primitive I/O helpers */

4544

4545 /* {{{ s_mp_tovalue(ch, r) */

4546

4547 /*

4548 Convert the given character to its digit value, in the given radix.

4549 If the given character is not understood in the given radix, -1 is

4550 returned. Otherwise the digit's numeric value is returned.

4551

4552 The results will be odd if you use a radix < 2 or > 62, you are

4553 expected to know what you're up to.

4554 */

4555 int s_mp_tovalue(char ch, int r)

4556 {

4557 int val, xch;

4558

4559 if(r > 36)

4560 xch = ch;

4561 else

4562 xch = toupper(ch);

4563

4564 if(isdigit(xch))

4565 val = xch - '0';

4566 else if(isupper(xch))

4567 val = xch - 'A' + 10;

4568 else if(islower(xch))

4569 val = xch - 'a' + 36;

4570 else if(xch == '+')

4571 val = 62;

4572 else if(xch == '/')

4573 val = 63;

4574 else

4575 return -1;

4576

4577 if(val < 0 \|\| val >= r)

4578 return -1;

4579

4580 return val;

4581

4582 } /* end s_mp_tovalue() */

4583

4584 /* }}} */

4585

4586 /* {{{ s_mp_todigit(val, r, low) */

4587

4588 /*

4589 Convert val to a radix-r digit, if possible. If val is out of range

4590 for r, returns zero. Otherwise, returns an ASCII character denoting

4591 the value in the given radix.

4592

4593 The results may be odd if you use a radix < 2 or > 64, you are

4594 expected to know what you're doing.

4595 */

4596

4597 char s_mp_todigit(mp_digit val, int r, int low)

4598 {

4599 char ch;

4600

4601 if(val >= r)

4602 return 0;

4603

4604 ch = s_dmap_1[val];

4605

4606 if(r <= 36 && low)

4607 ch = tolower(ch);

4608

4609 return ch;

4610

4611 } /* end s_mp_todigit() */

4612

4613 /* }}} */

4614

4615 /* {{{ s_mp_outlen(bits, radix) */

4616

4617 /*

4618 Return an estimate for how long a string is needed to hold a radix

4619 r representation of a number with 'bits' significant bits, plus an

4620 extra for a zero terminator (assuming C style strings here)

4621 */

4622 int s_mp_outlen(int bits, int r)

4623 {

4624 return (int)((double)bits * LOG_V_2(r) + 1.5) + 1;

4625

4626 } /* end s_mp_outlen() */

4627

4628 /* }}} */

4629

4630 /* }}} */

4631

4632 /* {{{ mp_read_unsigned_octets(mp, str, len) */

4633 /* mp_read_unsigned_octets(mp, str, len)

4634 Read in a raw value (base 256) into the given mp_int

4635 No sign bit, number is positive. Leading zeros ignored.

4636 */

4637

4638 mp_err

4639 mp_read_unsigned_octets(mp_int mp, const unsigned char str, mp_size len)

4640 {

4641 int count;

4642 mp_err res;

4643 mp_digit d;

4644

4645 ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);

4646

4647 mp_zero(mp);

4648

4649 count = len % sizeof(mp_digit);

4650 if (count) {

4651 for (d = 0; count-- > 0; --len) {

4652 d = (d << 8) \| *str++;

4653 }

4654 MP_DIGIT(mp, 0) = d;

4655 }

4656

4657 /* Read the rest of the digits */

4658 for(; len > 0; len -= sizeof(mp_digit)) {

4659 for (d = 0, count = sizeof(mp_digit); count > 0; --count) {

4660 d = (d << 8) \| *str++;

4661 }

4662 if (MP_EQ == mp_cmp_z(mp)) {

4663 if (!d)

4664 continue;

4665 } else {

4666 if((res = s_mp_lshd(mp, 1)) != MP_OKAY)

4667 return res;

4668 }

4669 MP_DIGIT(mp, 0) = d;

4670 }

4671 return MP_OKAY;

4672 } /* end mp_read_unsigned_octets() */

4673 /* }}} */

4674

4675 /* {{{ mp_unsigned_octet_size(mp) */

4676 int

4677 mp_unsigned_octet_size(const mp_int *mp)

4678 {

4679 int bytes;

4680 int ix;

4681 mp_digit d = 0;

4682

4683 ARGCHK(mp != NULL, MP_BADARG);

4684 ARGCHK(MP_ZPOS == SIGN(mp), MP_BADARG);

4685

4686 bytes = (USED(mp) * sizeof(mp_digit));

4687

4688 /* subtract leading zeros. */

4689 /* Iterate over each digit... */

4690 for(ix = USED(mp) - 1; ix >= 0; ix--) {

4691 d = DIGIT(mp, ix);

4692 if (d)

4693 break;

4694 bytes -= sizeof(d);

4695 }

4696 if (!bytes)

4697 return 1;

4698

4699 /* Have MSD, check digit bytes, high order first */

4700 for(ix = sizeof(mp_digit) - 1; ix >= 0; ix--) {

4701 unsigned char x = (unsigned char)(d >> (ix * CHAR_BIT));

4702 if (x)

4703 break;

4704 --bytes;

4705 }

4706 return bytes;

4707 } /* end mp_unsigned_octet_size() */

4708 /* }}} */

4709

4710 /* {{{ mp_to_unsigned_octets(mp, str) */

4711 /* output a buffer of big endian octets no longer than specified. */

4712 mp_err

4713 mp_to_unsigned_octets(const mp_int mp, unsigned char str, mp_size maxlen)

4714 {

4715 int ix, pos = 0;

4716 int bytes;

4717

4718 ARGCHK(mp != NULL && str != NULL && !SIGN(mp), MP_BADARG);

4719

4720 bytes = mp_unsigned_octet_size(mp);

4721 ARGCHK(bytes >= 0 && bytes <= maxlen, MP_BADARG);

4722

4723 /* Iterate over each digit... */

4724 for(ix = USED(mp) - 1; ix >= 0; ix--) {

4725 mp_digit d = DIGIT(mp, ix);

4726 int jx;

4727

4728 /* Unpack digit bytes, high order first */

4729 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {

4730 unsigned char x = (unsigned char)(d >> (jx * CHAR_BIT));

4731 if (!pos && !x) /* suppress leading zeros */

4732 continue;

4733 str[pos++] = x;

4734 }

4735 }

4736 if (!pos)

4737 str[pos++] = 0;

4738 return pos;

4739 } /* end mp_to_unsigned_octets() */

4740 /* }}} */

4741

4742 /* {{{ mp_to_signed_octets(mp, str) */

4743 /* output a buffer of big endian octets no longer than specified. */

4744 mp_err

4745 mp_to_signed_octets(const mp_int mp, unsigned char str, mp_size maxlen)

4746 {

4747 int ix, pos = 0;

4748 int bytes;

4749

4750 ARGCHK(mp != NULL && str != NULL && !SIGN(mp), MP_BADARG);

4751

4752 bytes = mp_unsigned_octet_size(mp);

4753 ARGCHK(bytes >= 0 && bytes <= maxlen, MP_BADARG);

4754

4755 /* Iterate over each digit... */

4756 for(ix = USED(mp) - 1; ix >= 0; ix--) {

4757 mp_digit d = DIGIT(mp, ix);

4758 int jx;

4759

4760 /* Unpack digit bytes, high order first */

4761 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {

4762 unsigned char x = (unsigned char)(d >> (jx * CHAR_BIT));

4763 if (!pos) {

4764 if (!x) /* suppress leading zeros */

4765 continue;

4766 if (x & 0x80) { /* add one leading zero to make output positive. */

4767 ARGCHK(bytes + 1 <= maxlen, MP_BADARG);

4768 if (bytes + 1 > maxlen)

4769 return MP_BADARG;

4770 str[pos++] = 0;

4771 }

4772 }

4773 str[pos++] = x;

4774 }

4775 }

4776 if (!pos)

4777 str[pos++] = 0;

4778 return pos;

4779 } /* end mp_to_signed_octets() */

4780 /* }}} */

4781

4782 /* {{{ mp_to_fixlen_octets(mp, str) */

4783 /* output a buffer of big endian octets exactly as long as requested. */

4784 mp_err

4785 mp_to_fixlen_octets(const mp_int mp, unsigned char str, mp_size length)

4786 {

4787 int ix, pos = 0;

4788 int bytes;

4789

4790 ARGCHK(mp != NULL && str != NULL && !SIGN(mp), MP_BADARG);

4791

4792 bytes = mp_unsigned_octet_size(mp);

4793 ARGCHK(bytes >= 0 && bytes <= length, MP_BADARG);

4794

4795 /* place any needed leading zeros */

4796 for (;length > bytes; --length) {

4797 *str++ = 0;

4798 }

4799

4800 /* Iterate over each digit... */

4801 for(ix = USED(mp) - 1; ix >= 0; ix--) {

4802 mp_digit d = DIGIT(mp, ix);

4803 int jx;

4804

4805 /* Unpack digit bytes, high order first */

4806 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {

4807 unsigned char x = (unsigned char)(d >> (jx * CHAR_BIT));

4808 if (!pos && !x) /* suppress leading zeros */

4809 continue;

4810 str[pos++] = x;

4811 }

4812 }

4813 if (!pos)

4814 str[pos++] = 0;

4815 return MP_OKAY;

4816 } /* end mp_to_fixlen_octets() */

4817 /* }}} */

4818

4819

4820 /------------------------------------------------------------------------/

4821 /* HERE THERE BE DRAGONS */

OLD	NEW

« no previous file with comments | « mozilla/security/nss/lib/freebl/mpi/mpi.h ('k') | mozilla/security/nss/lib/freebl/mpi/mpi-config.h » ('j') | no next file with comments »