1/* mpihelp-div.c - MPI helper functions 2 * Copyright (C) 1994, 1996 Free Software Foundation, Inc. 3 * Copyright (C) 1998, 1999 Free Software Foundation, Inc. 4 * 5 * This file is part of GnuPG. 6 * 7 * GnuPG is free software; you can redistribute it and/or modify 8 * it under the terms of the GNU General Public License as published by 9 * the Free Software Foundation; either version 2 of the License, or 10 * (at your option) any later version. 11 * 12 * GnuPG is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU General Public License for more details. 16 * 17 * You should have received a copy of the GNU General Public License 18 * along with this program; if not, write to the Free Software 19 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA 20 * 21 * Note: This code is heavily based on the GNU MP Library. 22 * Actually it's the same code with only minor changes in the 23 * way the data is stored; this is to support the abstraction 24 * of an optional secure memory allocation which may be used 25 * to avoid revealing of sensitive data due to paging etc. 26 * The GNU MP Library itself is published under the LGPL; 27 * however I decided to publish this code under the plain GPL. 28 */ 29 30#include "mpi-internal.h" 31#include "longlong.h" 32 33#ifndef UMUL_TIME 34#define UMUL_TIME 1 35#endif 36#ifndef UDIV_TIME 37#define UDIV_TIME UMUL_TIME 38#endif 39 40/* Divide num (NP/NSIZE) by den (DP/DSIZE) and write 41 * the NSIZE-DSIZE least significant quotient limbs at QP 42 * and the DSIZE long remainder at NP. If QEXTRA_LIMBS is 43 * non-zero, generate that many fraction bits and append them after the 44 * other quotient limbs. 45 * Return the most significant limb of the quotient, this is always 0 or 1. 46 * 47 * Preconditions: 48 * 0. NSIZE >= DSIZE. 49 * 1. The most significant bit of the divisor must be set. 50 * 2. QP must either not overlap with the input operands at all, or 51 * QP + DSIZE >= NP must hold true. (This means that it's 52 * possible to put the quotient in the high part of NUM, right after the 53 * remainder in NUM. 54 * 3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero. 55 */ 56 57mpi_limb_t 58mpihelp_divrem(mpi_ptr_t qp, mpi_size_t qextra_limbs, 59 mpi_ptr_t np, mpi_size_t nsize, mpi_ptr_t dp, mpi_size_t dsize) 60{ 61 mpi_limb_t most_significant_q_limb = 0; 62 63 switch (dsize) { 64 case 0: 65 /* We are asked to divide by zero, so go ahead and do it! (To make 66 the compiler not remove this statement, return the value.) */ 67 /* 68 * existing clients of this function have been modified 69 * not to call it with dsize == 0, so this should not happen 70 */ 71 return 1 / dsize; 72 73 case 1: 74 { 75 mpi_size_t i; 76 mpi_limb_t n1; 77 mpi_limb_t d; 78 79 d = dp[0]; 80 n1 = np[nsize - 1]; 81 82 if (n1 >= d) { 83 n1 -= d; 84 most_significant_q_limb = 1; 85 } 86 87 qp += qextra_limbs; 88 for (i = nsize - 2; i >= 0; i--) 89 udiv_qrnnd(qp[i], n1, n1, np[i], d); 90 qp -= qextra_limbs; 91 92 for (i = qextra_limbs - 1; i >= 0; i--) 93 udiv_qrnnd(qp[i], n1, n1, 0, d); 94 95 np[0] = n1; 96 } 97 break; 98 99 case 2: 100 { 101 mpi_size_t i; 102 mpi_limb_t n1, n0, n2; 103 mpi_limb_t d1, d0; 104 105 np += nsize - 2; 106 d1 = dp[1]; 107 d0 = dp[0]; 108 n1 = np[1]; 109 n0 = np[0]; 110 111 if (n1 >= d1 && (n1 > d1 || n0 >= d0)) { 112 sub_ddmmss(n1, n0, n1, n0, d1, d0); 113 most_significant_q_limb = 1; 114 } 115 116 for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--) { 117 mpi_limb_t q; 118 mpi_limb_t r; 119 120 if (i >= qextra_limbs) 121 np--; 122 else 123 np[0] = 0; 124 125 if (n1 == d1) { 126 /* Q should be either 111..111 or 111..110. Need special 127 * treatment of this rare case as normal division would 128 * give overflow. */ 129 q = ~(mpi_limb_t) 0; 130 131 r = n0 + d1; 132 if (r < d1) { /* Carry in the addition? */ 133 add_ssaaaa(n1, n0, r - d0, 134 np[0], 0, d0); 135 qp[i] = q; 136 continue; 137 } 138 n1 = d0 - (d0 != 0 ? 1 : 0); 139 n0 = -d0; 140 } else { 141 udiv_qrnnd(q, r, n1, n0, d1); 142 umul_ppmm(n1, n0, d0, q); 143 } 144 145 n2 = np[0]; 146q_test: 147 if (n1 > r || (n1 == r && n0 > n2)) { 148 /* The estimated Q was too large. */ 149 q--; 150 sub_ddmmss(n1, n0, n1, n0, 0, d0); 151 r += d1; 152 if (r >= d1) /* If not carry, test Q again. */ 153 goto q_test; 154 } 155 156 qp[i] = q; 157 sub_ddmmss(n1, n0, r, n2, n1, n0); 158 } 159 np[1] = n1; 160 np[0] = n0; 161 } 162 break; 163 164 default: 165 { 166 mpi_size_t i; 167 mpi_limb_t dX, d1, n0; 168 169 np += nsize - dsize; 170 dX = dp[dsize - 1]; 171 d1 = dp[dsize - 2]; 172 n0 = np[dsize - 1]; 173 174 if (n0 >= dX) { 175 if (n0 > dX 176 || mpihelp_cmp(np, dp, dsize - 1) >= 0) { 177 mpihelp_sub_n(np, np, dp, dsize); 178 n0 = np[dsize - 1]; 179 most_significant_q_limb = 1; 180 } 181 } 182 183 for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--) { 184 mpi_limb_t q; 185 mpi_limb_t n1, n2; 186 mpi_limb_t cy_limb; 187 188 if (i >= qextra_limbs) { 189 np--; 190 n2 = np[dsize]; 191 } else { 192 n2 = np[dsize - 1]; 193 MPN_COPY_DECR(np + 1, np, dsize - 1); 194 np[0] = 0; 195 } 196 197 if (n0 == dX) { 198 /* This might over-estimate q, but it's probably not worth 199 * the extra code here to find out. */ 200 q = ~(mpi_limb_t) 0; 201 } else { 202 mpi_limb_t r; 203 204 udiv_qrnnd(q, r, n0, np[dsize - 1], dX); 205 umul_ppmm(n1, n0, d1, q); 206 207 while (n1 > r 208 || (n1 == r 209 && n0 > np[dsize - 2])) { 210 q--; 211 r += dX; 212 if (r < dX) /* I.e. "carry in previous addition?" */ 213 break; 214 n1 -= n0 < d1; 215 n0 -= d1; 216 } 217 } 218 219 /* Possible optimization: We already have (q * n0) and (1 * n1) 220 * after the calculation of q. Taking advantage of that, we 221 * could make this loop make two iterations less. */ 222 cy_limb = mpihelp_submul_1(np, dp, dsize, q); 223 224 if (n2 != cy_limb) { 225 mpihelp_add_n(np, np, dp, dsize); 226 q--; 227 } 228 229 qp[i] = q; 230 n0 = np[dsize - 1]; 231 } 232 } 233 } 234 235 return most_significant_q_limb; 236} 237