← Back to branch summary

~ubuntu-branches/debian/sid/genius/sid

« back to all changes in this revision

Viewing changes to mpfr/gmp_op.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel Holbach
  • Date: 2006-08-21 12:57:45 UTC
  • Revision ID: james.westby@ubuntu.com-20060821125745-sl9ks8v7fq324bdf
Tags: upstream-0.7.6.1
Import upstream version 0.7.6.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
mpfr/gmp_op.c
 
1
/* Implementations of operations between mpfr and mpz/mpq data
 
2
 
 
3
Copyright 2001, 2003, 2004, 2005 Free Software Foundation, Inc.
 
4
 
 
5
This file is part of the MPFR Library.
 
6
 
 
7
The MPFR 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 MPFR 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 MPFR Library; see the file COPYING.LIB.  If not, write to
 
19
the Free Software Foundation, Inc., 51 Franklin Place, Fifth Floor, Boston,
 
20
MA 02110-1301, USA. */
 
21
 
 
22
#include <stddef.h>
 
23
 
 
24
#define MPFR_NEED_LONGLONG_H
 
25
#include "mpfr-impl.h"
 
26
 
 
27
/* Init and set a mpfr_t with enought precision to store a mpz */
 
28
static void
 
29
init_set_z (mpfr_ptr t, mpz_srcptr z)
 
30
{
 
31
  mp_prec_t p;
 
32
  int i;
 
33
 
 
34
  if (mpz_size (z) <= 1)
 
35
    p = BITS_PER_MP_LIMB;
 
36
  else
 
37
    MPFR_MPZ_SIZEINBASE2 (p, z);
 
38
  mpfr_init2 (t, p);
 
39
  i = mpfr_set_z (t, z, GMP_RNDN);
 
40
  MPFR_ASSERTD (i == 0);
 
41
}
 
42
 
 
43
/* Init, set a mpfr_t with enought precision to store a mpz_t without round,
 
44
   call the function, and clear the allocated mpfr_t  */
 
45
static int
 
46
foo (mpfr_ptr x, mpfr_srcptr y, mpz_srcptr z, mp_rnd_t r,
 
47
     int (*f)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t))
 
48
{
 
49
  mpfr_t t;
 
50
  int i;
 
51
  init_set_z (t, z);
 
52
  i = (*f) (x, y, t, r);
 
53
  mpfr_clear (t);
 
54
  return i;
 
55
}
 
56
 
 
57
int
 
58
mpfr_mul_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r)
 
59
{
 
60
  return foo (y, x, z, r, mpfr_mul);
 
61
}
 
62
 
 
63
int
 
64
mpfr_div_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r)
 
65
{
 
66
  return foo (y, x, z, r, mpfr_div);
 
67
}
 
68
 
 
69
int
 
70
mpfr_add_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r)
 
71
{
 
72
  /* Mpz 0 is unsigned */
 
73
  if (MPFR_UNLIKELY (mpz_sgn (z) == 0))
 
74
    return mpfr_set (y, x, r);
 
75
  else
 
76
    return foo (y, x, z, r, mpfr_add);
 
77
}
 
78
 
 
79
int
 
80
mpfr_sub_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r)
 
81
{
 
82
  /* Mpz 0 is unsigned */
 
83
  if (MPFR_UNLIKELY (mpz_sgn (z) == 0))
 
84
    return mpfr_set (y, x, r);
 
85
  else
 
86
    return foo (y, x, z, r, mpfr_sub);
 
87
}
 
88
 
 
89
int
 
90
mpfr_cmp_z (mpfr_srcptr x, mpz_srcptr z)
 
91
{
 
92
  mpfr_t t;
 
93
  int res;
 
94
  init_set_z (t, z);
 
95
  res = mpfr_cmp (x, t);
 
96
  mpfr_clear (t);
 
97
  return res;
 
98
}
 
99
 
 
100
int
 
101
mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode)
 
102
{
 
103
  mpfr_t tmp;
 
104
  int res;
 
105
  mp_prec_t p;
 
106
 
 
107
  if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
 
108
    return mpfr_mul_ui (y, x, 0, rnd_mode);
 
109
  else
 
110
    {
 
111
      MPFR_MPZ_SIZEINBASE2 (p, mpq_numref (z));
 
112
      mpfr_init2 (tmp, MPFR_PREC (x) + p);
 
113
      res = mpfr_mul_z (tmp, x, mpq_numref(z), GMP_RNDN );
 
114
      MPFR_ASSERTD (res == 0);
 
115
      res = mpfr_div_z (y, tmp, mpq_denref(z), rnd_mode);
 
116
      mpfr_clear (tmp);
 
117
      return res;
 
118
    }
 
119
}
 
120
 
 
121
int
 
122
mpfr_div_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode)
 
123
{
 
124
  mpfr_t tmp;
 
125
  int res;
 
126
  mp_prec_t p;
 
127
 
 
128
  if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
 
129
    return mpfr_div_ui (y, x, 0, rnd_mode);
 
130
  else if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z)) == 0))
 
131
    p = 0;
 
132
  else
 
133
    MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (z));
 
134
  mpfr_init2 (tmp, MPFR_PREC(x) + p);
 
135
  res = mpfr_mul_z (tmp, x, mpq_denref(z), GMP_RNDN );
 
136
  MPFR_ASSERTD( res == 0 );
 
137
  res = mpfr_div_z (y, tmp, mpq_numref(z), rnd_mode);
 
138
  mpfr_clear (tmp);
 
139
  return res;
 
140
}
 
141
 
 
142
int
 
143
mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode)
 
144
{
 
145
  mpfr_t    t,q;
 
146
  mp_prec_t p;
 
147
  mp_exp_t  err;
 
148
  int res;
 
149
  MPFR_ZIV_DECL (loop);
 
150
 
 
151
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 
152
    {
 
153
      if (MPFR_IS_NAN (x))
 
154
        {
 
155
          MPFR_SET_NAN (y);
 
156
          MPFR_RET_NAN;
 
157
        }
 
158
      else if (MPFR_IS_INF (x))
 
159
        {
 
160
          MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0);
 
161
          MPFR_SET_INF (y);
 
162
          MPFR_SET_SAME_SIGN (y, x);
 
163
          MPFR_RET (0);
 
164
        }
 
165
      else
 
166
        {
 
167
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
 
168
          if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
 
169
            return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */
 
170
          else
 
171
            return mpfr_set_q (y, z, rnd_mode);
 
172
        }
 
173
    }
 
174
 
 
175
  p = MPFR_PREC (y) + 10;
 
176
  mpfr_init2 (t, p);
 
177
  mpfr_init2 (q, p);
 
178
 
 
179
  MPFR_ZIV_INIT (loop, p);
 
180
  for (;;)
 
181
    {
 
182
      res = mpfr_set_q (q, z, GMP_RNDN);  /* Error <= 1/2 ulp(q) */
 
183
      /* If z if @INF@ (1/0), res = 0, so it quits immediately */
 
184
      if (MPFR_UNLIKELY (res == 0))
 
185
        /* Result is exact so we can add it directly! */
 
186
        {
 
187
          res = mpfr_add (y, x, q, rnd_mode);
 
188
          break;
 
189
        }
 
190
      mpfr_add (t, x, q, GMP_RNDN);       /* Error <= 1/2 ulp(t) */
 
191
      /* Error / ulp(t)      <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t))
 
192
         If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t)))
 
193
                             <= 2^(EXP(q)-EXP(t))
 
194
         If EXP(q)-EXP(t)<0, <= 2^0 */
 
195
      /* We can get 0, but we can't round since q is inexact */
 
196
      if (MPFR_LIKELY (!MPFR_IS_ZERO (t)))
 
197
        {
 
198
          err = (mp_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0);
 
199
          if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode)))
 
200
            {
 
201
              res = mpfr_set (y, t, rnd_mode);
 
202
              break;
 
203
            }
 
204
        }
 
205
      MPFR_ZIV_NEXT (loop, p);
 
206
      mpfr_set_prec (t, p);
 
207
      mpfr_set_prec (q, p);
 
208
    }
 
209
  MPFR_ZIV_FREE (loop);
 
210
  mpfr_clear (t);
 
211
  mpfr_clear (q);
 
212
  return res;
 
213
}
 
214
 
 
215
int
 
216
mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode)
 
217
{
 
218
  mpfr_t t,q;
 
219
  mp_prec_t p;
 
220
  int res;
 
221
  mp_exp_t err;
 
222
  MPFR_ZIV_DECL (loop);
 
223
 
 
224
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 
225
    {
 
226
      if (MPFR_IS_NAN (x))
 
227
        {
 
228
          MPFR_SET_NAN (y);
 
229
          MPFR_RET_NAN;
 
230
        }
 
231
      else if (MPFR_IS_INF (x))
 
232
        {
 
233
          MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0);
 
234
          MPFR_SET_INF (y);
 
235
          MPFR_SET_SAME_SIGN (y, x);
 
236
          MPFR_RET (0);
 
237
        }
 
238
      else
 
239
        {
 
240
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
 
241
 
 
242
          if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
 
243
            return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */
 
244
          else
 
245
            {
 
246
              res =  mpfr_set_q (y, z, MPFR_INVERT_RND (rnd_mode));
 
247
              MPFR_CHANGE_SIGN (y);
 
248
              return -res;
 
249
            }
 
250
        }
 
251
    }
 
252
 
 
253
  p = MPFR_PREC (y) + 10;
 
254
  mpfr_init2 (t, p);
 
255
  mpfr_init2 (q, p);
 
256
 
 
257
  MPFR_ZIV_INIT (loop, p);
 
258
  for(;;)
 
259
    {
 
260
      res = mpfr_set_q(q, z, GMP_RNDN);  /* Error <= 1/2 ulp(q) */
 
261
      /* If z if @INF@ (1/0), res = 0, so it quits immediately */
 
262
      if (MPFR_UNLIKELY (res == 0))
 
263
        /* Result is exact so we can add it directly!*/
 
264
        {
 
265
          res = mpfr_sub (y, x, q, rnd_mode);
 
266
          break;
 
267
        }
 
268
      mpfr_sub (t, x, q, GMP_RNDN);       /* Error <= 1/2 ulp(t) */
 
269
      /* Error / ulp(t)      <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t))
 
270
         If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t)))
 
271
                             <= 2^(EXP(q)-EXP(t))
 
272
                             If EXP(q)-EXP(t)<0, <= 2^0 */
 
273
      /* We can get 0, but we can't round since q is inexact */
 
274
      if (MPFR_LIKELY (!MPFR_IS_ZERO (t)))
 
275
        {
 
276
          err = (mp_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0);
 
277
          res = MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode);
 
278
          if (MPFR_LIKELY (res != 0))  /* We can round! */
 
279
            {
 
280
              res = mpfr_set (y, t, rnd_mode);
 
281
              break;
 
282
            }
 
283
        }
 
284
      MPFR_ZIV_NEXT (loop, p);
 
285
      mpfr_set_prec (t, p);
 
286
      mpfr_set_prec (q, p);
 
287
    }
 
288
  MPFR_ZIV_FREE (loop);
 
289
  mpfr_clear (t);
 
290
  mpfr_clear (q);
 
291
  return res;
 
292
}
 
293
 
 
294
int
 
295
mpfr_cmp_q (mpfr_srcptr x, mpq_srcptr z)
 
296
{
 
297
  mpfr_t t;
 
298
  int res;
 
299
  mp_prec_t p;
 
300
  /* x < a/b ? <=> x*b < a */
 
301
  MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0);
 
302
  MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (z));
 
303
  mpfr_init2 (t, MPFR_PREC(x) + p);
 
304
  res = mpfr_mul_z (t, x, mpq_denref (z), GMP_RNDN );
 
305
  MPFR_ASSERTD (res == 0);
 
306
  res = mpfr_cmp_z (t, mpq_numref (z) );
 
307
  mpfr_clear (t);
 
308
  return res;
 
309
}
 
310
 
 
311
int
 
312
mpfr_cmp_f (mpfr_srcptr x, mpf_srcptr z)
 
313
{
 
314
  mpfr_t t;
 
315
  int res;
 
316
 
 
317
  mpfr_init2 (t, MPFR_PREC_MIN + ABS(SIZ(z)) * BITS_PER_MP_LIMB );
 
318
  res = mpfr_set_f (t, z, GMP_RNDN);
 
319
  MPFR_ASSERTD (res == 0);
 
320
  res = mpfr_cmp (x, t);
 
321
  mpfr_clear (t);
 
322
  return res;
 
323
}