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