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     ! 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: }