github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mpf/sqrt_ui.c (about)

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