-
+ FB984A326A9BCDA1A4CBC05EE279E55CFEFCF157393D2F66405760B256395C3A73F1F41EBFC335722022EA04C79F6E02AB3179ECC9A66E037DD7A106572B4924
mpi/udiv-w-sdiv.c
(0 . 0)(1 . 132)
11213 /* mpihelp_udiv_w_sdiv -- implement udiv_qrnnd on machines with only signed
11214 * division.
11215 * Copyright (C) 1992, 1994, 1996, 1998 Free Software Foundation, Inc.
11216 * Contributed by Peter L. Montgomery.
11217 *
11218 * This file is part of GnuPG.
11219 *
11220 * GnuPG is free software; you can redistribute it and/or modify
11221 * it under the terms of the GNU General Public License as published by
11222 * the Free Software Foundation; either version 3 of the License, or
11223 * (at your option) any later version.
11224 *
11225 * GnuPG is distributed in the hope that it will be useful,
11226 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11227 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11228 * GNU General Public License for more details.
11229 *
11230 * You should have received a copy of the GNU General Public License
11231 * along with this program; if not, see <http://www.gnu.org/licenses/>.
11232 */
11233
11234 #include <config.h>
11235 #include <stdio.h>
11236 #include <stdlib.h>
11237 #include "mpi-internal.h"
11238 #include "longlong.h"
11239
11240
11241 #if 0 /* not yet ported to MPI */
11242
11243 mpi_limb_t
11244 mpihelp_udiv_w_sdiv( mpi_limp_t *rp,
11245 mpi_limp_t *a1,
11246 mpi_limp_t *a0,
11247 mpi_limp_t *d )
11248 {
11249 mp_limb_t q, r;
11250 mp_limb_t c0, c1, b1;
11251
11252 if ((mpi_limb_signed_t) d >= 0)
11253 {
11254 if (a1 < d - a1 - (a0 >> (BITS_PER_MP_LIMB - 1)))
11255 {
11256 /* dividend, divisor, and quotient are nonnegative */
11257 sdiv_qrnnd (q, r, a1, a0, d);
11258 }
11259 else
11260 {
11261 /* Compute c1*2^32 + c0 = a1*2^32 + a0 - 2^31*d */
11262 sub_ddmmss (c1, c0, a1, a0, d >> 1, d << (BITS_PER_MP_LIMB - 1));
11263 /* Divide (c1*2^32 + c0) by d */
11264 sdiv_qrnnd (q, r, c1, c0, d);
11265 /* Add 2^31 to quotient */
11266 q += (mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1);
11267 }
11268 }
11269 else
11270 {
11271 b1 = d >> 1; /* d/2, between 2^30 and 2^31 - 1 */
11272 c1 = a1 >> 1; /* A/2 */
11273 c0 = (a1 << (BITS_PER_MP_LIMB - 1)) + (a0 >> 1);
11274
11275 if (a1 < b1) /* A < 2^32*b1, so A/2 < 2^31*b1 */
11276 {
11277 sdiv_qrnnd (q, r, c1, c0, b1); /* (A/2) / (d/2) */
11278
11279 r = 2*r + (a0 & 1); /* Remainder from A/(2*b1) */
11280 if ((d & 1) != 0)
11281 {
11282 if (r >= q)
11283 r = r - q;
11284 else if (q - r <= d)
11285 {
11286 r = r - q + d;
11287 q--;
11288 }
11289 else
11290 {
11291 r = r - q + 2*d;
11292 q -= 2;
11293 }
11294 }
11295 }
11296 else if (c1 < b1) /* So 2^31 <= (A/2)/b1 < 2^32 */
11297 {
11298 c1 = (b1 - 1) - c1;
11299 c0 = ~c0; /* logical NOT */
11300
11301 sdiv_qrnnd (q, r, c1, c0, b1); /* (A/2) / (d/2) */
11302
11303 q = ~q; /* (A/2)/b1 */
11304 r = (b1 - 1) - r;
11305
11306 r = 2*r + (a0 & 1); /* A/(2*b1) */
11307
11308 if ((d & 1) != 0)
11309 {
11310 if (r >= q)
11311 r = r - q;
11312 else if (q - r <= d)
11313 {
11314 r = r - q + d;
11315 q--;
11316 }
11317 else
11318 {
11319 r = r - q + 2*d;
11320 q -= 2;
11321 }
11322 }
11323 }
11324 else /* Implies c1 = b1 */
11325 { /* Hence a1 = d - 1 = 2*b1 - 1 */
11326 if (a0 >= -d)
11327 {
11328 q = -1;
11329 r = a0 + d;
11330 }
11331 else
11332 {
11333 q = -2;
11334 r = a0 + 2*d;
11335 }
11336 }
11337 }
11338
11339 *rp = r;
11340 return q;
11341 }
11342
11343 #endif
11344