1356290Sjkim/* mpn_udiv_w_sdiv -- implement udiv_qrnnd on machines with only signed 2238405Sjkim division. 3238405Sjkim 4238405Sjkim Contributed by Peter L. Montgomery. 5238405Sjkim 6238405Sjkim THIS IS AN INTERNAL FUNCTION WITH A MUTABLE INTERFACE. IT IS ONLY SAFE 7238405Sjkim TO REACH THIS FUNCTION THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS 8238405Sjkim ALMOST GUARANTEED THAT THIS FUNCTION WILL CHANGE OR DISAPPEAR IN A FUTURE 9238405Sjkim GNU MP RELEASE. 10238405Sjkim 11238405Sjkim 12238405SjkimCopyright 1992, 1994, 1996, 2000, 2011, 2012 Free Software Foundation, Inc. 13238405Sjkim 14238405SjkimThis file is part of the GNU MP Library. 15238405Sjkim 16238405SjkimThe GNU MP Library is free software; you can redistribute it and/or modify 17238405Sjkimit under the terms of either: 18238405Sjkim 19238405Sjkim * the GNU Lesser General Public License as published by the Free 20238405Sjkim Software Foundation; either version 3 of the License, or (at your 21238405Sjkim option) any later version. 22238405Sjkim 23238405Sjkimor 24238405Sjkim 25238405Sjkim * the GNU General Public License as published by the Free Software 26238405Sjkim Foundation; either version 2 of the License, or (at your option) any 27238405Sjkim later version. 28238405Sjkim 29238405Sjkimor both in parallel, as here. 30238405Sjkim 31238405SjkimThe GNU MP Library is distributed in the hope that it will be useful, but 32238405SjkimWITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 33238405Sjkimor FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 34238405Sjkimfor more details. 35238405Sjkim 36238405SjkimYou should have received copies of the GNU General Public License and the 37238405SjkimGNU Lesser General Public License along with the GNU MP Library. If not, 38238405Sjkimsee https://www.gnu.org/licenses/. */ 39238405Sjkim 40238405Sjkim#include "gmp-impl.h" 41276861Sjkim#include "longlong.h" 42276861Sjkim 43238405Sjkimmp_limb_t 44238405Sjkimmpn_udiv_w_sdiv (mp_limb_t *rp, mp_limb_t a1, mp_limb_t a0, mp_limb_t d) 45238405Sjkim{ 46238405Sjkim mp_limb_t q, r; 47238405Sjkim mp_limb_t c0, c1, b1; 48238405Sjkim 49312826Sjkim ASSERT (d != 0); 50238405Sjkim ASSERT (a1 < d); 51238405Sjkim 52238405Sjkim if ((mp_limb_signed_t) d >= 0) 53276861Sjkim { 54276861Sjkim if (a1 < d - a1 - (a0 >> (GMP_LIMB_BITS - 1))) 55276861Sjkim { 56238405Sjkim /* dividend, divisor, and quotient are nonnegative */ 57344604Sjkim sdiv_qrnnd (q, r, a1, a0, d); 58344604Sjkim } 59344604Sjkim else 60344604Sjkim { 61344604Sjkim /* Compute c1*2^32 + c0 = a1*2^32 + a0 - 2^31*d */ 62344604Sjkim sub_ddmmss (c1, c0, a1, a0, d >> 1, d << (GMP_LIMB_BITS - 1)); 63238405Sjkim /* Divide (c1*2^32 + c0) by d */ 64344604Sjkim sdiv_qrnnd (q, r, c1, c0, d); 65344604Sjkim /* Add 2^31 to quotient */ 66344604Sjkim q += (mp_limb_t) 1 << (GMP_LIMB_BITS - 1); 67344604Sjkim } 68276861Sjkim } 69238405Sjkim else 70344604Sjkim { 71238405Sjkim b1 = d >> 1; /* d/2, between 2^30 and 2^31 - 1 */ 72238405Sjkim c1 = a1 >> 1; /* A/2 */ 73238405Sjkim c0 = (a1 << (GMP_LIMB_BITS - 1)) + (a0 >> 1); 74238405Sjkim 75238405Sjkim if (a1 < b1) /* A < 2^32*b1, so A/2 < 2^31*b1 */ 76238405Sjkim { 77238405Sjkim sdiv_qrnnd (q, r, c1, c0, b1); /* (A/2) / (d/2) */ 78238405Sjkim 79238405Sjkim r = 2*r + (a0 & 1); /* Remainder from A/(2*b1) */ 80238405Sjkim if ((d & 1) != 0) 81238405Sjkim { 82238405Sjkim if (r >= q) 83238405Sjkim r = r - q; 84238405Sjkim else if (q - r <= d) 85238405Sjkim { 86238405Sjkim r = r - q + d; 87238405Sjkim q--; 88238405Sjkim } 89238405Sjkim else 90238405Sjkim { 91238405Sjkim r = r - q + 2*d; 92238405Sjkim q -= 2; 93238405Sjkim } 94238405Sjkim } 95238405Sjkim } 96238405Sjkim else if (c1 < b1) /* So 2^31 <= (A/2)/b1 < 2^32 */ 97238405Sjkim { 98238405Sjkim c1 = (b1 - 1) - c1; 99238405Sjkim c0 = ~c0; /* logical NOT */ 100238405Sjkim 101238405Sjkim sdiv_qrnnd (q, r, c1, c0, b1); /* (A/2) / (d/2) */ 102238405Sjkim 103238405Sjkim q = ~q; /* (A/2)/b1 */ 104238405Sjkim r = (b1 - 1) - r; 105238405Sjkim 106238405Sjkim r = 2*r + (a0 & 1); /* A/(2*b1) */ 107238405Sjkim 108238405Sjkim if ((d & 1) != 0) 109238405Sjkim { 110238405Sjkim if (r >= q) 111238405Sjkim r = r - q; 112238405Sjkim else if (q - r <= d) 113238405Sjkim { 114238405Sjkim r = r - q + d; 115238405Sjkim q--; 116238405Sjkim } 117238405Sjkim else 118238405Sjkim { 119238405Sjkim r = r - q + 2*d; 120238405Sjkim q -= 2; 121238405Sjkim } 122238405Sjkim } 123238405Sjkim } 124238405Sjkim else /* Implies c1 = b1 */ 125238405Sjkim { /* Hence a1 = d - 1 = 2*b1 - 1 */ 126238405Sjkim if (a0 >= -d) 127238405Sjkim { 128238405Sjkim q = -CNST_LIMB(1); 129238405Sjkim r = a0 + d; 130238405Sjkim } 131238405Sjkim else 132238405Sjkim { 133238405Sjkim q = -CNST_LIMB(2); 134238405Sjkim r = a0 + 2*d; 135238405Sjkim } 136356290Sjkim } 137238405Sjkim } 138238405Sjkim 139238405Sjkim *rp = r; 140238405Sjkim return q; 141238405Sjkim} 142238405Sjkim