[BACK]Return to agm.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpfr

Annotation of OpenXM_contrib/gmp/mpfr/agm.c, Revision 1.1.1.2

1.1       maekawa     1: /* mpfr_agm -- arithmetic-geometric mean of two floating-point numbers
                      2:
1.1.1.2 ! ohara       3: Copyright 1999, 2000, 2001, 2002 Free Software Foundation.
1.1       maekawa     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
1.1.1.2 ! ohara       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
1.1       maekawa    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
1.1.1.2 ! ohara      14: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
1.1       maekawa    15: License for more details.
                     16:
1.1.1.2 ! ohara      17: You should have received a copy of the GNU Lesser General Public License
1.1       maekawa    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"
1.1.1.2 ! ohara      25: #include "mpfr-impl.h"
1.1       maekawa    26:
                     27: void
1.1.1.2 ! ohara      28: mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode)
1.1       maekawa    29: {
1.1.1.2 ! ohara      30:   int s, go_on, compare;
        !            31:   mp_prec_t p, q;
1.1       maekawa    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;
1.1.1.2 ! ohara      35:   TMP_DECL(marker);
1.1       maekawa    36:
                     37:   /* If a or b is NaN, the result is NaN */
1.1.1.2 ! ohara      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:     }
1.1       maekawa    44:
1.1.1.2 ! ohara      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:     }
1.1       maekawa    52:
1.1.1.2 ! ohara      53:   MPFR_CLEAR_NAN(r);
1.1       maekawa    54:
1.1.1.2 ! ohara      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:     }
1.1       maekawa    62:
1.1.1.2 ! ohara      63:   MPFR_CLEAR_INF(r);
1.1       maekawa    64:
                     65:   /* If a or b is 0, the result is 0 */
1.1.1.2 ! ohara      66:   if ((MPFR_NOTZERO(op1) && MPFR_NOTZERO(op2)) == 0)
        !            67:     {
        !            68:       MPFR_SET_ZERO(r);
        !            69:       return;
1.1       maekawa    70:     }
                     71:
                     72:  /* precision of the following calculus */
1.1.1.2 ! ohara      73:   q = MPFR_PREC(r);
1.1       maekawa    74:   p = q + 15;
                     75:
                     76:   /* Initialisations */
                     77:   go_on=1;
                     78:
1.1.1.2 ! ohara      79:   TMP_MARK(marker);
1.1       maekawa    80:   s=(p-1)/BITS_PER_MP_LIMB+1;
1.1.1.2 ! ohara      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);
1.1       maekawa    88:
                     89:
                     90:
1.1.1.2 ! ohara      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:     }
1.1       maekawa   108:
1.1.1.2 ! ohara     109:   vo = mpfr_get_d1 (b);
        !           110:   uo = mpfr_get_d1 (a);
1.1       maekawa   111:
                    112:   mpfr_set(u,a,GMP_RNDN);
                    113:   mpfr_set(v,b,GMP_RNDN);
                    114:
                    115:   /* Main loop */
                    116:
                    117:   while (go_on) {
1.1.1.2 ! ohara     118:     int err, can_round;
        !           119:     mp_prec_t eq;
1.1       maekawa   120:
1.1.1.2 ! ohara     121:     err=1 + (int) ((3.0/2.0*(double)_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p)
        !           122:             +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
1.1       maekawa   123:     if(p-err-3<=q) {
                    124:       p=q+err+4;
1.1.1.2 ! ohara     125:       err= 1 +
        !           126:        (int) ((3.0/2.0*_mpfr_ceil_log2((double)p)+1.0)*_mpfr_ceil_exp2(-(double)p)
        !           127:               +3.0*_mpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
1.1       maekawa   128:     }
                    129:
                    130:     /* Calculus of un and vn */
1.1.1.2 ! ohara     131:     do
        !           132:       {
        !           133:         mpfr_mul(tmp, u, v, GMP_RNDN);
        !           134:         mpfr_sqrt (tmpu, tmp, GMP_RNDN);
        !           135:         mpfr_add(tmp, u, v, GMP_RNDN);
        !           136:         mpfr_div_2ui(tmpv, tmp, 1, GMP_RNDN);
        !           137:         mpfr_set(u, tmpu, GMP_RNDN);
        !           138:         mpfr_set(v, tmpv, GMP_RNDN);
        !           139:       }
        !           140:     while (mpfr_cmp2(u, v, &eq) != 0 && eq <= p - 2);
1.1       maekawa   141:
                    142:     /* Roundability of the result */
                    143:       can_round=mpfr_can_round(v,p-err-3,GMP_RNDN,rnd_mode,q);
                    144:
                    145:       if (can_round)
                    146:        go_on=0;
                    147:
                    148:       else {
                    149:          go_on=1;
                    150:          p+=5;
                    151:          s=(p-1)/BITS_PER_MP_LIMB+1;
1.1.1.2 ! ohara     152:          MPFR_INIT(up, u, p, s);
        !           153:          MPFR_INIT(vp, v, p, s);
        !           154:          MPFR_INIT(tmpup, tmpu, p, s);
        !           155:          MPFR_INIT(tmpvp, tmpv, p, s);
        !           156:          MPFR_INIT(tmpp, tmp, p, s);
1.1       maekawa   157:          mpfr_set(u,a,GMP_RNDN);
                    158:          mpfr_set(v,b,GMP_RNDN);
                    159:       }
                    160:   }
                    161:   /* End of while */
                    162:
                    163:   /* Setting of the result */
                    164:
1.1.1.2 ! ohara     165:   mpfr_set(r,v,rnd_mode);
1.1       maekawa   166:
                    167:   /* Let's clean */
1.1.1.2 ! ohara     168:   TMP_FREE(marker);
        !           169:
1.1       maekawa   170:   return ;
                    171: }

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