← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/gmp/mpn/generic/rootrem.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/mpn/generic/rootrem.c
7
7
   FUTURE GNU MP RELEASE.
8
8
 
9
9
 
10
 
Copyright 2002 Free Software Foundation, Inc.
 
10
Copyright 2002, 2005 Free Software Foundation, Inc.
11
11
 
12
12
This file is part of the GNU MP Library.
13
13
 
22
22
License for more details.
23
23
 
24
24
You should have received a copy of the GNU Lesser General Public License
25
 
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
26
 
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
27
 
MA 02111-1307, USA. */
 
25
along with the GNU MP Library; see the file COPYING.LIB.  If not, write
 
26
to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 
27
Boston, MA 02110-1301, USA. */
28
28
 
29
29
/*
30
30
  We use Newton's method to compute the root of a:
62
62
  unsigned int cnt;
63
63
  mp_size_t i;
64
64
  unsigned long int n_valid_bits, adj;
65
 
  TMP_DECL (marker);
 
65
  TMP_DECL;
66
66
 
67
 
  TMP_MARK (marker);
 
67
  TMP_MARK;
68
68
 
69
69
  /* The extra factor 1.585 = log(3)/log(2) here is for the worst case
70
70
     overestimate of the root, i.e., when the code rounds a root that is
71
 
     2+epsilon to 3, and the powers this to a potentially huge power.  We
 
71
     2+epsilon to 3, and then powers this to a potentially huge power.  We
72
72
     could generalize the code for detecting root=1 a few lines below to deal
73
73
     with xnb <= k, for some small k.  For example, when xnb <= 2, meaning
74
74
     the root should be 1, 2, or 3, we could replace this factor by the much
88
88
      mpn_sub_1 (remp, up, un, (mp_limb_t) 1);
89
89
      MPN_NORMALIZE (remp, un);
90
90
      rootp[0] = 1;
91
 
      TMP_FREE (marker);
 
91
      TMP_FREE;
92
92
      return un;
93
93
    }
94
94
 
95
95
  xn = (xnb + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
96
96
 
97
 
  qp = TMP_ALLOC_LIMBS (un + 1);
 
97
  qp = TMP_ALLOC_LIMBS (PP_ALLOC);
98
98
  xp = TMP_ALLOC_LIMBS (xn + 1);
99
99
 
100
100
  /* Set initial root to only ones.  This is an overestimate of the actual root
125
125
  adj = n_valid_bits - 1;
126
126
 
127
127
  /* Newton loop.  Converges downwards towards root(U,nth).  Currently we use
128
 
   full precision from iteration 1.  Clearly, we should use just n_valid_bits
129
 
   of precision in each step, and thus save most of the computations.  */
 
128
     full precision from iteration 1.  Clearly, we should use just n_valid_bits
 
129
     of precision in each step, and thus save most of the computations.  */
130
130
  while (n_valid_bits <= xnb)
131
131
    {
132
132
      mp_limb_t cy;
172
172
  mpn_sub (remp, up, un, pp, pn);
173
173
  MPN_NORMALIZE (remp, un);
174
174
  MPN_COPY (rootp, xp, xn);
175
 
  TMP_FREE (marker);
 
175
  TMP_FREE;
176
176
  return un;
177
177
}