Return to agm.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr

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

1.1       maekawa     1: /* mpfr_agm -- arithmetic-geometric mean of two floating-point numbers
                      2:
                      3: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
                      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 Library General Public License as published by
                      9: the Free Software Foundation; either version 2 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 Library General Public
                     15: License for more details.
                     16:
                     17: You should have received a copy of the GNU Library 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 <math.h>
                     24: #include "gmp.h"
                     25: #include "gmp-impl.h"
                     26: #include "mpfr.h"
                     27:
                     28:
                     29: /*Memory gestion */
                     30: #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;
                     31:
                     32:
                     33:
                     34:
                     35: void
                     36: #ifdef __STDC__
                     37: mpfr_agm(mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, unsigned char rnd_mode)
                     38: #else
                     39: mpfr_agm(r, a, b, rnd_mode)
                     40:      mpfr_ptr r;
                     41:      mpfr_srcptr a;
                     42:      mpfr_srcptr b;
                     43:      unsigned char rnd_mode;
                     44: #endif
                     45: {
                     46:   int s, p, q, go_on;
                     47:   double uo, vo;
                     48:   mp_limb_t *up, *vp, *tmpp, *tmpup, *tmpvp, *ap, *bp;
                     49:   mpfr_t u, v, tmp, tmpu, tmpv, a, b;
                     50:   TMP_DECL(marker1);
                     51:   TMP_DECL(marker2);
                     52:
                     53:
                     54:   /* If a or b is NaN, the result is NaN */
                     55:   if (FLAG_NAN(op1) || FLAG_NAN(op2))
                     56:     { SET_NAN(r); return; }
                     57:
                     58:
                     59:   /* If a or b is negative, the result is NaN */
                     60:   if ((SIGN(op1)<0)||(SIGN(op2)<0))
                     61:     { SET_NAN(r); return; }
                     62:
                     63:
                     64:
                     65:   /* If a or b is 0, the result is 0 */
                     66:   if ((SIGN(op1)==0)||(SIGN(op2)==0))
                     67:     { SET_ZERO(r);
                     68:     return;
                     69:     }
                     70:
                     71:  /* precision of the following calculus */
                     72:   q = PREC(r);
                     73:   p = q + 15;
                     74:
                     75:
                     76:   /* Initialisations */
                     77:   go_on=1;
                     78:
                     79:   TMP_MARK(marker1);
                     80:   s=(p-1)/BITS_PER_MP_LIMB+1;
                     81:   MON_INIT(ap, a, p, s);
                     82:   MON_INIT(bp, b, p, s);
                     83:   TMP_MARK(marker2);
                     84:   MON_INIT(up, u, p, s);
                     85:   MON_INIT(vp, v, p, s);
                     86:   MON_INIT(tmpup, tmpu, p, s);
                     87:   MON_INIT(tmpvp, tmpv, p, s);
                     88:   MON_INIT(tmpp, tmp, p, s);
                     89:
                     90:
                     91:
                     92:   /* b and a will be the 2 operands but I want b>= a */
                     93:   if (mpfr_cmp(op1,op2) > 0) {
                     94:     mpfr_set(b,op1,GMP_RNDN); mpfr_set(a,op2,GMP_RNDN);
                     95:   }
                     96:   else {
                     97:     mpfr_set(b,op2,GMP_RNDN); mpfr_set(a,op1,GMP_RNDN);
                     98:   }
                     99:
                    100:   vo=mpfr_get_d(b);
                    101:   uo=mpfr_get_d(a);
                    102:
                    103:   mpfr_set(u,a,GMP_RNDN);
                    104:   mpfr_set(v,b,GMP_RNDN);
                    105:
                    106:
                    107:   /* Main loop */
                    108:
                    109:   while (go_on) {
                    110:     int err,  eq, can_round;
                    111:
                    112:     eq=0;
                    113:
                    114:     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)));
                    115:     if(p-err-3<=q) {
                    116:       p=q+err+4;
                    117:       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)));
                    118:     }
                    119:
                    120:     /* Calculus of un and vn */
                    121:     while (eq<=p-2) {
                    122:       mpfr_mul(tmp,u,v,GMP_RNDN);
                    123:       mpfr_sqrt(tmpu,tmp,GMP_RNDN);
                    124:       mpfr_add(tmp,u,v,GMP_RNDN);
                    125:       mpfr_div_2exp(tmpv,tmp,1,GMP_RNDN);
                    126:       mpfr_set(u,tmpu,GMP_RNDN);
                    127:       mpfr_set(v,tmpv,GMP_RNDN);
                    128:       if (mpfr_cmp(v,u)>=0)
                    129:        eq=mpfr_cmp2(v,u);
                    130:       else
                    131:        eq=mpfr_cmp2(u,v);
                    132:     }
                    133:
                    134:     /*  printf("avant can_round %i bits faux\n v : ",err+3);
                    135:        mpfr_print_raw(v); printf("\n u : ");
                    136:        mpfr_print_raw(u);printf("\n");*/
                    137:
                    138:
                    139:     /* Roundability of the result */
                    140:       can_round=mpfr_can_round(v,p-err-3,GMP_RNDN,rnd_mode,q);
                    141:
                    142:       if (can_round)
                    143:        go_on=0;
                    144:
                    145:       else {
                    146:          go_on=1;
                    147:          p+=5;
                    148:          TMP_FREE(marker2);
                    149:          TMP_MARK(marker2);
                    150:          s=(p-1)/BITS_PER_MP_LIMB+1;
                    151:          MON_INIT(up, u, p, s);
                    152:          MON_INIT(vp, v, p, s);
                    153:          MON_INIT(tmpup, tmpu, p, s);
                    154:          MON_INIT(tmpvp, tmpv, p, s);
                    155:          MON_INIT(tmpp, tmp, p, s);
                    156:          mpfr_set(u,a,GMP_RNDN);
                    157:          mpfr_set(v,b,GMP_RNDN);
                    158:       }
                    159:   }
                    160:   /* End of while */
                    161:
                    162:   /* Setting of the result */
                    163:
                    164:     mpfr_set(r,v,rnd_mode);
                    165:
                    166:
                    167:   /* Let's clean */
                    168:     TMP_FREE(marker1);
                    169:
                    170:   return ;
                    171: }
                    172: