source/libs/gmp/gmp-src/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, 2012, 2014 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 either:
   9
  10   * the GNU Lesser General Public License as published by the Free
  11     Software Foundation; either version 3 of the License, or (at your
  12     option) any later version.
  13
  14 or
  15
  16   * the GNU General Public License as published by the Free Software
  17     Foundation; either version 2 of the License, or (at your option) any
  18     later version.
  19
  20 or both in parallel, as here.
  21
  22 The GNU MP Library is distributed in the hope that it will be useful, but
  23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  24 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  25 for more details.
  26
  27 You should have received copies of the GNU General Public License and the
  28 GNU Lesser General Public License along with the GNU MP Library.  If not,
  29 see https://www.gnu.org/licenses/.  */
  30
  31 #include "config.h"
  32
  33 #if HAVE_FLOAT_H
  34 #include <float.h>  /* for DBL_MAX */
  35 #endif
  36
  37 #include "gmp.h"
  38 #include "gmp-impl.h"
  39 #include "longlong.h"
  40
  41 #if LIMBS_PER_DOUBLE > 4
  42   choke me
  43 #endif
  44
  45 void
  46 mpq_set_d (mpq_ptr dest, double d)
  47 {
  48   int negative;
  49   mp_exp_t exp;
  50   mp_limb_t tp[LIMBS_PER_DOUBLE];
  51   mp_ptr np, dp;
  52   mp_size_t nn, dn;
  53   int c;
  54
  55   DOUBLE_NAN_INF_ACTION (d,
  56                          __gmp_invalid_operation (),
  57                          __gmp_invalid_operation ());
  58
  59   negative = d < 0;
  60   d = ABS (d);
  61
  62   exp = __gmp_extract_double (tp, d);
  63
  64   /* There are two main version of the conversion.  The `then' arm handles
  65      numbers with a fractional part, while the `else' arm handles integers.  */
  66 #if LIMBS_PER_DOUBLE == 4
  67   if (exp <= 1 || (exp == 2 && (tp[0] | tp[1]) != 0))
  68 #endif
  69 #if LIMBS_PER_DOUBLE == 3
  70   if (exp <= 1 || (exp == 2 && tp[0] != 0))
  71 #endif
  72 #if LIMBS_PER_DOUBLE == 2
  73   if (exp <= 1)
  74 #endif
  75     {
  76       if (d == 0.0)
  77         {
  78           SIZ(NUM(dest)) = 0;
  79           SIZ(DEN(dest)) = 1;
  80           PTR(DEN(dest))[0] = 1;
  81           return;
  82         }
  83
  84 #if LIMBS_PER_DOUBLE == 4
  85       np = MPZ_NEWALLOC (NUM(dest), 4);
  86       if ((tp[0] | tp[1] | tp[2]) == 0)
  87         np[0] = tp[3], nn = 1;
  88       else if ((tp[0] | tp[1]) == 0)
  89         np[1] = tp[3], np[0] = tp[2], nn = 2;
  90       else if (tp[0] == 0)
  91         np[2] = tp[3], np[1] = tp[2], np[0] = tp[1], nn = 3;
  92       else
  93         np[3] = tp[3], np[2] = tp[2], np[1] = tp[1], np[0] = tp[0], nn = 4;
  94 #endif
  95 #if LIMBS_PER_DOUBLE == 3
  96       np = MPZ_NEWALLOC (NUM(dest), 3);
  97       if ((tp[0] | tp[1]) == 0)
  98         np[0] = tp[2], nn = 1;
  99       else if (tp[0] == 0)
 100         np[1] = tp[2], np[0] = tp[1], nn = 2;
 101       else
 102         np[2] = tp[2], np[1] = tp[1], np[0] = tp[0], nn = 3;
 103 #endif
 104 #if LIMBS_PER_DOUBLE == 2
 105       np = MPZ_NEWALLOC (NUM(dest), 2);
 106       if (tp[0] == 0)
 107         np[0] = tp[1], nn = 1;
 108       else
 109         np[1] = tp[1], np[0] = tp[0], nn = 2;
 110 #endif
 111       dn = nn + 1 - exp;
 112       ASSERT (dn > 0); /* -exp >= -1; nn >= 1*/
 113       dp = MPZ_NEWALLOC (DEN(dest), dn);
 114       MPN_ZERO (dp, dn - 1);
 115       dp[dn - 1] = 1;
 116       count_trailing_zeros (c, np[0] | dp[0]);
 117       if (c != 0)
 118         {
 119           mpn_rshift (np, np, nn, c);
 120           nn -= np[nn - 1] == 0;
 121           --dn;
 122           dp[dn - 1] = CNST_LIMB(1) << (GMP_LIMB_BITS - c);
 123         }
 124       SIZ(DEN(dest)) = dn;
 125     }
 126   else
 127     {
 128       nn = exp;
 129       np = MPZ_NEWALLOC (NUM(dest), nn);
 130       switch (nn)
 131         {
 132         default:
 133           MPN_ZERO (np, nn - LIMBS_PER_DOUBLE);
 134           np += nn - LIMBS_PER_DOUBLE;
 135           /* fall through */
 136 #if LIMBS_PER_DOUBLE == 2
 137         case 2:
 138           np[1] = tp[1], np[0] = tp[0];
 139           break;
 140 #endif
 141 #if LIMBS_PER_DOUBLE == 3
 142         case 3:
 143           np[2] = tp[2], np[1] = tp[1], np[0] = tp[0];
 144           break;
 145         case 2:
 146           np[1] = tp[2], np[0] = tp[1];
 147           break;
 148 #endif
 149 #if LIMBS_PER_DOUBLE == 4
 150         case 4:
 151           np[3] = tp[3], np[2] = tp[2], np[1] = tp[1], np[0] = tp[0];
 152           break;
 153         case 3:
 154           np[2] = tp[3], np[1] = tp[2], np[0] = tp[1];
 155           break;
 156         case 2:
 157           np[1] = tp[3], np[0] = tp[2];
 158           break;
 159 #endif
 160         }
 161       *PTR(DEN(dest)) = 1;
 162       SIZ(DEN(dest)) = 1;
 163     }
 164   SIZ(NUM(dest)) = negative ? -nn : nn;
 165 }