← Back to branch summary

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

« back to all changes in this revision

Viewing changes to gmp3/tests/refmpf.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/tests/refmpf.c
 
1
/* Reference floating point routines.
 
2
 
 
3
Copyright 1996, 2001 Free Software Foundation, Inc.
 
4
 
 
5
This file is part of the GNU MP Library.
 
6
 
 
7
The GNU MP Library is free software; you can redistribute it and/or modify
 
8
it under the terms of the GNU Lesser General Public License as published by
 
9
the Free Software Foundation; either version 2.1 of the License, or (at your
 
10
option) any later version.
 
11
 
 
12
The GNU MP Library is distributed in the hope that it will be useful, but
 
13
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 
14
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 
15
License for more details.
 
16
 
 
17
You should have received a copy of the GNU Lesser General Public License
 
18
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
 
19
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
20
MA 02111-1307, USA. */
 
21
 
 
22
#include "gmp.h"
 
23
#include "gmp-impl.h"
 
24
#include "tests.h"
 
25
 
 
26
 
 
27
void
 
28
refmpf_add (mpf_ptr w, mpf_srcptr u, mpf_srcptr v)
 
29
{
 
30
  mp_size_t hi, lo, size;
 
31
  mp_ptr ut, vt, wt;
 
32
  int neg;
 
33
  mp_exp_t exp;
 
34
  mp_limb_t cy;
 
35
  TMP_DECL (mark);
 
36
 
 
37
  TMP_MARK (mark);
 
38
 
 
39
  if (SIZ (u) == 0)
 
40
    {
 
41
      size = ABSIZ (v);
 
42
      wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
 
43
      MPN_COPY (wt, PTR (v), size);
 
44
      exp = EXP (v);
 
45
      neg = SIZ (v) < 0;
 
46
      goto done;
 
47
    }
 
48
  if (SIZ (v) == 0)
 
49
    {
 
50
      size = ABSIZ (u);
 
51
      wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
 
52
      MPN_COPY (wt, PTR (u), size);
 
53
      exp = EXP (u);
 
54
      neg = SIZ (u) < 0;
 
55
      goto done;
 
56
    }
 
57
  if ((SIZ (u) ^ SIZ (v)) < 0)
 
58
    {
 
59
      mpf_t tmp;
 
60
      SIZ (tmp) = -SIZ (v);
 
61
      EXP (tmp) = EXP (v);
 
62
      PTR (tmp) = PTR (v);
 
63
      refmpf_sub (w, u, tmp);
 
64
      return;
 
65
    }
 
66
  neg = SIZ (u) < 0;
 
67
 
 
68
  /* Compute the significance of the hi and lo end of the result.  */
 
69
  hi = MAX (EXP (u), EXP (v));
 
70
  lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v));
 
71
  size = hi - lo;
 
72
  ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
 
73
  vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
 
74
  wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
 
75
  MPN_ZERO (ut, size);
 
76
  MPN_ZERO (vt, size);
 
77
  {int off;
 
78
  off = size + (EXP (u) - hi) - ABSIZ (u);
 
79
  MPN_COPY (ut + off, PTR (u), ABSIZ (u));
 
80
  off = size + (EXP (v) - hi) - ABSIZ (v);
 
81
  MPN_COPY (vt + off, PTR (v), ABSIZ (v));
 
82
  }
 
83
 
 
84
  cy = mpn_add_n (wt, ut, vt, size);
 
85
  wt[size] = cy;
 
86
  size += cy;
 
87
  exp = hi + cy;
 
88
 
 
89
done:
 
90
  if (size > PREC (w))
 
91
    {
 
92
      wt += size - PREC (w);
 
93
      size = PREC (w);
 
94
    }
 
95
  MPN_COPY (PTR (w), wt, size);
 
96
  SIZ (w) = neg == 0 ? size : -size;
 
97
  EXP (w) = exp;
 
98
  TMP_FREE (mark);
 
99
}
 
100
 
 
101
void
 
102
refmpf_sub (mpf_ptr w, mpf_srcptr u, mpf_srcptr v)
 
103
{
 
104
  mp_size_t hi, lo, size;
 
105
  mp_ptr ut, vt, wt;
 
106
  int neg;
 
107
  mp_exp_t exp;
 
108
  TMP_DECL (mark);
 
109
 
 
110
  TMP_MARK (mark);
 
111
 
 
112
  if (SIZ (u) == 0)
 
113
    {
 
114
      size = ABSIZ (v);
 
115
      wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
 
116
      MPN_COPY (wt, PTR (v), size);
 
117
      exp = EXP (v);
 
118
      neg = SIZ (v) > 0;
 
119
      goto done;
 
120
    }
 
121
  if (SIZ (v) == 0)
 
122
    {
 
123
      size = ABSIZ (u);
 
124
      wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
 
125
      MPN_COPY (wt, PTR (u), size);
 
126
      exp = EXP (u);
 
127
      neg = SIZ (u) < 0;
 
128
      goto done;
 
129
    }
 
130
  if ((SIZ (u) ^ SIZ (v)) < 0)
 
131
    {
 
132
      mpf_t tmp;
 
133
      SIZ (tmp) = -SIZ (v);
 
134
      EXP (tmp) = EXP (v);
 
135
      PTR (tmp) = PTR (v);
 
136
      refmpf_add (w, u, tmp);
 
137
      if (SIZ (u) < 0)
 
138
        mpf_neg (w, w);
 
139
      return;
 
140
    }
 
141
  neg = SIZ (u) < 0;
 
142
 
 
143
  /* Compute the significance of the hi and lo end of the result.  */
 
144
  hi = MAX (EXP (u), EXP (v));
 
145
  lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v));
 
146
  size = hi - lo;
 
147
  ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
 
148
  vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
 
149
  wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
 
150
  MPN_ZERO (ut, size);
 
151
  MPN_ZERO (vt, size);
 
152
  {int off;
 
153
  off = size + (EXP (u) - hi) - ABSIZ (u);
 
154
  MPN_COPY (ut + off, PTR (u), ABSIZ (u));
 
155
  off = size + (EXP (v) - hi) - ABSIZ (v);
 
156
  MPN_COPY (vt + off, PTR (v), ABSIZ (v));
 
157
  }
 
158
 
 
159
  if (mpn_cmp (ut, vt, size) >= 0)
 
160
    mpn_sub_n (wt, ut, vt, size);
 
161
  else
 
162
    {
 
163
      mpn_sub_n (wt, vt, ut, size);
 
164
      neg ^= 1;
 
165
    }
 
166
  exp = hi;
 
167
  while (size != 0 && wt[size - 1] == 0)
 
168
    {
 
169
      size--;
 
170
      exp--;
 
171
    }
 
172
 
 
173
done:
 
174
  if (size > PREC (w))
 
175
    {
 
176
      wt += size - PREC (w);
 
177
      size = PREC (w);
 
178
    }
 
179
  MPN_COPY (PTR (w), wt, size);
 
180
  SIZ (w) = neg == 0 ? size : -size;
 
181
  EXP (w) = exp;
 
182
  TMP_FREE (mark);
 
183
}