sqrt_ui.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:3k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* mpf_sqrt_ui -- Compute the square root of an unsigned integer.
  2. Copyright 1993, 1994, 1996, 2000, 2001, 2004, 2005 Free Software Foundation,
  3. Inc.
  4. This file is part of the GNU MP Library.
  5. The GNU MP Library is free software; you can redistribute it and/or modify
  6. it under the terms of the GNU Lesser General Public License as published by
  7. the Free Software Foundation; either version 3 of the License, or (at your
  8. option) any later version.
  9. The GNU MP Library is distributed in the hope that it will be useful, but
  10. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  11. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  12. License for more details.
  13. You should have received a copy of the GNU Lesser General Public License
  14. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  15. #include <stdio.h> /* for NULL */
  16. #include "gmp.h"
  17. #include "gmp-impl.h"
  18. /* As usual the aim is to produce PREC(r) limbs of result with the high limb
  19.    non-zero.  That high limb will end up floor(sqrt(u)), and limbs below are
  20.    produced by padding the input with zeros, two for each desired result
  21.    limb, being 2*(prec-1) for a total 2*prec-1 limbs passed to mpn_sqrtrem.
  22.    The way mpn_sqrtrem calculates floor(sqrt(x)) ensures the root is correct
  23.    to the intended accuracy, ie. truncated to prec limbs.
  24.    With nails, u might be two limbs, in which case a total 2*prec limbs is
  25.    passed to mpn_sqrtrem (still giving a prec limb result).  If uhigh is
  26.    zero we adjust back to 2*prec-1, since mpn_sqrtrem requires the high
  27.    non-zero.  2*prec limbs are always allocated, even when uhigh is zero, so
  28.    the store of uhigh can be done without a conditional.
  29.    u==0 is a special case so the rest of the code can assume the result is
  30.    non-zero (ie. will have a non-zero high limb on the result).
  31.    Not done:
  32.    No attempt is made to identify perfect squares.  It's considered this can
  33.    be left to an application if it might occur with any frequency.  As it
  34.    stands, mpn_sqrtrem does its normal amount of work on a perfect square
  35.    followed by zero limbs, though of course only an mpn_sqrtrem1 would be
  36.    actually needed.  We also end up leaving our mpf result with lots of low
  37.    trailing zeros, slowing down subsequent operations.
  38.    We're not aware of any optimizations that can be made using the fact the
  39.    input has lots of trailing zeros (apart from the perfect square
  40.    case).  */
  41. /* 1 if we (might) need two limbs for u */
  42. #define U2   (GMP_NUMB_BITS < BITS_PER_ULONG)
  43. void
  44. mpf_sqrt_ui (mpf_ptr r, unsigned long int u)
  45. {
  46.   mp_size_t rsize, zeros;
  47.   mp_ptr tp;
  48.   mp_size_t prec;
  49.   TMP_DECL;
  50.   if (UNLIKELY (u == 0))
  51.     {
  52.       r->_mp_size = 0;
  53.       r->_mp_exp = 0;
  54.       return;
  55.     }
  56.   TMP_MARK;
  57.   prec = r->_mp_prec;
  58.   zeros = 2 * prec - 2;
  59.   rsize = zeros + 1 + U2;
  60.   tp = TMP_ALLOC_LIMBS (rsize);
  61.   MPN_ZERO (tp, zeros);
  62.   tp[zeros] = u & GMP_NUMB_MASK;
  63. #if U2
  64.   {
  65.     mp_limb_t uhigh = u >> GMP_NUMB_BITS;
  66.     tp[zeros + 1] = uhigh;
  67.     rsize -= (uhigh == 0);
  68.   }
  69. #endif
  70.   mpn_sqrtrem (r->_mp_d, NULL, tp, rsize);
  71.   r->_mp_size = prec;
  72.   r->_mp_exp = 1;
  73.   TMP_FREE;
  74. }