nss/lib/freebl/mpi/mpi.c - Issue 2078763002: Delete bundled copy of NSS and replace with README.

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

Issue 2078763002: Delete bundled copy of NSS and replace with README. (Closed) Base URL: https://chromium.googlesource.com/chromium/deps/nss@master

Patch Set: Delete bundled copy of NSS and replace with README. Created 4 years, 6 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

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

10 #include "mpi-priv.h"

11 #if defined(OSF1)

12 #include <c_asm.h>

13 #endif

14

15 #if defined(__arm__) && \

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

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

18 #undef MP_ASSEMBLY_MULTIPLY

19 #undef MP_ASSEMBLY_SQUARE

20 #endif

21

22 #if MP_LOGTAB

23 /*

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

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

26

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

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

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

30 greater than log_r(n), where:

31

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

33

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

35 which are the output bases supported.

36 */

37 #include "logtab.h"

38 #endif

39

40 /* {{{ Constant strings */

41

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

43 static const char *mp_err_string[] = {

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

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

46 "boolean false", /* MP_NO */

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

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

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

50 "result is undefined" /* MP_UNDEF */

51 };

52

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

54

55 /* s_dmap_1 - standard digits and letters */

56 static const char *s_dmap_1 =

57 "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";

58

59 /* }}} */

60

61 unsigned long mp_allocs;

62 unsigned long mp_frees;

63 unsigned long mp_copies;

64

65 /* {{{ Default precision manipulation */

66

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

68 static mp_size s_mp_defprec = MP_DEFPREC;

69

70 mp_size mp_get_prec(void)

71 {

72 return s_mp_defprec;

73

74 } /* end mp_get_prec() */

75

76 void mp_set_prec(mp_size prec)

77 {

78 if(prec == 0)

79 s_mp_defprec = MP_DEFPREC;

80 else

81 s_mp_defprec = prec;

82

83 } /* end mp_set_prec() */

84

85 /* }}} */

86

87 /------------------------------------------------------------------------/

88 /* {{{ mp_init(mp) */

89

90 /*

91 mp_init(mp)

92

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

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

95 */

96

97 mp_err mp_init(mp_int *mp)

98 {

99 return mp_init_size(mp, s_mp_defprec);

100

101 } /* end mp_init() */

102

103 /* }}} */

104

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

106

107 /*

108 mp_init_size(mp, prec)

109

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

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

112 not be allocated for the structure.

113 */

114

115 mp_err mp_init_size(mp_int *mp, mp_size prec)

116 {

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

118

119 prec = MP_ROUNDUP(prec, s_mp_defprec);

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

121 return MP_MEM;

122

123 SIGN(mp) = ZPOS;

124 USED(mp) = 1;

125 ALLOC(mp) = prec;

126

127 return MP_OKAY;

128

129 } /* end mp_init_size() */

130

131 /* }}} */

132

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

134

135 /*

136 mp_init_copy(mp, from)

137

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

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

140 structure.

141 */

142

143 mp_err mp_init_copy(mp_int mp, const mp_int from)

144 {

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

146

147 if(mp == from)

148 return MP_OKAY;

149

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

151 return MP_MEM;

152

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

154 USED(mp) = USED(from);

155 ALLOC(mp) = ALLOC(from);

156 SIGN(mp) = SIGN(from);

157

158 return MP_OKAY;

159

160 } /* end mp_init_copy() */

161

162 /* }}} */

163

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

165

166 /*

167 mp_copy(from, to)

168

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

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

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

172 */

173

174 mp_err mp_copy(const mp_int from, mp_int to)

175 {

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

177

178 if(from == to)

179 return MP_OKAY;

180

181 { /* copy */

182 mp_digit *tmp;

183

184 /*

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

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

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

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

189 usual

190 */

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

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

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

194

195 } else {

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

197 return MP_MEM;

198

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

200

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

202 #if MP_CRYPTO

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

204 #endif

205 s_mp_free(DIGITS(to));

206 }

207

208 DIGITS(to) = tmp;

209 ALLOC(to) = ALLOC(from);

210 }

211

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

213 USED(to) = USED(from);

214 SIGN(to) = SIGN(from);

215 } /* end copy */

216

217 return MP_OKAY;

218

219 } /* end mp_copy() */

220

221 /* }}} */

222

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

224

225 /*

226 mp_exch(mp1, mp2)

227

228 Exchange mp1 and mp2 without allocating any intermediate memory

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

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

231 */

232

233 void mp_exch(mp_int mp1, mp_int mp2)

234 {

235 #if MP_ARGCHK == 2

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

237 #else

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

239 return;

240 #endif

241

242 s_mp_exch(mp1, mp2);

243

244 } /* end mp_exch() */

245

246 /* }}} */

247

248 /* {{{ mp_clear(mp) */

249

250 /*

251 mp_clear(mp)

252

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

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

255 get tollchocked.

256 */

257

258 void mp_clear(mp_int *mp)

259 {

260 if(mp == NULL)

261 return;

262

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

264 #if MP_CRYPTO

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

266 #endif

267 s_mp_free(DIGITS(mp));

268 DIGITS(mp) = NULL;

269 }

270

271 USED(mp) = 0;

272 ALLOC(mp) = 0;

273

274 } /* end mp_clear() */

275

276 /* }}} */

277

278 /* {{{ mp_zero(mp) */

279

280 /*

281 mp_zero(mp)

282

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

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

285 */

286 void mp_zero(mp_int *mp)

287 {

288 if(mp == NULL)

289 return;

290

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

292 USED(mp) = 1;

293 SIGN(mp) = ZPOS;

294

295 } /* end mp_zero() */

296

297 /* }}} */

298

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

300

301 void mp_set(mp_int *mp, mp_digit d)

302 {

303 if(mp == NULL)

304 return;

305

306 mp_zero(mp);

307 DIGIT(mp, 0) = d;

308

309 } /* end mp_set() */

310

311 /* }}} */

312

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

314

315 mp_err mp_set_int(mp_int *mp, long z)

316 {

317 int ix;

318 unsigned long v = labs(z);

319 mp_err res;

320

321 ARGCHK(mp != NULL, MP_BADARG);

322

323 mp_zero(mp);

324 if(z == 0)

325 return MP_OKAY; /* shortcut for zero */

326

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

328 DIGIT(mp,0) = v;

329 } else {

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

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

332 return res;

333

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

335 if (res != MP_OKAY)

336 return res;

337 }

338 }

339 if(z < 0)

340 SIGN(mp) = NEG;

341

342 return MP_OKAY;

343

344 } /* end mp_set_int() */

345

346 /* }}} */

347

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

349

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

351 {

352 int ix;

353 mp_err res;

354

355 ARGCHK(mp != NULL, MP_BADARG);

356

357 mp_zero(mp);

358 if(z == 0)

359 return MP_OKAY; /* shortcut for zero */

360

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

362 DIGIT(mp,0) = z;

363 } else {

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

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

366 return res;

367

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

369 if (res != MP_OKAY)

370 return res;

371 }

372 }

373 return MP_OKAY;

374 } /* end mp_set_ulong() */

375

376 /* }}} */

377

378 /------------------------------------------------------------------------/

379 /* {{{ Digit arithmetic */

380

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

382

383 /*

384 mp_add_d(a, d, b)

385

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

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

388 */

389

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

391 {

392 mp_int tmp;

393 mp_err res;

394

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

396

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

398 return res;

399

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

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

402 goto CLEANUP;

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

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

405 goto CLEANUP;

406 } else {

407 mp_neg(&tmp, &tmp);

408

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

410 }

411

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

413 SIGN(&tmp) = ZPOS;

414

415 s_mp_exch(&tmp, b);

416

417 CLEANUP:

418 mp_clear(&tmp);

419 return res;

420

421 } /* end mp_add_d() */

422

423 /* }}} */

424

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

426

427 /*

428 mp_sub_d(a, d, b)

429

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

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

432 */

433

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

435 {

436 mp_int tmp;

437 mp_err res;

438

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

440

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

442 return res;

443

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

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

446 goto CLEANUP;

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

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

449 goto CLEANUP;

450 } else {

451 mp_neg(&tmp, &tmp);

452

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

454 SIGN(&tmp) = NEG;

455 }

456

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

458 SIGN(&tmp) = ZPOS;

459

460 s_mp_exch(&tmp, b);

461

462 CLEANUP:

463 mp_clear(&tmp);

464 return res;

465

466 } /* end mp_sub_d() */

467

468 /* }}} */

469

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

471

472 /*

473 mp_mul_d(a, d, b)

474

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

476 of its multiplicand (single digits are unsigned anyway)

477 */

478

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

480 {

481 mp_err res;

482

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

484

485 if(d == 0) {

486 mp_zero(b);

487 return MP_OKAY;

488 }

489

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

491 return res;

492

493 res = s_mp_mul_d(b, d);

494

495 return res;

496

497 } /* end mp_mul_d() */

498

499 /* }}} */

500

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

502

503 mp_err mp_mul_2(const mp_int a, mp_int c)

504 {

505 mp_err res;

506

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

508

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

510 return res;

511

512 return s_mp_mul_2(c);

513

514 } /* end mp_mul_2() */

515

516 /* }}} */

517

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

519

520 /*

521 mp_div_d(a, d, q, r)

522

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

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

525 unsigned anyway).

526 */

527

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

529 {

530 mp_err res;

531 mp_int qp;

532 mp_digit rem;

533 int pow;

534

535 ARGCHK(a != NULL, MP_BADARG);

536

537 if(d == 0)

538 return MP_RANGE;

539

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

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

542 mp_digit mask;

543

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

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

546

547 if(q) {

548 mp_copy(a, q);

549 s_mp_div_2d(q, pow);

550 }

551

552 if(r)

553 *r = rem;

554

555 return MP_OKAY;

556 }

557

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

559 return res;

560

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

562

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

564 SIGN(q) = ZPOS;

565

566 if(r)

567 *r = rem;

568

569 if(q)

570 s_mp_exch(&qp, q);

571

572 mp_clear(&qp);

573 return res;

574

575 } /* end mp_div_d() */

576

577 /* }}} */

578

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

580

581 /*

582 mp_div_2(a, c)

583

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

585 */

586

587 mp_err mp_div_2(const mp_int a, mp_int c)

588 {

589 mp_err res;

590

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

592

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

594 return res;

595

596 s_mp_div_2(c);

597

598 return MP_OKAY;

599

600 } /* end mp_div_2() */

601

602 /* }}} */

603

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

605

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

607 {

608 mp_int s, x;

609 mp_err res;

610

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

612

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

614 return res;

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

616 goto X;

617

618 DIGIT(&s, 0) = 1;

619

620 while(d != 0) {

621 if(d & 1) {

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

623 goto CLEANUP;

624 }

625

626 d /= 2;

627

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

629 goto CLEANUP;

630 }

631

632 s_mp_exch(&s, c);

633

634 CLEANUP:

635 mp_clear(&x);

636 X:

637 mp_clear(&s);

638

639 return res;

640

641 } /* end mp_expt_d() */

642

643 /* }}} */

644

645 /* }}} */

646

647 /------------------------------------------------------------------------/

648 /* {{{ Full arithmetic */

649

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

651

652 /*

653 mp_abs(a, b)

654

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

656 */

657

658 mp_err mp_abs(const mp_int a, mp_int b)

659 {

660 mp_err res;

661

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

663

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

665 return res;

666

667 SIGN(b) = ZPOS;

668

669 return MP_OKAY;

670

671 } /* end mp_abs() */

672

673 /* }}} */

674

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

676

677 /*

678 mp_neg(a, b)

679

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

681 */

682

683 mp_err mp_neg(const mp_int a, mp_int b)

684 {

685 mp_err res;

686

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

688

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

690 return res;

691

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

693 SIGN(b) = ZPOS;

694 else

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

696

697 return MP_OKAY;

698

699 } /* end mp_neg() */

700

701 /* }}} */

702

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

704

705 /*

706 mp_add(a, b, c)

707

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

709 */

710

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

712 {

713 mp_err res;

714

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

716

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

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

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

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

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

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

723 }

724

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

726 SIGN(c) = ZPOS;

727

728 CLEANUP:

729 return res;

730

731 } /* end mp_add() */

732

733 /* }}} */

734

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

736

737 /*

738 mp_sub(a, b, c)

739

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

741 */

742

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

744 {

745 mp_err res;

746 int magDiff;

747

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

749

750 if (a == b) {

751 mp_zero(c);

752 return MP_OKAY;

753 }

754

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

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

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

758 mp_zero(c);

759 res = MP_OKAY;

760 } else if (magDiff > 0) {

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

762 } else {

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

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

765 }

766

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

768 MP_SIGN(c) = MP_ZPOS;

769

770 CLEANUP:

771 return res;

772

773 } /* end mp_sub() */

774

775 /* }}} */

776

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

778

779 /*

780 mp_mul(a, b, c)

781

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

783 */

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

785 {

786 mp_digit *pb;

787 mp_int tmp;

788 mp_err res;

789 mp_size ib;

790 mp_size useda, usedb;

791

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

793

794 if (a == c) {

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

796 return res;

797 if (a == b)

798 b = &tmp;

799 a = &tmp;

800 } else if (b == c) {

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

802 return res;

803 b = &tmp;

804 } else {

805 MP_DIGITS(&tmp) = 0;

806 }

807

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

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

810 b = a;

811 a = xch;

812 }

813

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

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

816 goto CLEANUP;

817

818 #ifdef NSS_USE_COMBA

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

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

821 s_mp_mul_comba_4(a, b, c);

822 goto CLEANUP;

823 }

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

825 s_mp_mul_comba_8(a, b, c);

826 goto CLEANUP;

827 }

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

829 s_mp_mul_comba_16(a, b, c);

830 goto CLEANUP;

831 }

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

833 s_mp_mul_comba_32(a, b, c);

834 goto CLEANUP;

835 }

836 }

837 #endif

838

839 pb = MP_DIGITS(b);

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

841

842 /* Outer loop: Digits of b */

843 useda = MP_USED(a);

844 usedb = MP_USED(b);

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

846 mp_digit b_i = *pb++;

847

848 /* Inner product: Digits of a */

849 if (b_i)

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

851 else

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

853 }

854

855 s_mp_clamp(c);

856

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

858 SIGN(c) = ZPOS;

859 else

860 SIGN(c) = NEG;

861

862 CLEANUP:

863 mp_clear(&tmp);

864 return res;

865 } /* end mp_mul() */

866

867 /* }}} */

868

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

870

871 #if MP_SQUARE

872 /*

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

874 efficiently than a general multiplication, because many of the

875 computation steps are redundant when squaring. The inner product

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

877 iterations of the multiplication loop.

878 */

879

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

881 mp_err mp_sqr(const mp_int a, mp_int sqr)

882 {

883 mp_digit *pa;

884 mp_digit d;

885 mp_err res;

886 mp_size ix;

887 mp_int tmp;

888 int count;

889

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

891

892 if (a == sqr) {

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

894 return res;

895 a = &tmp;

896 } else {

897 DIGITS(&tmp) = 0;

898 res = MP_OKAY;

899 }

900

901 ix = 2 * MP_USED(a);

902 if (ix > MP_ALLOC(sqr)) {

903 MP_USED(sqr) = 1;

904 MP_CHECKOK( s_mp_grow(sqr, ix) );

905 }

906 MP_USED(sqr) = ix;

907 MP_DIGIT(sqr, 0) = 0;

908

909 #ifdef NSS_USE_COMBA

910 if (IS_POWER_OF_2(MP_USED(a))) {

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

912 s_mp_sqr_comba_4(a, sqr);

913 goto CLEANUP;

914 }

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

916 s_mp_sqr_comba_8(a, sqr);

917 goto CLEANUP;

918 }

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

920 s_mp_sqr_comba_16(a, sqr);

921 goto CLEANUP;

922 }

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

924 s_mp_sqr_comba_32(a, sqr);

925 goto CLEANUP;

926 }

927 }

928 #endif

929

930 pa = MP_DIGITS(a);

931 count = MP_USED(a) - 1;

932 if (count > 0) {

933 d = *pa++;

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

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

936 d = *pa++;

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

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

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

940

941 /* now sqr = 2 /

942 s_mp_mul_2(sqr);

943 } else {

944 MP_DIGIT(sqr, 1) = 0;

945 }

946

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

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

949

950 SIGN(sqr) = ZPOS;

951 s_mp_clamp(sqr);

952

953 CLEANUP:

954 mp_clear(&tmp);

955 return res;

956

957 } /* end mp_sqr() */

958 #endif

959

960 /* }}} */

961

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

963

964 /*

965 mp_div(a, b, q, r)

966

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

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

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

970 */

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

972 {

973 mp_err res;

974 mp_int pQ, pR;

975 mp_int qtmp, rtmp, btmp;

976 int cmp;

977 mp_sign signA;

978 mp_sign signB;

979

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

981

982 signA = MP_SIGN(a);

983 signB = MP_SIGN(b);

984

985 if(mp_cmp_z(b) == MP_EQ)

986 return MP_RANGE;

987

988 DIGITS(&qtmp) = 0;

989 DIGITS(&rtmp) = 0;

990 DIGITS(&btmp) = 0;

991

992 /* Set up some temporaries... */

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

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

995 pR = &rtmp;

996 } else {

997 MP_CHECKOK( mp_copy(a, r) );

998 pR = r;

999 }

1000

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

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

1003 pQ = &qtmp;

1004 } else {

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

1006 pQ = q;

1007 mp_zero(pQ);

1008 }

1009

1010 /*

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

1012 otherwise, we actually do the work required.

1013 */

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

1015 if (cmp) {

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

1017 mp_zero(pQ);

1018 } else {

1019 mp_set(pQ, 1);

1020 mp_zero(pR);

1021 }

1022 } else {

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

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

1025 }

1026

1027 /* Compute the signs for the output */

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

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

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

1031

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

1033 SIGN(pQ) = ZPOS;

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

1035 SIGN(pR) = ZPOS;

1036

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

1038 if(q && q != pQ)

1039 s_mp_exch(pQ, q);

1040

1041 if(r && r != pR)

1042 s_mp_exch(pR, r);

1043

1044 CLEANUP:

1045 mp_clear(&btmp);

1046 mp_clear(&rtmp);

1047 mp_clear(&qtmp);

1048

1049 return res;

1050

1051 } /* end mp_div() */

1052

1053 /* }}} */

1054

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

1056

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

1058 {

1059 mp_err res;

1060

1061 ARGCHK(a != NULL, MP_BADARG);

1062

1063 if(q) {

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

1065 return res;

1066 }

1067 if(r) {

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

1069 return res;

1070 }

1071 if(q) {

1072 s_mp_div_2d(q, d);

1073 }

1074 if(r) {

1075 s_mp_mod_2d(r, d);

1076 }

1077

1078 return MP_OKAY;

1079

1080 } /* end mp_div_2d() */

1081

1082 /* }}} */

1083

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

1085

1086 /*

1087 mp_expt(a, b, c)

1088

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

1090 standard iterative square-and-multiply technique.

1091 */

1092

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

1094 {

1095 mp_int s, x;

1096 mp_err res;

1097 mp_digit d;

1098 unsigned int dig, bit;

1099

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

1101

1102 if(mp_cmp_z(b) < 0)

1103 return MP_RANGE;

1104

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

1106 return res;

1107

1108 mp_set(&s, 1);

1109

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

1111 goto X;

1112

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

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

1115 d = DIGIT(b, dig);

1116

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

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

1119 if(d & 1) {

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

1121 goto CLEANUP;

1122 }

1123

1124 d >>= 1;

1125

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

1127 goto CLEANUP;

1128 }

1129 }

1130

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

1132 d = DIGIT(b, dig);

1133

1134 while(d) {

1135 if(d & 1) {

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

1137 goto CLEANUP;

1138 }

1139

1140 d >>= 1;

1141

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

1143 goto CLEANUP;

1144 }

1145

1146 if(mp_iseven(b))

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

1148

1149 res = mp_copy(&s, c);

1150

1151 CLEANUP:

1152 mp_clear(&x);

1153 X:

1154 mp_clear(&s);

1155

1156 return res;

1157

1158 } /* end mp_expt() */

1159

1160 /* }}} */

1161

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

1163

1164 /* Compute a = 2^k */

1165

1166 mp_err mp_2expt(mp_int *a, mp_digit k)

1167 {

1168 ARGCHK(a != NULL, MP_BADARG);

1169

1170 return s_mp_2expt(a, k);

1171

1172 } /* end mp_2expt() */

1173

1174 /* }}} */

1175

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

1177

1178 /*

1179 mp_mod(a, m, c)

1180

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

1182 */

1183

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

1185 {

1186 mp_err res;

1187 int mag;

1188

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

1190

1191 if(SIGN(m) == NEG)

1192 return MP_RANGE;

1193

1194 /*

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

1196 absolute value.

1197

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

1199 the sign (if a is negative).

1200

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

1202

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

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

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

1206 */

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

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

1209 return res;

1210

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

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

1213 return res;

1214 }

1215

1216 } else if(mag < 0) {

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

1218 return res;

1219

1220 if(mp_cmp_z(a) < 0) {

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

1222 return res;

1223

1224 }

1225

1226 } else {

1227 mp_zero(c);

1228

1229 }

1230

1231 return MP_OKAY;

1232

1233 } /* end mp_mod() */

1234

1235 /* }}} */

1236

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

1238

1239 /*

1240 mp_mod_d(a, d, c)

1241

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

1243 */

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

1245 {

1246 mp_err res;

1247 mp_digit rem;

1248

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

1250

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

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

1253 return res;

1254

1255 } else {

1256 if(SIGN(a) == NEG)

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

1258 else

1259 rem = DIGIT(a, 0);

1260 }

1261

1262 if(c)

1263 *c = rem;

1264

1265 return MP_OKAY;

1266

1267 } /* end mp_mod_d() */

1268

1269 /* }}} */

1270

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

1272

1273 /*

1274 mp_sqrt(a, b)

1275

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

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

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

1279 following two properties:

1280

1281 b^2 <= a

1282 (b+1)^2 >= a

1283

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

1285 */

1286 mp_err mp_sqrt(const mp_int a, mp_int b)

1287 {

1288 mp_int x, t;

1289 mp_err res;

1290 mp_size used;

1291

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

1293

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

1295 if(SIGN(a) == NEG)

1296 return MP_RANGE;

1297

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

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

1300 return mp_copy(a, b);

1301

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

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

1304 return res;

1305

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

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

1308 goto X;

1309

1310 used = MP_USED(&x);

1311 if (used > 1) {

1312 s_mp_rshd(&x, used / 2);

1313 }

1314

1315 for(;;) {

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

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

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

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

1320 goto CLEANUP;

1321

1322 /* t = t / 2x */

1323 s_mp_mul_2(&x);

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

1325 goto CLEANUP;

1326 s_mp_div_2(&x);

1327

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

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

1330 break;

1331

1332 /* x = x - t */

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

1334 goto CLEANUP;

1335

1336 }

1337

1338 /* Copy result to output parameter */

1339 MP_CHECKOK(mp_sub_d(&x, 1, &x));

1340 s_mp_exch(&x, b);

1341

1342 CLEANUP:

1343 mp_clear(&x);

1344 X:

1345 mp_clear(&t);

1346

1347 return res;

1348

1349 } /* end mp_sqrt() */

1350

1351 /* }}} */

1352

1353 /* }}} */

1354

1355 /------------------------------------------------------------------------/

1356 /* {{{ Modular arithmetic */

1357

1358 #if MP_MODARITH

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

1360

1361 /*

1362 mp_addmod(a, b, m, c)

1363

1364 Compute c = (a + b) mod m

1365 */

1366

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

1368 {

1369 mp_err res;

1370

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

1372

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

1374 return res;

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

1376 return res;

1377

1378 return MP_OKAY;

1379

1380 }

1381

1382 /* }}} */

1383

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

1385

1386 /*

1387 mp_submod(a, b, m, c)

1388

1389 Compute c = (a - b) mod m

1390 */

1391

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

1393 {

1394 mp_err res;

1395

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

1397

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

1399 return res;

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

1401 return res;

1402

1403 return MP_OKAY;

1404

1405 }

1406

1407 /* }}} */

1408

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

1410

1411 /*

1412 mp_mulmod(a, b, m, c)

1413

1414 Compute c = (a * b) mod m

1415 */

1416

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

1418 {

1419 mp_err res;

1420

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

1422

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

1424 return res;

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

1426 return res;

1427

1428 return MP_OKAY;

1429

1430 }

1431

1432 /* }}} */

1433

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

1435

1436 #if MP_SQUARE

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

1438 {

1439 mp_err res;

1440

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

1442

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

1444 return res;

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

1446 return res;

1447

1448 return MP_OKAY;

1449

1450 } /* end mp_sqrmod() */

1451 #endif

1452

1453 /* }}} */

1454

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

1456

1457 /*

1458 s_mp_exptmod(a, b, m, c)

1459

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

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

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

1463

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

1465 s_mp_reduce() below for details)

1466 */

1467

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

1469 {

1470 mp_int s, x, mu;

1471 mp_err res;

1472 mp_digit d;

1473 unsigned int dig, bit;

1474

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

1476

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

1478 return MP_RANGE;

1479

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

1481 return res;

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

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

1484 goto X;

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

1486 goto MU;

1487

1488 mp_set(&s, 1);

1489

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

1491 s_mp_add_d(&mu, 1);

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

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

1494 goto CLEANUP;

1495

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

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

1498 d = DIGIT(b, dig);

1499

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

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

1502 if(d & 1) {

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

1504 goto CLEANUP;

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

1506 goto CLEANUP;

1507 }

1508

1509 d >>= 1;

1510

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

1512 goto CLEANUP;

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

1514 goto CLEANUP;

1515 }

1516 }

1517

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

1519 d = DIGIT(b, dig);

1520

1521 while(d) {

1522 if(d & 1) {

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

1524 goto CLEANUP;

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

1526 goto CLEANUP;

1527 }

1528

1529 d >>= 1;

1530

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

1532 goto CLEANUP;

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

1534 goto CLEANUP;

1535 }

1536

1537 s_mp_exch(&s, c);

1538

1539 CLEANUP:

1540 mp_clear(&mu);

1541 MU:

1542 mp_clear(&x);

1543 X:

1544 mp_clear(&s);

1545

1546 return res;

1547

1548 } /* end s_mp_exptmod() */

1549

1550 /* }}} */

1551

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

1553

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

1555 {

1556 mp_int s, x;

1557 mp_err res;

1558

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

1560

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

1562 return res;

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

1564 goto X;

1565

1566 mp_set(&s, 1);

1567

1568 while(d != 0) {

1569 if(d & 1) {

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

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

1572 goto CLEANUP;

1573 }

1574

1575 d /= 2;

1576

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

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

1579 goto CLEANUP;

1580 }

1581

1582 s_mp_exch(&s, c);

1583

1584 CLEANUP:

1585 mp_clear(&x);

1586 X:

1587 mp_clear(&s);

1588

1589 return res;

1590

1591 } /* end mp_exptmod_d() */

1592

1593 /* }}} */

1594 #endif /* if MP_MODARITH */

1595

1596 /* }}} */

1597

1598 /------------------------------------------------------------------------/

1599 /* {{{ Comparison functions */

1600

1601 /* {{{ mp_cmp_z(a) */

1602

1603 /*

1604 mp_cmp_z(a)

1605

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

1607 */

1608

1609 int mp_cmp_z(const mp_int *a)

1610 {

1611 if(SIGN(a) == NEG)

1612 return MP_LT;

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

1614 return MP_EQ;

1615 else

1616 return MP_GT;

1617

1618 } /* end mp_cmp_z() */

1619

1620 /* }}} */

1621

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

1623

1624 /*

1625 mp_cmp_d(a, d)

1626

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

1628 */

1629

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

1631 {

1632 ARGCHK(a != NULL, MP_EQ);

1633

1634 if(SIGN(a) == NEG)

1635 return MP_LT;

1636

1637 return s_mp_cmp_d(a, d);

1638

1639 } /* end mp_cmp_d() */

1640

1641 /* }}} */

1642

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

1644

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

1646 {

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

1648

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

1650 int mag;

1651

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

1653 return MP_EQ;

1654

1655 if(SIGN(a) == ZPOS)

1656 return mag;

1657 else

1658 return -mag;

1659

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

1661 return MP_GT;

1662 } else {

1663 return MP_LT;

1664 }

1665

1666 } /* end mp_cmp() */

1667

1668 /* }}} */

1669

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

1671

1672 /*

1673 mp_cmp_mag(a, b)

1674

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

1676 */

1677

1678 int mp_cmp_mag(mp_int a, mp_int b)

1679 {

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

1681

1682 return s_mp_cmp(a, b);

1683

1684 } /* end mp_cmp_mag() */

1685

1686 /* }}} */

1687

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

1689

1690 /*

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

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

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

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

1695 */

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

1697 {

1698 mp_int tmp;

1699 int out;

1700

1701 ARGCHK(a != NULL, MP_EQ);

1702

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

1704 out = mp_cmp(a, &tmp);

1705 mp_clear(&tmp);

1706

1707 return out;

1708

1709 } /* end mp_cmp_int() */

1710

1711 /* }}} */

1712

1713 /* {{{ mp_isodd(a) */

1714

1715 /*

1716 mp_isodd(a)

1717

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

1719 */

1720 int mp_isodd(const mp_int *a)

1721 {

1722 ARGCHK(a != NULL, 0);

1723

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

1725

1726 } /* end mp_isodd() */

1727

1728 /* }}} */

1729

1730 /* {{{ mp_iseven(a) */

1731

1732 int mp_iseven(const mp_int *a)

1733 {

1734 return !mp_isodd(a);

1735

1736 } /* end mp_iseven() */

1737

1738 /* }}} */

1739

1740 /* }}} */

1741

1742 /------------------------------------------------------------------------/

1743 /* {{{ Number theoretic functions */

1744

1745 #if MP_NUMTH

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

1747

1748 /*

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

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

1751 */

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

1753 {

1754 mp_err res;

1755 mp_int u, v, t;

1756 mp_size k = 0;

1757

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

1759

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

1761 return MP_RANGE;

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

1763 return mp_copy(b, c);

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

1765 return mp_copy(a, c);

1766 }

1767

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

1769 return res;

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

1771 goto U;

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

1773 goto V;

1774

1775 SIGN(&u) = ZPOS;

1776 SIGN(&v) = ZPOS;

1777

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

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

1780 s_mp_div_2(&u);

1781 s_mp_div_2(&v);

1782 ++k;

1783 }

1784

1785 /* Initialize t */

1786 if(mp_isodd(&u)) {

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

1788 goto CLEANUP;

1789

1790 /* t = -v */

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

1792 SIGN(&t) = NEG;

1793 else

1794 SIGN(&t) = ZPOS;

1795

1796 } else {

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

1798 goto CLEANUP;

1799

1800 }

1801

1802 for(;;) {

1803 while(mp_iseven(&t)) {

1804 s_mp_div_2(&t);

1805 }

1806

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

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

1809 goto CLEANUP;

1810

1811 } else {

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

1813 goto CLEANUP;

1814

1815 /* v = -t */

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

1817 SIGN(&v) = NEG;

1818 else

1819 SIGN(&v) = ZPOS;

1820 }

1821

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

1823 goto CLEANUP;

1824

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

1826 break;

1827 }

1828

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

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

1831

1832 CLEANUP:

1833 mp_clear(&v);

1834 V:

1835 mp_clear(&u);

1836 U:

1837 mp_clear(&t);

1838

1839 return res;

1840

1841 } /* end mp_gcd() */

1842

1843 /* }}} */

1844

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

1846

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

1848

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

1850

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

1852 */

1853

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

1855 {

1856 mp_int gcd, prod;

1857 mp_err res;

1858

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

1860

1861 /* Set up temporaries */

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

1863 return res;

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

1865 goto GCD;

1866

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

1868 goto CLEANUP;

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

1870 goto CLEANUP;

1871

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

1873

1874 CLEANUP:

1875 mp_clear(&prod);

1876 GCD:

1877 mp_clear(&gcd);

1878

1879 return res;

1880

1881 } /* end mp_lcm() */

1882

1883 /* }}} */

1884

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

1886

1887 /*

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

1889

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

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

1892 based on the Stein algorithm used for mp_gcd()

1893 See algorithm 14.61 in Handbook of Applied Cryptogrpahy.

1894 */

1895

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

1897 {

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

1899 mp_int *clean[9];

1900 mp_err res;

1901 int last = -1;

1902

1903 if(mp_cmp_z(b) == 0)

1904 return MP_RANGE;

1905

1906 /* Initialize all these variables we need */

1907 MP_CHECKOK( mp_init(&u) );

1908 clean[++last] = &u;

1909 MP_CHECKOK( mp_init(&v) );

1910 clean[++last] = &v;

1911 MP_CHECKOK( mp_init(&gx) );

1912 clean[++last] = &gx;

1913 MP_CHECKOK( mp_init(&A) );

1914 clean[++last] = &A;

1915 MP_CHECKOK( mp_init(&B) );

1916 clean[++last] = &B;

1917 MP_CHECKOK( mp_init(&C) );

1918 clean[++last] = &C;

1919 MP_CHECKOK( mp_init(&D) );

1920 clean[++last] = &D;

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

1922 clean[++last] = &xc;

1923 mp_abs(&xc, &xc);

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

1925 clean[++last] = &yc;

1926 mp_abs(&yc, &yc);

1927

1928 mp_set(&gx, 1);

1929

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

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

1932 mp_size nx = mp_trailing_zeros(&xc);

1933 mp_size ny = mp_trailing_zeros(&yc);

1934 mp_size n = MP_MIN(nx, ny);

1935 s_mp_div_2d(&xc,n);

1936 s_mp_div_2d(&yc,n);

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

1938 }

1939

1940 mp_copy(&xc, &u);

1941 mp_copy(&yc, &v);

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

1943

1944 /* Loop through binary GCD algorithm */

1945 do {

1946 while(mp_iseven(&u)) {

1947 s_mp_div_2(&u);

1948

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

1950 s_mp_div_2(&A); s_mp_div_2(&B);

1951 } else {

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

1953 s_mp_div_2(&A);

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

1955 s_mp_div_2(&B);

1956 }

1957 }

1958

1959 while(mp_iseven(&v)) {

1960 s_mp_div_2(&v);

1961

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

1963 s_mp_div_2(&C); s_mp_div_2(&D);

1964 } else {

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

1966 s_mp_div_2(&C);

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

1968 s_mp_div_2(&D);

1969 }

1970 }

1971

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

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

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

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

1976 } else {

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

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

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

1980 }

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

1982

1983 /* copy results to output */

1984 if(x)

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

1986

1987 if(y)

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

1989

1990 if(g)

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

1992

1993 CLEANUP:

1994 while(last >= 0)

1995 mp_clear(clean[last--]);

1996

1997 return res;

1998

1999 } /* end mp_xgcd() */

2000

2001 /* }}} */

2002

2003 mp_size mp_trailing_zeros(const mp_int *mp)

2004 {

2005 mp_digit d;

2006 mp_size n = 0;

2007 unsigned int ix;

2008

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

2010 return n;

2011

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

2013 n += MP_DIGIT_BIT;

2014 if (!d)

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

2016 #if !defined(MP_USE_UINT_DIGIT)

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

2018 d >>= 32;

2019 n += 32;

2020 }

2021 #endif

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

2023 d >>= 16;

2024 n += 16;

2025 }

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

2027 d >>= 8;

2028 n += 8;

2029 }

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

2031 d >>= 4;

2032 n += 4;

2033 }

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

2035 d >>= 2;

2036 n += 2;

2037 }

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

2039 d >>= 1;

2040 n += 1;

2041 }

2042 #if MP_ARGCHK == 2

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

2044 #endif

2045 return n;

2046 }

2047

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

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

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

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

2052 */

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

2054 {

2055 mp_err res;

2056 mp_err k = 0;

2057 mp_int d, f, g;

2058

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

2060

2061 MP_DIGITS(&d) = 0;

2062 MP_DIGITS(&f) = 0;

2063 MP_DIGITS(&g) = 0;

2064 MP_CHECKOK( mp_init(&d) );

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

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

2067

2068 mp_set(c, 1);

2069 mp_zero(&d);

2070

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

2072 res = MP_UNDEF;

2073 } else

2074 for (;;) {

2075 int diff_sign;

2076 while (mp_iseven(&f)) {

2077 mp_size n = mp_trailing_zeros(&f);

2078 if (!n) {

2079 res = MP_UNDEF;

2080 goto CLEANUP;

2081 }

2082 s_mp_div_2d(&f, n);

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

2084 k += n;

2085 }

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

2087 res = k;

2088 break;

2089 }

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

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

2092 s_mp_exch(&f, &g);

2093 s_mp_exch(c, &d);

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

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

2096 break;

2097 }

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

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

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

2101 } else {

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

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

2104 }

2105 }

2106 if (res >= 0) {

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

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

2109 }

2110 res = k;

2111 }

2112

2113 CLEANUP:

2114 mp_clear(&d);

2115 mp_clear(&f);

2116 mp_clear(&g);

2117 return res;

2118 }

2119

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

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

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

2123 */

2124 mp_digit s_mp_invmod_radix(mp_digit P)

2125 {

2126 mp_digit T = P;

2127 T = 2 - (P T);

2128 T = 2 - (P T);

2129 T = 2 - (P T);

2130 T = 2 - (P T);

2131 #if !defined(MP_USE_UINT_DIGIT)

2132 T = 2 - (P T);

2133 T = 2 - (P T);

2134 #endif

2135 return T;

2136 }

2137

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

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

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

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

2142 */

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

2144 {

2145 int k_orig = k;

2146 mp_digit r;

2147 mp_size ix;

2148 mp_err res;

2149

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

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

2152 } else {

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

2154 }

2155

2156 /* make sure x is large enough */

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

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

2159 MP_CHECKOK( s_mp_pad(x, ix) );

2160

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

2162

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

2164 int j = MP_MIN(k, MP_DIGIT_BIT);

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

2166 if (j < MP_DIGIT_BIT) {

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

2168 }

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

2170 k -= j;

2171 }

2172 s_mp_clamp(x);

2173 s_mp_div_2d(x, k_orig);

2174 res = MP_OKAY;

2175

2176 CLEANUP:

2177 return res;

2178 }

2179

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

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

2182 {

2183 int k;

2184 mp_err res;

2185 mp_int x;

2186

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

2188

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

2190 return MP_RANGE;

2191 if (mp_iseven(m))

2192 return MP_UNDEF;

2193

2194 MP_DIGITS(&x) = 0;

2195

2196 if (a == c) {

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

2198 return res;

2199 if (a == m)

2200 m = &x;

2201 a = &x;

2202 } else if (m == c) {

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

2204 return res;

2205 m = &x;

2206 } else {

2207 MP_DIGITS(&x) = 0;

2208 }

2209

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

2211 k = res;

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

2213 CLEANUP:

2214 mp_clear(&x);

2215 return res;

2216 }

2217

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

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

2220 {

2221 mp_int g, x;

2222 mp_err res;

2223

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

2225

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

2227 return MP_RANGE;

2228

2229 MP_DIGITS(&g) = 0;

2230 MP_DIGITS(&x) = 0;

2231 MP_CHECKOK( mp_init(&x) );

2232 MP_CHECKOK( mp_init(&g) );

2233

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

2235

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

2237 res = MP_UNDEF;

2238 goto CLEANUP;

2239 }

2240

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

2242 SIGN(c) = SIGN(a);

2243

2244 CLEANUP:

2245 mp_clear(&x);

2246 mp_clear(&g);

2247

2248 return res;

2249 }

2250

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

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

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

2254 {

2255 mp_err res;

2256 mp_size ix = k + 4;

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

2258

2259 static const mp_digit d2 = 2;

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

2261

2262 if (mp_iseven(a))

2263 return MP_UNDEF;

2264 if (k <= MP_DIGIT_BIT) {

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

2266 if (k < MP_DIGIT_BIT)

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

2268 mp_set(c, i);

2269 return MP_OKAY;

2270 }

2271 MP_DIGITS(&t0) = 0;

2272 MP_DIGITS(&t1) = 0;

2273 MP_DIGITS(&val) = 0;

2274 MP_DIGITS(&tmp) = 0;

2275 MP_DIGITS(&two2k) = 0;

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

2277 s_mp_mod_2d(&val, k);

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

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

2280 MP_CHECKOK( mp_init(&tmp) );

2281 MP_CHECKOK( mp_init(&two2k) );

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

2283 do {

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

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

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

2287 s_mp_mod_2d(&t1, k);

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

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

2290 }

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

2292 break;

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

2294 } while (--ix > 0);

2295 if (!ix) {

2296 res = MP_UNDEF;

2297 } else {

2298 mp_exch(c, &t1);

2299 }

2300

2301 CLEANUP:

2302 mp_clear(&t0);

2303 mp_clear(&t1);

2304 mp_clear(&val);

2305 mp_clear(&tmp);

2306 mp_clear(&two2k);

2307 return res;

2308 }

2309

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

2311 {

2312 mp_err res;

2313 mp_size k;

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

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

2316 mp_int C2, tmp1, tmp2;

2317

2318 /static const mp_digit d1 = 1; /

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

2320

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

2322 k = res;

2323 return s_mp_invmod_2d(a, k, c);

2324 }

2325 MP_DIGITS(&oddFactor) = 0;

2326 MP_DIGITS(&evenFactor) = 0;

2327 MP_DIGITS(&oddPart) = 0;

2328 MP_DIGITS(&evenPart) = 0;

2329 MP_DIGITS(&C2) = 0;

2330 MP_DIGITS(&tmp1) = 0;

2331 MP_DIGITS(&tmp2) = 0;

2332

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

2334 MP_CHECKOK( mp_init(&evenFactor) );

2335 MP_CHECKOK( mp_init(&oddPart) );

2336 MP_CHECKOK( mp_init(&evenPart) );

2337 MP_CHECKOK( mp_init(&C2) );

2338 MP_CHECKOK( mp_init(&tmp1) );

2339 MP_CHECKOK( mp_init(&tmp2) );

2340

2341 k = mp_trailing_zeros(m);

2342 s_mp_div_2d(&oddFactor, k);

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

2344

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

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

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

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

2349

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

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

2352 * let m2 = evenFactor, v2 = evenPart.

2353 */

2354

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

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

2357

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

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

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

2361 s_mp_mod_2d(&tmp2, k);

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

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

2364 }

2365

2366 /* compute answer = v1 + um1 /

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

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

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

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

2371

2372 CLEANUP:

2373 mp_clear(&oddFactor);

2374 mp_clear(&evenFactor);

2375 mp_clear(&oddPart);

2376 mp_clear(&evenPart);

2377 mp_clear(&C2);

2378 mp_clear(&tmp1);

2379 mp_clear(&tmp2);

2380 return res;

2381 }

2382

2383

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

2385

2386 /*

2387 mp_invmod(a, m, c)

2388

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

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

2391 MP_UNDEF is returned, and there is no inverse.

2392 */

2393

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

2395 {

2396

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

2398

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

2400 return MP_RANGE;

2401

2402 if (mp_isodd(m)) {

2403 return s_mp_invmod_odd_m(a, m, c);

2404 }

2405 if (mp_iseven(a))

2406 return MP_UNDEF; /* not invertable */

2407

2408 return s_mp_invmod_even_m(a, m, c);

2409

2410 } /* end mp_invmod() */

2411

2412 /* }}} */

2413 #endif /* if MP_NUMTH */

2414

2415 /* }}} */

2416

2417 /------------------------------------------------------------------------/

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

2419

2420 #if MP_IOFUNC

2421 /*

2422 mp_print(mp, ofp)

2423

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

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

2426 */

2427

2428 void mp_print(mp_int mp, FILE ofp)

2429 {

2430 int ix;

2431

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

2433 return;

2434

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

2436

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

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

2439 }

2440

2441 } /* end mp_print() */

2442

2443 #endif /* if MP_IOFUNC */

2444

2445 /* }}} */

2446

2447 /------------------------------------------------------------------------/

2448 /* {{{ More I/O Functions */

2449

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

2451

2452 /*

2453 mp_read_raw(mp, str, len)

2454

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

2456 */

2457

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

2459 {

2460 int ix;

2461 mp_err res;

2462 unsigned char ustr = (unsigned char )str;

2463

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

2465

2466 mp_zero(mp);

2467

2468 /* Get sign from first byte */

2469 if(ustr[0])

2470 SIGN(mp) = NEG;

2471 else

2472 SIGN(mp) = ZPOS;

2473

2474 /* Read the rest of the digits */

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

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

2477 return res;

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

2479 return res;

2480 }

2481

2482 return MP_OKAY;

2483

2484 } /* end mp_read_raw() */

2485

2486 /* }}} */

2487

2488 /* {{{ mp_raw_size(mp) */

2489

2490 int mp_raw_size(mp_int *mp)

2491 {

2492 ARGCHK(mp != NULL, 0);

2493

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

2495

2496 } /* end mp_raw_size() */

2497

2498 /* }}} */

2499

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

2501

2502 mp_err mp_toraw(mp_int mp, char str)

2503 {

2504 int ix, jx, pos = 1;

2505

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

2507

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

2509

2510 /* Iterate over each digit... */

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

2512 mp_digit d = DIGIT(mp, ix);

2513

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

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

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

2517 }

2518 }

2519

2520 return MP_OKAY;

2521

2522 } /* end mp_toraw() */

2523

2524 /* }}} */

2525

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

2527

2528 /*

2529 mp_read_radix(mp, str, radix)

2530

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

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

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

2534 character or the end of the string.

2535 */

2536

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

2538 {

2539 int ix = 0, val = 0;

2540 mp_err res;

2541 mp_sign sig = ZPOS;

2542

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

2544 MP_BADARG);

2545

2546 mp_zero(mp);

2547

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

2549 while(str[ix] &&

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

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

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

2553 ++ix;

2554 }

2555

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

2557 sig = NEG;

2558 ++ix;

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

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

2561 ++ix;

2562 }

2563

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

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

2566 return res;

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

2568 return res;

2569 ++ix;

2570 }

2571

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

2573 SIGN(mp) = ZPOS;

2574 else

2575 SIGN(mp) = sig;

2576

2577 return MP_OKAY;

2578

2579 } /* end mp_read_radix() */

2580

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

2582 {

2583 int radix = default_radix;

2584 int cx;

2585 mp_sign sig = ZPOS;

2586 mp_err res;

2587

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

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

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

2591 cx != '-' &&

2592 cx != '+') {

2593 ++str;

2594 }

2595

2596 if (cx == '-') {

2597 sig = NEG;

2598 ++str;

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

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

2601 ++str;

2602 }

2603

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

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

2606 radix = 16;

2607 str += 2;

2608 } else {

2609 radix = 8;

2610 str++;

2611 }

2612 }

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

2614 if (res == MP_OKAY) {

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

2616 }

2617 return res;

2618 }

2619

2620 /* }}} */

2621

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

2623

2624 int mp_radix_size(mp_int *mp, int radix)

2625 {

2626 int bits;

2627

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

2629 return 0;

2630

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

2632

2633 return s_mp_outlen(bits, radix);

2634

2635 } /* end mp_radix_size() */

2636

2637 /* }}} */

2638

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

2640

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

2642 {

2643 int ix, pos = 0;

2644

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

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

2647

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

2649 str[0] = '0';

2650 str[1] = '\0';

2651 } else {

2652 mp_err res;

2653 mp_int tmp;

2654 mp_sign sgn;

2655 mp_digit rem, rdx = (mp_digit)radix;

2656 char ch;

2657

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

2659 return res;

2660

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

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

2663

2664 /* Generate output digits in reverse order */

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

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

2667 mp_clear(&tmp);

2668 return res;

2669 }

2670

2671 /* Generate digits, use capital letters */

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

2673

2674 str[pos++] = ch;

2675 }

2676

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

2678 if(sgn == NEG)

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

2680

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

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

2683

2684 /* Reverse the digits and sign indicator */

2685 ix = 0;

2686 while(ix < pos) {

2687 char tmp = str[ix];

2688

2689 str[ix] = str[pos];

2690 str[pos] = tmp;

2691 ++ix;

2692 --pos;

2693 }

2694

2695 mp_clear(&tmp);

2696 }

2697

2698 return MP_OKAY;

2699

2700 } /* end mp_toradix() */

2701

2702 /* }}} */

2703

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

2705

2706 int mp_tovalue(char ch, int r)

2707 {

2708 return s_mp_tovalue(ch, r);

2709

2710 } /* end mp_tovalue() */

2711

2712 /* }}} */

2713

2714 /* }}} */

2715

2716 /* {{{ mp_strerror(ec) */

2717

2718 /*

2719 mp_strerror(ec)

2720

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

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

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

2724 string.

2725 */

2726 const char *mp_strerror(mp_err ec)

2727 {

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

2729

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

2731 are accurate */

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

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

2734 } else {

2735 return mp_err_string[aec + 1];

2736 }

2737

2738 } /* end mp_strerror() */

2739

2740 /* }}} */

2741

2742 /========================================================================/

2743 /------------------------------------------------------------------------/

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

2745

2746 /* {{{ Memory management */

2747

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

2749

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

2751 mp_err s_mp_grow(mp_int *mp, mp_size min)

2752 {

2753 if(min > ALLOC(mp)) {

2754 mp_digit *tmp;

2755

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

2757 min = MP_ROUNDUP(min, s_mp_defprec);

2758

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

2760 return MP_MEM;

2761

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

2763

2764 #if MP_CRYPTO

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

2766 #endif

2767 s_mp_free(DIGITS(mp));

2768 DIGITS(mp) = tmp;

2769 ALLOC(mp) = min;

2770 }

2771

2772 return MP_OKAY;

2773

2774 } /* end s_mp_grow() */

2775

2776 /* }}} */

2777

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

2779

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

2781 mp_err s_mp_pad(mp_int *mp, mp_size min)

2782 {

2783 if(min > USED(mp)) {

2784 mp_err res;

2785

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

2787 if (min > ALLOC(mp)) {

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

2789 return res;

2790 } else {

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

2792 }

2793

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

2795 USED(mp) = min;

2796 }

2797

2798 return MP_OKAY;

2799

2800 } /* end s_mp_pad() */

2801

2802 /* }}} */

2803

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

2805

2806 #if MP_MACRO == 0

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

2808 void s_mp_setz(mp_digit *dp, mp_size count)

2809 {

2810 #if MP_MEMSET == 0

2811 int ix;

2812

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

2814 dp[ix] = 0;

2815 #else

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

2817 #endif

2818

2819 } /* end s_mp_setz() */

2820 #endif

2821

2822 /* }}} */

2823

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

2825

2826 #if MP_MACRO == 0

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

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

2829 {

2830 #if MP_MEMCPY == 0

2831 int ix;

2832

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

2834 dp[ix] = sp[ix];

2835 #else

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

2837 #endif

2838 ++mp_copies;

2839

2840 } /* end s_mp_copy() */

2841 #endif

2842

2843 /* }}} */

2844

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

2846

2847 #if MP_MACRO == 0

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

2849 void *s_mp_alloc(size_t nb, size_t ni)

2850 {

2851 ++mp_allocs;

2852 return calloc(nb, ni);

2853

2854 } /* end s_mp_alloc() */

2855 #endif

2856

2857 /* }}} */

2858

2859 /* {{{ s_mp_free(ptr) */

2860

2861 #if MP_MACRO == 0

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

2863 void s_mp_free(void *ptr)

2864 {

2865 if(ptr) {

2866 ++mp_frees;

2867 free(ptr);

2868 }

2869 } /* end s_mp_free() */

2870 #endif

2871

2872 /* }}} */

2873

2874 /* {{{ s_mp_clamp(mp) */

2875

2876 #if MP_MACRO == 0

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

2878 void s_mp_clamp(mp_int *mp)

2879 {

2880 mp_size used = MP_USED(mp);

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

2882 --used;

2883 MP_USED(mp) = used;

2884 } /* end s_mp_clamp() */

2885 #endif

2886

2887 /* }}} */

2888

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

2890

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

2892 void s_mp_exch(mp_int a, mp_int b)

2893 {

2894 mp_int tmp;

2895

2896 tmp = *a;

2897 a = b;

2898 *b = tmp;

2899

2900 } /* end s_mp_exch() */

2901

2902 /* }}} */

2903

2904 /* }}} */

2905

2906 /* {{{ Arithmetic helpers */

2907

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

2909

2910 /*

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

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

2913 alternative to multiplication by powers of the radix

2914 */

2915

2916 mp_err s_mp_lshd(mp_int *mp, mp_size p)

2917 {

2918 mp_err res;

2919 unsigned int ix;

2920

2921 if(p == 0)

2922 return MP_OKAY;

2923

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

2925 return MP_OKAY;

2926

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

2928 return res;

2929

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

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

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

2933 }

2934

2935 /* Fill the bottom digits with zeroes */

2936 for(ix = 0; (mp_size)ix < p; ix++)

2937 DIGIT(mp, ix) = 0;

2938

2939 return MP_OKAY;

2940

2941 } /* end s_mp_lshd() */

2942

2943 /* }}} */

2944

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

2946

2947 /*

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

2949 amounts to a bitwise shift of the value.

2950 */

2951 mp_err s_mp_mul_2d(mp_int *mp, mp_digit d)

2952 {

2953 mp_err res;

2954 mp_digit dshift, bshift;

2955 mp_digit mask;

2956

2957 ARGCHK(mp != NULL, MP_BADARG);

2958

2959 dshift = d / MP_DIGIT_BIT;

2960 bshift = d % MP_DIGIT_BIT;

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

2962 if (bshift) {

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

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

2965 } else {

2966 mask = 0;

2967 }

2968

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

2970 return res;

2971

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

2973 return res;

2974

2975 if (bshift) {

2976 mp_digit *pa = MP_DIGITS(mp);

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

2978 mp_digit prev = 0;

2979

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

2981 mp_digit x = *pa;

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

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

2984 }

2985 }

2986

2987 s_mp_clamp(mp);

2988 return MP_OKAY;

2989 } /* end s_mp_mul_2d() */

2990

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

2992

2993 /*

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

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

2996 end are lost. Cannot fail.

2997 */

2998

2999 void s_mp_rshd(mp_int *mp, mp_size p)

3000 {

3001 mp_size ix;

3002 mp_digit src, dst;

3003

3004 if(p == 0)

3005 return;

3006

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

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

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

3010 USED(mp) = 1;

3011 SIGN(mp) = ZPOS;

3012 return;

3013 }

3014

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

3016 dst = MP_DIGITS(mp);

3017 src = dst + p;

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

3019 dst++ = src++;

3020

3021 MP_USED(mp) -= p;

3022 /* Fill the top digits with zeroes */

3023 while (p-- > 0)

3024 *dst++ = 0;

3025

3026 #if 0

3027 /* Strip off any leading zeroes */

3028 s_mp_clamp(mp);

3029 #endif

3030

3031 } /* end s_mp_rshd() */

3032

3033 /* }}} */

3034

3035 /* {{{ s_mp_div_2(mp) */

3036

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

3038 void s_mp_div_2(mp_int *mp)

3039 {

3040 s_mp_div_2d(mp, 1);

3041

3042 } /* end s_mp_div_2() */

3043

3044 /* }}} */

3045

3046 /* {{{ s_mp_mul_2(mp) */

3047

3048 mp_err s_mp_mul_2(mp_int *mp)

3049 {

3050 mp_digit *pd;

3051 unsigned int ix, used;

3052 mp_digit kin = 0;

3053

3054 /* Shift digits leftward by 1 bit */

3055 used = MP_USED(mp);

3056 pd = MP_DIGITS(mp);

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

3058 mp_digit d = *pd;

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

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

3061 }

3062

3063 /* Deal with rollover from last digit */

3064 if (kin) {

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

3066 mp_err res;

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

3068 return res;

3069 }

3070

3071 DIGIT(mp, ix) = kin;

3072 USED(mp) += 1;

3073 }

3074

3075 return MP_OKAY;

3076

3077 } /* end s_mp_mul_2() */

3078

3079 /* }}} */

3080

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

3082

3083 /*

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

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

3086 division code

3087 */

3088 void s_mp_mod_2d(mp_int *mp, mp_digit d)

3089 {

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

3091 mp_size ix;

3092 mp_digit dmask;

3093

3094 if(ndig >= USED(mp))

3095 return;

3096

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

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

3099 DIGIT(mp, ndig) &= dmask;

3100

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

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

3103 DIGIT(mp, ix) = 0;

3104

3105 s_mp_clamp(mp);

3106

3107 } /* end s_mp_mod_2d() */

3108

3109 /* }}} */

3110

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

3112

3113 /*

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

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

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

3117 */

3118 void s_mp_div_2d(mp_int *mp, mp_digit d)

3119 {

3120 int ix;

3121 mp_digit save, next, mask;

3122

3123 s_mp_rshd(mp, d / DIGIT_BIT);

3124 d %= DIGIT_BIT;

3125 if (d) {

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

3127 save = 0;

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

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

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

3131 save = next;

3132 }

3133 }

3134 s_mp_clamp(mp);

3135

3136 } /* end s_mp_div_2d() */

3137

3138 /* }}} */

3139

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

3141

3142 /*

3143 s_mp_norm(a, b, *d)

3144

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

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

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

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

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

3150 back at the end of the division process.

3151 */

3152

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

3154 {

3155 mp_digit d;

3156 mp_digit mask;

3157 mp_digit b_msd;

3158 mp_err res = MP_OKAY;

3159

3160 d = 0;

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

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

3163 while (!(b_msd & mask)) {

3164 b_msd <<= 1;

3165 ++d;

3166 }

3167

3168 if (d) {

3169 MP_CHECKOK( s_mp_mul_2d(a, d) );

3170 MP_CHECKOK( s_mp_mul_2d(b, d) );

3171 }

3172

3173 *pd = d;

3174 CLEANUP:

3175 return res;

3176

3177 } /* end s_mp_norm() */

3178

3179 /* }}} */

3180

3181 /* }}} */

3182

3183 /* {{{ Primitive digit arithmetic */

3184

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

3186

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

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

3189 {

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

3191 mp_word w, k = 0;

3192 mp_size ix = 1;

3193

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

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

3196 k = CARRYOUT(w);

3197

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

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

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

3201 k = CARRYOUT(w);

3202 ++ix;

3203 }

3204

3205 if(k != 0) {

3206 mp_err res;

3207

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

3209 return res;

3210

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

3212 }

3213

3214 return MP_OKAY;

3215 #else

3216 mp_digit * pmp = MP_DIGITS(mp);

3217 mp_digit sum, mp_i, carry = 0;

3218 mp_err res = MP_OKAY;

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

3220

3221 mp_i = *pmp;

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

3223 carry = (sum < d);

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

3225 mp_i = *pmp;

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

3227 carry = !sum;

3228 }

3229 if (carry && !used) {

3230 /* mp is growing */

3231 used = MP_USED(mp);

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

3233 MP_DIGIT(mp, used) = carry;

3234 }

3235 CLEANUP:

3236 return res;

3237 #endif

3238 } /* end s_mp_add_d() */

3239

3240 /* }}} */

3241

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

3243

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

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

3246 {

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

3248 mp_word w, b = 0;

3249 mp_size ix = 1;

3250

3251 /* Compute initial subtraction */

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

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

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

3255

3256 /* Propagate borrows leftward */

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

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

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

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

3261 ++ix;

3262 }

3263

3264 /* Remove leading zeroes */

3265 s_mp_clamp(mp);

3266

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

3268 if(b)

3269 return MP_RANGE;

3270 else

3271 return MP_OKAY;

3272 #else

3273 mp_digit *pmp = MP_DIGITS(mp);

3274 mp_digit mp_i, diff, borrow;

3275 mp_size used = MP_USED(mp);

3276

3277 mp_i = *pmp;

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

3279 borrow = (diff > mp_i);

3280 while (borrow && --used) {

3281 mp_i = *pmp;

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

3283 borrow = (diff > mp_i);

3284 }

3285 s_mp_clamp(mp);

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

3287 #endif

3288 } /* end s_mp_sub_d() */

3289

3290 /* }}} */

3291

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

3293

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

3295 mp_err s_mp_mul_d(mp_int *a, mp_digit d)

3296 {

3297 mp_err res;

3298 mp_size used;

3299 int pow;

3300

3301 if (!d) {

3302 mp_zero(a);

3303 return MP_OKAY;

3304 }

3305 if (d == 1)

3306 return MP_OKAY;

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

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

3309 }

3310

3311 used = MP_USED(a);

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

3313

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

3315

3316 s_mp_clamp(a);

3317

3318 CLEANUP:

3319 return res;

3320

3321 } /* end s_mp_mul_d() */

3322

3323 /* }}} */

3324

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

3326

3327 /*

3328 s_mp_div_d(mp, d, r)

3329

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

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

3332 */

3333

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

3335 {

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

3337 mp_word w = 0, q;

3338 #else

3339 mp_digit w, q;

3340 #endif

3341 int ix;

3342 mp_err res;

3343 mp_int quot;

3344 mp_int rem;

3345

3346 if(d == 0)

3347 return MP_RANGE;

3348 if (d == 1) {

3349 if (r)

3350 *r = 0;

3351 return MP_OKAY;

3352 }

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

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

3355 mp_digit n = MP_DIGIT(mp,0);

3356 mp_digit rem;

3357

3358 q = n / d;

3359 rem = n % d;

3360 MP_DIGIT(mp,0) = q;

3361 if (r)

3362 *r = rem;

3363 return MP_OKAY;

3364 }

3365

3366 MP_DIGITS(&rem) = 0;

3367 MP_DIGITS(&quot) = 0;

3368 /* Make room for the quotient */

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

3370

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

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

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

3374

3375 if(w >= d) {

3376 q = w / d;

3377 w = w % d;

3378 } else {

3379 q = 0;

3380 }

3381

3382 s_mp_lshd(&quot, 1);

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

3384 }

3385 #else

3386 {

3387 mp_digit p;

3388 #if !defined(MP_ASSEMBLY_DIV_2DX1D)

3389 mp_digit norm;

3390 #endif

3391

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

3393

3394 #if !defined(MP_ASSEMBLY_DIV_2DX1D)

3395 MP_DIGIT(&quot, 0) = d;

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

3397 if (norm)

3398 d <<= norm;

3399 MP_DIGIT(&quot, 0) = 0;

3400 #endif

3401

3402 p = 0;

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

3404 w = DIGIT(&rem, ix);

3405

3406 if (p) {

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

3408 } else if (w >= d) {

3409 q = w / d;

3410 w = w % d;

3411 } else {

3412 q = 0;

3413 }

3414

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

3416 DIGIT(&quot, 0) = q;

3417 p = w;

3418 }

3419 #if !defined(MP_ASSEMBLY_DIV_2DX1D)

3420 if (norm)

3421 w >>= norm;

3422 #endif

3423 }

3424 #endif

3425

3426 /* Deliver the remainder, if desired */

3427 if(r)

3428 *r = (mp_digit)w;

3429

3430 s_mp_clamp(&quot);

3431 mp_exch(&quot, mp);

3432 CLEANUP:

3433 mp_clear(&quot);

3434 mp_clear(&rem);

3435

3436 return res;

3437 } /* end s_mp_div_d() */

3438

3439 /* }}} */

3440

3441

3442 /* }}} */

3443

3444 /* {{{ Primitive full arithmetic */

3445

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

3447

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

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

3450 {

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

3452 mp_word w = 0;

3453 #else

3454 mp_digit d, sum, carry = 0;

3455 #endif

3456 mp_digit pa, pb;

3457 mp_size ix;

3458 mp_size used;

3459 mp_err res;

3460

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

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

3463 return res;

3464

3465 /*

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

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

3468 padding step above, so there is no problem. If b had initially

3469 less precision, we'll have to make sure the carry out is duly

3470 propagated upward among the higher-order digits of the sum.

3471 */

3472 pa = MP_DIGITS(a);

3473 pb = MP_DIGITS(b);

3474 used = MP_USED(b);

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

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

3477 w = w + pa + pb++;

3478 *pa++ = ACCUM(w);

3479 w = CARRYOUT(w);

3480 #else

3481 d = *pa;

3482 sum = d + *pb++;

3483 d = (sum < d); /* detect overflow */

3484 *pa++ = sum += carry;

3485 carry = d + (sum < carry); /* detect overflow */

3486 #endif

3487 }

3488

3489 /* If we run out of 'b' digits before we're actually done, make

3490 sure the carries get propagated upward...

3491 */

3492 used = MP_USED(a);

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

3494 while (w && ix < used) {

3495 w = w + *pa;

3496 *pa++ = ACCUM(w);

3497 w = CARRYOUT(w);

3498 ++ix;

3499 }

3500 #else

3501 while (carry && ix < used) {

3502 sum = carry + *pa;

3503 *pa++ = sum;

3504 carry = !sum;

3505 ++ix;

3506 }

3507 #endif

3508

3509 /* If there's an overall carry out, increase precision and include

3510 it. We could have done this initially, but why touch the memory

3511 allocator unless we're sure we have to?

3512 */

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

3514 if (w) {

3515 if((res = s_mp_pad(a, used + 1)) != MP_OKAY)

3516 return res;

3517

3518 DIGIT(a, ix) = (mp_digit)w;

3519 }

3520 #else

3521 if (carry) {

3522 if((res = s_mp_pad(a, used + 1)) != MP_OKAY)

3523 return res;

3524

3525 DIGIT(a, used) = carry;

3526 }

3527 #endif

3528

3529 return MP_OKAY;

3530 } /* end s_mp_add() */

3531

3532 /* }}} */

3533

3534 /* Compute c = \|a\| + \|b\| / / magnitude addition */

3535 mp_err s_mp_add_3arg(const mp_int a, const mp_int b, mp_int *c)

3536 {

3537 mp_digit pa, pb, *pc;

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

3539 mp_word w = 0;

3540 #else

3541 mp_digit sum, carry = 0, d;

3542 #endif

3543 mp_size ix;

3544 mp_size used;

3545 mp_err res;

3546

3547 MP_SIGN(c) = MP_SIGN(a);

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

3549 const mp_int *xch = a;

3550 a = b;

3551 b = xch;

3552 }

3553

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

3555 if (MP_OKAY != (res = s_mp_pad(c, MP_USED(a))))

3556 return res;

3557

3558 /*

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

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

3561 exchange step above, so there is no problem. If b had initially

3562 less precision, we'll have to make sure the carry out is duly

3563 propagated upward among the higher-order digits of the sum.

3564 */

3565 pa = MP_DIGITS(a);

3566 pb = MP_DIGITS(b);

3567 pc = MP_DIGITS(c);

3568 used = MP_USED(b);

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

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

3571 w = w + pa++ + pb++;

3572 *pc++ = ACCUM(w);

3573 w = CARRYOUT(w);

3574 #else

3575 d = *pa++;

3576 sum = d + *pb++;

3577 d = (sum < d); /* detect overflow */

3578 *pc++ = sum += carry;

3579 carry = d + (sum < carry); /* detect overflow */

3580 #endif

3581 }

3582

3583 /* If we run out of 'b' digits before we're actually done, make

3584 sure the carries get propagated upward...

3585 */

3586 for (used = MP_USED(a); ix < used; ++ix) {

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

3588 w = w + *pa++;

3589 *pc++ = ACCUM(w);

3590 w = CARRYOUT(w);

3591 #else

3592 pc++ = sum = carry + pa++;

3593 carry = (sum < carry);

3594 #endif

3595 }

3596

3597 /* If there's an overall carry out, increase precision and include

3598 it. We could have done this initially, but why touch the memory

3599 allocator unless we're sure we have to?

3600 */

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

3602 if (w) {

3603 if((res = s_mp_pad(c, used + 1)) != MP_OKAY)

3604 return res;

3605

3606 DIGIT(c, used) = (mp_digit)w;

3607 ++used;

3608 }

3609 #else

3610 if (carry) {

3611 if((res = s_mp_pad(c, used + 1)) != MP_OKAY)

3612 return res;

3613

3614 DIGIT(c, used) = carry;

3615 ++used;

3616 }

3617 #endif

3618 MP_USED(c) = used;

3619 return MP_OKAY;

3620 }

3621 /* {{{ s_mp_add_offset(a, b, offset) */

3622

3623 /* Compute a = \|a\| + ( \|b\| * (RADIX ** offset) ) */

3624 mp_err s_mp_add_offset(mp_int a, mp_int b, mp_size offset)

3625 {

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

3627 mp_word w, k = 0;

3628 #else

3629 mp_digit d, sum, carry = 0;

3630 #endif

3631 mp_size ib;

3632 mp_size ia;

3633 mp_size lim;

3634 mp_err res;

3635

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

3637 lim = MP_USED(b) + offset;

3638 if((lim > USED(a)) && (res = s_mp_pad(a, lim)) != MP_OKAY)

3639 return res;

3640

3641 /*

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

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

3644 padding step above, so there is no problem. If b had initially

3645 less precision, we'll have to make sure the carry out is duly

3646 propagated upward among the higher-order digits of the sum.

3647 */

3648 lim = USED(b);

3649 for(ib = 0, ia = offset; ib < lim; ib++, ia++) {

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

3651 w = (mp_word)DIGIT(a, ia) + DIGIT(b, ib) + k;

3652 DIGIT(a, ia) = ACCUM(w);

3653 k = CARRYOUT(w);

3654 #else

3655 d = MP_DIGIT(a, ia);

3656 sum = d + MP_DIGIT(b, ib);

3657 d = (sum < d);

3658 MP_DIGIT(a,ia) = sum += carry;

3659 carry = d + (sum < carry);

3660 #endif

3661 }

3662

3663 /* If we run out of 'b' digits before we're actually done, make

3664 sure the carries get propagated upward...

3665 */

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

3667 for (lim = MP_USED(a); k && (ia < lim); ++ia) {

3668 w = (mp_word)DIGIT(a, ia) + k;

3669 DIGIT(a, ia) = ACCUM(w);

3670 k = CARRYOUT(w);

3671 }

3672 #else

3673 for (lim = MP_USED(a); carry && (ia < lim); ++ia) {

3674 d = MP_DIGIT(a, ia);

3675 MP_DIGIT(a,ia) = sum = d + carry;

3676 carry = (sum < d);

3677 }

3678 #endif

3679

3680 /* If there's an overall carry out, increase precision and include

3681 it. We could have done this initially, but why touch the memory

3682 allocator unless we're sure we have to?

3683 */

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

3685 if(k) {

3686 if((res = s_mp_pad(a, USED(a) + 1)) != MP_OKAY)

3687 return res;

3688

3689 DIGIT(a, ia) = (mp_digit)k;

3690 }

3691 #else

3692 if (carry) {

3693 if((res = s_mp_pad(a, lim + 1)) != MP_OKAY)

3694 return res;

3695

3696 DIGIT(a, lim) = carry;

3697 }

3698 #endif

3699 s_mp_clamp(a);

3700

3701 return MP_OKAY;

3702

3703 } /* end s_mp_add_offset() */

3704

3705 /* }}} */

3706

3707 /* {{{ s_mp_sub(a, b) */

3708

3709 /* Compute a = \|a\| - \|b\|, assumes \|a\| >= \|b\| */

3710 mp_err s_mp_sub(mp_int a, const mp_int b) /* magnitude subtract */

3711 {

3712 mp_digit pa, pb, *limit;

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

3714 mp_sword w = 0;

3715 #else

3716 mp_digit d, diff, borrow = 0;

3717 #endif

3718

3719 /*

3720 Subtract and propagate borrow. Up to the precision of b, this

3721 accounts for the digits of b; after that, we just make sure the

3722 carries get to the right place. This saves having to pad b out to

3723 the precision of a just to make the loops work right...

3724 */

3725 pa = MP_DIGITS(a);

3726 pb = MP_DIGITS(b);

3727 limit = pb + MP_USED(b);

3728 while (pb < limit) {

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

3730 w = w + pa - pb++;

3731 *pa++ = ACCUM(w);

3732 w >>= MP_DIGIT_BIT;

3733 #else

3734 d = *pa;

3735 diff = d - *pb++;

3736 d = (diff > d); /* detect borrow */

3737 if (borrow && --diff == MP_DIGIT_MAX)

3738 ++d;

3739 *pa++ = diff;

3740 borrow = d;

3741 #endif

3742 }

3743 limit = MP_DIGITS(a) + MP_USED(a);

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

3745 while (w && pa < limit) {

3746 w = w + *pa;

3747 *pa++ = ACCUM(w);

3748 w >>= MP_DIGIT_BIT;

3749 }

3750 #else

3751 while (borrow && pa < limit) {

3752 d = *pa;

3753 *pa++ = diff = d - borrow;

3754 borrow = (diff > d);

3755 }

3756 #endif

3757

3758 /* Clobber any leading zeroes we created */

3759 s_mp_clamp(a);

3760

3761 /*

3762 If there was a borrow out, then \|b\| > \|a\| in violation

3763 of our input invariant. We've already done the work,

3764 but we'll at least complain about it...

3765 */

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

3767 return w ? MP_RANGE : MP_OKAY;

3768 #else

3769 return borrow ? MP_RANGE : MP_OKAY;

3770 #endif

3771 } /* end s_mp_sub() */

3772

3773 /* }}} */

3774

3775 /* Compute c = \|a\| - \|b\|, assumes \|a\| >= \|b\| / / magnitude subtract */

3776 mp_err s_mp_sub_3arg(const mp_int a, const mp_int b, mp_int *c)

3777 {

3778 mp_digit pa, pb, *pc;

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

3780 mp_sword w = 0;

3781 #else

3782 mp_digit d, diff, borrow = 0;

3783 #endif

3784 int ix, limit;

3785 mp_err res;

3786

3787 MP_SIGN(c) = MP_SIGN(a);

3788

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

3790 if (MP_OKAY != (res = s_mp_pad(c, MP_USED(a))))

3791 return res;

3792

3793 /*

3794 Subtract and propagate borrow. Up to the precision of b, this

3795 accounts for the digits of b; after that, we just make sure the

3796 carries get to the right place. This saves having to pad b out to

3797 the precision of a just to make the loops work right...

3798 */

3799 pa = MP_DIGITS(a);

3800 pb = MP_DIGITS(b);

3801 pc = MP_DIGITS(c);

3802 limit = MP_USED(b);

3803 for (ix = 0; ix < limit; ++ix) {

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

3805 w = w + pa++ - pb++;

3806 *pc++ = ACCUM(w);

3807 w >>= MP_DIGIT_BIT;

3808 #else

3809 d = *pa++;

3810 diff = d - *pb++;

3811 d = (diff > d);

3812 if (borrow && --diff == MP_DIGIT_MAX)

3813 ++d;

3814 *pc++ = diff;

3815 borrow = d;

3816 #endif

3817 }

3818 for (limit = MP_USED(a); ix < limit; ++ix) {

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

3820 w = w + *pa++;

3821 *pc++ = ACCUM(w);

3822 w >>= MP_DIGIT_BIT;

3823 #else

3824 d = *pa++;

3825 *pc++ = diff = d - borrow;

3826 borrow = (diff > d);

3827 #endif

3828 }

3829

3830 /* Clobber any leading zeroes we created */

3831 MP_USED(c) = ix;

3832 s_mp_clamp(c);

3833

3834 /*

3835 If there was a borrow out, then \|b\| > \|a\| in violation

3836 of our input invariant. We've already done the work,

3837 but we'll at least complain about it...

3838 */

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

3840 return w ? MP_RANGE : MP_OKAY;

3841 #else

3842 return borrow ? MP_RANGE : MP_OKAY;

3843 #endif

3844 }

3845 /* {{{ s_mp_mul(a, b) */

3846

3847 /* Compute a = \|a\| * \|b\| */

3848 mp_err s_mp_mul(mp_int a, const mp_int b)

3849 {

3850 return mp_mul(a, b, a);

3851 } /* end s_mp_mul() */

3852

3853 /* }}} */

3854

3855 #if defined(MP_USE_UINT_DIGIT) && defined(MP_USE_LONG_LONG_MULTIPLY)

3856 /* This trick works on Sparc V8 CPUs with the Workshop compilers. */

3857 #define MP_MUL_DxD(a, b, Phi, Plo) \

3858 { unsigned long long product = (unsigned long long)a * b; \

3859 Plo = (mp_digit)product; \

3860 Phi = (mp_digit)(product >> MP_DIGIT_BIT); }

3861 #elif defined(OSF1)

3862 #define MP_MUL_DxD(a, b, Phi, Plo) \

3863 { Plo = asm ("mulq %a0, %a1, %v0", a, b);\

3864 Phi = asm ("umulh %a0, %a1, %v0", a, b); }

3865 #else

3866 #define MP_MUL_DxD(a, b, Phi, Plo) \

3867 { mp_digit a0b1, a1b0; \

3868 Plo = (a & MP_HALF_DIGIT_MAX) * (b & MP_HALF_DIGIT_MAX); \

3869 Phi = (a >> MP_HALF_DIGIT_BIT) * (b >> MP_HALF_DIGIT_BIT); \

3870 a0b1 = (a & MP_HALF_DIGIT_MAX) * (b >> MP_HALF_DIGIT_BIT); \

3871 a1b0 = (a >> MP_HALF_DIGIT_BIT) * (b & MP_HALF_DIGIT_MAX); \

3872 a1b0 += a0b1; \

3873 Phi += a1b0 >> MP_HALF_DIGIT_BIT; \

3874 if (a1b0 < a0b1) \

3875 Phi += MP_HALF_RADIX; \

3876 a1b0 <<= MP_HALF_DIGIT_BIT; \

3877 Plo += a1b0; \

3878 if (Plo < a1b0) \

3879 ++Phi; \

3880 }

3881 #endif

3882

3883 #if !defined(MP_ASSEMBLY_MULTIPLY)

3884 /* c = a * b */

3885 void s_mpv_mul_d(const mp_digit a, mp_size a_len, mp_digit b, mp_digit c)

3886 {

3887 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)

3888 mp_digit d = 0;

3889

3890 /* Inner product: Digits of a */

3891 while (a_len--) {

3892 mp_word w = ((mp_word)b * *a++) + d;

3893 *c++ = ACCUM(w);

3894 d = CARRYOUT(w);

3895 }

3896 *c = d;

3897 #else

3898 mp_digit carry = 0;

3899 while (a_len--) {

3900 mp_digit a_i = *a++;

3901 mp_digit a0b0, a1b1;

3902

3903 MP_MUL_DxD(a_i, b, a1b1, a0b0);

3904

3905 a0b0 += carry;

3906 if (a0b0 < carry)

3907 ++a1b1;

3908 *c++ = a0b0;

3909 carry = a1b1;

3910 }

3911 *c = carry;

3912 #endif

3913 }

3914

3915 /* c += a * b */

3916 void s_mpv_mul_d_add(const mp_digit *a, mp_size a_len, mp_digit b,

3917 mp_digit *c)

3918 {

3919 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)

3920 mp_digit d = 0;

3921

3922 /* Inner product: Digits of a */

3923 while (a_len--) {

3924 mp_word w = ((mp_word)b * a++) + c + d;

3925 *c++ = ACCUM(w);

3926 d = CARRYOUT(w);

3927 }

3928 *c = d;

3929 #else

3930 mp_digit carry = 0;

3931 while (a_len--) {

3932 mp_digit a_i = *a++;

3933 mp_digit a0b0, a1b1;

3934

3935 MP_MUL_DxD(a_i, b, a1b1, a0b0);

3936

3937 a0b0 += carry;

3938 if (a0b0 < carry)

3939 ++a1b1;

3940 a0b0 += a_i = *c;

3941 if (a0b0 < a_i)

3942 ++a1b1;

3943 *c++ = a0b0;

3944 carry = a1b1;

3945 }

3946 *c = carry;

3947 #endif

3948 }

3949

3950 /* Presently, this is only used by the Montgomery arithmetic code. */

3951 /* c += a * b */

3952 void s_mpv_mul_d_add_prop(const mp_digit a, mp_size a_len, mp_digit b, mp_digit c)

3953 {

3954 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)

3955 mp_digit d = 0;

3956

3957 /* Inner product: Digits of a */

3958 while (a_len--) {

3959 mp_word w = ((mp_word)b * a++) + c + d;

3960 *c++ = ACCUM(w);

3961 d = CARRYOUT(w);

3962 }

3963

3964 while (d) {

3965 mp_word w = (mp_word)*c + d;

3966 *c++ = ACCUM(w);

3967 d = CARRYOUT(w);

3968 }

3969 #else

3970 mp_digit carry = 0;

3971 while (a_len--) {

3972 mp_digit a_i = *a++;

3973 mp_digit a0b0, a1b1;

3974

3975 MP_MUL_DxD(a_i, b, a1b1, a0b0);

3976

3977 a0b0 += carry;

3978 if (a0b0 < carry)

3979 ++a1b1;

3980

3981 a0b0 += a_i = *c;

3982 if (a0b0 < a_i)

3983 ++a1b1;

3984

3985 *c++ = a0b0;

3986 carry = a1b1;

3987 }

3988 while (carry) {

3989 mp_digit c_i = *c;

3990 carry += c_i;

3991 *c++ = carry;

3992 carry = carry < c_i;

3993 }

3994 #endif

3995 }

3996 #endif

3997

3998 #if defined(MP_USE_UINT_DIGIT) && defined(MP_USE_LONG_LONG_MULTIPLY)

3999 /* This trick works on Sparc V8 CPUs with the Workshop compilers. */

4000 #define MP_SQR_D(a, Phi, Plo) \

4001 { unsigned long long square = (unsigned long long)a * a; \

4002 Plo = (mp_digit)square; \

4003 Phi = (mp_digit)(square >> MP_DIGIT_BIT); }

4004 #elif defined(OSF1)

4005 #define MP_SQR_D(a, Phi, Plo) \

4006 { Plo = asm ("mulq %a0, %a0, %v0", a);\

4007 Phi = asm ("umulh %a0, %a0, %v0", a); }

4008 #else

4009 #define MP_SQR_D(a, Phi, Plo) \

4010 { mp_digit Pmid; \

4011 Plo = (a & MP_HALF_DIGIT_MAX) * (a & MP_HALF_DIGIT_MAX); \

4012 Phi = (a >> MP_HALF_DIGIT_BIT) * (a >> MP_HALF_DIGIT_BIT); \

4013 Pmid = (a & MP_HALF_DIGIT_MAX) * (a >> MP_HALF_DIGIT_BIT); \

4014 Phi += Pmid >> (MP_HALF_DIGIT_BIT - 1); \

4015 Pmid <<= (MP_HALF_DIGIT_BIT + 1); \

4016 Plo += Pmid; \

4017 if (Plo < Pmid) \

4018 ++Phi; \

4019 }

4020 #endif

4021

4022 #if !defined(MP_ASSEMBLY_SQUARE)

4023 /* Add the squares of the digits of a to the digits of b. */

4024 void s_mpv_sqr_add_prop(const mp_digit pa, mp_size a_len, mp_digit ps)

4025 {

4026 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)

4027 mp_word w;

4028 mp_digit d;

4029 mp_size ix;

4030

4031 w = 0;

4032 #define ADD_SQUARE(n) \

4033 d = pa[n]; \

4034 w += (d * (mp_word)d) + ps[2*n]; \

4035 ps[2*n] = ACCUM(w); \

4036 w = (w >> DIGIT_BIT) + ps[2*n+1]; \

4037 ps[2*n+1] = ACCUM(w); \

4038 w = (w >> DIGIT_BIT)

4039

4040 for (ix = a_len; ix >= 4; ix -= 4) {

4041 ADD_SQUARE(0);

4042 ADD_SQUARE(1);

4043 ADD_SQUARE(2);

4044 ADD_SQUARE(3);

4045 pa += 4;

4046 ps += 8;

4047 }

4048 if (ix) {

4049 ps += 2*ix;

4050 pa += ix;

4051 switch (ix) {

4052 case 3: ADD_SQUARE(-3); /* FALLTHRU */

4053 case 2: ADD_SQUARE(-2); /* FALLTHRU */

4054 case 1: ADD_SQUARE(-1); /* FALLTHRU */

4055 case 0: break;

4056 }

4057 }

4058 while (w) {

4059 w += *ps;

4060 *ps++ = ACCUM(w);

4061 w = (w >> DIGIT_BIT);

4062 }

4063 #else

4064 mp_digit carry = 0;

4065 while (a_len--) {

4066 mp_digit a_i = *pa++;

4067 mp_digit a0a0, a1a1;

4068

4069 MP_SQR_D(a_i, a1a1, a0a0);

4070

4071 /* here a1a1 and a0a0 constitute a_i ** 2 */

4072 a0a0 += carry;

4073 if (a0a0 < carry)

4074 ++a1a1;

4075

4076 /* now add to ps */

4077 a0a0 += a_i = *ps;

4078 if (a0a0 < a_i)

4079 ++a1a1;

4080 *ps++ = a0a0;

4081 a1a1 += a_i = *ps;

4082 carry = (a1a1 < a_i);

4083 *ps++ = a1a1;

4084 }

4085 while (carry) {

4086 mp_digit s_i = *ps;

4087 carry += s_i;

4088 *ps++ = carry;

4089 carry = carry < s_i;

4090 }

4091 #endif

4092 }

4093 #endif

4094

4095 #if (defined(MP_NO_MP_WORD) \|\| defined(MP_NO_DIV_WORD)) \

4096 && !defined(MP_ASSEMBLY_DIV_2DX1D)

4097 /*

4098 ** Divide 64-bit (Nhi,Nlo) by 32-bit divisor, which must be normalized

4099 ** so its high bit is 1. This code is from NSPR.

4100 */

4101 mp_err s_mpv_div_2dx1d(mp_digit Nhi, mp_digit Nlo, mp_digit divisor,

4102 mp_digit qp, mp_digit rp)

4103 {

4104 mp_digit d1, d0, q1, q0;

4105 mp_digit r1, r0, m;

4106

4107 d1 = divisor >> MP_HALF_DIGIT_BIT;

4108 d0 = divisor & MP_HALF_DIGIT_MAX;

4109 r1 = Nhi % d1;

4110 q1 = Nhi / d1;

4111 m = q1 * d0;

4112 r1 = (r1 << MP_HALF_DIGIT_BIT) \| (Nlo >> MP_HALF_DIGIT_BIT);

4113 if (r1 < m) {

4114 q1--, r1 += divisor;

4115 if (r1 >= divisor && r1 < m) {

4116 q1--, r1 += divisor;

4117 }

4118 }

4119 r1 -= m;

4120 r0 = r1 % d1;

4121 q0 = r1 / d1;

4122 m = q0 * d0;

4123 r0 = (r0 << MP_HALF_DIGIT_BIT) \| (Nlo & MP_HALF_DIGIT_MAX);

4124 if (r0 < m) {

4125 q0--, r0 += divisor;

4126 if (r0 >= divisor && r0 < m) {

4127 q0--, r0 += divisor;

4128 }

4129 }

4130 if (qp)

4131 *qp = (q1 << MP_HALF_DIGIT_BIT) \| q0;

4132 if (rp)

4133 *rp = r0 - m;

4134 return MP_OKAY;

4135 }

4136 #endif

4137

4138 #if MP_SQUARE

4139 /* {{{ s_mp_sqr(a) */

4140

4141 mp_err s_mp_sqr(mp_int *a)

4142 {

4143 mp_err res;

4144 mp_int tmp;

4145

4146 if((res = mp_init_size(&tmp, 2 * USED(a))) != MP_OKAY)

4147 return res;

4148 res = mp_sqr(a, &tmp);

4149 if (res == MP_OKAY) {

4150 s_mp_exch(&tmp, a);

4151 }

4152 mp_clear(&tmp);

4153 return res;

4154 }

4155

4156 /* }}} */

4157 #endif

4158

4159 /* {{{ s_mp_div(a, b) */

4160

4161 /*

4162 s_mp_div(a, b)

4163

4164 Compute a = a / b and b = a mod b. Assumes b > a.

4165 */

4166

4167 mp_err s_mp_div(mp_int rem, / i: dividend, o: remainder */

4168 mp_int div, / i: divisor */

4169 mp_int quot) / i: 0; o: quotient */

4170 {

4171 mp_int part, t;

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

4173 mp_word q_msd;

4174 #else

4175 mp_digit q_msd;

4176 #endif

4177 mp_err res;

4178 mp_digit d;

4179 mp_digit div_msd;

4180 int ix;

4181

4182 if(mp_cmp_z(div) == 0)

4183 return MP_RANGE;

4184

4185 DIGITS(&t) = 0;

4186 /* Shortcut if divisor is power of two */

4187 if((ix = s_mp_ispow2(div)) >= 0) {

4188 MP_CHECKOK( mp_copy(rem, quot) );

4189 s_mp_div_2d(quot, (mp_digit)ix);

4190 s_mp_mod_2d(rem, (mp_digit)ix);

4191

4192 return MP_OKAY;

4193 }

4194

4195 MP_SIGN(rem) = ZPOS;

4196 MP_SIGN(div) = ZPOS;

4197 MP_SIGN(&part) = ZPOS;

4198

4199 /* A working temporary for division */

4200 MP_CHECKOK( mp_init_size(&t, MP_ALLOC(rem)));

4201

4202 /* Normalize to optimize guessing */

4203 MP_CHECKOK( s_mp_norm(rem, div, &d) );

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 int partExtended = 0; /* set to true if we need to extend part */

4214

4215 unusedRem = MP_USED(rem) - MP_USED(div);

4216 MP_DIGITS(&part) = MP_DIGITS(rem) + unusedRem;

4217 MP_ALLOC(&part) = MP_ALLOC(rem) - unusedRem;

4218 MP_USED(&part) = MP_USED(div);

4219

4220 /* We have now truncated the part of the remainder to the same length as

4221 * the divisor. If part is smaller than div, extend part by one digit. */

4222 if (s_mp_cmp(&part, div) < 0) {

4223 -- unusedRem;

4224 #if MP_ARGCHK == 2

4225 assert(unusedRem >= 0);

4226 #endif

4227 -- MP_DIGITS(&part);

4228 ++ MP_USED(&part);

4229 ++ MP_ALLOC(&part);

4230 partExtended = 1;

4231 }

4232

4233 /* Compute a guess for the next quotient digit */

4234 q_msd = MP_DIGIT(&part, MP_USED(&part) - 1);

4235 div_msd = MP_DIGIT(div, MP_USED(div) - 1);

4236 if (!partExtended) {

4237 /* In this case, q_msd /= div_msd is always 1. First, since div_msd is

4238 * normalized to have the high bit set, 2*div_msd > MP_DIGIT_MAX. Since

4239 * we didn't extend part, q_msd >= div_msd. Therefore we know that

4240 * div_msd <= q_msd <= MP_DIGIT_MAX < 2*div_msd. Dividing by div_msd we

4241 * get 1 <= q_msd/div_msd < 2. So q_msd /= div_msd must be 1. */

4242 q_msd = 1;

4243 } else {

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

4245 q_msd = (q_msd << MP_DIGIT_BIT) \| MP_DIGIT(&part, MP_USED(&part) - 2);

4246 q_msd /= div_msd;

4247 if (q_msd == RADIX)

4248 --q_msd;

4249 #else

4250 if (q_msd == div_msd) {

4251 q_msd = MP_DIGIT_MAX;

4252 } else {

4253 mp_digit r;

4254 MP_CHECKOK( s_mpv_div_2dx1d(q_msd, MP_DIGIT(&part, MP_USED(&part) - 2),

4255 div_msd, &q_msd, &r) );

4256 }

4257 #endif

4258 }

4259 #if MP_ARGCHK == 2

4260 assert(q_msd > 0); /* This case should never occur any more. */

4261 #endif

4262 if (q_msd <= 0)

4263 break;

4264

4265 /* See what that multiplies out to */

4266 mp_copy(div, &t);

4267 MP_CHECKOK( s_mp_mul_d(&t, (mp_digit)q_msd) );

4268

4269 /*

4270 If it's too big, back it off. We should not have to do this

4271 more than once, or, in rare cases, twice. Knuth describes a

4272 method by which this could be reduced to a maximum of once, but

4273 I didn't implement that here.

4274 * When using s_mpv_div_2dx1d, we may have to do this 3 times.

4275 */

4276 for (i = 4; s_mp_cmp(&t, &part) > 0 && i > 0; --i) {

4277 --q_msd;

4278 s_mp_sub(&t, div); /* t -= div */

4279 }

4280 if (i < 0) {

4281 res = MP_RANGE;

4282 goto CLEANUP;

4283 }

4284

4285 /* At this point, q_msd should be the right next digit */

4286 MP_CHECKOK( s_mp_sub(&part, &t) ); /* part -= t */

4287 s_mp_clamp(rem);

4288

4289 /*

4290 Include the digit in the quotient. We allocated enough memory

4291 for any quotient we could ever possibly get, so we should not

4292 have to check for failures here

4293 */

4294 MP_DIGIT(quot, unusedRem) = (mp_digit)q_msd;

4295 }

4296

4297 /* Denormalize remainder */

4298 if (d) {

4299 s_mp_div_2d(rem, d);

4300 }

4301

4302 s_mp_clamp(quot);

4303

4304 CLEANUP:

4305 mp_clear(&t);

4306

4307 return res;

4308

4309 } /* end s_mp_div() */

4310

4311

4312 /* }}} */

4313

4314 /* {{{ s_mp_2expt(a, k) */

4315

4316 mp_err s_mp_2expt(mp_int *a, mp_digit k)

4317 {

4318 mp_err res;

4319 mp_size dig, bit;

4320

4321 dig = k / DIGIT_BIT;

4322 bit = k % DIGIT_BIT;

4323

4324 mp_zero(a);

4325 if((res = s_mp_pad(a, dig + 1)) != MP_OKAY)

4326 return res;

4327

4328 DIGIT(a, dig) \|= ((mp_digit)1 << bit);

4329

4330 return MP_OKAY;

4331

4332 } /* end s_mp_2expt() */

4333

4334 /* }}} */

4335

4336 /* {{{ s_mp_reduce(x, m, mu) */

4337

4338 /*

4339 Compute Barrett reduction, x (mod m), given a precomputed value for

4340 mu = b^2k / m, where b = RADIX and k = #digits(m). This should be

4341 faster than straight division, when many reductions by the same

4342 value of m are required (such as in modular exponentiation). This

4343 can nearly halve the time required to do modular exponentiation,

4344 as compared to using the full integer divide to reduce.

4345

4346 This algorithm was derived from the _Handbook of Applied

4347 Cryptography_ by Menezes, Oorschot and VanStone, Ch. 14,

4348 pp. 603-604.

4349 */

4350

4351 mp_err s_mp_reduce(mp_int x, const mp_int m, const mp_int *mu)

4352 {

4353 mp_int q;

4354 mp_err res;

4355

4356 if((res = mp_init_copy(&q, x)) != MP_OKAY)

4357 return res;

4358

4359 s_mp_rshd(&q, USED(m) - 1); /* q1 = x / b^(k-1) */

4360 s_mp_mul(&q, mu); /* q2 = q1 * mu */

4361 s_mp_rshd(&q, USED(m) + 1); /* q3 = q2 / b^(k+1) */

4362

4363 /* x = x mod b^(k+1), quick (no division) */

4364 s_mp_mod_2d(x, DIGIT_BIT * (USED(m) + 1));

4365

4366 /* q = q * m mod b^(k+1), quick (no division) */

4367 s_mp_mul(&q, m);

4368 s_mp_mod_2d(&q, DIGIT_BIT * (USED(m) + 1));

4369

4370 /* x = x - q */

4371 if((res = mp_sub(x, &q, x)) != MP_OKAY)

4372 goto CLEANUP;

4373

4374 /* If x < 0, add b^(k+1) to it */

4375 if(mp_cmp_z(x) < 0) {

4376 mp_set(&q, 1);

4377 if((res = s_mp_lshd(&q, USED(m) + 1)) != MP_OKAY)

4378 goto CLEANUP;

4379 if((res = mp_add(x, &q, x)) != MP_OKAY)

4380 goto CLEANUP;

4381 }

4382

4383 /* Back off if it's too big */

4384 while(mp_cmp(x, m) >= 0) {

4385 if((res = s_mp_sub(x, m)) != MP_OKAY)

4386 break;

4387 }

4388

4389 CLEANUP:

4390 mp_clear(&q);

4391

4392 return res;

4393

4394 } /* end s_mp_reduce() */

4395

4396 /* }}} */

4397

4398 /* }}} */

4399

4400 /* {{{ Primitive comparisons */

4401

4402 /* {{{ s_mp_cmp(a, b) */

4403

4404 /* Compare \|a\| <=> \|b\|, return 0 if equal, <0 if a<b, >0 if a>b */

4405 int s_mp_cmp(const mp_int a, const mp_int b)

4406 {

4407 mp_size used_a = MP_USED(a);

4408 {

4409 mp_size used_b = MP_USED(b);

4410

4411 if (used_a > used_b)

4412 goto IS_GT;

4413 if (used_a < used_b)

4414 goto IS_LT;

4415 }

4416 {

4417 mp_digit pa, pb;

4418 mp_digit da = 0, db = 0;

4419

4420 #define CMP_AB(n) if ((da = pa[n]) != (db = pb[n])) goto done

4421

4422 pa = MP_DIGITS(a) + used_a;

4423 pb = MP_DIGITS(b) + used_a;

4424 while (used_a >= 4) {

4425 pa -= 4;

4426 pb -= 4;

4427 used_a -= 4;

4428 CMP_AB(3);

4429 CMP_AB(2);

4430 CMP_AB(1);

4431 CMP_AB(0);

4432 }

4433 while (used_a-- > 0 && ((da = --pa) == (db = --pb)))

4434 /* do nothing */;

4435 done:

4436 if (da > db)

4437 goto IS_GT;

4438 if (da < db)

4439 goto IS_LT;

4440 }

4441 return MP_EQ;

4442 IS_LT:

4443 return MP_LT;

4444 IS_GT:

4445 return MP_GT;

4446 } /* end s_mp_cmp() */

4447

4448 /* }}} */

4449

4450 /* {{{ s_mp_cmp_d(a, d) */

4451

4452 /* Compare \|a\| <=> d, return 0 if equal, <0 if a<d, >0 if a>d */

4453 int s_mp_cmp_d(const mp_int *a, mp_digit d)

4454 {

4455 if(USED(a) > 1)

4456 return MP_GT;

4457

4458 if(DIGIT(a, 0) < d)

4459 return MP_LT;

4460 else if(DIGIT(a, 0) > d)

4461 return MP_GT;

4462 else

4463 return MP_EQ;

4464

4465 } /* end s_mp_cmp_d() */

4466

4467 /* }}} */

4468

4469 /* {{{ s_mp_ispow2(v) */

4470

4471 /*

4472 Returns -1 if the value is not a power of two; otherwise, it returns

4473 k such that v = 2^k, i.e. lg(v).

4474 */

4475 int s_mp_ispow2(const mp_int *v)

4476 {

4477 mp_digit d;

4478 int extra = 0, ix;

4479

4480 ix = MP_USED(v) - 1;

4481 d = MP_DIGIT(v, ix); /* most significant digit of v */

4482

4483 extra = s_mp_ispow2d(d);

4484 if (extra < 0 \|\| ix == 0)

4485 return extra;

4486

4487 while (--ix >= 0) {

4488 if (DIGIT(v, ix) != 0)

4489 return -1; /* not a power of two */

4490 extra += MP_DIGIT_BIT;

4491 }

4492

4493 return extra;

4494

4495 } /* end s_mp_ispow2() */

4496

4497 /* }}} */

4498

4499 /* {{{ s_mp_ispow2d(d) */

4500

4501 int s_mp_ispow2d(mp_digit d)

4502 {

4503 if ((d != 0) && ((d & (d-1)) == 0)) { /* d is a power of 2 */

4504 int pow = 0;

4505 #if defined (MP_USE_UINT_DIGIT)

4506 if (d & 0xffff0000U)

4507 pow += 16;

4508 if (d & 0xff00ff00U)

4509 pow += 8;

4510 if (d & 0xf0f0f0f0U)

4511 pow += 4;

4512 if (d & 0xccccccccU)

4513 pow += 2;

4514 if (d & 0xaaaaaaaaU)

4515 pow += 1;

4516 #elif defined(MP_USE_LONG_LONG_DIGIT)

4517 if (d & 0xffffffff00000000ULL)

4518 pow += 32;

4519 if (d & 0xffff0000ffff0000ULL)

4520 pow += 16;

4521 if (d & 0xff00ff00ff00ff00ULL)

4522 pow += 8;

4523 if (d & 0xf0f0f0f0f0f0f0f0ULL)

4524 pow += 4;

4525 if (d & 0xccccccccccccccccULL)

4526 pow += 2;

4527 if (d & 0xaaaaaaaaaaaaaaaaULL)

4528 pow += 1;

4529 #elif defined(MP_USE_LONG_DIGIT)

4530 if (d & 0xffffffff00000000UL)

4531 pow += 32;

4532 if (d & 0xffff0000ffff0000UL)

4533 pow += 16;

4534 if (d & 0xff00ff00ff00ff00UL)

4535 pow += 8;

4536 if (d & 0xf0f0f0f0f0f0f0f0UL)

4537 pow += 4;

4538 if (d & 0xccccccccccccccccUL)

4539 pow += 2;

4540 if (d & 0xaaaaaaaaaaaaaaaaUL)

4541 pow += 1;

4542 #else

4543 #error "unknown type for mp_digit"

4544 #endif

4545 return pow;

4546 }

4547 return -1;

4548

4549 } /* end s_mp_ispow2d() */

4550

4551 /* }}} */

4552

4553 /* }}} */

4554

4555 /* {{{ Primitive I/O helpers */

4556

4557 /* {{{ s_mp_tovalue(ch, r) */

4558

4559 /*

4560 Convert the given character to its digit value, in the given radix.

4561 If the given character is not understood in the given radix, -1 is

4562 returned. Otherwise the digit's numeric value is returned.

4563

4564 The results will be odd if you use a radix < 2 or > 62, you are

4565 expected to know what you're up to.

4566 */

4567 int s_mp_tovalue(char ch, int r)

4568 {

4569 int val, xch;

4570

4571 if(r > 36)

4572 xch = ch;

4573 else

4574 xch = toupper(ch);

4575

4576 if(isdigit(xch))

4577 val = xch - '0';

4578 else if(isupper(xch))

4579 val = xch - 'A' + 10;

4580 else if(islower(xch))

4581 val = xch - 'a' + 36;

4582 else if(xch == '+')

4583 val = 62;

4584 else if(xch == '/')

4585 val = 63;

4586 else

4587 return -1;

4588

4589 if(val < 0 \|\| val >= r)

4590 return -1;

4591

4592 return val;

4593

4594 } /* end s_mp_tovalue() */

4595

4596 /* }}} */

4597

4598 /* {{{ s_mp_todigit(val, r, low) */

4599

4600 /*

4601 Convert val to a radix-r digit, if possible. If val is out of range

4602 for r, returns zero. Otherwise, returns an ASCII character denoting

4603 the value in the given radix.

4604

4605 The results may be odd if you use a radix < 2 or > 64, you are

4606 expected to know what you're doing.

4607 */

4608

4609 char s_mp_todigit(mp_digit val, int r, int low)

4610 {

4611 char ch;

4612

4613 if(val >= r)

4614 return 0;

4615

4616 ch = s_dmap_1[val];

4617

4618 if(r <= 36 && low)

4619 ch = tolower(ch);

4620

4621 return ch;

4622

4623 } /* end s_mp_todigit() */

4624

4625 /* }}} */

4626

4627 /* {{{ s_mp_outlen(bits, radix) */

4628

4629 /*

4630 Return an estimate for how long a string is needed to hold a radix

4631 r representation of a number with 'bits' significant bits, plus an

4632 extra for a zero terminator (assuming C style strings here)

4633 */

4634 int s_mp_outlen(int bits, int r)

4635 {

4636 return (int)((double)bits * LOG_V_2(r) + 1.5) + 1;

4637

4638 } /* end s_mp_outlen() */

4639

4640 /* }}} */

4641

4642 /* }}} */

4643

4644 /* {{{ mp_read_unsigned_octets(mp, str, len) */

4645 /* mp_read_unsigned_octets(mp, str, len)

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

4647 No sign bit, number is positive. Leading zeros ignored.

4648 */

4649

4650 mp_err

4651 mp_read_unsigned_octets(mp_int mp, const unsigned char str, mp_size len)

4652 {

4653 int count;

4654 mp_err res;

4655 mp_digit d;

4656

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

4658

4659 mp_zero(mp);

4660

4661 count = len % sizeof(mp_digit);

4662 if (count) {

4663 for (d = 0; count-- > 0; --len) {

4664 d = (d << 8) \| *str++;

4665 }

4666 MP_DIGIT(mp, 0) = d;

4667 }

4668

4669 /* Read the rest of the digits */

4670 for(; len > 0; len -= sizeof(mp_digit)) {

4671 for (d = 0, count = sizeof(mp_digit); count > 0; --count) {

4672 d = (d << 8) \| *str++;

4673 }

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

4675 if (!d)

4676 continue;

4677 } else {

4678 if((res = s_mp_lshd(mp, 1)) != MP_OKAY)

4679 return res;

4680 }

4681 MP_DIGIT(mp, 0) = d;

4682 }

4683 return MP_OKAY;

4684 } /* end mp_read_unsigned_octets() */

4685 /* }}} */

4686

4687 /* {{{ mp_unsigned_octet_size(mp) */

4688 unsigned int

4689 mp_unsigned_octet_size(const mp_int *mp)

4690 {

4691 unsigned int bytes;

4692 int ix;

4693 mp_digit d = 0;

4694

4695 ARGCHK(mp != NULL, MP_BADARG);

4696 ARGCHK(MP_ZPOS == SIGN(mp), MP_BADARG);

4697

4698 bytes = (USED(mp) * sizeof(mp_digit));

4699

4700 /* subtract leading zeros. */

4701 /* Iterate over each digit... */

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

4703 d = DIGIT(mp, ix);

4704 if (d)

4705 break;

4706 bytes -= sizeof(d);

4707 }

4708 if (!bytes)

4709 return 1;

4710

4711 /* Have MSD, check digit bytes, high order first */

4712 for(ix = sizeof(mp_digit) - 1; ix >= 0; ix--) {

4713 unsigned char x = (unsigned char)(d >> (ix * CHAR_BIT));

4714 if (x)

4715 break;

4716 --bytes;

4717 }

4718 return bytes;

4719 } /* end mp_unsigned_octet_size() */

4720 /* }}} */

4721

4722 /* {{{ mp_to_unsigned_octets(mp, str) */

4723 /* output a buffer of big endian octets no longer than specified. */

4724 mp_err

4725 mp_to_unsigned_octets(const mp_int mp, unsigned char str, mp_size maxlen)

4726 {

4727 int ix, pos = 0;

4728 unsigned int bytes;

4729

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

4731

4732 bytes = mp_unsigned_octet_size(mp);

4733 ARGCHK(bytes <= maxlen, MP_BADARG);

4734

4735 /* Iterate over each digit... */

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

4737 mp_digit d = DIGIT(mp, ix);

4738 int jx;

4739

4740 /* Unpack digit bytes, high order first */

4741 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {

4742 unsigned char x = (unsigned char)(d >> (jx * CHAR_BIT));

4743 if (!pos && !x) /* suppress leading zeros */

4744 continue;

4745 str[pos++] = x;

4746 }

4747 }

4748 if (!pos)

4749 str[pos++] = 0;

4750 return pos;

4751 } /* end mp_to_unsigned_octets() */

4752 /* }}} */

4753

4754 /* {{{ mp_to_signed_octets(mp, str) */

4755 /* output a buffer of big endian octets no longer than specified. */

4756 mp_err

4757 mp_to_signed_octets(const mp_int mp, unsigned char str, mp_size maxlen)

4758 {

4759 int ix, pos = 0;

4760 unsigned int bytes;

4761

4762 ARGCHK(mp != NULL && str != NULL && !SIGN(mp), MP_BADARG);

4763

4764 bytes = mp_unsigned_octet_size(mp);

4765 ARGCHK(bytes <= maxlen, MP_BADARG);

4766

4767 /* Iterate over each digit... */

4768 for(ix = USED(mp) - 1; ix >= 0; ix--) {

4769 mp_digit d = DIGIT(mp, ix);

4770 int jx;

4771

4772 /* Unpack digit bytes, high order first */

4773 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {

4774 unsigned char x = (unsigned char)(d >> (jx * CHAR_BIT));

4775 if (!pos) {

4776 if (!x) /* suppress leading zeros */

4777 continue;

4778 if (x & 0x80) { /* add one leading zero to make output positive. */

4779 ARGCHK(bytes + 1 <= maxlen, MP_BADARG);

4780 if (bytes + 1 > maxlen)

4781 return MP_BADARG;

4782 str[pos++] = 0;

4783 }

4784 }

4785 str[pos++] = x;

4786 }

4787 }

4788 if (!pos)

4789 str[pos++] = 0;

4790 return pos;

4791 } /* end mp_to_signed_octets() */

4792 /* }}} */

4793

4794 /* {{{ mp_to_fixlen_octets(mp, str) */

4795 /* output a buffer of big endian octets exactly as long as requested. */

4796 mp_err

4797 mp_to_fixlen_octets(const mp_int mp, unsigned char str, mp_size length)

4798 {

4799 int ix, pos = 0;

4800 unsigned int bytes;

4801

4802 ARGCHK(mp != NULL && str != NULL && !SIGN(mp), MP_BADARG);

4803

4804 bytes = mp_unsigned_octet_size(mp);

4805 ARGCHK(bytes <= length, MP_BADARG);

4806

4807 /* place any needed leading zeros */

4808 for (;length > bytes; --length) {

4809 *str++ = 0;

4810 }

4811

4812 /* Iterate over each digit... */

4813 for(ix = USED(mp) - 1; ix >= 0; ix--) {

4814 mp_digit d = DIGIT(mp, ix);

4815 int jx;

4816

4817 /* Unpack digit bytes, high order first */

4818 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {

4819 unsigned char x = (unsigned char)(d >> (jx * CHAR_BIT));

4820 if (!pos && !x) /* suppress leading zeros */

4821 continue;

4822 str[pos++] = x;

4823 }

4824 }

4825 if (!pos)

4826 str[pos++] = 0;

4827 return MP_OKAY;

4828 } /* end mp_to_fixlen_octets() */

4829 /* }}} */

4830

4831

4832 /------------------------------------------------------------------------/

4833 /* HERE THERE BE DRAGONS */

OLD	NEW

« no previous file with comments | « nss/lib/freebl/mpi/mpi.h ('k') | nss/lib/freebl/mpi/mpi-config.h » ('j') | no next file with comments »