contrib/gmp/mpq/set_d.c

   1 /* mpq_set_d(mpq_t q, double d) -- Set q to d without rounding.
   2
   3 Copyright 2000, 2002, 2003 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP 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 Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  19
  20 #include "config.h"
  21
  22 #if HAVE_FLOAT_H
  23 #include <float.h>  /* for DBL_MAX */
  24 #endif
  25
  26 #include "gmp.h"
  27 #include "gmp-impl.h"
  28 #include "longlong.h"
  29
  30 #if LIMBS_PER_DOUBLE > 4
  31   choke me
  32 #endif
  33
  34 void
  35 mpq_set_d (mpq_ptr dest, double d)
  36 {
  37   int negative;
  38   mp_exp_t exp;
  39   mp_limb_t tp[LIMBS_PER_DOUBLE];
  40   mp_ptr np, dp;
  41   mp_size_t nn, dn;
  42   int c;
  43
  44   DOUBLE_NAN_INF_ACTION (d,
  45                          __gmp_invalid_operation (),
  46                          __gmp_invalid_operation ());
  47
  48   negative = d < 0;
  49   d = ABS (d);
  50
  51   exp = __gmp_extract_double (tp, d);
  52
  53   /* There are two main version of the conversion.  The `then' arm handles
  54      numbers with a fractional part, while the `else' arm handles integers.  */
  55 #if LIMBS_PER_DOUBLE == 4
  56   if (exp <= 1 || (exp == 2 && (tp[0] | tp[1]) != 0))
  57 #endif
  58 #if LIMBS_PER_DOUBLE == 3
  59   if (exp <= 1 || (exp == 2 && tp[0] != 0))
  60 #endif
  61 #if LIMBS_PER_DOUBLE == 2
  62   if (exp <= 1)
  63 #endif
  64     {
  65       if (d == 0.0)
  66         {
  67           SIZ(&(dest->_mp_num)) = 0;
  68           SIZ(&(dest->_mp_den)) = 1;
  69           PTR(&(dest->_mp_den))[0] = 1;
  70           return;
  71         }
  72
  73       dn = -exp;
  74       MPZ_REALLOC (&(dest->_mp_num), 3);
  75       np = PTR(&(dest->_mp_num));
  76 #if LIMBS_PER_DOUBLE == 4
  77       if ((tp[0] | tp[1] | tp[2]) == 0)
  78         np[0] = tp[3], nn = 1;
  79       else if ((tp[0] | tp[1]) == 0)
  80         np[1] = tp[3], np[0] = tp[2], nn = 2;
  81       else if (tp[0] == 0)
  82         np[2] = tp[3], np[1] = tp[2], np[0] = tp[1], nn = 3;
  83       else
  84         np[3] = tp[3], np[2] = tp[2], np[1] = tp[1], np[0] = tp[0], nn = 4;
  85 #endif
  86 #if LIMBS_PER_DOUBLE == 3
  87       if ((tp[0] | tp[1]) == 0)
  88         np[0] = tp[2], nn = 1;
  89       else if (tp[0] == 0)
  90         np[1] = tp[2], np[0] = tp[1], nn = 2;
  91       else
  92         np[2] = tp[2], np[1] = tp[1], np[0] = tp[0], nn = 3;
  93 #endif
  94 #if LIMBS_PER_DOUBLE == 2
  95       if (tp[0] == 0)
  96         np[0] = tp[1], nn = 1;
  97       else
  98         np[1] = tp[1], np[0] = tp[0], nn = 2;
  99 #endif
 100       dn += nn + 1;
 101       ASSERT_ALWAYS (dn > 0);
 102       MPZ_REALLOC (&(dest->_mp_den), dn);
 103       dp = PTR(&(dest->_mp_den));
 104       MPN_ZERO (dp, dn - 1);
 105       dp[dn - 1] = 1;
 106       count_trailing_zeros (c, np[0] | dp[0]);
 107       if (c != 0)
 108         {
 109           mpn_rshift (np, np, nn, c);
 110           nn -= np[nn - 1] == 0;
 111           mpn_rshift (dp, dp, dn, c);
 112           dn -= dp[dn - 1] == 0;
 113         }
 114       SIZ(&(dest->_mp_den)) = dn;
 115       SIZ(&(dest->_mp_num)) = negative ? -nn : nn;
 116     }
 117   else
 118     {
 119       nn = exp;
 120       MPZ_REALLOC (&(dest->_mp_num), nn);
 121       np = PTR(&(dest->_mp_num));
 122       switch (nn)
 123         {
 124         default:
 125           MPN_ZERO (np, nn - LIMBS_PER_DOUBLE);
 126           np += nn - LIMBS_PER_DOUBLE;
 127           /* fall through */
 128 #if LIMBS_PER_DOUBLE == 2
 129         case 2:
 130           np[1] = tp[1], np[0] = tp[0];
 131           break;
 132 #endif
 133 #if LIMBS_PER_DOUBLE == 3
 134         case 3:
 135           np[2] = tp[2], np[1] = tp[1], np[0] = tp[0];
 136           break;
 137         case 2:
 138           np[1] = tp[2], np[0] = tp[1];
 139           break;
 140 #endif
 141 #if LIMBS_PER_DOUBLE == 4
 142         case 4:
 143           np[3] = tp[3], np[2] = tp[2], np[1] = tp[1], np[0] = tp[0];
 144           break;
 145         case 3:
 146           np[2] = tp[3], np[1] = tp[2], np[0] = tp[1];
 147           break;
 148         case 2:
 149           np[1] = tp[3], np[0] = tp[2];
 150           break;
 151 #endif
 152         }
 153       dp = PTR(&(dest->_mp_den));
 154       dp[0] = 1;
 155       SIZ(&(dest->_mp_den)) = 1;
 156       SIZ(&(dest->_mp_num)) = negative ? -nn : nn;
 157     }
 158 }