1 /* mpn_preinv_divrem_1 -- mpn by limb division with pre-inverted divisor.
2
3 THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
4 CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
5 FUTURE GNU MP RELEASES.
6
7 Copyright 2000-2003 Free Software Foundation, Inc.
8
9 This file is part of the GNU MP Library.
10
11 The GNU MP Library is free software; you can redistribute it and/or modify
12 it under the terms of either:
13
14 * the GNU Lesser General Public License as published by the Free
15 Software Foundation; either version 3 of the License, or (at your
16 option) any later version.
17
18 or
19
20 * the GNU General Public License as published by the Free Software
21 Foundation; either version 2 of the License, or (at your option) any
22 later version.
23
24 or both in parallel, as here.
25
26 The GNU MP Library is distributed in the hope that it will be useful, but
27 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
28 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
29 for more details.
30
31 You should have received copies of the GNU General Public License and the
32 GNU Lesser General Public License along with the GNU MP Library. If not,
33 see https://www.gnu.org/licenses/. */
34
35 #include "gmp-impl.h"
36 #include "longlong.h"
37
38
39 /* Don't bloat a shared library with unused code. */
40 #if USE_PREINV_DIVREM_1
41
42 /* Same test here for skipping one divide step as in mpn_divrem_1.
43
44 The main reason for a separate shift==0 case is that not all CPUs give
45 zero for "n0 >> GMP_LIMB_BITS" which would arise in the general case
46 code used on shift==0. shift==0 is also reasonably common in mp_bases
47 big_base, for instance base==10 on a 64-bit limb.
48
49 Under shift!=0 it would be possible to call mpn_lshift to adjust the
50 dividend all in one go (into the quotient space say), rather than
51 limb-by-limb in the loop. This might help if mpn_lshift is a lot faster
52 than what the compiler can generate for EXTRACT. But this is left to CPU
53 specific implementations to consider, especially since EXTRACT isn't on
54 the dependent chain.
55
56 If size==0 then the result is simply xsize limbs of zeros, but nothing
57 special is done for that, since it wouldn't be a usual call, and
58 certainly never arises from mpn_get_str which is our main caller. */
59
60 mp_limb_t
mpn_preinv_divrem_1(mp_ptr qp,mp_size_t xsize,mp_srcptr ap,mp_size_t size,mp_limb_t d_unnorm,mp_limb_t dinv,int shift)61 mpn_preinv_divrem_1 (mp_ptr qp, mp_size_t xsize,
62 mp_srcptr ap, mp_size_t size, mp_limb_t d_unnorm,
63 mp_limb_t dinv, int shift)
64 {
65 mp_limb_t ahigh, qhigh, r;
66 mp_size_t i;
67 mp_limb_t n1, n0;
68 mp_limb_t d;
69
70 ASSERT (xsize >= 0);
71 ASSERT (size >= 1);
72 ASSERT (d_unnorm != 0);
73 #if WANT_ASSERT
74 {
75 int want_shift;
76 mp_limb_t want_dinv;
77 count_leading_zeros (want_shift, d_unnorm);
78 ASSERT (shift == want_shift);
79 invert_limb (want_dinv, d_unnorm << shift);
80 ASSERT (dinv == want_dinv);
81 }
82 #endif
83 /* FIXME: What's the correct overlap rule when xsize!=0? */
84 ASSERT (MPN_SAME_OR_SEPARATE_P (qp+xsize, ap, size));
85
86 ahigh = ap[size-1];
87 d = d_unnorm << shift;
88 qp += (size + xsize - 1); /* dest high limb */
89
90 if (shift == 0)
91 {
92 /* High quotient limb is 0 or 1, and skip a divide step. */
93 r = ahigh;
94 qhigh = (r >= d);
95 r = (qhigh ? r-d : r);
96 *qp-- = qhigh;
97 size--;
98
99 for (i = size-1; i >= 0; i--)
100 {
101 n0 = ap[i];
102 udiv_qrnnd_preinv (*qp, r, r, n0, d, dinv);
103 qp--;
104 }
105 }
106 else
107 {
108 r = 0;
109 if (ahigh < d_unnorm)
110 {
111 r = ahigh << shift;
112 *qp-- = 0;
113 size--;
114 if (size == 0)
115 goto done_integer;
116 }
117
118 n1 = ap[size-1];
119 r |= n1 >> (GMP_LIMB_BITS - shift);
120
121 for (i = size-2; i >= 0; i--)
122 {
123 ASSERT (r < d);
124 n0 = ap[i];
125 udiv_qrnnd_preinv (*qp, r, r,
126 ((n1 << shift) | (n0 >> (GMP_LIMB_BITS - shift))),
127 d, dinv);
128 qp--;
129 n1 = n0;
130 }
131 udiv_qrnnd_preinv (*qp, r, r, n1 << shift, d, dinv);
132 qp--;
133 }
134
135 done_integer:
136 for (i = 0; i < xsize; i++)
137 {
138 udiv_qrnnd_preinv (*qp, r, r, CNST_LIMB(0), d, dinv);
139 qp--;
140 }
141
142 return r >> shift;
143 }
144
145 #endif /* USE_PREINV_DIVREM_1 */
146