← Back to branch summary

~ubuntu-branches/ubuntu/quantal/gclcvs/quantal

« back to all changes in this revision

Viewing changes to gmp3/mpn/generic/sb_divrem_mn.c

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2004-06-24 15:13:46 UTC
  • Revision ID: james.westby@ubuntu.com-20040624151346-xh0xaaktyyp7aorc
Tags: 2.7.0-26
C_GC_OFFSET is 2 on m68k-linux

Show diffs side-by-side

added added

removed removed

Lines of Context:
gmp3/mpn/generic/sb_divrem_mn.c
 
1
/* mpn_sb_divrem_mn -- Divide natural numbers, producing both remainder and
 
2
   quotient.
 
3
 
 
4
   THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE
 
5
   INTERFACES.  IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.
 
6
   IN FACT, IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A
 
7
   FUTURE GNU MP RELEASE.
 
8
 
 
9
 
 
10
Copyright 1993, 1994, 1995, 1996, 2000, 2001 Free Software Foundation, Inc.
 
11
 
 
12
This file is part of the GNU MP Library.
 
13
 
 
14
The GNU MP Library is free software; you can redistribute it and/or modify
 
15
it under the terms of the GNU Lesser General Public License as published by
 
16
the Free Software Foundation; either version 2.1 of the License, or (at your
 
17
option) any later version.
 
18
 
 
19
The GNU MP Library is distributed in the hope that it will be useful, but
 
20
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 
21
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 
22
License for more details.
 
23
 
 
24
You should have received a copy of the GNU Lesser General Public License
 
25
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
 
26
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
27
MA 02111-1307, USA. */
 
28
 
 
29
#include "gmp.h"
 
30
#include "gmp-impl.h"
 
31
#include "longlong.h"
 
32
 
 
33
 
 
34
/* The size where udiv_qrnnd_preinv should be used rather than udiv_qrnnd,
 
35
   meaning the quotient size where that should happen, the quotient size
 
36
   being how many udiv divisions will be done.  */
 
37
 
 
38
#ifndef SB_PREINV_THRESHOLD
 
39
# if UDIV_PREINV_ALWAYS
 
40
#  define SB_PREINV_THRESHOLD     0
 
41
# else
 
42
#  ifdef DIVREM_1_NORM_THRESHOLD
 
43
#   define SB_PREINV_THRESHOLD    DIVREM_1_NORM_THRESHOLD
 
44
#  else
 
45
#   if UDIV_TIME <= UDIV_NORM_PREINV_TIME
 
46
#    define SB_PREINV_THRESHOLD   MP_LIMB_T_MAX
 
47
#   else
 
48
#    define SB_PREINV_THRESHOLD \
 
49
       (1 + UDIV_TIME / (UDIV_TIME - UDIV_NORM_PREINV_TIME))
 
50
#   endif
 
51
#  endif
 
52
# endif
 
53
#endif
 
54
 
 
55
 
 
56
/* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
 
57
   the NSIZE-DSIZE least significant quotient limbs at QP
 
58
   and the DSIZE long remainder at NP.
 
59
   Return the most significant limb of the quotient, this is always 0 or 1.
 
60
 
 
61
   Preconditions:
 
62
   0. NSIZE >= DSIZE.
 
63
   1. The most significant bit of the divisor must be set.
 
64
   2. QP must either not overlap with the input operands at all, or
 
65
      QP + DSIZE >= NP must hold true.  (This means that it's
 
66
      possible to put the quotient in the high part of NUM, right after the
 
67
      remainder in NUM.
 
68
   3. NSIZE >= DSIZE.
 
69
   4. DSIZE >= 2.  */
 
70
 
 
71
 
 
72
mp_limb_t
 
73
mpn_sb_divrem_mn (mp_ptr qp,
 
74
                  mp_ptr np, mp_size_t nsize,
 
75
                  mp_srcptr dp, mp_size_t dsize)
 
76
{
 
77
  mp_limb_t most_significant_q_limb = 0;
 
78
  mp_size_t qsize = nsize - dsize;
 
79
  mp_size_t i;
 
80
  mp_limb_t dx, d1, n0;
 
81
  mp_limb_t dxinv;
 
82
  int use_preinv;
 
83
 
 
84
  ASSERT (dsize > 2);
 
85
  ASSERT (nsize >= dsize);
 
86
  ASSERT (dp[dsize-1] & MP_LIMB_T_HIGHBIT);
 
87
  ASSERT (! MPN_OVERLAP_P (np, nsize, dp, dsize));
 
88
  ASSERT (! MPN_OVERLAP_P (qp, nsize-dsize, dp, dsize));
 
89
  ASSERT (! MPN_OVERLAP_P (qp, nsize-dsize, np, nsize) || qp+dsize >= np);
 
90
 
 
91
  np += qsize;
 
92
  dx = dp[dsize - 1];
 
93
  d1 = dp[dsize - 2];
 
94
  n0 = np[dsize - 1];
 
95
 
 
96
  if (n0 >= dx)
 
97
    {
 
98
      if (n0 > dx || mpn_cmp (np, dp, dsize - 1) >= 0)
 
99
        {
 
100
          mpn_sub_n (np, np, dp, dsize);
 
101
          most_significant_q_limb = 1;
 
102
        }
 
103
    }
 
104
 
 
105
  /* use_preinv is possibly a constant, but it's left to the compiler to
 
106
     optimize away the unused code in that case.  */
 
107
  use_preinv = ABOVE_THRESHOLD (qsize, SB_PREINV_THRESHOLD);
 
108
  if (use_preinv)
 
109
    invert_limb (dxinv, dx);
 
110
 
 
111
  for (i = qsize-1; i >= 0; i--)
 
112
    {
 
113
      mp_limb_t q;
 
114
      mp_limb_t nx;
 
115
      mp_limb_t cy_limb;
 
116
 
 
117
      nx = np[dsize - 1];
 
118
      np--;
 
119
 
 
120
      if (nx == dx)
 
121
        {
 
122
          /* This might over-estimate q, but it's probably not worth
 
123
             the extra code here to find out.  */
 
124
          q = ~(mp_limb_t) 0;
 
125
 
 
126
#if 1
 
127
          cy_limb = mpn_submul_1 (np, dp, dsize, q);
 
128
#else
 
129
          /* This should be faster on many machines */
 
130
          cy_limb = mpn_sub_n (np + 1, np + 1, dp, dsize);
 
131
          cy = mpn_add_n (np, np, dp, dsize);
 
132
          np[dsize] += cy;
 
133
#endif
 
134
 
 
135
          if (nx != cy_limb)
 
136
            {
 
137
              mpn_add_n (np, np, dp, dsize);
 
138
              q--;
 
139
            }
 
140
 
 
141
          qp[i] = q;
 
142
        }
 
143
      else
 
144
        {
 
145
          mp_limb_t rx, r1, r0, p1, p0;
 
146
 
 
147
          /* "workaround" avoids a problem with gcc 2.7.2.3 i386 register
 
148
             usage when np[dsize-1] is used in an asm statement like
 
149
             umul_ppmm in udiv_qrnnd_preinv.  The symptom is seg faults due
 
150
             to registers being clobbered.  gcc 2.95 i386 doesn't have the
 
151
             problem. */
 
152
          {
 
153
            mp_limb_t  workaround = np[dsize - 1];
 
154
            if (use_preinv)
 
155
              udiv_qrnnd_preinv (q, r1, nx, workaround, dx, dxinv);
 
156
            else
 
157
              udiv_qrnnd (q, r1, nx, workaround, dx);
 
158
          }
 
159
          umul_ppmm (p1, p0, d1, q);
 
160
 
 
161
          r0 = np[dsize - 2];
 
162
          rx = 0;
 
163
          if (r1 < p1 || (r1 == p1 && r0 < p0))
 
164
            {
 
165
              p1 -= p0 < d1;
 
166
              p0 -= d1;
 
167
              q--;
 
168
              r1 += dx;
 
169
              rx = r1 < dx;
 
170
            }
 
171
 
 
172
          p1 += r0 < p0;        /* cannot carry! */
 
173
          rx -= r1 < p1;        /* may become 11..1 if q is still too large */
 
174
          r1 -= p1;
 
175
          r0 -= p0;
 
176
 
 
177
          cy_limb = mpn_submul_1 (np, dp, dsize - 2, q);
 
178
 
 
179
          {
 
180
            mp_limb_t cy1, cy2;
 
181
            cy1 = r0 < cy_limb;
 
182
            r0 -= cy_limb;
 
183
            cy2 = r1 < cy1;
 
184
            r1 -= cy1;
 
185
            np[dsize - 1] = r1;
 
186
            np[dsize - 2] = r0;
 
187
            if (cy2 != rx)
 
188
              {
 
189
                mpn_add_n (np, np, dp, dsize);
 
190
                q--;
 
191
              }
 
192
          }
 
193
          qp[i] = q;
 
194
        }
 
195
    }
 
196
 
 
197
  /* ______ ______ ______
 
198
    |__rx__|__r1__|__r0__|              partial remainder
 
199
            ______ ______
 
200
         - |__p1__|__p0__|              partial product to subtract
 
201
            ______ ______
 
202
         - |______|cylimb|              
 
203
 
 
204
     rx is -1, 0 or 1.  If rx=1, then q is correct (it should match
 
205
     carry out).  If rx=-1 then q is too large.  If rx=0, then q might
 
206
     be too large, but it is most likely correct.
 
207
  */
 
208
 
 
209
  return most_significant_q_limb;
 
210
}