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.h"
24 #include "gmp-impl.h"
25 #include "tests.h"
26
27 #ifndef SIZE
28 #define SIZE 16
29 #endif
30
31 void
check_rand1(int argc,char ** argv)32 check_rand1 (int argc, char **argv)
33 {
34 mp_size_t size;
35 mp_exp_t exp;
36 int reps = 20000;
37 int i;
38 mpf_t x, y, y2;
39 mp_size_t bprec = 100;
40 mpf_t rerr, max_rerr, limit_rerr;
41
42 if (argc > 1)
43 {
44 reps = strtol (argv[1], 0, 0);
45 if (argc > 2)
46 bprec = strtol (argv[2], 0, 0);
47 }
48
49 mpf_set_default_prec (bprec);
50
51 mpf_init_set_ui (limit_rerr, 1);
52 mpf_div_2exp (limit_rerr, limit_rerr, bprec);
53 #if VERBOSE
54 mpf_dump (limit_rerr);
55 #endif
56 mpf_init (rerr);
57 mpf_init_set_ui (max_rerr, 0);
58
59 mpf_init (x);
60 mpf_init (y);
61 mpf_init (y2);
62 for (i = 0; i < reps; i++)
63 {
64 size = urandom () % SIZE;
65 exp = urandom () % SIZE;
66 mpf_random2 (x, size, exp);
67
68 mpf_sqrt (y, x);
69 MPF_CHECK_FORMAT (y);
70 mpf_mul (y2, y, y);
71
72 mpf_reldiff (rerr, x, y2);
73 if (mpf_cmp (rerr, max_rerr) > 0)
74 {
75 mpf_set (max_rerr, rerr);
76 #if VERBOSE
77 mpf_dump (max_rerr);
78 #endif
79 if (mpf_cmp (rerr, limit_rerr) > 0)
80 {
81 printf ("ERROR after %d tests\n", i);
82 printf (" x = "); mpf_dump (x);
83 printf (" y = "); mpf_dump (y);
84 printf (" y2 = "); mpf_dump (y2);
85 printf (" rerr = "); mpf_dump (rerr);
86 printf (" limit_rerr = "); mpf_dump (limit_rerr);
87 printf ("in hex:\n");
88 mp_trace_base = 16;
89 mpf_trace (" x ", x);
90 mpf_trace (" y ", y);
91 mpf_trace (" y2 ", y2);
92 mpf_trace (" rerr ", rerr);
93 mpf_trace (" limit_rerr", limit_rerr);
94 abort ();
95 }
96 }
97 }
98
99 mpf_clear (limit_rerr);
100 mpf_clear (rerr);
101 mpf_clear (max_rerr);
102
103 mpf_clear (x);
104 mpf_clear (y);
105 mpf_clear (y2);
106 }
107
108 void
check_rand2(void)109 check_rand2 (void)
110 {
111 unsigned long max_prec = 20;
112 unsigned long min_prec = __GMPF_BITS_TO_PREC (1);
113 gmp_randstate_ptr rands = RANDS;
114 unsigned long x_prec, r_prec;
115 mpf_t x, r, s;
116 int i;
117
118 mpf_init (x);
119 mpf_init (r);
120 mpf_init (s);
121 refmpf_set_prec_limbs (s, 2*max_prec+10);
122
123 for (i = 0; i < 500; i++)
124 {
125 /* input precision */
126 x_prec = gmp_urandomm_ui (rands, max_prec-min_prec) + min_prec;
127 refmpf_set_prec_limbs (x, x_prec);
128
129 /* result precision */
130 r_prec = gmp_urandomm_ui (rands, max_prec-min_prec) + min_prec;
131 refmpf_set_prec_limbs (r, r_prec);
132
133 mpf_random2 (x, x_prec, 1000);
134
135 mpf_sqrt (r, x);
136 MPF_CHECK_FORMAT (r);
137
138 /* Expect to prec limbs of result.
139 In the current implementation there's no stripping of low zero
140 limbs in mpf_sqrt, so size should be exactly prec. */
141 if (SIZ(r) != r_prec)
142 {
143 printf ("mpf_sqrt wrong number of result limbs\n");
144 mpf_trace (" x", x);
145 mpf_trace (" r", r);
146 printf (" r_prec=%lu\n", r_prec);
147 printf (" SIZ(r) %ld\n", (long) SIZ(r));
148 printf (" PREC(r) %ld\n", (long) PREC(r));
149 abort ();
150 }
151
152 /* Must have r^2 <= x, since r has been truncated. */
153 mpf_mul (s, r, r);
154 if (! (mpf_cmp (s, x) <= 0))
155 {
156 printf ("mpf_sqrt result too big\n");
157 mpf_trace (" x", x);
158 printf (" r_prec=%lu\n", r_prec);
159 mpf_trace (" r", r);
160 mpf_trace (" s", s);
161 abort ();
162 }
163
164 /* Must have (r+ulp)^2 > x, or else r is too small. */
165 refmpf_add_ulp (r);
166 mpf_mul (s, r, r);
167 if (! (mpf_cmp (s, x) > 0))
168 {
169 printf ("mpf_sqrt result too small\n");
170 mpf_trace (" x", x);
171 printf (" r_prec=%lu\n", r_prec);
172 mpf_trace (" r+ulp", r);
173 mpf_trace (" s", s);
174 abort ();
175 }
176 }
177
178 mpf_clear (x);
179 mpf_clear (r);
180 mpf_clear (s);
181 }
182
183 int
main(int argc,char ** argv)184 main (int argc, char **argv)
185 {
186 tests_start ();
187 mp_trace_base = -16;
188
189 check_rand1 (argc, argv);
190 check_rand2 ();
191
192 tests_end ();
193 exit (0);
194 }
195