← Back to branch summary

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

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2006-05-17 02:46:26 UTC
  • Revision ID: james.westby@ubuntu.com-20060517024626-lljr08ftv9g9vefl
Tags: upstream-0.9h-20060510
Import upstream version 0.9h-20060510

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gmp/mpfr/exp3.c
 
1
/* mpfr_exp -- exponential of a floating-point number
 
2
 
 
3
Copyright 1999, 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 <stdio.h>
 
23
#include "gmp.h"
 
24
#include "gmp-impl.h"
 
25
#include "mpfr.h"
 
26
#include "mpfr-impl.h"
 
27
 
 
28
static int mpfr_exp_rational _PROTO ((mpfr_ptr, mpz_srcptr, int, int));
 
29
 
 
30
static int
 
31
mpfr_exp_rational (mpfr_ptr y, mpz_srcptr p, int r, int m)
 
32
{
 
33
  int n,i,k,j,l;
 
34
  mpz_t* P,*S;
 
35
  mpz_t* ptoj;
 
36
  int diff,expo;
 
37
  int precy = MPFR_PREC(y);
 
38
  int * mult;
 
39
  int prec_i_have;
 
40
  int *nb_terms;
 
41
  int accu;
 
42
  TMP_DECL (marker);
 
43
 
 
44
  TMP_MARK (marker);
 
45
  n = 1 << m;
 
46
  P = (mpz_t*) TMP_ALLOC((m+1) * sizeof(mpz_t));
 
47
  S = (mpz_t*) TMP_ALLOC((m+1) * sizeof(mpz_t));
 
48
  ptoj = (mpz_t*) TMP_ALLOC((m+1) * sizeof(mpz_t)); /* ptoj[i] = mantissa^(2^i) */
 
49
  mult = (int*) TMP_ALLOC((m+1) * sizeof(int)); 
 
50
  nb_terms = (int*) TMP_ALLOC((m+1) * sizeof(int)); 
 
51
  mult[0] = 0;
 
52
  for (i=0;i<=m;i++) { mpz_init(P[i]); mpz_init(S[i]); mpz_init(ptoj[i]); }
 
53
  mpz_set(ptoj[0], p);
 
54
  for (i=1;i<m;i++) mpz_mul(ptoj[i], ptoj[i-1], ptoj[i-1]);
 
55
  mpz_set_ui(P[0], 1);
 
56
  mpz_set_ui(S[0], 1);
 
57
  k = 0;
 
58
  nb_terms[0] = 1;
 
59
   prec_i_have = 0; 
 
60
   for (i=1;(prec_i_have < precy) && (i < n) ;i++) {
 
61
    k++;
 
62
    nb_terms[k] = 1;
 
63
    mpz_set_ui(P[k], i+1);
 
64
    mpz_set(S[k], P[k]);;
 
65
    j=i+1; l=0; while ((j & 1) == 0) {      
 
66
      mpz_mul(S[k], S[k], ptoj[l]);
 
67
      mpz_mul(S[k-1], S[k-1], P[k]);
 
68
      mpz_mul_2exp(S[k-1], S[k-1], r*(1<<l));
 
69
      mpz_add(S[k-1], S[k-1], S[k]);
 
70
      mpz_mul(P[k-1], P[k-1], P[k]);
 
71
      nb_terms[k-1] = nb_terms[k-1]+ nb_terms[k];
 
72
      mult[k] = mult[k-1] + (1 << l)*(r >> 2) + mpz_sizeinbase(P[k],2) - 1;
 
73
      prec_i_have = mult[k];
 
74
      l++; j>>=1; k--;
 
75
    }
 
76
  }
 
77
   l = 0;
 
78
   accu = 0;
 
79
   while (k > 0){
 
80
     mpz_mul(S[k], S[k], ptoj[_mpfr_ceil_log2((double) nb_terms[k])]);
 
81
     mpz_mul(S[k-1], S[k-1], P[k]);
 
82
     accu += nb_terms[k];
 
83
     mpz_mul_2exp(S[k-1], S[k-1], r* accu);
 
84
     mpz_add(S[k-1], S[k-1], S[k]);
 
85
     mpz_mul(P[k-1], P[k-1], P[k]);     
 
86
     l++; k--;
 
87
   }
 
88
   
 
89
  diff = mpz_sizeinbase(S[0],2) - 2*precy;
 
90
  expo = diff;
 
91
  if (diff >=0)
 
92
    {
 
93
      mpz_div_2exp(S[0],S[0],diff);
 
94
    } else 
 
95
      {
 
96
        mpz_mul_2exp(S[0],S[0],-diff);
 
97
      }
 
98
  diff = mpz_sizeinbase(P[0],2) - precy;
 
99
  expo -= diff;
 
100
  if (diff >=0)
 
101
    {
 
102
      mpz_div_2exp(P[0],P[0],diff);
 
103
    } else
 
104
      {
 
105
        mpz_mul_2exp(P[0],P[0],-diff);
 
106
        }
 
107
 
 
108
  mpz_tdiv_q(S[0], S[0], P[0]);
 
109
  mpfr_set_z(y,S[0], GMP_RNDD);
 
110
  MPFR_EXP(y) += expo;
 
111
 
 
112
  mpfr_div_2ui(y, y, r*(i-1),GMP_RNDN); 
 
113
  for (i=0;i<=m;i++) { mpz_clear(P[i]); mpz_clear(S[i]); mpz_clear(ptoj[i]); }
 
114
  TMP_FREE (marker);
 
115
  return 0;
 
116
}
 
117
 
 
118
#define shift (BITS_PER_MP_LIMB/2)
 
119
 
 
120
int
 
121
mpfr_exp3 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
 
122
{
 
123
  mpfr_t t;
 
124
  mpfr_t x_copy;
 
125
  int i,k;
 
126
  mpz_t uk;
 
127
  mpfr_t tmp;
 
128
  int ttt;
 
129
  int twopoweri;
 
130
  int Prec;
 
131
  int loop;
 
132
  int prec_x;
 
133
  int shift_x = 0;
 
134
  int good = 0;
 
135
  int realprec = 0;
 
136
  int iter;
 
137
  int logn, inexact = 0;
 
138
 
 
139
  /* decompose x */
 
140
  /* we first write x = 1.xxxxxxxxxxxxx
 
141
     ----- k bits -- */
 
142
  prec_x = _mpfr_ceil_log2 ((double) (MPFR_PREC(x)) / BITS_PER_MP_LIMB);
 
143
  if (prec_x < 0) prec_x = 0;
 
144
  logn =  _mpfr_ceil_log2 ((double) prec_x + MPFR_PREC(y));
 
145
  if (logn < 2) logn = 2;
 
146
  ttt = MPFR_EXP(x);
 
147
  mpfr_init2(x_copy,MPFR_PREC(x));
 
148
  mpfr_set(x_copy,x,GMP_RNDD);
 
149
  /* we shift to get a number less than 1 */
 
150
  if (ttt > 0) 
 
151
    {
 
152
      shift_x = ttt;
 
153
      mpfr_div_2ui(x_copy, x, ttt, GMP_RNDN);
 
154
      ttt = MPFR_EXP(x_copy);
 
155
    }
 
156
  realprec = MPFR_PREC(y)+logn;
 
157
  mpz_init (uk);
 
158
  while (!good){      
 
159
    Prec = realprec+shift+2+shift_x;
 
160
    k = _mpfr_ceil_log2 ((double) Prec / BITS_PER_MP_LIMB);
 
161
 
 
162
    /* now we have to extract */
 
163
    mpfr_init2 (t, Prec);
 
164
    mpfr_init2 (tmp, Prec);
 
165
    mpfr_set_ui(tmp,1,GMP_RNDN);
 
166
    twopoweri = BITS_PER_MP_LIMB;
 
167
    if (k <= prec_x) iter = k; else iter= prec_x;
 
168
    for(i = 0; i <= iter; i++){
 
169
      mpfr_extract (uk, x_copy, i);
 
170
        if (i)
 
171
            mpfr_exp_rational (t, uk, twopoweri - ttt, k  - i + 1);
 
172
        else
 
173
          {
 
174
            /* particular case: we have to compute with x/2^., then
 
175
               do squarings (this is faster) */    
 
176
              mpfr_exp_rational (t, uk, shift + twopoweri - ttt, k+1);
 
177
            for (loop= 0 ; loop < shift; loop++)
 
178
              mpfr_mul(t,t,t,GMP_RNDD);
 
179
 
 
180
          }
 
181
        mpfr_mul(tmp,tmp,t,GMP_RNDD); 
 
182
        twopoweri <<= 1;
 
183
    }
 
184
      mpfr_clear (t);
 
185
      for (loop= 0 ; loop < shift_x; loop++)
 
186
        mpfr_mul (tmp, tmp, tmp, GMP_RNDD);
 
187
      if (mpfr_can_round (tmp, realprec, GMP_RNDD, rnd_mode, MPFR_PREC(y)))
 
188
        {
 
189
          inexact = mpfr_set (y, tmp, rnd_mode);
 
190
          good = 1;
 
191
        }
 
192
      else
 
193
        realprec += 3*logn;
 
194
      mpfr_clear (tmp);
 
195
  }
 
196
  mpz_clear (uk);
 
197
  mpfr_clear(x_copy);
 
198
  return inexact;
 
199
}