raw
ch1_mpi                 1 /* mpi-div.c  -  MPI functions
ch1_mpi 2 * Modified by No Such Labs. (C) 2015. See README.
ch1_mpi 3 *
ch1_mpi 4 * This file was originally part of Gnu Privacy Guard (GPG), ver. 1.4.10,
ch1_mpi 5 * SHA256(gnupg-1.4.10.tar.gz):
ch1_mpi 6 * 0bfd74660a2f6cedcf7d8256db4a63c996ffebbcdc2cf54397bfb72878c5a85a
ch1_mpi 7 * (C) 1994-2005 Free Software Foundation, Inc.
ch1_mpi 8 *
ch1_mpi 9 * This program is free software: you can redistribute it and/or modify
ch1_mpi 10 * it under the terms of the GNU General Public License as published by
ch1_mpi 11 * the Free Software Foundation, either version 3 of the License, or
ch1_mpi 12 * (at your option) any later version.
ch1_mpi 13 *
ch1_mpi 14 * This program is distributed in the hope that it will be useful,
ch1_mpi 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
ch1_mpi 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
ch1_mpi 17 * GNU General Public License for more details.
ch1_mpi 18 *
ch1_mpi 19 * You should have received a copy of the GNU General Public License
ch1_mpi 20 * along with this program. If not, see <http://www.gnu.org/licenses/>.
ch1_mpi 21 */
ch1_mpi 22
ch1_mpi 23 #include <stdio.h>
ch1_mpi 24 #include <stdlib.h>
ch1_mpi 25
ch1_mpi 26 #include "knobs.h"
ch1_mpi 27 #include "mpi-internal.h"
ch1_mpi 28 #include "longlong.h"
ch1_mpi 29
ch1_mpi 30
ch1_mpi 31
ch1_mpi 32 void
ch1_mpi 33 mpi_fdiv_r( MPI rem, MPI dividend, MPI divisor )
ch1_mpi 34 {
ch1_mpi 35 int divisor_sign = divisor->sign;
ch1_mpi 36 MPI temp_divisor = NULL;
ch1_mpi 37
ch1_mpi 38 /* We need the original value of the divisor after the remainder has been
ch1_mpi 39 * preliminary calculated. We have to copy it to temporary space if it's
ch1_mpi 40 * the same variable as REM. */
ch1_mpi 41 if( rem == divisor ) {
ch1_mpi 42 temp_divisor = mpi_copy( divisor );
ch1_mpi 43 divisor = temp_divisor;
ch1_mpi 44 }
ch1_mpi 45
ch1_mpi 46 mpi_tdiv_r( rem, dividend, divisor );
ch1_mpi 47
ch1_mpi 48 if( ((divisor_sign?1:0) ^ (dividend->sign?1:0)) && rem->nlimbs )
ch1_mpi 49 mpi_add( rem, rem, divisor);
ch1_mpi 50
ch1_mpi 51 if( temp_divisor )
ch1_mpi 52 mpi_free(temp_divisor);
ch1_mpi 53 }
ch1_mpi 54
ch1_mpi 55
ch1_mpi 56
ch1_mpi 57 /****************
ch1_mpi 58 * Division rounding the quotient towards -infinity.
ch1_mpi 59 * The remainder gets the same sign as the denominator.
ch1_mpi 60 * rem is optional
ch1_mpi 61 */
ch1_mpi 62
ch1_mpi 63 ulong
ch1_mpi 64 mpi_fdiv_r_ui( MPI rem, MPI dividend, ulong divisor )
ch1_mpi 65 {
ch1_mpi 66 mpi_limb_t rlimb;
ch1_mpi 67
ch1_mpi 68 rlimb = mpihelp_mod_1( dividend->d, dividend->nlimbs, divisor );
ch1_mpi 69 if( rlimb && dividend->sign )
ch1_mpi 70 rlimb = divisor - rlimb;
ch1_mpi 71
ch1_mpi 72 if( rem ) {
ch1_mpi 73 rem->d[0] = rlimb;
ch1_mpi 74 rem->nlimbs = rlimb? 1:0;
ch1_mpi 75 }
ch1_mpi 76 return rlimb;
ch1_mpi 77 }
ch1_mpi 78
ch1_mpi 79
ch1_mpi 80 void
ch1_mpi 81 mpi_fdiv_q( MPI quot, MPI dividend, MPI divisor )
ch1_mpi 82 {
ch1_mpi 83 MPI tmp = mpi_alloc( mpi_get_nlimbs(quot) );
ch1_mpi 84 mpi_fdiv_qr( quot, tmp, dividend, divisor);
ch1_mpi 85 mpi_free(tmp);
ch1_mpi 86 }
ch1_mpi 87
ch1_mpi 88 void
ch1_mpi 89 mpi_fdiv_qr( MPI quot, MPI rem, MPI dividend, MPI divisor )
ch1_mpi 90 {
ch1_mpi 91 int divisor_sign = divisor->sign;
ch1_mpi 92 MPI temp_divisor = NULL;
ch1_mpi 93
ch1_mpi 94 if( quot == divisor || rem == divisor ) {
ch1_mpi 95 temp_divisor = mpi_copy( divisor );
ch1_mpi 96 divisor = temp_divisor;
ch1_mpi 97 }
ch1_mpi 98
ch1_mpi 99 mpi_tdiv_qr( quot, rem, dividend, divisor );
ch1_mpi 100
ch1_mpi 101 if( (divisor_sign ^ dividend->sign) && rem->nlimbs ) {
ch1_mpi 102 mpi_sub_ui( quot, quot, 1 );
ch1_mpi 103 mpi_add( rem, rem, divisor);
ch1_mpi 104 }
ch1_mpi 105
ch1_mpi 106 if( temp_divisor )
ch1_mpi 107 mpi_free(temp_divisor);
ch1_mpi 108 }
ch1_mpi 109
ch1_mpi 110
ch1_mpi 111 /* If den == quot, den needs temporary storage.
ch1_mpi 112 * If den == rem, den needs temporary storage.
ch1_mpi 113 * If num == quot, num needs temporary storage.
ch1_mpi 114 * If den has temporary storage, it can be normalized while being copied,
ch1_mpi 115 * i.e no extra storage should be allocated.
ch1_mpi 116 */
ch1_mpi 117
ch1_mpi 118 void
ch1_mpi 119 mpi_tdiv_r( MPI rem, MPI num, MPI den)
ch1_mpi 120 {
ch1_mpi 121 mpi_tdiv_qr(NULL, rem, num, den );
ch1_mpi 122 }
ch1_mpi 123
ch1_mpi 124 void
ch1_mpi 125 mpi_tdiv_qr( MPI quot, MPI rem, MPI num, MPI den)
ch1_mpi 126 {
ch1_mpi 127 mpi_ptr_t np, dp;
ch1_mpi 128 mpi_ptr_t qp, rp;
ch1_mpi 129 mpi_size_t nsize = num->nlimbs;
ch1_mpi 130 mpi_size_t dsize = den->nlimbs;
ch1_mpi 131 mpi_size_t qsize, rsize;
ch1_mpi 132 mpi_size_t sign_remainder = num->sign;
ch1_mpi 133 mpi_size_t sign_quotient = num->sign ^ den->sign;
ch1_mpi 134 unsigned normalization_steps;
ch1_mpi 135 mpi_limb_t q_limb;
ch1_mpi 136 mpi_ptr_t marker[5];
ch1_mpi 137 int markidx=0;
ch1_mpi 138
ch1_mpi 139 /* Ensure space is enough for quotient and remainder.
ch1_mpi 140 * We need space for an extra limb in the remainder, because it's
ch1_mpi 141 * up-shifted (normalized) below. */
ch1_mpi 142 rsize = nsize + 1;
ch1_mpi 143 mpi_resize( rem, rsize);
ch1_mpi 144
ch1_mpi 145 qsize = rsize - dsize; /* qsize cannot be bigger than this. */
ch1_mpi 146 if( qsize <= 0 ) {
ch1_mpi 147 if( num != rem ) {
ch1_mpi 148 rem->nlimbs = num->nlimbs;
ch1_mpi 149 rem->sign = num->sign;
ch1_mpi 150 MPN_COPY(rem->d, num->d, nsize);
ch1_mpi 151 }
ch1_mpi 152 if( quot ) {
ch1_mpi 153 /* This needs to follow the assignment to rem, in case the
ch1_mpi 154 * numerator and quotient are the same. */
ch1_mpi 155 quot->nlimbs = 0;
ch1_mpi 156 quot->sign = 0;
ch1_mpi 157 }
ch1_mpi 158 return;
ch1_mpi 159 }
ch1_mpi 160
ch1_mpi 161 if( quot )
ch1_mpi 162 mpi_resize( quot, qsize);
ch1_mpi 163
ch1_mpi 164 /* Read pointers here, when reallocation is finished. */
ch1_mpi 165 np = num->d;
ch1_mpi 166 dp = den->d;
ch1_mpi 167 rp = rem->d;
ch1_mpi 168
ch1_mpi 169 /* Optimize division by a single-limb divisor. */
ch1_mpi 170 if( dsize == 1 ) {
ch1_mpi 171 mpi_limb_t rlimb;
ch1_mpi 172 if( quot ) {
ch1_mpi 173 qp = quot->d;
ch1_mpi 174 rlimb = mpihelp_divmod_1( qp, np, nsize, dp[0] );
ch1_mpi 175 qsize -= qp[qsize - 1] == 0;
ch1_mpi 176 quot->nlimbs = qsize;
ch1_mpi 177 quot->sign = sign_quotient;
ch1_mpi 178 }
ch1_mpi 179 else
ch1_mpi 180 rlimb = mpihelp_mod_1( np, nsize, dp[0] );
ch1_mpi 181 rp[0] = rlimb;
ch1_mpi 182 rsize = rlimb != 0?1:0;
ch1_mpi 183 rem->nlimbs = rsize;
ch1_mpi 184 rem->sign = sign_remainder;
ch1_mpi 185 return;
ch1_mpi 186 }
ch1_mpi 187
ch1_mpi 188
ch1_mpi 189 if( quot ) {
ch1_mpi 190 qp = quot->d;
ch1_mpi 191 /* Make sure QP and NP point to different objects. Otherwise the
ch1_mpi 192 * numerator would be gradually overwritten by the quotient limbs. */
ch1_mpi 193 if(qp == np) { /* Copy NP object to temporary space. */
ch1_mpi 194 np = marker[markidx++] = mpi_alloc_limb_space(nsize,
ch1_mpi 195 mpi_is_secure(quot));
ch1_mpi 196 MPN_COPY(np, qp, nsize);
ch1_mpi 197 }
ch1_mpi 198 }
ch1_mpi 199 else /* Put quotient at top of remainder. */
ch1_mpi 200 qp = rp + dsize;
ch1_mpi 201
ch1_mpi 202 count_leading_zeros( normalization_steps, dp[dsize - 1] );
ch1_mpi 203
ch1_mpi 204 /* Normalize the denominator, i.e. make its most significant bit set by
ch1_mpi 205 * shifting it NORMALIZATION_STEPS bits to the left. Also shift the
ch1_mpi 206 * numerator the same number of steps (to keep the quotient the same!).
ch1_mpi 207 */
ch1_mpi 208 if( normalization_steps ) {
ch1_mpi 209 mpi_ptr_t tp;
ch1_mpi 210 mpi_limb_t nlimb;
ch1_mpi 211
ch1_mpi 212 /* Shift up the denominator setting the most significant bit of
ch1_mpi 213 * the most significant word. Use temporary storage not to clobber
ch1_mpi 214 * the original contents of the denominator. */
ch1_mpi 215 tp = marker[markidx++] = mpi_alloc_limb_space(dsize,mpi_is_secure(den));
ch1_mpi 216 mpihelp_lshift( tp, dp, dsize, normalization_steps );
ch1_mpi 217 dp = tp;
ch1_mpi 218
ch1_mpi 219 /* Shift up the numerator, possibly introducing a new most
ch1_mpi 220 * significant word. Move the shifted numerator in the remainder
ch1_mpi 221 * meanwhile. */
ch1_mpi 222 nlimb = mpihelp_lshift(rp, np, nsize, normalization_steps);
ch1_mpi 223 if( nlimb ) {
ch1_mpi 224 rp[nsize] = nlimb;
ch1_mpi 225 rsize = nsize + 1;
ch1_mpi 226 }
ch1_mpi 227 else
ch1_mpi 228 rsize = nsize;
ch1_mpi 229 }
ch1_mpi 230 else {
ch1_mpi 231 /* The denominator is already normalized, as required. Copy it to
ch1_mpi 232 * temporary space if it overlaps with the quotient or remainder. */
ch1_mpi 233 if( dp == rp || (quot && (dp == qp))) {
ch1_mpi 234 mpi_ptr_t tp;
ch1_mpi 235
ch1_mpi 236 tp = marker[markidx++] = mpi_alloc_limb_space(dsize, mpi_is_secure(den));
ch1_mpi 237 MPN_COPY( tp, dp, dsize );
ch1_mpi 238 dp = tp;
ch1_mpi 239 }
ch1_mpi 240
ch1_mpi 241 /* Move the numerator to the remainder. */
ch1_mpi 242 if( rp != np )
ch1_mpi 243 MPN_COPY(rp, np, nsize);
ch1_mpi 244
ch1_mpi 245 rsize = nsize;
ch1_mpi 246 }
ch1_mpi 247
ch1_mpi 248 q_limb = mpihelp_divrem( qp, 0, rp, rsize, dp, dsize );
ch1_mpi 249
ch1_mpi 250 if( quot ) {
ch1_mpi 251 qsize = rsize - dsize;
ch1_mpi 252 if(q_limb) {
ch1_mpi 253 qp[qsize] = q_limb;
ch1_mpi 254 qsize += 1;
ch1_mpi 255 }
ch1_mpi 256
ch1_mpi 257 quot->nlimbs = qsize;
ch1_mpi 258 quot->sign = sign_quotient;
ch1_mpi 259 }
ch1_mpi 260
ch1_mpi 261 rsize = dsize;
ch1_mpi 262 MPN_NORMALIZE (rp, rsize);
ch1_mpi 263
ch1_mpi 264 if( normalization_steps && rsize ) {
ch1_mpi 265 mpihelp_rshift(rp, rp, rsize, normalization_steps);
ch1_mpi 266 rsize -= rp[rsize - 1] == 0?1:0;
ch1_mpi 267 }
ch1_mpi 268
ch1_mpi 269 rem->nlimbs = rsize;
ch1_mpi 270 rem->sign = sign_remainder;
ch1_mpi 271 while( markidx )
ch1_mpi 272 mpi_free_limb_space(marker[--markidx]);
ch1_mpi 273 }
ch1_mpi 274
ch1_mpi 275 void
ch1_mpi 276 mpi_tdiv_q_2exp( MPI w, MPI u, unsigned count )
ch1_mpi 277 {
ch1_mpi 278 mpi_size_t usize, wsize;
ch1_mpi 279 mpi_size_t limb_cnt;
ch1_mpi 280
ch1_mpi 281 usize = u->nlimbs;
ch1_mpi 282 limb_cnt = count / BITS_PER_MPI_LIMB;
ch1_mpi 283 wsize = usize - limb_cnt;
ch1_mpi 284 if( limb_cnt >= usize )
ch1_mpi 285 w->nlimbs = 0;
ch1_mpi 286 else {
ch1_mpi 287 mpi_ptr_t wp;
ch1_mpi 288 mpi_ptr_t up;
ch1_mpi 289
ch1_mpi 290 RESIZE_IF_NEEDED( w, wsize );
ch1_mpi 291 wp = w->d;
ch1_mpi 292 up = u->d;
ch1_mpi 293
ch1_mpi 294 count %= BITS_PER_MPI_LIMB;
ch1_mpi 295 if( count ) {
ch1_mpi 296 mpihelp_rshift( wp, up + limb_cnt, wsize, count );
ch1_mpi 297 wsize -= !wp[wsize - 1];
ch1_mpi 298 }
ch1_mpi 299 else {
ch1_mpi 300 MPN_COPY_INCR( wp, up + limb_cnt, wsize);
ch1_mpi 301 }
ch1_mpi 302
ch1_mpi 303 w->nlimbs = wsize;
ch1_mpi 304 }
ch1_mpi 305 }
ch1_mpi 306
ch1_mpi 307 /****************
ch1_mpi 308 * Check whether dividend is divisible by divisor
ch1_mpi 309 * (note: divisor must fit into a limb)
ch1_mpi 310 */
ch1_mpi 311 int
ch1_mpi 312 mpi_divisible_ui(MPI dividend, ulong divisor )
ch1_mpi 313 {
ch1_mpi 314 return !mpihelp_mod_1( dividend->d, dividend->nlimbs, divisor );
ch1_mpi 315 }
ch1_mpi 316