← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/gmp/mpf/sqrt_ui.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2007-04-09 11:51:51 UTC
  • mfrom: (1.1.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20070409115151-ql8cr0kalzx1jmla
Tags: 0.9i-20070324-2
Upload to unstable. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gmp/mpf/sqrt_ui.c
1
1
/* mpf_sqrt_ui -- Compute the square root of an unsigned integer.
2
2
 
3
 
Copyright 1993, 1994, 1996, 2000, 2001 Free Software Foundation, Inc.
 
3
Copyright 1993, 1994, 1996, 2000, 2001, 2004, 2005 Free Software Foundation,
 
4
Inc.
4
5
 
5
6
This file is part of the GNU MP Library.
6
7
 
16
17
 
17
18
You should have received a copy of the GNU Lesser General Public License
18
19
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
19
 
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20
 
MA 02111-1307, USA. */
 
20
the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 
21
MA 02110-1301, USA. */
21
22
 
22
23
#include <stdio.h> /* for NULL */
23
24
#include "gmp.h"
24
25
#include "gmp-impl.h"
25
26
 
 
27
 
 
28
/* As usual the aim is to produce PREC(r) limbs of result with the high limb
 
29
   non-zero.  That high limb will end up floor(sqrt(u)), and limbs below are
 
30
   produced by padding the input with zeros, two for each desired result
 
31
   limb, being 2*(prec-1) for a total 2*prec-1 limbs passed to mpn_sqrtrem.
 
32
   The way mpn_sqrtrem calculates floor(sqrt(x)) ensures the root is correct
 
33
   to the intended accuracy, ie. truncated to prec limbs.
 
34
 
 
35
   With nails, u might be two limbs, in which case a total 2*prec limbs is
 
36
   passed to mpn_sqrtrem (still giving a prec limb result).  If uhigh is
 
37
   zero we adjust back to 2*prec-1, since mpn_sqrtrem requires the high
 
38
   non-zero.  2*prec limbs are always allocated, even when uhigh is zero, so
 
39
   the store of uhigh can be done without a conditional.
 
40
 
 
41
   u==0 is a special case so the rest of the code can assume the result is
 
42
   non-zero (ie. will have a non-zero high limb on the result).
 
43
 
 
44
   Not done:
 
45
 
 
46
   No attempt is made to identify perfect squares.  It's considered this can
 
47
   be left to an application if it might occur with any frequency.  As it
 
48
   stands, mpn_sqrtrem does its normal amount of work on a perfect square
 
49
   followed by zero limbs, though of course only an mpn_sqrtrem1 would be
 
50
   actually needed.  We also end up leaving our mpf result with lots of low
 
51
   trailing zeros, slowing down subsequent operations.
 
52
 
 
53
   We're not aware of any optimizations that can be made using the fact the
 
54
   input has lots of trailing zeros (apart from the perfect square
 
55
   case).  */
 
56
 
 
57
 
 
58
/* 1 if we (might) need two limbs for u */
 
59
#define U2   (GMP_NUMB_BITS < BITS_PER_ULONG)
 
60
 
26
61
void
27
62
mpf_sqrt_ui (mpf_ptr r, unsigned long int u)
28
63
{
29
 
  mp_size_t rsize;
 
64
  mp_size_t rsize, zeros;
30
65
  mp_ptr tp;
31
66
  mp_size_t prec;
32
 
  TMP_DECL (marker);
 
67
  TMP_DECL;
33
68
 
34
 
  if (u == 0)
 
69
  if (UNLIKELY (u == 0))
35
70
    {
36
71
      r->_mp_size = 0;
37
72
      r->_mp_exp = 0;
38
73
      return;
39
74
    }
40
75
 
41
 
  TMP_MARK (marker);
 
76
  TMP_MARK;
42
77
 
43
78
  prec = r->_mp_prec;
44
 
  rsize = 2 * prec + 1;
 
79
  zeros = 2 * prec - 2;
 
80
  rsize = zeros + 1 + U2;
45
81
 
46
82
  tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
47
83
 
48
 
  MPN_ZERO (tp, rsize - 1);
49
 
  tp[rsize - 1] = u;
 
84
  MPN_ZERO (tp, zeros);
 
85
  tp[zeros] = u & GMP_NUMB_MASK;
 
86
 
 
87
#if U2
 
88
  {
 
89
    mp_limb_t uhigh = u >> GMP_NUMB_BITS;
 
90
    tp[zeros + 1] = uhigh;
 
91
    rsize -= (uhigh == 0);
 
92
  }
 
93
#endif
50
94
 
51
95
  mpn_sqrtrem (r->_mp_d, NULL, tp, rsize);
52
96
 
53
 
  r->_mp_size = prec + 1;
 
97
  r->_mp_size = prec;
54
98
  r->_mp_exp = 1;
55
 
  TMP_FREE (marker);
 
99
  TMP_FREE;
56
100
}