1 /* Test mpf_sqrt, mpf_mul.
2
3 Copyright 1996, 2001, 2004 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library test suite.
6
7 The GNU MP Library test suite is free software; you can redistribute it
8 and/or modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11
12 The GNU MP Library test suite is distributed in the hope that it will be
13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15 Public License for more details.
16
17 You should have received a copy of the GNU General Public License along with
18 the 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 #ifndef SIZE
27 #define SIZE 16
28 #endif
29
30 void
31 check_rand1 (int argc, char **argv)
32 {
33 mp_size_t size;
34 mp_exp_t exp;
35 int reps = 20000;
36 int i;
37 mpf_t x, y, y2;
38 mp_size_t bprec = 100;
39 mpf_t rerr, max_rerr, limit_rerr;
40
41 if (argc > 1)
42 {
43 reps = strtol (argv[1], 0, 0);
44 if (argc > 2)
45 bprec = strtol (argv[2], 0, 0);
46 }
47
48 mpf_set_default_prec (bprec);
49
50 mpf_init_set_ui (limit_rerr, 1);
51 mpf_div_2exp (limit_rerr, limit_rerr, bprec);
52 #if VERBOSE
53 mpf_dump (limit_rerr);
54 #endif
55 mpf_init (rerr);
56 mpf_init_set_ui (max_rerr, 0);
57
58 mpf_init (x);
59 mpf_init (y);
60 mpf_init (y2);
61 for (i = 0; i < reps; i++)
62 {
63 size = urandom () % SIZE;
64 exp = urandom () % SIZE;
65 mpf_random2 (x, size, exp);
66
67 mpf_sqrt (y, x);
68 MPF_CHECK_FORMAT (y);
69 mpf_mul (y2, y, y);
70
71 mpf_reldiff (rerr, x, y2);
72 if (mpf_cmp (rerr, max_rerr) > 0)
73 {
74 mpf_set (max_rerr, rerr);
75 #if VERBOSE
76 mpf_dump (max_rerr);
77 #endif
78 if (mpf_cmp (rerr, limit_rerr) > 0)
79 {
80 printf ("ERROR after %d tests\n", i);
81 printf (" x = "); mpf_dump (x);
82 printf (" y = "); mpf_dump (y);
83 printf (" y2 = "); mpf_dump (y2);
84 printf (" rerr = "); mpf_dump (rerr);
85 printf (" limit_rerr = "); mpf_dump (limit_rerr);
86 printf ("in hex:\n");
87 mp_trace_base = 16;
88 mpf_trace (" x ", x);
89 mpf_trace (" y ", y);
90 mpf_trace (" y2 ", y2);
91 mpf_trace (" rerr ", rerr);
92 mpf_trace (" limit_rerr", limit_rerr);
93 abort ();
94 }
95 }
96 }
97
98 mpf_clear (limit_rerr);
99 mpf_clear (rerr);
100 mpf_clear (max_rerr);
101
102 mpf_clear (x);
103 mpf_clear (y);
104 mpf_clear (y2);
105 }
106
107 void
108 check_rand2 (void)
109 {
110 unsigned long max_prec = 20;
111 unsigned long min_prec = __GMPF_BITS_TO_PREC (1);
112 gmp_randstate_ptr rands = RANDS;
113 unsigned long x_prec, r_prec;
114 mpf_t x, r, s;
115 int i;
116
117 mpf_init (x);
118 mpf_init (r);
119 mpf_init (s);
120 refmpf_set_prec_limbs (s, 2*max_prec+10);
121
122 for (i = 0; i < 500; i++)
123 {
124 /* input precision */
125 x_prec = gmp_urandomm_ui (rands, max_prec-min_prec) + min_prec;
126 refmpf_set_prec_limbs (x, x_prec);
127
128 /* result precision */
129 r_prec = gmp_urandomm_ui (rands, max_prec-min_prec) + min_prec;
130 refmpf_set_prec_limbs (r, r_prec);
131
132 mpf_random2 (x, x_prec, 1000);
133
134 mpf_sqrt (r, x);
135 MPF_CHECK_FORMAT (r);
136
137 /* Expect to prec limbs of result.
138 In the current implementation there's no stripping of low zero
139 limbs in mpf_sqrt, so size should be exactly prec. */
140 if (SIZ(r) != r_prec)
141 {
142 printf ("mpf_sqrt wrong number of result limbs\n");
143 mpf_trace (" x", x);
144 mpf_trace (" r", r);
145 printf (" r_prec=%lu\n", r_prec);
146 printf (" SIZ(r) %ld\n", (long) SIZ(r));
147 printf (" PREC(r) %ld\n", (long) PREC(r));
148 abort ();
149 }
150
151 /* Must have r^2 <= x, since r has been truncated. */
152 mpf_mul (s, r, r);
153 if (! (mpf_cmp (s, x) <= 0))
154 {
155 printf ("mpf_sqrt result too big\n");
156 mpf_trace (" x", x);
157 printf (" r_prec=%lu\n", r_prec);
158 mpf_trace (" r", r);
159 mpf_trace (" s", s);
160 abort ();
161 }
162
163 /* Must have (r+ulp)^2 > x, or else r is too small. */
164 refmpf_add_ulp (r);
165 mpf_mul (s, r, r);
166 if (! (mpf_cmp (s, x) > 0))
167 {
168 printf ("mpf_sqrt result too small\n");
169 mpf_trace (" x", x);
170 printf (" r_prec=%lu\n", r_prec);
171 mpf_trace (" r+ulp", r);
172 mpf_trace (" s", s);
173 abort ();
174 }
175 }
176
177 mpf_clear (x);
178 mpf_clear (r);
179 mpf_clear (s);
180 }
181
182 int
183 main (int argc, char **argv)
184 {
185 tests_start ();
186 mp_trace_base = -16;
187
188 check_rand1 (argc, argv);
189 check_rand2 ();
190
191 tests_end ();
192 exit (0);
193 }