← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/gmp/mpfr/agm.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/agm.c
1
 
/* mpfr_agm -- arithmetic-geometric mean of two floating-point numbers
2
 
 
3
 
Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation.
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
 
void 
28
 
mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode)
29
 
{
30
 
  int s, go_on, compare;
31
 
  mp_prec_t p, q;
32
 
  double uo, vo;
33
 
  mp_limb_t *up, *vp, *tmpp, *tmpup, *tmpvp, *ap, *bp;
34
 
  mpfr_t u, v, tmp, tmpu, tmpv, a, b;
35
 
  TMP_DECL(marker);
36
 
 
37
 
  /* If a or b is NaN, the result is NaN */
38
 
  if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2))
39
 
    {
40
 
      MPFR_SET_NAN(r);
41
 
      __mpfr_flags |= MPFR_FLAGS_NAN;
42
 
      return;
43
 
    }
44
 
 
45
 
  /* If a or b is negative (including -Infinity), the result is NaN */
46
 
  if ((MPFR_SIGN(op1) < 0) || (MPFR_SIGN(op2) < 0))
47
 
    {
48
 
      MPFR_SET_NAN(r);
49
 
      __mpfr_flags |= MPFR_FLAGS_NAN;
50
 
      return;
51
 
    }
52
 
 
53
 
  MPFR_CLEAR_NAN(r);
54
 
 
55
 
  /* If a or b is +Infinity, the result is +Infinity */
56
 
  if (MPFR_IS_INF(op1) || MPFR_IS_INF(op2))
57
 
    {
58
 
      MPFR_SET_INF(r);
59
 
      MPFR_SET_SAME_SIGN(r, op1);
60
 
      return;
61
 
    }
62
 
 
63
 
  MPFR_CLEAR_INF(r);
64
 
  
65
 
  /* If a or b is 0, the result is 0 */
66
 
  if ((MPFR_NOTZERO(op1) && MPFR_NOTZERO(op2)) == 0)
67
 
    {
68
 
      MPFR_SET_ZERO(r);
69
 
      return;
70
 
    }
71
 
 
72
 
 /* precision of the following calculus */
73
 
  q = MPFR_PREC(r);
74
 
  p = q + 15;
75
 
 
76
 
  /* Initialisations */
77
 
  go_on=1;
78
 
 
79
 
  TMP_MARK(marker);
80
 
  s=(p-1)/BITS_PER_MP_LIMB+1;
81
 
  MPFR_INIT(ap, a, p, s);  
82
 
  MPFR_INIT(bp, b, p, s);
83
 
  MPFR_INIT(up, u, p, s);
84
 
  MPFR_INIT(vp, v, p, s);   
85
 
  MPFR_INIT(tmpup, tmpu, p, s);  
86
 
  MPFR_INIT(tmpvp, tmpv, p, s);  
87
 
  MPFR_INIT(tmpp, tmp, p, s);  
88
 
 
89
 
 
90
 
 
91
 
  /* b and a are the 2 operands but we want b >= a */
92
 
  if ((compare = mpfr_cmp (op1,op2)) > 0)
93
 
    {
94
 
      mpfr_set (b,op1,GMP_RNDN);
95
 
      mpfr_set (a,op2,GMP_RNDN);  
96
 
    }
97
 
  else if (compare == 0)
98
 
    {
99
 
      mpfr_set (r, op1, rnd_mode);
100
 
      TMP_FREE(marker);
101
 
      return;
102
 
    }
103
 
  else
104
 
    {
105
 
      mpfr_set (b,op2,GMP_RNDN);
106
 
      mpfr_set (a,op1,GMP_RNDN);  
107
 
    }
108
 
    
109
 
  vo = mpfr_get_d1 (b);
110
 
  uo = mpfr_get_d1 (a);
111
 
 
112
 
  mpfr_set(u,a,GMP_RNDN);
113
 
  mpfr_set(v,b,GMP_RNDN);
114
 
 
115
 
  /* Main loop */
116
 
 
117
 
  while (go_on) {
118
 
    int err, can_round;
119
 
    mp_prec_t eq;
120
 
    double erraux;
121
 
    
122
 
    erraux = (vo == uo) ? 0.0 : _mpfr_ceil_exp2 (-2.0 * (double) p * uo
123
 
                                                   / (vo - uo));
124
 
    err = 1 + (int) ((3.0 / 2.0 * (double) _mpfr_ceil_log2 ((double) p)
125
 
                      + 1.0) * _mpfr_ceil_exp2 (- (double) p)
126
 
                     + 3.0 * erraux);
127
 
    if(p-err-3<=q) {
128
 
      p=q+err+4;
129
 
      err= 1 + 
130
 
        (int) ((3.0/2.0*_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p)
131
 
               +3.0 * erraux);
132
 
    }
133
 
 
134
 
    /* Calculus of un and vn */
135
 
    do
136
 
      {
137
 
        mpfr_mul(tmp, u, v, GMP_RNDN);
138
 
        mpfr_sqrt (tmpu, tmp, GMP_RNDN); 
139
 
        mpfr_add(tmp, u, v, GMP_RNDN);
140
 
        mpfr_div_2ui(tmpv, tmp, 1, GMP_RNDN);
141
 
        mpfr_set(u, tmpu, GMP_RNDN);
142
 
        mpfr_set(v, tmpv, GMP_RNDN);
143
 
      }
144
 
    while (mpfr_cmp2(u, v, &eq) != 0 && eq <= p - 2);
145
 
 
146
 
    /* Roundability of the result */
147
 
      can_round=mpfr_can_round(v,p-err-3,GMP_RNDN,rnd_mode,q);
148
 
    
149
 
      if (can_round)
150
 
        go_on=0;
151
 
 
152
 
      else {
153
 
          go_on=1;
154
 
          p+=5;
155
 
          s=(p-1)/BITS_PER_MP_LIMB+1;
156
 
          MPFR_INIT(up, u, p, s);
157
 
          MPFR_INIT(vp, v, p, s);   
158
 
          MPFR_INIT(tmpup, tmpu, p, s);  
159
 
          MPFR_INIT(tmpvp, tmpv, p, s);  
160
 
          MPFR_INIT(tmpp, tmp, p, s);
161
 
          mpfr_set(u,a,GMP_RNDN);
162
 
          mpfr_set(v,b,GMP_RNDN);
163
 
      }
164
 
  }
165
 
  /* End of while */
166
 
 
167
 
  /* Setting of the result */
168
 
 
169
 
  mpfr_set(r,v,rnd_mode);
170
 
 
171
 
  /* Let's clean */
172
 
  TMP_FREE(marker); 
173
 
 
174
 
  return ;
175
 
}