Return to fma.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr

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

1.1       ohara       1: /* mpfr_fma -- Floating multiply-add
                      2:
                      3: Copyright 2001, 2002 Free Software Foundation, Inc.
                      4:
                      5: This file is part of the MPFR Library.
                      6:
                      7: The MPFR Library is free software; you can redistribute
                      8: it and/or modify it under the terms of the GNU Lesser
                      9: General Public License as published by the Free Software
                     10: Foundation; either version 2.1 of the License, or (at your
                     11: option) any later version.
                     12:
                     13: The MPFR Library is distributed in the hope that it will
                     14: be useful, but WITHOUT ANY WARRANTY; without even the
                     15: implied warranty of MERCHANTABILITY or FITNESS FOR A
                     16: PARTICULAR PURPOSE.  See the GNU Lesser General Public
                     17: License for more details.
                     18:
                     19: You should have received a copy of the GNU Lesser
                     20: General Public License along with the MPFR Library; see
                     21: the file COPYING.LIB.  If not, write to the Free Software
                     22: Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     23: MA 02111-1307, USA. */
                     24:
                     25: #include "gmp.h"
                     26: #include "gmp-impl.h"
                     27: #include "mpfr.h"
                     28: #include "mpfr-impl.h"
                     29:
                     30:
                     31: /* The computation of fma of x y and u is done by
                     32:     fma(s,x,y,z)= z + x*y = s
                     33: */
                     34:
                     35: int
                     36: mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
                     37:           mp_rnd_t rnd_mode)
                     38: {
                     39:   int inexact;
                     40:   mpfr_t u;
                     41:
                     42:   /* particular cases */
                     43:   if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y) || MPFR_IS_NAN(z))
                     44:     {
                     45:       MPFR_SET_NAN(s);
                     46:       MPFR_RET_NAN;
                     47:     }
                     48:
                     49:   if (MPFR_IS_INF(x) || MPFR_IS_INF(y))
                     50:     {
                     51:       /* cases Inf*0+z, 0*Inf+z, Inf-Inf */
                     52:       if ((MPFR_IS_FP(y) && MPFR_IS_ZERO(y)) ||
                     53:           (MPFR_IS_FP(x) && MPFR_IS_ZERO(x)) ||
                     54:           (MPFR_IS_INF(z) && ((MPFR_SIGN(x) * MPFR_SIGN(y)) != MPFR_SIGN(z))))
                     55:         {
                     56:           MPFR_SET_NAN(s);
                     57:           MPFR_RET_NAN;
                     58:         }
                     59:
                     60:       MPFR_CLEAR_NAN(s);
                     61:
                     62:       if (MPFR_IS_INF(z)) /* case Inf-Inf already checked above */
                     63:         {
                     64:           MPFR_SET_INF(s);
                     65:           MPFR_SET_SAME_SIGN(s, z);
                     66:           MPFR_RET(0);
                     67:         }
                     68:       else /* z is finite */
                     69:         {
                     70:           MPFR_SET_INF(s);
                     71:           if (MPFR_SIGN(s) != (MPFR_SIGN(x) * MPFR_SIGN(y)))
                     72:             MPFR_CHANGE_SIGN(s);
                     73:           MPFR_RET(0);
                     74:         }
                     75:     }
                     76:
                     77:   MPFR_CLEAR_NAN(s);
                     78:
                     79:   /* now x and y are finite */
                     80:   if (MPFR_IS_INF(z))
                     81:     {
                     82:       MPFR_SET_INF(s);
                     83:       MPFR_SET_SAME_SIGN(s, z);
                     84:       MPFR_RET(0);
                     85:     }
                     86:
                     87:   MPFR_CLEAR_INF(s);
                     88:
                     89:   if (MPFR_IS_ZERO(x) || MPFR_IS_ZERO(y))
                     90:     {
                     91:       if (MPFR_IS_ZERO(z))
                     92:         {
                     93:           int sign_p, sign_z;
                     94:           sign_p = MPFR_SIGN(x) * MPFR_SIGN(y);
                     95:           sign_z = MPFR_SIGN(z);
                     96:           if (MPFR_SIGN(s) !=
                     97:               (rnd_mode != GMP_RNDD ?
                     98:                ((sign_p < 0 && sign_z < 0) ? -1 : 1) :
                     99:                ((sign_p > 0 && sign_z > 0) ? 1 : -1)))
                    100:             {
                    101:               MPFR_CHANGE_SIGN(s);
                    102:             }
                    103:           MPFR_SET_ZERO(s);
                    104:           MPFR_RET(0);
                    105:         }
                    106:       else
                    107:         return mpfr_set (s, z, rnd_mode);
                    108:     }
                    109:
                    110:   if (MPFR_IS_ZERO(z))
                    111:     return mpfr_mul (s, x, y, rnd_mode);
                    112:
                    113:   /* if we take prec(u) >= prec(x) + prec(y), the product
                    114:      u <- x*y is always exact */
                    115:   mpfr_init2 (u, MPFR_PREC(x) + MPFR_PREC(y));
                    116:   mpfr_mul (u, x, y, GMP_RNDN); /* always exact */
                    117:   inexact = mpfr_add (s, z, u, rnd_mode);
                    118:   mpfr_clear(u);
                    119:
                    120:   return inexact;
                    121: }