← Back to branch summary

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

« back to all changes in this revision

Viewing changes to mpfr/root.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/root.c
 
1
/* mpfr_root -- kth root.
 
2
 
 
3
Copyright 2005 Free Software Foundation.
 
4
Contributed by the Spaces project, INRIA Lorraine.
 
5
 
 
6
This file is part of the MPFR Library.
 
7
 
 
8
The MPFR 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 MPFR 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 MPFR Library; see the file COPYING.LIB.  If not, write to
 
20
the Free Software Foundation, Inc., 51 Franklin Place, Fifth Floor, Boston,
 
21
MA 02110-1301, USA. */
 
22
 
 
23
#define MPFR_NEED_LONGLONG_H
 
24
#include "mpfr-impl.h"
 
25
 
 
26
 /* The computation of y = x^(1/k) is done as follows:
 
27
 
 
28
    Let x = sign * m * 2^(k*e) where m is an integer
 
29
 
 
30
    with 2^(k*(n-1)) <= m < 2^(k*n) where n = PREC(y)
 
31
 
 
32
    and m = s^k + r where 0 <= r and m < (s+1)^k
 
33
 
 
34
    we want that s has n bits i.e. s >= 2^(n-1), or m >= 2^(k*(n-1))
 
35
    i.e. m must have at least k*(n-1)+1 bits
 
36
 
 
37
    then, not taking into account the sign, the result will be
 
38
    x^(1/k) = s * 2^e or (s+1) * 2^e according to the rounding mode.
 
39
 */
 
40
 
 
41
int
 
42
mpfr_root (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mp_rnd_t rnd_mode)
 
43
{
 
44
  mpz_t m;
 
45
  mp_exp_t e, r, sh;
 
46
  mp_prec_t n, size_m, tmp;
 
47
  int inexact, negative;
 
48
  MPFR_SAVE_EXPO_DECL (expo);
 
49
 
 
50
  if (MPFR_UNLIKELY (k <= 1))
 
51
    {
 
52
      if (k < 1) /* k==0 => y=x^(1/0)=x^(+Inf) */
 
53
#if 0
 
54
        /* For 0 <= x < 1 => +0.
 
55
           For x = 1      => 1.
 
56
           For x > 1,     => +Inf.
 
57
           For x < 0      => NaN.
 
58
        */
 
59
        {
 
60
          if (MPFR_IS_NEG (x) && !MPFR_IS_ZERO (x))
 
61
            {
 
62
              MPFR_SET_NAN (y);
 
63
              MPFR_RET_NAN;
 
64
            }
 
65
          inexact = mpfr_cmp (x, __gmpfr_one);
 
66
          if (inexact == 0)
 
67
            return mpfr_set_ui (y, 1, rnd_mode); /* 1 may be Out of Range */
 
68
          else if (inexact < 0)
 
69
            return mpfr_set_ui (y, 0, rnd_mode); /* 0+ */
 
70
          else
 
71
            {
 
72
              mpfr_set_inf (y, 1);
 
73
              return 0;
 
74
            }
 
75
        }
 
76
#endif
 
77
      {
 
78
        MPFR_SET_NAN (y);
 
79
        MPFR_RET_NAN;
 
80
      }
 
81
      else /* y =x^(1/1)=x */
 
82
        return mpfr_set (y, x, rnd_mode);
 
83
    }
 
84
 
 
85
  /* Singular values */
 
86
  else if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 
87
    {
 
88
      if (MPFR_IS_NAN (x))
 
89
        {
 
90
          MPFR_SET_NAN (y); /* NaN^(1/k) = NaN */
 
91
          MPFR_RET_NAN;
 
92
        }
 
93
      else if (MPFR_IS_INF (x)) /* +Inf^(1/k) = +Inf
 
94
                                   -Inf^(1/k) = -Inf if k odd
 
95
                                   -Inf^(1/k) = NaN if k even */
 
96
        {
 
97
          if (MPFR_IS_NEG(x) && (k % 2 == 0))
 
98
            {
 
99
              MPFR_SET_NAN (y);
 
100
              MPFR_RET_NAN;
 
101
            }
 
102
          MPFR_SET_INF (y);
 
103
          MPFR_SET_SAME_SIGN (y, x);
 
104
          MPFR_RET (0);
 
105
        }
 
106
      else /* x is necessarily 0: (+0)^(1/k) = +0
 
107
                                  (-0)^(1/k) = -0 */
 
108
        {
 
109
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
 
110
          MPFR_SET_ZERO (y);
 
111
          MPFR_SET_SAME_SIGN (y, x);
 
112
          MPFR_RET (0);
 
113
        }
 
114
    }
 
115
 
 
116
  /* Returns NAN for x < 0 and k even */
 
117
  else if (MPFR_IS_NEG (x) && (k % 2 == 0))
 
118
    {
 
119
      MPFR_SET_NAN (y);
 
120
      MPFR_RET_NAN;
 
121
    }
 
122
 
 
123
  /* General case */
 
124
  MPFR_SAVE_EXPO_MARK (expo);
 
125
  mpz_init (m);
 
126
 
 
127
  e = mpfr_get_z_exp (m, x);                /* x = m * 2^e */
 
128
  if ((negative = MPFR_IS_NEG(x)))
 
129
    mpz_neg (m, m);
 
130
  r = e % (mp_exp_t) k;
 
131
  if (r < 0)
 
132
    r += k; /* now r = e (mod k) with 0 <= e < r */
 
133
  /* x = (m*2^r) * 2^(e-r) where e-r is a multiple of k */
 
134
 
 
135
  MPFR_MPZ_SIZEINBASE2 (size_m, m);
 
136
  /* for rounding to nearest, we want the round bit to be in the root */
 
137
  n = MPFR_PREC (y) + (rnd_mode == GMP_RNDN);
 
138
 
 
139
  /* we now multiply m by 2^(r+k*sh) so that root(m,k) will give
 
140
     exactly n bits: we want k*(n-1)+1 <= size_m + k*sh + r <= k*n
 
141
     i.e. sh = floor ((kn-size_m-r)/k) */
 
142
  if ((mp_exp_t) size_m + r > k * (mp_exp_t) n)
 
143
    sh = 0; /* we already have too many bits */
 
144
  else
 
145
    sh = (k * (mp_exp_t) n - (mp_exp_t) size_m - r) / k;
 
146
  sh = k * sh + r;
 
147
  if (sh >= 0)
 
148
    {
 
149
      mpz_mul_2exp (m, m, sh);
 
150
      e = e - sh;
 
151
    }
 
152
  else if (r > 0)
 
153
    {
 
154
      mpz_mul_2exp (m, m, r);
 
155
      e = e - r;
 
156
    }
 
157
 
 
158
  /* invariant: x = m*2^e, with e divisible by k */
 
159
 
 
160
  /* we reuse the variable m to store the cube root, since it is not needed
 
161
     any more: we just need to know if the root is exact */
 
162
  inexact = mpz_root (m, m, k) == 0;
 
163
 
 
164
  MPFR_MPZ_SIZEINBASE2 (tmp, m);
 
165
  sh = tmp - n;
 
166
  if (sh > 0) /* we have to flush to 0 the last sh bits from m */
 
167
    {
 
168
      inexact = inexact || ((mp_exp_t) mpz_scan1 (m, 0) < sh);
 
169
      mpz_div_2exp (m, m, sh);
 
170
      e += k * sh;
 
171
    }
 
172
 
 
173
  if (inexact)
 
174
    {
 
175
      if (negative)
 
176
        rnd_mode = MPFR_INVERT_RND (rnd_mode);
 
177
      if (rnd_mode == GMP_RNDU
 
178
          || (rnd_mode == GMP_RNDN && mpz_tstbit (m, 0)))
 
179
        inexact = 1, mpz_add_ui (m, m, 1);
 
180
      else
 
181
        inexact = -1;
 
182
    }
 
183
 
 
184
  /* either inexact is not zero, and the conversion is exact, i.e. inexact
 
185
     is not changed; or inexact=0, and inexact is set only when
 
186
     rnd_mode=GMP_RNDN and bit (n+1) from m is 1 */
 
187
  inexact += mpfr_set_z (y, m, GMP_RNDN);
 
188
  MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / (mp_exp_t) k);
 
189
 
 
190
  if (negative)
 
191
    {
 
192
      MPFR_CHANGE_SIGN (y);
 
193
      inexact = -inexact;
 
194
    }
 
195
 
 
196
  mpz_clear (m);
 
197
  MPFR_SAVE_EXPO_FREE (expo);
 
198
  return mpfr_check_range (y, inexact, rnd_mode);
 
199
}