OpenXM_contrib/gmp/mpfr/agm.c - diff

Return to agm.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr

Diff for /OpenXM_contrib/gmp/mpfr/Attic/agm.c between version 1.1 and 1.1.1.2

-version 1.1, 2000/09/09 14:12:19
+version 1.1.1.2, 2003/08/25 16:06:06
 Line 1
 Line 1
 Line 1
  /* mpfr_agm -- arithmetic-geometric mean of two floating-point numbers
- Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+ Copyright 1999, 2000, 2001, 2002 Free Software Foundation.
  This file is part of the MPFR Library.
  The MPFR Library is free software; you can redistribute it and/or modify
- it under the terms of the GNU Library General Public License as published by
+ it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation; either version 2 of the License, or (at your
+ the Free Software Foundation; either version 2.1 of the License, or (at your
  option) any later version.
  The MPFR Library is distributed in the hope that it will be useful, but
  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
+ or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  License for more details.
- You should have received a copy of the GNU Library General Public License
+ You should have received a copy of the GNU Lesser General Public License
  along with the MPFR Library; see the file COPYING.LIB.  If not, write to
  the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  MA 02111-1307, USA. */
- #include <stdio.h>
- #include <math.h>
  #include "gmp.h"
  #include "gmp-impl.h"
  #include "mpfr.h"
+ #include "mpfr-impl.h"
- /*Memory gestion */
- #define MON_INIT(xp, x, p, s) xp = (mp_ptr) TMP_ALLOC(s*BYTES_PER_MP_LIMB);    x -> _mp_prec = p; x -> _mp_d = xp; x -> _mp_size = s; x -> _mp_exp = 0;
  void
- #ifdef __STDC__
+ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode)
- mpfr_agm(mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, unsigned char rnd_mode)
- #else
- mpfr_agm(r, a, b, rnd_mode)
-      mpfr_ptr r;
-      mpfr_srcptr a;
-      mpfr_srcptr b;
-      unsigned char rnd_mode;
- #endif
  {
-   int s, p, q, go_on;
+   int s, go_on, compare;
+   mp_prec_t p, q;
    double uo, vo;
    mp_limb_t *up, *vp, *tmpp, *tmpup, *tmpvp, *ap, *bp;
    mpfr_t u, v, tmp, tmpu, tmpv, a, b;
-   TMP_DECL(marker1);
+   TMP_DECL(marker);
-   TMP_DECL(marker2);
    /* If a or b is NaN, the result is NaN */
-   if (FLAG_NAN(op1) || FLAG_NAN(op2))
+   if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2))
-     { SET_NAN(r); return; }
+     {
+       MPFR_SET_NAN(r);
+       __mpfr_flags |= MPFR_FLAGS_NAN;
+       return;
+     }
+   /* If a or b is negative (including -Infinity), the result is NaN */
+   if ((MPFR_SIGN(op1) < 0) || (MPFR_SIGN(op2) < 0))
+     {
+       MPFR_SET_NAN(r);
+       __mpfr_flags |= MPFR_FLAGS_NAN;
+       return;
+     }
-   /* If a or b is negative, the result is NaN */
+   MPFR_CLEAR_NAN(r);
-   if ((SIGN(op1)<0)||(SIGN(op2)<0))
-     { SET_NAN(r); return; }
+   /* If a or b is +Infinity, the result is +Infinity */
+   if (MPFR_IS_INF(op1) || MPFR_IS_INF(op2))
+     {
+       MPFR_SET_INF(r);
+       MPFR_SET_SAME_SIGN(r, op1);
+       return;
+     }
+   MPFR_CLEAR_INF(r);
    /* If a or b is 0, the result is 0 */
-   if ((SIGN(op1)==0)||(SIGN(op2)==0))
+   if ((MPFR_NOTZERO(op1) && MPFR_NOTZERO(op2)) == 0)
-     { SET_ZERO(r);
+     {
-     return;
+       MPFR_SET_ZERO(r);
+       return;
      }
   /* precision of the following calculus */
-   q = PREC(r);
+   q = MPFR_PREC(r);
    p = q + 15;
    /* Initialisations */
    go_on=1;
-   TMP_MARK(marker1);
+   TMP_MARK(marker);
    s=(p-1)/BITS_PER_MP_LIMB+1;
-   MON_INIT(ap, a, p, s);
+   MPFR_INIT(ap, a, p, s);
-   MON_INIT(bp, b, p, s);
+   MPFR_INIT(bp, b, p, s);
-   TMP_MARK(marker2);
+   MPFR_INIT(up, u, p, s);
-   MON_INIT(up, u, p, s);
+   MPFR_INIT(vp, v, p, s);
-   MON_INIT(vp, v, p, s);
+   MPFR_INIT(tmpup, tmpu, p, s);
-   MON_INIT(tmpup, tmpu, p, s);
+   MPFR_INIT(tmpvp, tmpv, p, s);
-   MON_INIT(tmpvp, tmpv, p, s);
+   MPFR_INIT(tmpp, tmp, p, s);
-   MON_INIT(tmpp, tmp, p, s);
-   /* b and a will be the 2 operands but I want b>= a */
+   /* b and a are the 2 operands but we want b >= a */
-   if (mpfr_cmp(op1,op2) > 0) {
+   if ((compare = mpfr_cmp (op1,op2)) > 0)
-     mpfr_set(b,op1,GMP_RNDN); mpfr_set(a,op2,GMP_RNDN);
+     {
-   }
+       mpfr_set (b,op1,GMP_RNDN);
-   else {
+       mpfr_set (a,op2,GMP_RNDN);
-     mpfr_set(b,op2,GMP_RNDN); mpfr_set(a,op1,GMP_RNDN);
+     }
-   }
+   else if (compare == 0)
+     {
+       mpfr_set (r, op1, rnd_mode);
+       TMP_FREE(marker);
+       return;
+     }
+   else
+     {
+       mpfr_set (b,op2,GMP_RNDN);
+       mpfr_set (a,op1,GMP_RNDN);
+     }
-   vo=mpfr_get_d(b);
+   vo = mpfr_get_d1 (b);
-   uo=mpfr_get_d(a);
+   uo = mpfr_get_d1 (a);
    mpfr_set(u,a,GMP_RNDN);
    mpfr_set(v,b,GMP_RNDN);
    /* Main loop */
    while (go_on) {
-     int err,  eq, can_round;
+     int err, can_round;
+     mp_prec_t eq;
-     eq=0;
+     err=1 + (int) ((3.0/2.0*(double)_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p)
+              +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
-     err=ceil((3.0/2.0*log((double)p)/log(2.0)+1.0)*exp(-(double)p*log(2.0))+3.0*exp(-2.0*(double)p*uo*log(2.0)/(vo-uo)));
      if(p-err-3<=q) {
        p=q+err+4;
-       err=ceil((3.0/2.0*log((double)p)/log(2.0)+1.0)*exp(-(double)p*log(2.0))+3.0*exp(-2.0*(double)p*uo*log(2.0)/(vo-uo)));
+       err= 1 +
+         (int) ((3.0/2.0*_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p)
+                +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
      }
      /* Calculus of un and vn */
-     while (eq<=p-2) {
+     do
-       mpfr_mul(tmp,u,v,GMP_RNDN);
+       {
-       mpfr_sqrt(tmpu,tmp,GMP_RNDN);
+         mpfr_mul(tmp, u, v, GMP_RNDN);
-       mpfr_add(tmp,u,v,GMP_RNDN);
+         mpfr_sqrt (tmpu, tmp, GMP_RNDN);
-       mpfr_div_2exp(tmpv,tmp,1,GMP_RNDN);
+         mpfr_add(tmp, u, v, GMP_RNDN);
-       mpfr_set(u,tmpu,GMP_RNDN);
+         mpfr_div_2ui(tmpv, tmp, 1, GMP_RNDN);
-       mpfr_set(v,tmpv,GMP_RNDN);
+         mpfr_set(u, tmpu, GMP_RNDN);
-       if (mpfr_cmp(v,u)>=0)
+         mpfr_set(v, tmpv, GMP_RNDN);
-         eq=mpfr_cmp2(v,u);
+       }
-       else
+     while (mpfr_cmp2(u, v, &eq) != 0 && eq <= p - 2);
-         eq=mpfr_cmp2(u,v);
-     }
-     /*  printf("avant can_round %i bits faux\n v : ",err+3);
-         mpfr_print_raw(v); printf("\n u : ");
-         mpfr_print_raw(u);printf("\n");*/
      /* Roundability of the result */
        can_round=mpfr_can_round(v,p-err-3,GMP_RNDN,rnd_mode,q);
-Line 145  mpfr_agm(r, a, b, rnd_mode)
+Line 148  mpfr_agm(r, a, b, rnd_mode)
 Line 145  mpfr_agm(r, a, b, rnd_mode)
 Line 148  mpfr_agm(r, a, b, rnd_mode)
        else {
            go_on=1;
            p+=5;
-           TMP_FREE(marker2);
-           TMP_MARK(marker2);
            s=(p-1)/BITS_PER_MP_LIMB+1;
-           MON_INIT(up, u, p, s);
+           MPFR_INIT(up, u, p, s);
-           MON_INIT(vp, v, p, s);
+           MPFR_INIT(vp, v, p, s);
-           MON_INIT(tmpup, tmpu, p, s);
+           MPFR_INIT(tmpup, tmpu, p, s);
-           MON_INIT(tmpvp, tmpv, p, s);
+           MPFR_INIT(tmpvp, tmpv, p, s);
-           MON_INIT(tmpp, tmp, p, s);
+           MPFR_INIT(tmpp, tmp, p, s);
            mpfr_set(u,a,GMP_RNDN);
            mpfr_set(v,b,GMP_RNDN);
        }
-Line 161  mpfr_agm(r, a, b, rnd_mode)
+Line 162  mpfr_agm(r, a, b, rnd_mode)
 Line 161  mpfr_agm(r, a, b, rnd_mode)
 Line 162  mpfr_agm(r, a, b, rnd_mode)
    /* Setting of the result */
-     mpfr_set(r,v,rnd_mode);
+   mpfr_set(r,v,rnd_mode);
    /* Let's clean */
-     TMP_FREE(marker1);
+   TMP_FREE(marker);
    return ;
  }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>