contrib/mpfr/agm.c

   1 /* mpfr_agm -- arithmetic-geometric mean of two floating-point numbers
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU 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 GNU 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 GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* agm(x,y) is between x and y, so we don't need to save exponent range */
  27 int
  28 mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode)
  29 {
  30   int compare, inexact;
  31   mp_size_t s;
  32   mp_prec_t p, q;
  33   mp_limb_t *up, *vp, *tmpp;
  34   mpfr_t u, v, tmp;
  35   unsigned long n; /* number of iterations */
  36   unsigned long err = 0;
  37   MPFR_ZIV_DECL (loop);
  38   MPFR_TMP_DECL(marker);
  39
  40   MPFR_LOG_FUNC (("op2[%#R]=%R op1[%#R]=%R rnd=%d", op2,op2,op1,op1,rnd_mode),
  41                  ("r[%#R]=%R inexact=%d", r, r, inexact));
  42
  43   /* Deal with special values */
  44   if (MPFR_ARE_SINGULAR (op1, op2))
  45     {
  46       /* If a or b is NaN, the result is NaN */
  47       if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2))
  48         {
  49           MPFR_SET_NAN(r);
  50           MPFR_RET_NAN;
  51         }
  52       /* now one of a or b is Inf or 0 */
  53       /* If a and b is +Inf, the result is +Inf.
  54          Otherwise if a or b is -Inf or 0, the result is NaN */
  55       else if (MPFR_IS_INF(op1) || MPFR_IS_INF(op2))
  56         {
  57           if (MPFR_IS_STRICTPOS(op1) && MPFR_IS_STRICTPOS(op2))
  58             {
  59               MPFR_SET_INF(r);
  60               MPFR_SET_SAME_SIGN(r, op1);
  61               MPFR_RET(0); /* exact */
  62             }
  63           else
  64             {
  65               MPFR_SET_NAN(r);
  66               MPFR_RET_NAN;
  67             }
  68         }
  69       else /* a and b are neither NaN nor Inf, and one is zero */
  70         {  /* If a or b is 0, the result is +0 since a sqrt is positive */
  71           MPFR_ASSERTD (MPFR_IS_ZERO (op1) || MPFR_IS_ZERO (op2));
  72           MPFR_SET_POS (r);
  73           MPFR_SET_ZERO (r);
  74           MPFR_RET (0); /* exact */
  75         }
  76     }
  77   MPFR_CLEAR_FLAGS (r);
  78
  79   /* If a or b is negative (excluding -Infinity), the result is NaN */
  80   if (MPFR_UNLIKELY(MPFR_IS_NEG(op1) || MPFR_IS_NEG(op2)))
  81     {
  82       MPFR_SET_NAN(r);
  83       MPFR_RET_NAN;
  84     }
  85
  86   /* Precision of the following calculus */
  87   q = MPFR_PREC(r);
  88   p = q + MPFR_INT_CEIL_LOG2(q) + 15;
  89   MPFR_ASSERTD (p >= 7); /* see algorithms.tex */
  90   s = (p - 1) / BITS_PER_MP_LIMB + 1;
  91
  92   /* b (op2) and a (op1) are the 2 operands but we want b >= a */
  93   compare = mpfr_cmp (op1, op2);
  94   if (MPFR_UNLIKELY( compare == 0 ))
  95     {
  96       mpfr_set (r, op1, rnd_mode);
  97       MPFR_RET (0); /* exact */
  98     }
  99   else if (compare > 0)
 100     {
 101       mpfr_srcptr t = op1;
 102       op1 = op2;
 103       op2 = t;
 104     }
 105   /* Now b(=op2) >= a (=op1) */
 106
 107   MPFR_TMP_MARK(marker);
 108
 109   /* Main loop */
 110   MPFR_ZIV_INIT (loop, p);
 111   for (;;)
 112     {
 113       mp_prec_t eq;
 114
 115       /* Init temporary vars */
 116       MPFR_TMP_INIT (up, u, p, s);
 117       MPFR_TMP_INIT (vp, v, p, s);
 118       MPFR_TMP_INIT (tmpp, tmp, p, s);
 119
 120       /* Calculus of un and vn */
 121       mpfr_mul (u, op1, op2, GMP_RNDN); /* Faster since PREC(op) < PREC(u) */
 122       mpfr_sqrt (u, u, GMP_RNDN);
 123       mpfr_add (v, op1, op2, GMP_RNDN); /* add with !=prec is still good*/
 124       mpfr_div_2ui (v, v, 1, GMP_RNDN);
 125       n = 1;
 126       while (mpfr_cmp2 (u, v, &eq) != 0 && eq <= p - 2)
 127         {
 128           mpfr_add (tmp, u, v, GMP_RNDN);
 129           mpfr_div_2ui (tmp, tmp, 1, GMP_RNDN);
 130           /* See proof in algorithms.tex */
 131           if (4*eq > p)
 132             {
 133               mpfr_t w;
 134               /* tmp = U(k) */
 135               mpfr_init2 (w, (p + 1) / 2);
 136               mpfr_sub (w, v, u, GMP_RNDN);         /* e = V(k-1)-U(k-1) */
 137               mpfr_sqr (w, w, GMP_RNDN);            /* e = e^2 */
 138               mpfr_div_2ui (w, w, 4, GMP_RNDN);     /* e*= (1/2)^2*1/4  */
 139               mpfr_div (w, w, tmp, GMP_RNDN);       /* 1/4*e^2/U(k) */
 140               mpfr_sub (v, tmp, w, GMP_RNDN);
 141               err = MPFR_GET_EXP (tmp) - MPFR_GET_EXP (v); /* 0 or 1 */
 142               mpfr_clear (w);
 143               break;
 144             }
 145           mpfr_mul (u, u, v, GMP_RNDN);
 146           mpfr_sqrt (u, u, GMP_RNDN);
 147           mpfr_swap (v, tmp);
 148           n ++;
 149         }
 150       /* the error on v is bounded by (18n+51) ulps, or twice if there
 151          was an exponent loss in the final subtraction */
 152       err += MPFR_INT_CEIL_LOG2(18 * n + 51); /* 18n+51 should not overflow
 153                                                  since n is about log(p) */
 154       /* we should have n+2 <= 2^(p/4) [see algorithms.tex] */
 155       if (MPFR_LIKELY (MPFR_INT_CEIL_LOG2(n + 2) <= p / 4 &&
 156                        MPFR_CAN_ROUND (v, p - err, q, rnd_mode)))
 157         break; /* Stop the loop */
 158
 159       /* Next iteration */
 160       MPFR_ZIV_NEXT (loop, p);
 161       s = (p - 1) / BITS_PER_MP_LIMB + 1;
 162     }
 163   MPFR_ZIV_FREE (loop);
 164
 165   /* Setting of the result */
 166   inexact = mpfr_set (r, v, rnd_mode);
 167
 168   /* Let's clean */
 169   MPFR_TMP_FREE(marker);
 170
 171   return inexact; /* agm(u,v) can be exact for u, v rational only for u=v.
 172                      Proof (due to Nicolas Brisebarre): it suffices to consider
 173                      u=1 and v<1. Then 1/AGM(1,v) = 2F1(1/2,1/2,1;1-v^2),
 174                      and a theorem due to G.V. Chudnovsky states that for x a
 175                      non-zero algebraic number with |x|<1, then
 176                      2F1(1/2,1/2,1;x) and 2F1(-1/2,1/2,1;x) are algebraically
 177                      independent over Q. */
 178 }