← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/gmp/mpfr/mul.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2007-04-09 11:51:51 UTC
  • mfrom: (1.1.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20070409115151-ql8cr0kalzx1jmla
Tags: 0.9i-20070324-2
Upload to unstable. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gmp/mpfr/mul.c
1
 
/* mpfr_mul -- multiply two floating-point numbers
2
 
 
3
 
Copyright 1999, 2000, 2001, 2002 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., 59 Temple Place - Suite 330, Boston,
20
 
MA 02111-1307, USA. */
21
 
 
22
 
#include "gmp.h"
23
 
#include "gmp-impl.h"
24
 
#include "mpfr.h"
25
 
#include "mpfr-impl.h"
26
 
 
27
 
int
28
 
mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) 
29
 
{
30
 
  int sign_product, cc, inexact, ec, em = 0;
31
 
  mp_exp_t bx, cx;
32
 
  mp_limb_t *ap, *bp, *cp, *tmp;
33
 
  mp_limb_t b1;
34
 
  mp_prec_t aq, bq, cq;
35
 
  mp_size_t an, bn, cn, tn, k;
36
 
  TMP_DECL(marker);
37
 
 
38
 
  /* deal with NaN and zero */
39
 
  if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c))
40
 
    {
41
 
      MPFR_SET_NAN(a);
42
 
      MPFR_RET_NAN;
43
 
    }
44
 
 
45
 
  MPFR_CLEAR_NAN(a);
46
 
 
47
 
  sign_product = MPFR_SIGN(b) * MPFR_SIGN(c);
48
 
 
49
 
  if (MPFR_IS_INF(b))
50
 
    {
51
 
      if (MPFR_IS_INF(c) || MPFR_NOTZERO(c))
52
 
        {
53
 
          if (MPFR_SIGN(a) != sign_product)
54
 
            MPFR_CHANGE_SIGN(a);
55
 
          MPFR_SET_INF(a);
56
 
          MPFR_RET(0); /* exact */
57
 
        }
58
 
      else
59
 
        {
60
 
          MPFR_SET_NAN(a);
61
 
          MPFR_RET_NAN;
62
 
        }
63
 
    }
64
 
  else if (MPFR_IS_INF(c))
65
 
    {
66
 
      if (MPFR_NOTZERO(b))
67
 
        {
68
 
          if (MPFR_SIGN(a) != sign_product)
69
 
            MPFR_CHANGE_SIGN(a);
70
 
          MPFR_SET_INF(a);
71
 
          MPFR_RET(0); /* exact */
72
 
        }
73
 
      else
74
 
        {
75
 
          MPFR_SET_NAN(a);
76
 
          MPFR_RET_NAN;
77
 
        }
78
 
    }
79
 
 
80
 
  MPFR_ASSERTN(MPFR_IS_FP(b) && MPFR_IS_FP(c));
81
 
  MPFR_CLEAR_INF(a); /* clear Inf flag */
82
 
 
83
 
  if (MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c))
84
 
    {
85
 
      if (MPFR_SIGN(a) != sign_product)
86
 
        MPFR_CHANGE_SIGN(a);
87
 
      MPFR_SET_ZERO(a);
88
 
      MPFR_RET(0); /* 0 * 0 is exact */
89
 
    }
90
 
 
91
 
  bx = MPFR_EXP(b);
92
 
  cx = MPFR_EXP(c);
93
 
  /* Note: exponent of the result will be bx + cx + ec with ec in {-1,0,1} */
94
 
  if (bx >= 0 && cx > 0)
95
 
    { /* bx + cx > 0 */
96
 
      if (__mpfr_emax < 0 ||
97
 
          (mp_exp_unsigned_t) bx + cx > (mp_exp_unsigned_t) __mpfr_emax + 1)
98
 
        return mpfr_set_overflow(a, rnd_mode, sign_product);
99
 
 
100
 
      if ((mp_exp_unsigned_t) bx + cx == (mp_exp_unsigned_t) __mpfr_emax + 1)
101
 
        em = 1;
102
 
    }
103
 
  else if (bx <= 0 && cx < 0)
104
 
    { /* bx + cx < 0 */
105
 
      if (__mpfr_emin > 0 ||
106
 
          (mp_exp_unsigned_t) bx + cx < (mp_exp_unsigned_t) __mpfr_emin - 1)
107
 
        return mpfr_set_underflow(a, rnd_mode, sign_product);
108
 
 
109
 
      if ((mp_exp_unsigned_t) bx + cx == (mp_exp_unsigned_t) __mpfr_emin - 1)
110
 
        em = -1;
111
 
    }
112
 
  else
113
 
    { /* bx != 0 and cx doesn't have the same sign */
114
 
      if ((bx + cx) - 1 > __mpfr_emax)
115
 
        return mpfr_set_overflow(a, rnd_mode, sign_product);
116
 
 
117
 
      if ((bx + cx) - 1 == __mpfr_emax)
118
 
        em = 1;
119
 
 
120
 
      if ((bx + cx) + 1 < __mpfr_emin)
121
 
        return mpfr_set_underflow(a, rnd_mode, sign_product);
122
 
 
123
 
      if ((bx + cx) + 1 == __mpfr_emin)
124
 
        em = -1;
125
 
    }
126
 
 
127
 
  ap = MPFR_MANT(a);
128
 
  bp = MPFR_MANT(b);
129
 
  cp = MPFR_MANT(c);
130
 
 
131
 
  aq = MPFR_PREC(a);
132
 
  bq = MPFR_PREC(b);
133
 
  cq = MPFR_PREC(c);
134
 
 
135
 
  an = (aq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of a */
136
 
  bn = (bq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */
137
 
  cn = (cq-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of c */
138
 
 
139
 
  MPFR_ASSERTN((mp_size_unsigned_t) bn + cn <= MP_SIZE_T_MAX);
140
 
  k = bn + cn; /* effective nb of limbs used by b*c (= tn or tn+1) below */
141
 
 
142
 
  MPFR_ASSERTN(bq + cq >= bq); /* no integer overflow */
143
 
  tn = (bq + cq - 1) / BITS_PER_MP_LIMB + 1; /* <= k, thus no int overflow */
144
 
 
145
 
  MPFR_ASSERTN(k <= ((size_t) -1) / BYTES_PER_MP_LIMB);
146
 
  TMP_MARK(marker); 
147
 
  tmp = (mp_limb_t *) TMP_ALLOC((size_t) k * BYTES_PER_MP_LIMB);
148
 
 
149
 
  /* multiplies two mantissa in temporary allocated space */
150
 
  b1 = (bn >= cn) ? mpn_mul (tmp, bp, bn, cp, cn)
151
 
    : mpn_mul (tmp, cp, cn, bp, bn);
152
 
 
153
 
  /* now tmp[0]..tmp[k-1] contains the product of both mantissa,
154
 
     with tmp[k-1]>=2^(BITS_PER_MP_LIMB-2) */
155
 
  b1 >>= BITS_PER_MP_LIMB - 1; /* msb from the product */
156
 
 
157
 
  tmp += k - tn;
158
 
  if (b1 == 0)
159
 
    mpn_lshift (tmp, tmp, tn, 1);
160
 
  cc = mpfr_round_raw (ap, tmp, bq + cq, sign_product < 0, aq,
161
 
                       rnd_mode, &inexact);
162
 
  if (cc) /* cc = 1 ==> result is a power of two */
163
 
    ap[an-1] = GMP_LIMB_HIGHBIT;
164
 
 
165
 
  TMP_FREE(marker);
166
 
 
167
 
  ec = b1 - 1 + cc;
168
 
 
169
 
  if (em == 0)
170
 
    {
171
 
      mp_exp_t ax = bx + cx;
172
 
 
173
 
      if (ax == __mpfr_emax && ec > 0)
174
 
        return mpfr_set_overflow(a, rnd_mode, sign_product);
175
 
 
176
 
      if (ax == __mpfr_emin && ec < 0)
177
 
        return mpfr_set_underflow(a, rnd_mode, sign_product);
178
 
 
179
 
      MPFR_EXP(a) = ax + ec;
180
 
    }
181
 
  else if (em > 0)
182
 
    {
183
 
      if (ec >= 0)
184
 
        return mpfr_set_overflow(a, rnd_mode, sign_product);
185
 
 
186
 
      MPFR_EXP(a) = __mpfr_emax;
187
 
    }
188
 
  else
189
 
    {
190
 
      if (ec <= 0)
191
 
        return mpfr_set_underflow(a, rnd_mode, sign_product);
192
 
 
193
 
      MPFR_EXP(a) = __mpfr_emin;
194
 
    }
195
 
 
196
 
  if (MPFR_SIGN(a) != sign_product)
197
 
    MPFR_CHANGE_SIGN(a);
198
 
 
199
 
  return inexact;
200
 
}