contrib/mpfr/gmp_op.c

   1 /* Implementations of operations between mpfr and mpz/mpq data
   2
   3 Copyright 2001, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 2.1 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* Init and set a mpfr_t with enough precision to store a mpz */
  27 static void
  28 init_set_z (mpfr_ptr t, mpz_srcptr z)
  29 {
  30   mp_prec_t p;
  31   int i;
  32
  33   if (mpz_size (z) <= 1)
  34     p = BITS_PER_MP_LIMB;
  35   else
  36     MPFR_MPZ_SIZEINBASE2 (p, z);
  37   mpfr_init2 (t, p);
  38   i = mpfr_set_z (t, z, GMP_RNDN);
  39   MPFR_ASSERTD (i == 0);
  40 }
  41
  42 /* Init, set a mpfr_t with enough precision to store a mpz_t without round,
  43    call the function, and clear the allocated mpfr_t  */
  44 static int
  45 foo (mpfr_ptr x, mpfr_srcptr y, mpz_srcptr z, mp_rnd_t r,
  46      int (*f)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t))
  47 {
  48   mpfr_t t;
  49   int i;
  50   init_set_z (t, z);
  51   i = (*f) (x, y, t, r);
  52   mpfr_clear (t);
  53   return i;
  54 }
  55
  56 int
  57 mpfr_mul_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r)
  58 {
  59   return foo (y, x, z, r, mpfr_mul);
  60 }
  61
  62 int
  63 mpfr_div_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r)
  64 {
  65   return foo (y, x, z, r, mpfr_div);
  66 }
  67
  68 int
  69 mpfr_add_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r)
  70 {
  71   /* Mpz 0 is unsigned */
  72   if (MPFR_UNLIKELY (mpz_sgn (z) == 0))
  73     return mpfr_set (y, x, r);
  74   else
  75     return foo (y, x, z, r, mpfr_add);
  76 }
  77
  78 int
  79 mpfr_sub_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t r)
  80 {
  81   /* Mpz 0 is unsigned */
  82   if (MPFR_UNLIKELY (mpz_sgn (z) == 0))
  83     return mpfr_set (y, x, r);
  84   else
  85     return foo (y, x, z, r, mpfr_sub);
  86 }
  87
  88 int
  89 mpfr_cmp_z (mpfr_srcptr x, mpz_srcptr z)
  90 {
  91   mpfr_t t;
  92   int res;
  93   init_set_z (t, z);
  94   res = mpfr_cmp (x, t);
  95   mpfr_clear (t);
  96   return res;
  97 }
  98
  99 int
 100 mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode)
 101 {
 102   mpfr_t tmp;
 103   int res;
 104   mp_prec_t p;
 105
 106   if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
 107     return mpfr_mul_ui (y, x, 0, rnd_mode);
 108   else
 109     {
 110       MPFR_MPZ_SIZEINBASE2 (p, mpq_numref (z));
 111       mpfr_init2 (tmp, MPFR_PREC (x) + p);
 112       res = mpfr_mul_z (tmp, x, mpq_numref(z), GMP_RNDN );
 113       MPFR_ASSERTD (res == 0);
 114       res = mpfr_div_z (y, tmp, mpq_denref(z), rnd_mode);
 115       mpfr_clear (tmp);
 116       return res;
 117     }
 118 }
 119
 120 int
 121 mpfr_div_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode)
 122 {
 123   mpfr_t tmp;
 124   int res;
 125   mp_prec_t p;
 126
 127   if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
 128     return mpfr_div_ui (y, x, 0, rnd_mode);
 129   else if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z)) == 0))
 130     p = 0;
 131   else
 132     MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (z));
 133   mpfr_init2 (tmp, MPFR_PREC(x) + p);
 134   res = mpfr_mul_z (tmp, x, mpq_denref(z), GMP_RNDN );
 135   MPFR_ASSERTD( res == 0 );
 136   res = mpfr_div_z (y, tmp, mpq_numref(z), rnd_mode);
 137   mpfr_clear (tmp);
 138   return res;
 139 }
 140
 141 int
 142 mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode)
 143 {
 144   mpfr_t    t,q;
 145   mp_prec_t p;
 146   mp_exp_t  err;
 147   int res;
 148   MPFR_ZIV_DECL (loop);
 149
 150   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 151     {
 152       if (MPFR_IS_NAN (x))
 153         {
 154           MPFR_SET_NAN (y);
 155           MPFR_RET_NAN;
 156         }
 157       else if (MPFR_IS_INF (x))
 158         {
 159           MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0);
 160           MPFR_SET_INF (y);
 161           MPFR_SET_SAME_SIGN (y, x);
 162           MPFR_RET (0);
 163         }
 164       else
 165         {
 166           MPFR_ASSERTD (MPFR_IS_ZERO (x));
 167           if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
 168             return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */
 169           else
 170             return mpfr_set_q (y, z, rnd_mode);
 171         }
 172     }
 173
 174   p = MPFR_PREC (y) + 10;
 175   mpfr_init2 (t, p);
 176   mpfr_init2 (q, p);
 177
 178   MPFR_ZIV_INIT (loop, p);
 179   for (;;)
 180     {
 181       res = mpfr_set_q (q, z, GMP_RNDN);  /* Error <= 1/2 ulp(q) */
 182       /* If z if @INF@ (1/0), res = 0, so it quits immediately */
 183       if (MPFR_UNLIKELY (res == 0))
 184         /* Result is exact so we can add it directly! */
 185         {
 186           res = mpfr_add (y, x, q, rnd_mode);
 187           break;
 188         }
 189       mpfr_add (t, x, q, GMP_RNDN);       /* Error <= 1/2 ulp(t) */
 190       /* Error / ulp(t)      <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t))
 191          If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t)))
 192                              <= 2^(EXP(q)-EXP(t))
 193          If EXP(q)-EXP(t)<0, <= 2^0 */
 194       /* We can get 0, but we can't round since q is inexact */
 195       if (MPFR_LIKELY (!MPFR_IS_ZERO (t)))
 196         {
 197           err = (mp_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0);
 198           if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode)))
 199             {
 200               res = mpfr_set (y, t, rnd_mode);
 201               break;
 202             }
 203         }
 204       MPFR_ZIV_NEXT (loop, p);
 205       mpfr_set_prec (t, p);
 206       mpfr_set_prec (q, p);
 207     }
 208   MPFR_ZIV_FREE (loop);
 209   mpfr_clear (t);
 210   mpfr_clear (q);
 211   return res;
 212 }
 213
 214 int
 215 mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode)
 216 {
 217   mpfr_t t,q;
 218   mp_prec_t p;
 219   int res;
 220   mp_exp_t err;
 221   MPFR_ZIV_DECL (loop);
 222
 223   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 224     {
 225       if (MPFR_IS_NAN (x))
 226         {
 227           MPFR_SET_NAN (y);
 228           MPFR_RET_NAN;
 229         }
 230       else if (MPFR_IS_INF (x))
 231         {
 232           MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0);
 233           MPFR_SET_INF (y);
 234           MPFR_SET_SAME_SIGN (y, x);
 235           MPFR_RET (0);
 236         }
 237       else
 238         {
 239           MPFR_ASSERTD (MPFR_IS_ZERO (x));
 240
 241           if (MPFR_UNLIKELY (mpq_sgn (z) == 0))
 242             return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */
 243           else
 244             {
 245               res =  mpfr_set_q (y, z, MPFR_INVERT_RND (rnd_mode));
 246               MPFR_CHANGE_SIGN (y);
 247               return -res;
 248             }
 249         }
 250     }
 251
 252   p = MPFR_PREC (y) + 10;
 253   mpfr_init2 (t, p);
 254   mpfr_init2 (q, p);
 255
 256   MPFR_ZIV_INIT (loop, p);
 257   for(;;)
 258     {
 259       res = mpfr_set_q(q, z, GMP_RNDN);  /* Error <= 1/2 ulp(q) */
 260       /* If z if @INF@ (1/0), res = 0, so it quits immediately */
 261       if (MPFR_UNLIKELY (res == 0))
 262         /* Result is exact so we can add it directly!*/
 263         {
 264           res = mpfr_sub (y, x, q, rnd_mode);
 265           break;
 266         }
 267       mpfr_sub (t, x, q, GMP_RNDN);       /* Error <= 1/2 ulp(t) */
 268       /* Error / ulp(t)      <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t))
 269          If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t)))
 270                              <= 2^(EXP(q)-EXP(t))
 271                              If EXP(q)-EXP(t)<0, <= 2^0 */
 272       /* We can get 0, but we can't round since q is inexact */
 273       if (MPFR_LIKELY (!MPFR_IS_ZERO (t)))
 274         {
 275           err = (mp_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0);
 276           res = MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode);
 277           if (MPFR_LIKELY (res != 0))  /* We can round! */
 278             {
 279               res = mpfr_set (y, t, rnd_mode);
 280               break;
 281             }
 282         }
 283       MPFR_ZIV_NEXT (loop, p);
 284       mpfr_set_prec (t, p);
 285       mpfr_set_prec (q, p);
 286     }
 287   MPFR_ZIV_FREE (loop);
 288   mpfr_clear (t);
 289   mpfr_clear (q);
 290   return res;
 291 }
 292
 293 int
 294 mpfr_cmp_q (mpfr_srcptr x, mpq_srcptr z)
 295 {
 296   mpfr_t t;
 297   int res;
 298   mp_prec_t p;
 299   /* x < a/b ? <=> x*b < a */
 300   MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0);
 301   MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (z));
 302   mpfr_init2 (t, MPFR_PREC(x) + p);
 303   res = mpfr_mul_z (t, x, mpq_denref (z), GMP_RNDN );
 304   MPFR_ASSERTD (res == 0);
 305   res = mpfr_cmp_z (t, mpq_numref (z) );
 306   mpfr_clear (t);
 307   return res;
 308 }
 309
 310 int
 311 mpfr_cmp_f (mpfr_srcptr x, mpf_srcptr z)
 312 {
 313   mpfr_t t;
 314   int res;
 315
 316   mpfr_init2 (t, MPFR_PREC_MIN + ABS(SIZ(z)) * BITS_PER_MP_LIMB );
 317   res = mpfr_set_f (t, z, GMP_RNDN);
 318   MPFR_ASSERTD (res == 0);
 319   res = mpfr_cmp (x, t);
 320   mpfr_clear (t);
 321   return res;
 322 }