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 } 194