Blame mpf/set_q.c

Packit 5c3484
/* mpf_set_q (mpf_t rop, mpq_t op) -- Convert the rational op to the float rop.
Packit 5c3484
Packit 5c3484
Copyright 1996, 1999, 2001, 2002, 2004, 2005 Free Software Foundation, Inc.
Packit 5c3484
Packit 5c3484
This file is part of the GNU MP Library.
Packit 5c3484
Packit 5c3484
The GNU MP Library is free software; you can redistribute it and/or modify
Packit 5c3484
it under the terms of either:
Packit 5c3484
Packit 5c3484
  * the GNU Lesser General Public License as published by the Free
Packit 5c3484
    Software Foundation; either version 3 of the License, or (at your
Packit 5c3484
    option) any later version.
Packit 5c3484
Packit 5c3484
or
Packit 5c3484
Packit 5c3484
  * the GNU General Public License as published by the Free Software
Packit 5c3484
    Foundation; either version 2 of the License, or (at your option) any
Packit 5c3484
    later version.
Packit 5c3484
Packit 5c3484
or both in parallel, as here.
Packit 5c3484
Packit 5c3484
The GNU MP Library is distributed in the hope that it will be useful, but
Packit 5c3484
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
Packit 5c3484
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
Packit 5c3484
for more details.
Packit 5c3484
Packit 5c3484
You should have received copies of the GNU General Public License and the
Packit 5c3484
GNU Lesser General Public License along with the GNU MP Library.  If not,
Packit 5c3484
see https://www.gnu.org/licenses/.  */
Packit 5c3484
Packit 5c3484
#include <stdio.h>  /* for NULL */
Packit 5c3484
#include "gmp.h"
Packit 5c3484
#include "gmp-impl.h"
Packit 5c3484
#include "longlong.h"
Packit 5c3484
Packit 5c3484
Packit 5c3484
/* As usual the aim is to produce PREC(r) limbs, with the high non-zero.
Packit 5c3484
   The basic mpn_tdiv_qr produces a quotient of nsize-dsize+1 limbs, with
Packit 5c3484
   either the high or second highest limb non-zero.  We arrange for
Packit 5c3484
   nsize-dsize+1 to equal prec+1, hence giving either prec or prec+1 result
Packit 5c3484
   limbs at PTR(r).
Packit 5c3484
Packit 5c3484
   nsize-dsize+1 == prec+1 is achieved by adjusting num(q), either dropping
Packit 5c3484
   low limbs if it's too big, or padding with low zeros if it's too small.
Packit 5c3484
   The full given den(q) is always used.
Packit 5c3484
Packit 5c3484
   We cannot truncate den(q), because even when it's much bigger than prec
Packit 5c3484
   the last limbs can still influence the final quotient.  Often they don't,
Packit 5c3484
   but we leave optimization of that to a prospective quotient-only mpn
Packit 5c3484
   division.
Packit 5c3484
Packit 5c3484
   Not done:
Packit 5c3484
Packit 5c3484
   If den(q) is a power of 2 then we may end up with low zero limbs on the
Packit 5c3484
   result.  But nothing is done about this, since it should be unlikely on
Packit 5c3484
   random data, and can be left to an application to call mpf_div_2exp if it
Packit 5c3484
   might occur with any frequency.
Packit 5c3484
Packit 5c3484
   Enhancements:
Packit 5c3484
Packit 5c3484
   The high quotient limb is non-zero when high{np,dsize} >= {dp,dsize}.  We
Packit 5c3484
   could make that comparison and use qsize==prec instead of qsize==prec+1,
Packit 5c3484
   to save one limb in the division.
Packit 5c3484
Packit 5c3484
   Future:
Packit 5c3484
Packit 5c3484
   If/when mpn_tdiv_qr supports its qxn parameter we can use that instead of
Packit 5c3484
   padding n with zeros in temporary space.
Packit 5c3484
Packit 5c3484
   If/when a quotient-only division exists it can be used here immediately.
Packit 5c3484
   remp is only to satisfy mpn_tdiv_qr, the remainder is not used.  */
Packit 5c3484
Packit 5c3484
void
Packit 5c3484
mpf_set_q (mpf_t r, mpq_srcptr q)
Packit 5c3484
{
Packit 5c3484
  mp_srcptr np, dp;
Packit 5c3484
  mp_size_t prec, nsize, dsize, qsize, prospective_qsize, tsize, zeros;
Packit 5c3484
  mp_size_t sign_quotient, high_zero;
Packit 5c3484
  mp_ptr qp, tp, remp;
Packit 5c3484
  mp_exp_t exp;
Packit 5c3484
  TMP_DECL;
Packit 5c3484
Packit 5c3484
  ASSERT (SIZ(&q->_mp_den) > 0);  /* canonical q */
Packit 5c3484
Packit 5c3484
  nsize = SIZ (&q->_mp_num);
Packit 5c3484
  dsize = SIZ (&q->_mp_den);
Packit 5c3484
Packit 5c3484
  if (UNLIKELY (nsize == 0))
Packit 5c3484
    {
Packit 5c3484
      SIZ (r) = 0;
Packit 5c3484
      EXP (r) = 0;
Packit 5c3484
      return;
Packit 5c3484
    }
Packit 5c3484
Packit 5c3484
  TMP_MARK;
Packit 5c3484
Packit 5c3484
  prec = PREC (r);
Packit 5c3484
  qp = PTR (r);
Packit 5c3484
Packit 5c3484
  sign_quotient = nsize;
Packit 5c3484
  nsize = ABS (nsize);
Packit 5c3484
  np = PTR (&q->_mp_num);
Packit 5c3484
  dp = PTR (&q->_mp_den);
Packit 5c3484
Packit 5c3484
  prospective_qsize = nsize - dsize + 1;  /* q from using given n,d sizes */
Packit 5c3484
  exp = prospective_qsize;                /* ie. number of integer limbs */
Packit 5c3484
  qsize = prec + 1;                       /* desired q */
Packit 5c3484
Packit 5c3484
  zeros = qsize - prospective_qsize;   /* n zeros to get desired qsize */
Packit 5c3484
  tsize = nsize + zeros;               /* possible copy of n */
Packit 5c3484
Packit 5c3484
  if (WANT_TMP_DEBUG)
Packit 5c3484
    {
Packit 5c3484
      /* separate alloc blocks, for malloc debugging */
Packit 5c3484
      remp = TMP_ALLOC_LIMBS (dsize);
Packit 5c3484
      tp = NULL;
Packit 5c3484
      if (zeros > 0)
Packit 5c3484
        tp = TMP_ALLOC_LIMBS (tsize);
Packit 5c3484
    }
Packit 5c3484
  else
Packit 5c3484
    {
Packit 5c3484
      /* one alloc with a conditionalized size, for efficiency */
Packit 5c3484
      mp_size_t size = dsize + (zeros > 0 ? tsize : 0);
Packit 5c3484
      remp = TMP_ALLOC_LIMBS (size);
Packit 5c3484
      tp = remp + dsize;
Packit 5c3484
    }
Packit 5c3484
Packit 5c3484
  if (zeros > 0)
Packit 5c3484
    {
Packit 5c3484
      /* pad n with zeros into temporary space */
Packit 5c3484
      MPN_ZERO (tp, zeros);
Packit 5c3484
      MPN_COPY (tp+zeros, np, nsize);
Packit 5c3484
      np = tp;
Packit 5c3484
      nsize = tsize;
Packit 5c3484
    }
Packit 5c3484
  else
Packit 5c3484
    {
Packit 5c3484
      /* shorten n to get desired qsize */
Packit 5c3484
      nsize += zeros;
Packit 5c3484
      np -= zeros;
Packit 5c3484
    }
Packit 5c3484
Packit 5c3484
  ASSERT (nsize-dsize+1 == qsize);
Packit 5c3484
  mpn_tdiv_qr (qp, remp, (mp_size_t) 0, np, nsize, dp, dsize);
Packit 5c3484
Packit 5c3484
  /* strip possible zero high limb */
Packit 5c3484
  high_zero = (qp[qsize-1] == 0);
Packit 5c3484
  qsize -= high_zero;
Packit 5c3484
  exp -= high_zero;
Packit 5c3484
Packit 5c3484
  EXP (r) = exp;
Packit 5c3484
  SIZ (r) = sign_quotient >= 0 ? qsize : -qsize;
Packit 5c3484
Packit 5c3484
  TMP_FREE;
Packit 5c3484
}