lib/mpi/mpi-div.c

0001 /* mpi-div.c  -  MPI functions
0002  * Copyright (C) 1994, 1996, 1998, 2001, 2002,
0003  *               2003 Free Software Foundation, Inc.
0004  *
0005  * This file is part of Libgcrypt.
0006  *
0007  * Note: This code is heavily based on the GNU MP Library.
0008  *   Actually it's the same code with only minor changes in the
0009  *   way the data is stored; this is to support the abstraction
0010  *   of an optional secure memory allocation which may be used
0011  *   to avoid revealing of sensitive data due to paging etc.
0012  */
0013
0014 #include "mpi-internal.h"
0015 #include "longlong.h"
0016
0017 void mpi_tdiv_qr(MPI quot, MPI rem, MPI num, MPI den);
0018 void mpi_fdiv_qr(MPI quot, MPI rem, MPI dividend, MPI divisor);
0019
0020 void mpi_fdiv_r(MPI rem, MPI dividend, MPI divisor)
0021 {
0022     int divisor_sign = divisor->sign;
0023     MPI temp_divisor = NULL;
0024
0025     /* We need the original value of the divisor after the remainder has been
0026      * preliminary calculated.  We have to copy it to temporary space if it's
0027      * the same variable as REM.
0028      */
0029     if (rem == divisor) {
0030         temp_divisor = mpi_copy(divisor);
0031         divisor = temp_divisor;
0032     }
0033
0034     mpi_tdiv_r(rem, dividend, divisor);
0035
0036     if (((divisor_sign?1:0) ^ (dividend->sign?1:0)) && rem->nlimbs)
0037         mpi_add(rem, rem, divisor);
0038
0039     if (temp_divisor)
0040         mpi_free(temp_divisor);
0041 }
0042
0043 void mpi_fdiv_q(MPI quot, MPI dividend, MPI divisor)
0044 {
0045     MPI tmp = mpi_alloc(mpi_get_nlimbs(quot));
0046     mpi_fdiv_qr(quot, tmp, dividend, divisor);
0047     mpi_free(tmp);
0048 }
0049
0050 void mpi_fdiv_qr(MPI quot, MPI rem, MPI dividend, MPI divisor)
0051 {
0052     int divisor_sign = divisor->sign;
0053     MPI temp_divisor = NULL;
0054
0055     if (quot == divisor || rem == divisor) {
0056         temp_divisor = mpi_copy(divisor);
0057         divisor = temp_divisor;
0058     }
0059
0060     mpi_tdiv_qr(quot, rem, dividend, divisor);
0061
0062     if ((divisor_sign ^ dividend->sign) && rem->nlimbs) {
0063         mpi_sub_ui(quot, quot, 1);
0064         mpi_add(rem, rem, divisor);
0065     }
0066
0067     if (temp_divisor)
0068         mpi_free(temp_divisor);
0069 }
0070
0071 /* If den == quot, den needs temporary storage.
0072  * If den == rem, den needs temporary storage.
0073  * If num == quot, num needs temporary storage.
0074  * If den has temporary storage, it can be normalized while being copied,
0075  *   i.e no extra storage should be allocated.
0076  */
0077
0078 void mpi_tdiv_r(MPI rem, MPI num, MPI den)
0079 {
0080     mpi_tdiv_qr(NULL, rem, num, den);
0081 }
0082
0083 void mpi_tdiv_qr(MPI quot, MPI rem, MPI num, MPI den)
0084 {
0085     mpi_ptr_t np, dp;
0086     mpi_ptr_t qp, rp;
0087     mpi_size_t nsize = num->nlimbs;
0088     mpi_size_t dsize = den->nlimbs;
0089     mpi_size_t qsize, rsize;
0090     mpi_size_t sign_remainder = num->sign;
0091     mpi_size_t sign_quotient = num->sign ^ den->sign;
0092     unsigned int normalization_steps;
0093     mpi_limb_t q_limb;
0094     mpi_ptr_t marker[5];
0095     int markidx = 0;
0096
0097     /* Ensure space is enough for quotient and remainder.
0098      * We need space for an extra limb in the remainder, because it's
0099      * up-shifted (normalized) below.
0100      */
0101     rsize = nsize + 1;
0102     mpi_resize(rem, rsize);
0103
0104     qsize = rsize - dsize;    /* qsize cannot be bigger than this.  */
0105     if (qsize <= 0) {
0106         if (num != rem) {
0107             rem->nlimbs = num->nlimbs;
0108             rem->sign = num->sign;
0109             MPN_COPY(rem->d, num->d, nsize);
0110         }
0111         if (quot) {
0112             /* This needs to follow the assignment to rem, in case the
0113              * numerator and quotient are the same.
0114              */
0115             quot->nlimbs = 0;
0116             quot->sign = 0;
0117         }
0118         return;
0119     }
0120
0121     if (quot)
0122         mpi_resize(quot, qsize);
0123
0124     /* Read pointers here, when reallocation is finished.  */
0125     np = num->d;
0126     dp = den->d;
0127     rp = rem->d;
0128
0129     /* Optimize division by a single-limb divisor.  */
0130     if (dsize == 1) {
0131         mpi_limb_t rlimb;
0132         if (quot) {
0133             qp = quot->d;
0134             rlimb = mpihelp_divmod_1(qp, np, nsize, dp[0]);
0135             qsize -= qp[qsize - 1] == 0;
0136             quot->nlimbs = qsize;
0137             quot->sign = sign_quotient;
0138         } else
0139             rlimb = mpihelp_mod_1(np, nsize, dp[0]);
0140         rp[0] = rlimb;
0141         rsize = rlimb != 0?1:0;
0142         rem->nlimbs = rsize;
0143         rem->sign = sign_remainder;
0144         return;
0145     }
0146
0147
0148     if (quot) {
0149         qp = quot->d;
0150         /* Make sure QP and NP point to different objects.  Otherwise the
0151          * numerator would be gradually overwritten by the quotient limbs.
0152          */
0153         if (qp == np) { /* Copy NP object to temporary space.  */
0154             np = marker[markidx++] = mpi_alloc_limb_space(nsize);
0155             MPN_COPY(np, qp, nsize);
0156         }
0157     } else /* Put quotient at top of remainder. */
0158         qp = rp + dsize;
0159
0160     normalization_steps = count_leading_zeros(dp[dsize - 1]);
0161
0162     /* Normalize the denominator, i.e. make its most significant bit set by
0163      * shifting it NORMALIZATION_STEPS bits to the left.  Also shift the
0164      * numerator the same number of steps (to keep the quotient the same!).
0165      */
0166     if (normalization_steps) {
0167         mpi_ptr_t tp;
0168         mpi_limb_t nlimb;
0169
0170         /* Shift up the denominator setting the most significant bit of
0171          * the most significant word.  Use temporary storage not to clobber
0172          * the original contents of the denominator.
0173          */
0174         tp = marker[markidx++] = mpi_alloc_limb_space(dsize);
0175         mpihelp_lshift(tp, dp, dsize, normalization_steps);
0176         dp = tp;
0177
0178         /* Shift up the numerator, possibly introducing a new most
0179          * significant word.  Move the shifted numerator in the remainder
0180          * meanwhile.
0181          */
0182         nlimb = mpihelp_lshift(rp, np, nsize, normalization_steps);
0183         if (nlimb) {
0184             rp[nsize] = nlimb;
0185             rsize = nsize + 1;
0186         } else
0187             rsize = nsize;
0188     } else {
0189         /* The denominator is already normalized, as required.  Copy it to
0190          * temporary space if it overlaps with the quotient or remainder.
0191          */
0192         if (dp == rp || (quot && (dp == qp))) {
0193             mpi_ptr_t tp;
0194
0195             tp = marker[markidx++] = mpi_alloc_limb_space(dsize);
0196             MPN_COPY(tp, dp, dsize);
0197             dp = tp;
0198         }
0199
0200         /* Move the numerator to the remainder.  */
0201         if (rp != np)
0202             MPN_COPY(rp, np, nsize);
0203
0204         rsize = nsize;
0205     }
0206
0207     q_limb = mpihelp_divrem(qp, 0, rp, rsize, dp, dsize);
0208
0209     if (quot) {
0210         qsize = rsize - dsize;
0211         if (q_limb) {
0212             qp[qsize] = q_limb;
0213             qsize += 1;
0214         }
0215
0216         quot->nlimbs = qsize;
0217         quot->sign = sign_quotient;
0218     }
0219
0220     rsize = dsize;
0221     MPN_NORMALIZE(rp, rsize);
0222
0223     if (normalization_steps && rsize) {
0224         mpihelp_rshift(rp, rp, rsize, normalization_steps);
0225         rsize -= rp[rsize - 1] == 0?1:0;
0226     }
0227
0228     rem->nlimbs = rsize;
0229     rem->sign   = sign_remainder;
0230     while (markidx) {
0231         markidx--;
0232         mpi_free_limb_space(marker[markidx]);
0233     }
0234 }