1/* mpn_udiv_w_sdiv -- implement udiv_qrnnd on machines with only signed 2 division. 3 4 Contributed by Peter L. Montgomery. 5 6 THIS IS AN INTERNAL FUNCTION WITH A MUTABLE INTERFACE. IT IS ONLY SAFE 7 TO REACH THIS FUNCTION THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS 8 ALMOST GUARANTEED THAT THIS FUNCTION WILL CHANGE OR DISAPPEAR IN A FUTURE 9 GNU MP RELEASE. 10 11 12Copyright 1992, 1994, 1996, 2000, 2011, 2012 Free Software Foundation, Inc. 13 14This file is part of the GNU MP Library. 15 16The GNU MP Library is free software; you can redistribute it and/or modify 17it under the terms of the GNU Lesser General Public License as published by 18the Free Software Foundation; either version 3 of the License, or (at your 19option) any later version. 20 21The GNU MP Library is distributed in the hope that it will be useful, but 22WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 23or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 24License for more details. 25 26You should have received a copy of the GNU Lesser General Public License 27along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 28 29#include "gmp.h" 30#include "gmp-impl.h" 31#include "longlong.h" 32 33mp_limb_t 34mpn_udiv_w_sdiv (mp_limb_t *rp, mp_limb_t a1, mp_limb_t a0, mp_limb_t d) 35{ 36 mp_limb_t q, r; 37 mp_limb_t c0, c1, b1; 38 39 ASSERT (d != 0); 40 ASSERT (a1 < d); 41 42 if ((mp_limb_signed_t) d >= 0) 43 { 44 if (a1 < d - a1 - (a0 >> (GMP_LIMB_BITS - 1))) 45 { 46 /* dividend, divisor, and quotient are nonnegative */ 47 sdiv_qrnnd (q, r, a1, a0, d); 48 } 49 else 50 { 51 /* Compute c1*2^32 + c0 = a1*2^32 + a0 - 2^31*d */ 52 sub_ddmmss (c1, c0, a1, a0, d >> 1, d << (GMP_LIMB_BITS - 1)); 53 /* Divide (c1*2^32 + c0) by d */ 54 sdiv_qrnnd (q, r, c1, c0, d); 55 /* Add 2^31 to quotient */ 56 q += (mp_limb_t) 1 << (GMP_LIMB_BITS - 1); 57 } 58 } 59 else 60 { 61 b1 = d >> 1; /* d/2, between 2^30 and 2^31 - 1 */ 62 c1 = a1 >> 1; /* A/2 */ 63 c0 = (a1 << (GMP_LIMB_BITS - 1)) + (a0 >> 1); 64 65 if (a1 < b1) /* A < 2^32*b1, so A/2 < 2^31*b1 */ 66 { 67 sdiv_qrnnd (q, r, c1, c0, b1); /* (A/2) / (d/2) */ 68 69 r = 2*r + (a0 & 1); /* Remainder from A/(2*b1) */ 70 if ((d & 1) != 0) 71 { 72 if (r >= q) 73 r = r - q; 74 else if (q - r <= d) 75 { 76 r = r - q + d; 77 q--; 78 } 79 else 80 { 81 r = r - q + 2*d; 82 q -= 2; 83 } 84 } 85 } 86 else if (c1 < b1) /* So 2^31 <= (A/2)/b1 < 2^32 */ 87 { 88 c1 = (b1 - 1) - c1; 89 c0 = ~c0; /* logical NOT */ 90 91 sdiv_qrnnd (q, r, c1, c0, b1); /* (A/2) / (d/2) */ 92 93 q = ~q; /* (A/2)/b1 */ 94 r = (b1 - 1) - r; 95 96 r = 2*r + (a0 & 1); /* A/(2*b1) */ 97 98 if ((d & 1) != 0) 99 { 100 if (r >= q) 101 r = r - q; 102 else if (q - r <= d) 103 { 104 r = r - q + d; 105 q--; 106 } 107 else 108 { 109 r = r - q + 2*d; 110 q -= 2; 111 } 112 } 113 } 114 else /* Implies c1 = b1 */ 115 { /* Hence a1 = d - 1 = 2*b1 - 1 */ 116 if (a0 >= -d) 117 { 118 q = -CNST_LIMB(1); 119 r = a0 + d; 120 } 121 else 122 { 123 q = -CNST_LIMB(2); 124 r = a0 + 2*d; 125 } 126 } 127 } 128 129 *rp = r; 130 return q; 131} 132