1 /*
2 * Copyright (c) 2007, 2014, Oracle and/or its affiliates. All rights reserved.
3 * Use is subject to license terms.
4 *
5 * This library is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 2.1 of the License, or (at your option) any later version.
9 *
10 * This library is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with this library; if not, write to the Free Software Foundation,
17 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 *
19 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
20 * or visit www.oracle.com if you need additional information or have any
21 * questions.
22 */
23
24 /* *********************************************************************
25 *
26 * The Original Code is the MPI Arbitrary Precision Integer Arithmetic library.
27 *
28 * The Initial Developer of the Original Code is
29 * Michael J. Fromberger.
30 * Portions created by the Initial Developer are Copyright (C) 1998
31 * the Initial Developer. All Rights Reserved.
32 *
33 * Contributor(s):
34 * Netscape Communications Corporation
35 * Douglas Stebila <douglas@stebila.ca> of Sun Laboratories.
36 *
37 * Last Modified Date from the Original Code: June 2014
38 *********************************************************************** */
39
40 /* Arbitrary precision integer arithmetic library */
41
42 #include "mpi-priv.h"
43 #if defined(OSF1)
44 #include <c_asm.h>
45 #endif
46
47 #if MP_LOGTAB
48 /*
49 A table of the logs of 2 for various bases (the 0 and 1 entries of
50 this table are meaningless and should not be referenced).
51
52 This table is used to compute output lengths for the mp_toradix()
53 function. Since a number n in radix r takes up about log_r(n)
54 digits, we estimate the output size by taking the least integer
55 greater than log_r(n), where:
56
57 log_r(n) = log_2(n) * log_r(2)
642
643 if((res = mp_init(&s, FLAG(a))) != MP_OKAY)
644 return res;
645 if((res = mp_init_copy(&x, a)) != MP_OKAY)
646 goto X;
647
648 DIGIT(&s, 0) = 1;
649
650 while(d != 0) {
651 if(d & 1) {
652 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
653 goto CLEANUP;
654 }
655
656 d /= 2;
657
658 if((res = s_mp_sqr(&x)) != MP_OKAY)
659 goto CLEANUP;
660 }
661
662 s_mp_exch(&s, c);
663
664 CLEANUP:
665 mp_clear(&x);
666 X:
667 mp_clear(&s);
668
669 return res;
670
671 } /* end mp_expt_d() */
672
673 /* }}} */
674
675 /* }}} */
676
677 /*------------------------------------------------------------------------*/
678 /* {{{ Full arithmetic */
679
680 /* {{{ mp_abs(a, b) */
681
1592 return res;
1593 if((res = mp_init_copy(&x, a)) != MP_OKAY)
1594 goto X;
1595
1596 mp_set(&s, 1);
1597
1598 while(d != 0) {
1599 if(d & 1) {
1600 if((res = s_mp_mul(&s, &x)) != MP_OKAY ||
1601 (res = mp_mod(&s, m, &s)) != MP_OKAY)
1602 goto CLEANUP;
1603 }
1604
1605 d /= 2;
1606
1607 if((res = s_mp_sqr(&x)) != MP_OKAY ||
1608 (res = mp_mod(&x, m, &x)) != MP_OKAY)
1609 goto CLEANUP;
1610 }
1611
1612 s_mp_exch(&s, c);
1613
1614 CLEANUP:
1615 mp_clear(&x);
1616 X:
1617 mp_clear(&s);
1618
1619 return res;
1620
1621 } /* end mp_exptmod_d() */
1622
1623 /* }}} */
1624 #endif /* if MP_MODARITH */
1625
1626 /* }}} */
1627
1628 /*------------------------------------------------------------------------*/
1629 /* {{{ Comparison functions */
1630
1631 /* {{{ mp_cmp_z(a) */
4166 q0--, r0 += divisor;
4167 if (r0 >= divisor && r0 < m) {
4168 q0--, r0 += divisor;
4169 }
4170 }
4171 if (qp)
4172 *qp = (q1 << MP_HALF_DIGIT_BIT) | q0;
4173 if (rp)
4174 *rp = r0 - m;
4175 return MP_OKAY;
4176 }
4177 #endif
4178
4179 #if MP_SQUARE
4180 /* {{{ s_mp_sqr(a) */
4181
4182 mp_err s_mp_sqr(mp_int *a)
4183 {
4184 mp_err res;
4185 mp_int tmp;
4186
4187 if((res = mp_init_size(&tmp, 2 * USED(a), FLAG(a))) != MP_OKAY)
4188 return res;
4189 res = mp_sqr(a, &tmp);
4190 if (res == MP_OKAY) {
4191 s_mp_exch(&tmp, a);
4192 }
4193 mp_clear(&tmp);
4194 return res;
4195 }
4196
4197 /* }}} */
4198 #endif
4199
4200 /* {{{ s_mp_div(a, b) */
4201
4202 /*
4203 s_mp_div(a, b)
4204
4205 Compute a = a / b and b = a mod b. Assumes b > a.
4206 */
4207
4208 mp_err s_mp_div(mp_int *rem, /* i: dividend, o: remainder */
4209 mp_int *div, /* i: divisor */
4210 mp_int *quot) /* i: 0; o: quotient */
4211 {
4212 mp_int part, t;
4213 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)
4214 mp_word q_msd;
4215 #else
4216 mp_digit q_msd;
4217 #endif
4218 mp_err res;
4219 mp_digit d;
4220 mp_digit div_msd;
4221 int ix;
4222
4223 if(mp_cmp_z(div) == 0)
4224 return MP_RANGE;
4225
4226 /* Shortcut if divisor is power of two */
4227 if((ix = s_mp_ispow2(div)) >= 0) {
4228 MP_CHECKOK( mp_copy(rem, quot) );
4229 s_mp_div_2d(quot, (mp_digit)ix);
4230 s_mp_mod_2d(rem, (mp_digit)ix);
4231
4232 return MP_OKAY;
4233 }
4234
4235 DIGITS(&t) = 0;
4236 MP_SIGN(rem) = ZPOS;
4237 MP_SIGN(div) = ZPOS;
4238
4239 /* A working temporary for division */
4240 MP_CHECKOK( mp_init_size(&t, MP_ALLOC(rem), FLAG(rem)));
4241
|
1 /*
2 * Copyright (c) 2007, 2016, Oracle and/or its affiliates. All rights reserved.
3 * Use is subject to license terms.
4 *
5 * This library is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 2.1 of the License, or (at your option) any later version.
9 *
10 * This library is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with this library; if not, write to the Free Software Foundation,
17 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 *
19 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
20 * or visit www.oracle.com if you need additional information or have any
21 * questions.
22 */
23
24 /* *********************************************************************
25 *
26 * The Original Code is the MPI Arbitrary Precision Integer Arithmetic library.
27 *
28 * The Initial Developer of the Original Code is
29 * Michael J. Fromberger.
30 * Portions created by the Initial Developer are Copyright (C) 1998
31 * the Initial Developer. All Rights Reserved.
32 *
33 * Contributor(s):
34 * Netscape Communications Corporation
35 * Douglas Stebila <douglas@stebila.ca> of Sun Laboratories.
36 *
37 * Last Modified Date from the Original Code: Nov 2016
38 *********************************************************************** */
39
40 /* Arbitrary precision integer arithmetic library */
41
42 #include "mpi-priv.h"
43 #if defined(OSF1)
44 #include <c_asm.h>
45 #endif
46
47 #if MP_LOGTAB
48 /*
49 A table of the logs of 2 for various bases (the 0 and 1 entries of
50 this table are meaningless and should not be referenced).
51
52 This table is used to compute output lengths for the mp_toradix()
53 function. Since a number n in radix r takes up about log_r(n)
54 digits, we estimate the output size by taking the least integer
55 greater than log_r(n), where:
56
57 log_r(n) = log_2(n) * log_r(2)
642
643 if((res = mp_init(&s, FLAG(a))) != MP_OKAY)
644 return res;
645 if((res = mp_init_copy(&x, a)) != MP_OKAY)
646 goto X;
647
648 DIGIT(&s, 0) = 1;
649
650 while(d != 0) {
651 if(d & 1) {
652 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
653 goto CLEANUP;
654 }
655
656 d /= 2;
657
658 if((res = s_mp_sqr(&x)) != MP_OKAY)
659 goto CLEANUP;
660 }
661
662 s.flag = (mp_flag)0;
663 s_mp_exch(&s, c);
664
665 CLEANUP:
666 mp_clear(&x);
667 X:
668 mp_clear(&s);
669
670 return res;
671
672 } /* end mp_expt_d() */
673
674 /* }}} */
675
676 /* }}} */
677
678 /*------------------------------------------------------------------------*/
679 /* {{{ Full arithmetic */
680
681 /* {{{ mp_abs(a, b) */
682
1593 return res;
1594 if((res = mp_init_copy(&x, a)) != MP_OKAY)
1595 goto X;
1596
1597 mp_set(&s, 1);
1598
1599 while(d != 0) {
1600 if(d & 1) {
1601 if((res = s_mp_mul(&s, &x)) != MP_OKAY ||
1602 (res = mp_mod(&s, m, &s)) != MP_OKAY)
1603 goto CLEANUP;
1604 }
1605
1606 d /= 2;
1607
1608 if((res = s_mp_sqr(&x)) != MP_OKAY ||
1609 (res = mp_mod(&x, m, &x)) != MP_OKAY)
1610 goto CLEANUP;
1611 }
1612
1613 s.flag = (mp_flag)0;
1614 s_mp_exch(&s, c);
1615
1616 CLEANUP:
1617 mp_clear(&x);
1618 X:
1619 mp_clear(&s);
1620
1621 return res;
1622
1623 } /* end mp_exptmod_d() */
1624
1625 /* }}} */
1626 #endif /* if MP_MODARITH */
1627
1628 /* }}} */
1629
1630 /*------------------------------------------------------------------------*/
1631 /* {{{ Comparison functions */
1632
1633 /* {{{ mp_cmp_z(a) */
4168 q0--, r0 += divisor;
4169 if (r0 >= divisor && r0 < m) {
4170 q0--, r0 += divisor;
4171 }
4172 }
4173 if (qp)
4174 *qp = (q1 << MP_HALF_DIGIT_BIT) | q0;
4175 if (rp)
4176 *rp = r0 - m;
4177 return MP_OKAY;
4178 }
4179 #endif
4180
4181 #if MP_SQUARE
4182 /* {{{ s_mp_sqr(a) */
4183
4184 mp_err s_mp_sqr(mp_int *a)
4185 {
4186 mp_err res;
4187 mp_int tmp;
4188 tmp.flag = (mp_flag)0;
4189
4190 if((res = mp_init_size(&tmp, 2 * USED(a), FLAG(a))) != MP_OKAY)
4191 return res;
4192 res = mp_sqr(a, &tmp);
4193 if (res == MP_OKAY) {
4194 s_mp_exch(&tmp, a);
4195 }
4196 mp_clear(&tmp);
4197 return res;
4198 }
4199
4200 /* }}} */
4201 #endif
4202
4203 /* {{{ s_mp_div(a, b) */
4204
4205 /*
4206 s_mp_div(a, b)
4207
4208 Compute a = a / b and b = a mod b. Assumes b > a.
4209 */
4210
4211 mp_err s_mp_div(mp_int *rem, /* i: dividend, o: remainder */
4212 mp_int *div, /* i: divisor */
4213 mp_int *quot) /* i: 0; o: quotient */
4214 {
4215 mp_int part, t;
4216 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)
4217 mp_word q_msd;
4218 #else
4219 mp_digit q_msd;
4220 #endif
4221 mp_err res;
4222 mp_digit d;
4223 mp_digit div_msd;
4224 int ix;
4225
4226 t.dp = (mp_digit)0;
4227
4228 if(mp_cmp_z(div) == 0)
4229 return MP_RANGE;
4230
4231 /* Shortcut if divisor is power of two */
4232 if((ix = s_mp_ispow2(div)) >= 0) {
4233 MP_CHECKOK( mp_copy(rem, quot) );
4234 s_mp_div_2d(quot, (mp_digit)ix);
4235 s_mp_mod_2d(rem, (mp_digit)ix);
4236
4237 return MP_OKAY;
4238 }
4239
4240 DIGITS(&t) = 0;
4241 MP_SIGN(rem) = ZPOS;
4242 MP_SIGN(div) = ZPOS;
4243
4244 /* A working temporary for division */
4245 MP_CHECKOK( mp_init_size(&t, MP_ALLOC(rem), FLAG(rem)));
4246
|