1/* Test mpf_sqrt_ui.
2
3Copyright 2004, 2015 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library test suite.
6
7The GNU MP Library test suite is free software; you can redistribute it
8and/or modify it under the terms of the GNU General Public License as
9published by the Free Software Foundation; either version 3 of the License,
10or (at your option) any later version.
11
12The GNU MP Library test suite is distributed in the hope that it will be
13useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15Public License for more details.
16
17You should have received a copy of the GNU General Public License along with
18the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
19
20#include <stdio.h>
21#include <stdlib.h>
22
23#include "gmp-impl.h"
24#include "tests.h"
25
26
27void
28check_rand (void)
29{
30  unsigned long      max_prec = 15;
31  unsigned long      min_prec = __GMPF_BITS_TO_PREC (1);
32  gmp_randstate_ptr  rands = RANDS;
33  unsigned long      x, prec;
34  mpf_t              r, s;
35  int                i;
36
37  mpf_init (r);
38  mpf_init (s);
39  refmpf_set_prec_limbs (s, 2*max_prec+10);
40
41  for (x = 0; x < 2; x++)
42    {
43      mpf_sqrt_ui (r, x);
44      MPF_CHECK_FORMAT (r);
45      if (mpf_cmp_ui (r, x) != 0)
46	{
47	  printf    ("mpf_sqrt_ui wrong for special case:\n");
48          printf    ("  x=%lu\n", x);
49          mpf_trace ("  r", r);
50	  abort ();
51	}
52    }
53
54  for (i = 0; i < 50; i++)
55    {
56      /* input, a random non-zero ulong, exponentially distributed */
57      do {
58        x = gmp_urandomb_ui (rands,
59                             gmp_urandomm_ui (rands, BITS_PER_ULONG) + 1);
60      } while (x <= 1);
61
62      /* result precision */
63      prec = gmp_urandomm_ui (rands, max_prec-min_prec) + min_prec;
64      refmpf_set_prec_limbs (r, prec);
65
66      mpf_sqrt_ui (r, x);
67      MPF_CHECK_FORMAT (r);
68
69      /* Expect to prec limbs of result.
70         In the current implementation there's no stripping of low zero
71         limbs in mpf_sqrt_ui, not even on perfect squares, so size should
72         be exactly prec.  */
73      if (SIZ(r) != prec)
74        {
75          printf ("mpf_sqrt_ui result not enough result limbs\n");
76          printf    ("  x=%lu\n", x);
77          printf    ("  want prec=%lu\n", prec);
78          mpf_trace ("  r", r);
79          printf    ("  r size %ld\n", (long) SIZ(r));
80          printf    ("  r prec %ld\n", (long) PREC(r));
81          abort ();
82        }
83
84      /* Must have r^2 <= x, since r has been truncated. */
85      mpf_mul (s, r, r);
86      if (! (mpf_cmp_ui (s, x) <= 0))
87        {
88          printf    ("mpf_sqrt_ui result too big\n");
89          printf    ("  x=%lu\n", x);
90          printf    ("  want prec=%lu\n", prec);
91          mpf_trace ("  r", r);
92          mpf_trace ("  s", s);
93          abort ();
94        }
95
96      /* Must have (r+ulp)^2 > x.
97         No overflow from refmpf_add_ulp since r is only prec limbs. */
98      refmpf_add_ulp (r);
99      mpf_mul (s, r, r);
100      if (! (mpf_cmp_ui (s, x) > 0))
101        {
102          printf    ("mpf_sqrt_ui result too small\n");
103          printf    ("  x=%lu\n", x);
104          printf    ("  want prec=%lu\n", prec);
105          mpf_trace ("  r+ulp", r);
106          mpf_trace ("  s", s);
107          abort ();
108        }
109    }
110
111  mpf_clear (r);
112  mpf_clear (s);
113}
114
115int
116main (int argc, char **argv)
117{
118  tests_start ();
119  mp_trace_base = -16;
120
121  check_rand ();
122
123  tests_end ();
124  exit (0);
125}
126