322 lines
6.5 KiB
C
322 lines
6.5 KiB
C
/* mpn_div -- Divide natural numbers, producing both remainder and
|
|
quotient.
|
|
|
|
Copyright (C) 1991 Free Software Foundation, Inc.
|
|
|
|
This file is part of the GNU MP Library.
|
|
|
|
The GNU MP Library is free software; you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation; either version 2, or (at your option)
|
|
any later version.
|
|
|
|
The GNU MP Library is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with the GNU MP Library; see the file COPYING. If not, write to
|
|
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
|
|
|
|
#include "gmp.h"
|
|
#include "gmp-impl.h"
|
|
#include "longlong.h"
|
|
|
|
/* Divide num (NUM_PTR/NUM_SIZE) by den (DEN_PTR/DEN_SIZE) and write
|
|
the quotient at QUOT_PTR and the remainder at NUM_PTR.
|
|
|
|
Return 0 or 1, depending on if the quotient size is (NSIZE - DSIZE)
|
|
or (NSIZE - DSIZE + 1).
|
|
|
|
Argument constraints:
|
|
1. The most significant bit of d must be set.
|
|
2. QUOT_PTR != DEN_PTR and QUOT_PTR != NUM_PTR, i.e. the quotient storage
|
|
area must be distinct from either input operands.
|
|
|
|
The exact sizes of the quotient and remainder must be determined
|
|
by the caller, in spite of the return value. The return value just
|
|
informs the caller about if the highest digit is written or not, and
|
|
it may very well be 0. */
|
|
|
|
/* THIS WILL BE IMPROVED SOON. MORE COMMENTS AND FASTER CODE. */
|
|
|
|
mp_size
|
|
#ifdef __STDC__
|
|
mpn_div (mp_ptr quot_ptr,
|
|
mp_ptr num_ptr, mp_size num_size,
|
|
mp_srcptr den_ptr, mp_size den_size)
|
|
#else
|
|
mpn_div (quot_ptr, num_ptr, num_size, den_ptr, den_size)
|
|
mp_ptr quot_ptr;
|
|
mp_ptr num_ptr;
|
|
mp_size num_size;
|
|
mp_srcptr den_ptr;
|
|
mp_size den_size;
|
|
#endif
|
|
{
|
|
mp_size q_is_long = 0;
|
|
|
|
switch (den_size)
|
|
{
|
|
case 0:
|
|
/* We are asked to divide by zero, so go ahead and do it!
|
|
(To make the compiler not remove this statement, assign NUM_SIZE
|
|
and fall through.) */
|
|
num_size = 1 / den_size;
|
|
|
|
case 1:
|
|
{
|
|
mp_size i;
|
|
mp_limb n1, n0;
|
|
mp_limb d;
|
|
|
|
d = den_ptr[0];
|
|
i = num_size - 1;
|
|
n1 = num_ptr[i];
|
|
i--;
|
|
|
|
if (n1 >= d)
|
|
{
|
|
q_is_long = 1;
|
|
n1 = 0;
|
|
i++;
|
|
}
|
|
|
|
for (; i >= 0; i--)
|
|
{
|
|
n0 = num_ptr[i];
|
|
udiv_qrnnd (quot_ptr[i], n1, n1, n0, d);
|
|
}
|
|
|
|
num_ptr[0] = n1;
|
|
}
|
|
break;
|
|
|
|
case 2:
|
|
{
|
|
mp_size i;
|
|
mp_limb n0, n1, n2;
|
|
mp_limb d0, d1;
|
|
|
|
num_ptr += num_size - 2;
|
|
d0 = den_ptr[1];
|
|
d1 = den_ptr[0];
|
|
n0 = num_ptr[1];
|
|
n1 = num_ptr[0];
|
|
|
|
if (n0 >= d0)
|
|
{
|
|
q_is_long = 1;
|
|
n1 = n0;
|
|
n0 = 0;
|
|
num_ptr++;
|
|
num_size++;
|
|
}
|
|
|
|
for (i = num_size - den_size - 1; i >= 0; i--)
|
|
{
|
|
mp_limb q;
|
|
mp_limb r;
|
|
|
|
num_ptr--;
|
|
if (n0 == d0)
|
|
{
|
|
/* Q should be either 111..111 or 111..110. Need special
|
|
treatment of this rare case as normal division would
|
|
give overflow. */
|
|
q = ~0;
|
|
|
|
r = n1 + d0;
|
|
if (r < d0) /* Carry in the addition? */
|
|
{
|
|
n2 = num_ptr[0];
|
|
|
|
add_ssaaaa (n0, n1, r - d1, n2, 0, d1);
|
|
quot_ptr[i] = q;
|
|
continue;
|
|
}
|
|
n0 = d1 - (d1 != 0);
|
|
n1 = -d1;
|
|
}
|
|
else
|
|
{
|
|
udiv_qrnnd (q, r, n0, n1, d0);
|
|
umul_ppmm (n0, n1, d1, q);
|
|
}
|
|
|
|
n2 = num_ptr[0];
|
|
q_test:
|
|
if (n0 > r || (n0 == r && n1 > n2))
|
|
{
|
|
/* The estimated Q was too large. */
|
|
q--;
|
|
|
|
sub_ddmmss (n0, n1, n0, n1, 0, d1);
|
|
r += d0;
|
|
if (r >= d0) /* If not carry, test q again. */
|
|
goto q_test;
|
|
}
|
|
|
|
quot_ptr[i] = q;
|
|
sub_ddmmss (n0, n1, r, n2, n0, n1);
|
|
}
|
|
num_ptr[1] = n0;
|
|
num_ptr[0] = n1;
|
|
}
|
|
break;
|
|
|
|
default:
|
|
{
|
|
mp_size i;
|
|
mp_limb d0 = den_ptr[den_size - 1];
|
|
mp_limb d1 = den_ptr[den_size - 2];
|
|
mp_limb n0 = num_ptr[num_size - 1];
|
|
int ugly_hack_flag = 0;
|
|
|
|
if (n0 >= d0)
|
|
{
|
|
|
|
/* There's a problem with this case, which shows up later in the
|
|
code. q becomes 1 (and sometimes 0) the first time when
|
|
we've been here, and w_cy == 0 after the main do-loops below.
|
|
But c = num_ptr[j] reads rubbish outside the num_ptr vector!
|
|
Maybe I can solve this cleanly when I fix the early-end
|
|
optimization here in the default case. For now, I change the
|
|
add_back entering condition, to kludge. Leaving the stray
|
|
memref behind!
|
|
|
|
HACK: Added ugly_hack_flag to make it work. */
|
|
|
|
q_is_long = 1;
|
|
n0 = 0;
|
|
num_size++;
|
|
ugly_hack_flag = 1;
|
|
}
|
|
|
|
num_ptr += num_size;
|
|
den_ptr += den_size;
|
|
for (i = num_size - den_size - 1; i >= 0; i--)
|
|
{
|
|
mp_limb q;
|
|
mp_limb n1;
|
|
mp_limb w_cy;
|
|
mp_limb d, c;
|
|
mp_size j;
|
|
|
|
num_ptr--;
|
|
if (n0 == d0)
|
|
/* This might over-estimate q, but it's probably not worth
|
|
the extra code here to find out. */
|
|
q = ~0;
|
|
else
|
|
{
|
|
mp_limb r;
|
|
|
|
udiv_qrnnd (q, r, n0, num_ptr[-1], d0);
|
|
umul_ppmm (n1, n0, d1, q);
|
|
|
|
while (n1 > r || (n1 == r && n0 > num_ptr[-2]))
|
|
{
|
|
q--;
|
|
r += d0;
|
|
if (r < d0) /* I.e. "carry in previous addition?" */
|
|
break;
|
|
n1 -= n0 < d1;
|
|
n0 -= d1;
|
|
}
|
|
}
|
|
|
|
w_cy = 0;
|
|
j = -den_size;
|
|
do
|
|
{
|
|
d = den_ptr[j];
|
|
c = num_ptr[j];
|
|
umul_ppmm (n1, n0, d, q);
|
|
n0 += w_cy;
|
|
w_cy = (n0 < w_cy) + n1;
|
|
n0 = c - n0;
|
|
num_ptr[j] = n0;
|
|
if (n0 > c)
|
|
goto cy_loop;
|
|
ncy_loop:
|
|
j++;
|
|
}
|
|
while (j < 0);
|
|
|
|
if (ugly_hack_flag)
|
|
{
|
|
c = 0;
|
|
ugly_hack_flag = 0;
|
|
}
|
|
else
|
|
c = num_ptr[j];
|
|
if (c >= w_cy)
|
|
goto store_q;
|
|
goto add_back;
|
|
|
|
do
|
|
{
|
|
d = den_ptr[j];
|
|
c = num_ptr[j];
|
|
umul_ppmm (n1, n0, d, q);
|
|
n0 += w_cy;
|
|
w_cy = (n0 < w_cy) + n1;
|
|
n0 = c - n0 - 1;
|
|
num_ptr[j] = n0;
|
|
if (n0 < c)
|
|
goto ncy_loop;
|
|
cy_loop:
|
|
j++;
|
|
}
|
|
while (j < 0);
|
|
|
|
if (ugly_hack_flag)
|
|
{
|
|
c = 0;
|
|
ugly_hack_flag = 0;
|
|
}
|
|
else
|
|
c = num_ptr[j];
|
|
w_cy++;
|
|
if (c >= w_cy)
|
|
goto store_q;
|
|
|
|
add_back:
|
|
j = -den_size;
|
|
do
|
|
{
|
|
d = den_ptr[j];
|
|
n0 = num_ptr[j] + d;
|
|
num_ptr[j] = n0;
|
|
if (n0 < d)
|
|
goto ab_cy_loop;
|
|
ab_ncy_loop:
|
|
j++;
|
|
}
|
|
while (j < 0);
|
|
abort (); /* We should always have a carry out! */
|
|
|
|
do
|
|
{
|
|
d = den_ptr[j];
|
|
n0 = num_ptr[j] + d + 1;
|
|
num_ptr[j] = n0;
|
|
if (n0 > d)
|
|
goto ab_ncy_loop;
|
|
ab_cy_loop:
|
|
j++;
|
|
}
|
|
while (j < 0);
|
|
q--;
|
|
|
|
store_q:
|
|
quot_ptr[i] = q;
|
|
}
|
|
}
|
|
}
|
|
|
|
return q_is_long;
|
|
}
|