← Back to branch summary

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

« back to all changes in this revision

Viewing changes to gmp3/mpf/sub.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/mpf/sub.c
 
1
/* mpf_sub -- Subtract two floats.
 
2
 
 
3
Copyright 1993, 1994, 1995, 1996, 1999, 2000, 2001 Free Software Foundation,
 
4
Inc.
 
5
 
 
6
This file is part of the GNU MP Library.
 
7
 
 
8
The GNU MP Library is free software; you can redistribute it and/or modify
 
9
it under the terms of the GNU Lesser General Public License as published by
 
10
the Free Software Foundation; either version 2.1 of the License, or (at your
 
11
option) any later version.
 
12
 
 
13
The GNU MP Library is distributed in the hope that it will be useful, but
 
14
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 
15
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 
16
License for more details.
 
17
 
 
18
You should have received a copy of the GNU Lesser General Public License
 
19
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
 
20
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
21
MA 02111-1307, USA. */
 
22
 
 
23
#include "gmp.h"
 
24
#include "gmp-impl.h"
 
25
 
 
26
void
 
27
mpf_sub (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
 
28
{
 
29
  mp_srcptr up, vp;
 
30
  mp_ptr rp, tp;
 
31
  mp_size_t usize, vsize, rsize;
 
32
  mp_size_t prec;
 
33
  mp_exp_t exp;
 
34
  mp_size_t ediff;
 
35
  int negate;
 
36
  TMP_DECL (marker);
 
37
 
 
38
  usize = u->_mp_size;
 
39
  vsize = v->_mp_size;
 
40
 
 
41
  /* Handle special cases that don't work in generic code below.  */
 
42
  if (usize == 0)
 
43
    {
 
44
      mpf_neg (r, v);
 
45
      return;
 
46
    }
 
47
  if (vsize == 0)
 
48
    {
 
49
      if (r != u)
 
50
        mpf_set (r, u);
 
51
      return;
 
52
    }
 
53
 
 
54
  /* If signs of U and V are different, perform addition.  */
 
55
  if ((usize ^ vsize) < 0)
 
56
    {
 
57
      __mpf_struct v_negated;
 
58
      v_negated._mp_size = -vsize;
 
59
      v_negated._mp_exp = v->_mp_exp;
 
60
      v_negated._mp_d = v->_mp_d;
 
61
      mpf_add (r, u, &v_negated);
 
62
      return;
 
63
    }
 
64
 
 
65
  TMP_MARK (marker);
 
66
 
 
67
  /* Signs are now known to be the same.  */
 
68
  negate = usize < 0;
 
69
 
 
70
  /* Make U be the operand with the largest exponent.  */
 
71
  if (u->_mp_exp < v->_mp_exp)
 
72
    {
 
73
      mpf_srcptr t;
 
74
      t = u; u = v; v = t;
 
75
      negate ^= 1;
 
76
      usize = u->_mp_size;
 
77
      vsize = v->_mp_size;
 
78
    }
 
79
 
 
80
  usize = ABS (usize);
 
81
  vsize = ABS (vsize);
 
82
  up = u->_mp_d;
 
83
  vp = v->_mp_d;
 
84
  rp = r->_mp_d;
 
85
  prec = r->_mp_prec + 1;
 
86
  exp = u->_mp_exp;
 
87
  ediff = u->_mp_exp - v->_mp_exp;
 
88
 
 
89
  /* If ediff is 0 or 1, we might have a situation where the operands are
 
90
     extremely close.  We need to scan the operands from the most significant
 
91
     end ignore the initial parts that are equal.  */
 
92
  if (ediff <= 1)
 
93
    {
 
94
      if (ediff == 0)
 
95
        {
 
96
          /* Skip leading limbs in U and V that are equal.  */
 
97
          if (up[usize - 1] == vp[vsize - 1])
 
98
            {
 
99
              /* This loop normally exits immediately.  Optimize for that.  */
 
100
              do
 
101
                {
 
102
                  usize--;
 
103
                  vsize--;
 
104
                  exp--;
 
105
 
 
106
                  if (usize == 0)
 
107
                    {
 
108
                      if (vsize > prec)
 
109
                        {
 
110
                          vp += vsize - prec;
 
111
                          vsize = prec;
 
112
                        }
 
113
                      rsize = vsize;
 
114
                      tp = (mp_ptr) vp;
 
115
                      negate ^= 1;
 
116
                      goto normalize;
 
117
                    }
 
118
                  if (vsize == 0)
 
119
                    {
 
120
                      if (usize > prec)
 
121
                        {
 
122
                          up += usize - prec;
 
123
                          usize = prec;
 
124
                        }
 
125
                      rsize = usize;
 
126
                      tp = (mp_ptr) up;
 
127
                      goto normalize;
 
128
                    }
 
129
                }
 
130
              while (up[usize - 1] == vp[vsize - 1]);
 
131
            }
 
132
 
 
133
          if (up[usize - 1] < vp[vsize - 1])
 
134
            {
 
135
              /* For simplicity, swap U and V.  Note that since the loop above
 
136
                 wouldn't have exited unless up[usize - 1] and vp[vsize - 1]
 
137
                 were non-equal, this if-statement catches all cases where U
 
138
                 is smaller than V.  */
 
139
              MPN_SRCPTR_SWAP (up,usize, vp,vsize);
 
140
              negate ^= 1;
 
141
              /* negating ediff not necessary since it is 0.  */
 
142
            }
 
143
 
 
144
          /* Check for
 
145
             x+1 00000000 ...
 
146
              x  ffffffff ... */
 
147
          if (up[usize - 1] != vp[vsize - 1] + 1)
 
148
            goto general_case;
 
149
          usize--;
 
150
          vsize--;
 
151
          exp--;
 
152
        }
 
153
      else /* ediff == 1 */
 
154
        {
 
155
          /* Check for
 
156
             1 00000000 ...
 
157
             0 ffffffff ... */
 
158
 
 
159
          if (up[usize - 1] != 1 || vp[vsize - 1] != ~(mp_limb_t) 0
 
160
              || (usize >= 2 && up[usize - 2] != 0))
 
161
            goto general_case;
 
162
 
 
163
          usize--;
 
164
          exp--;
 
165
        }
 
166
 
 
167
      /* Skip sequences of 00000000/ffffffff */
 
168
      while (vsize != 0 && usize != 0 && up[usize - 1] == 0
 
169
             && vp[vsize - 1] == ~(mp_limb_t) 0)
 
170
        {
 
171
          usize--;
 
172
          vsize--;
 
173
          exp--;
 
174
        }
 
175
 
 
176
      if (usize == 0)
 
177
        {
 
178
          while (vsize != 0 && vp[vsize - 1] == ~(mp_limb_t) 0)
 
179
            {
 
180
              vsize--;
 
181
              exp--;
 
182
            }
 
183
        }
 
184
 
 
185
      if (usize > prec - 1)
 
186
        {
 
187
          up += usize - (prec - 1);
 
188
          usize = prec - 1;
 
189
        }
 
190
      if (vsize > prec - 1)
 
191
        {
 
192
          vp += vsize - (prec - 1);
 
193
          vsize = prec - 1;
 
194
        }
 
195
 
 
196
      tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
 
197
      {
 
198
        mp_limb_t cy_limb;
 
199
        if (vsize == 0)
 
200
          {
 
201
            mp_size_t size, i;
 
202
            size = usize;
 
203
            for (i = 0; i < size; i++)
 
204
              tp[i] = up[i];
 
205
            tp[size] = 1;
 
206
            rsize = size + 1;
 
207
            exp++;
 
208
            goto normalize;
 
209
          }
 
210
        if (usize == 0)
 
211
          {
 
212
            mp_size_t size, i;
 
213
            size = vsize;
 
214
            for (i = 0; i < size; i++)
 
215
              tp[i] = ~vp[i];
 
216
            cy_limb = 1 - mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1);
 
217
            rsize = vsize;
 
218
            if (cy_limb == 0)
 
219
              {
 
220
                tp[rsize] = 1;
 
221
                rsize++;
 
222
                exp++;
 
223
              }
 
224
            goto normalize;
 
225
          }
 
226
        if (usize >= vsize)
 
227
          {
 
228
            /* uuuu     */
 
229
            /* vv       */
 
230
            mp_size_t size;
 
231
            size = usize - vsize;
 
232
            MPN_COPY (tp, up, size);
 
233
            cy_limb = mpn_sub_n (tp + size, up + size, vp, vsize);
 
234
            rsize = usize;
 
235
          }
 
236
        else /* (usize < vsize) */
 
237
          {
 
238
            /* uuuu     */
 
239
            /* vvvvvvv  */
 
240
            mp_size_t size, i;
 
241
            size = vsize - usize;
 
242
            for (i = 0; i < size; i++)
 
243
              tp[i] = ~vp[i];
 
244
            cy_limb = mpn_sub_n (tp + size, up, vp + size, usize);
 
245
            cy_limb+= mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
 
246
            cy_limb-= mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1);
 
247
            rsize = vsize;
 
248
          }
 
249
        if (cy_limb == 0)
 
250
          {
 
251
            tp[rsize] = 1;
 
252
            rsize++;
 
253
            exp++;
 
254
          }
 
255
        goto normalize;
 
256
      }
 
257
    }
 
258
 
 
259
general_case:
 
260
  /* If U extends beyond PREC, ignore the part that does.  */
 
261
  if (usize > prec)
 
262
    {
 
263
      up += usize - prec;
 
264
      usize = prec;
 
265
    }
 
266
 
 
267
  /* If V extends beyond PREC, ignore the part that does.
 
268
     Note that this may make vsize negative.  */
 
269
  if (vsize + ediff > prec)
 
270
    {
 
271
      vp += vsize + ediff - prec;
 
272
      vsize = prec - ediff;
 
273
    }
 
274
 
 
275
  /* Allocate temp space for the result.  Allocate
 
276
     just vsize + ediff later???  */
 
277
  tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB);
 
278
 
 
279
  if (ediff >= prec)
 
280
    {
 
281
      /* V completely cancelled.  */
 
282
      if (tp != up)
 
283
        MPN_COPY (rp, up, usize);
 
284
      rsize = usize;
 
285
    }
 
286
  else
 
287
    {
 
288
      /* Locate the least significant non-zero limb in (the needed
 
289
         parts of) U and V, to simplify the code below.  */
 
290
      for (;;)
 
291
        {
 
292
          if (vsize == 0)
 
293
            {
 
294
              MPN_COPY (rp, up, usize);
 
295
              rsize = usize;
 
296
              goto done;
 
297
            }
 
298
          if (vp[0] != 0)
 
299
            break;
 
300
          vp++, vsize--;
 
301
        }
 
302
      for (;;)
 
303
        {
 
304
          if (usize == 0)
 
305
            {
 
306
              MPN_COPY (rp, vp, vsize);
 
307
              rsize = vsize;
 
308
              negate ^= 1;
 
309
              goto done;
 
310
            }
 
311
          if (up[0] != 0)
 
312
            break;
 
313
          up++, usize--;
 
314
        }
 
315
 
 
316
      /* uuuu     |  uuuu     |  uuuu     |  uuuu     |  uuuu    */
 
317
      /* vvvvvvv  |  vv       |    vvvvv  |    v      |       vv */
 
318
 
 
319
      if (usize > ediff)
 
320
        {
 
321
          /* U and V partially overlaps.  */
 
322
          if (ediff == 0)
 
323
            {
 
324
              /* Have to compare the leading limbs of u and v
 
325
                 to determine whether to compute u - v or v - u.  */
 
326
              if (usize >= vsize)
 
327
                {
 
328
                  /* uuuu     */
 
329
                  /* vv       */
 
330
                  mp_size_t size;
 
331
                  size = usize - vsize;
 
332
                  MPN_COPY (tp, up, size);
 
333
                  mpn_sub_n (tp + size, up + size, vp, vsize);
 
334
                  rsize = usize;
 
335
                }
 
336
              else /* (usize < vsize) */
 
337
                {
 
338
                  /* uuuu     */
 
339
                  /* vvvvvvv  */
 
340
                  mp_size_t size, i;
 
341
                  size = vsize - usize;
 
342
                  tp[0] = -vp[0];
 
343
                  for (i = 1; i < size; i++)
 
344
                    tp[i] = ~vp[i];
 
345
                  mpn_sub_n (tp + size, up, vp + size, usize);
 
346
                  mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
 
347
                  rsize = vsize;
 
348
                }
 
349
            }
 
350
          else
 
351
            {
 
352
              if (vsize + ediff <= usize)
 
353
                {
 
354
                  /* uuuu     */
 
355
                  /*   v      */
 
356
                  mp_size_t size;
 
357
                  size = usize - ediff - vsize;
 
358
                  MPN_COPY (tp, up, size);
 
359
                  mpn_sub (tp + size, up + size, usize - size, vp, vsize);
 
360
                  rsize = usize;
 
361
                }
 
362
              else
 
363
                {
 
364
                  /* uuuu     */
 
365
                  /*   vvvvv  */
 
366
                  mp_size_t size, i;
 
367
                  size = vsize + ediff - usize;
 
368
                  tp[0] = -vp[0];
 
369
                  for (i = 1; i < size; i++)
 
370
                    tp[i] = ~vp[i];
 
371
                  mpn_sub (tp + size, up, usize, vp + size, usize - ediff);
 
372
                  mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1);
 
373
                  rsize = vsize + ediff;
 
374
                }
 
375
            }
 
376
        }
 
377
      else
 
378
        {
 
379
          /* uuuu     */
 
380
          /*      vv  */
 
381
          mp_size_t size, i;
 
382
          size = vsize + ediff - usize;
 
383
          tp[0] = -vp[0];
 
384
          for (i = 1; i < vsize; i++)
 
385
            tp[i] = ~vp[i];
 
386
          for (i = vsize; i < size; i++)
 
387
            tp[i] = ~(mp_limb_t) 0;
 
388
          mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1);
 
389
          rsize = size + usize;
 
390
        }
 
391
 
 
392
    normalize:
 
393
      /* Full normalize.  Optimize later.  */
 
394
      while (rsize != 0 && tp[rsize - 1] == 0)
 
395
        {
 
396
          rsize--;
 
397
          exp--;
 
398
        }
 
399
      MPN_COPY (rp, tp, rsize);
 
400
    }
 
401
 
 
402
 done:
 
403
  r->_mp_size = negate ? -rsize : rsize;
 
404
  r->_mp_exp = exp;
 
405
  TMP_FREE (marker);
 
406
}