github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/tests/refmpf.c (about) 1 /* Reference floating point routines. 2 3 Copyright 1996, 2001, 2004, 2005 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 28 void 29 refmpf_add (mpf_ptr w, mpf_srcptr u, mpf_srcptr v) 30 { 31 mp_size_t hi, lo, size; 32 mp_ptr ut, vt, wt; 33 int neg; 34 mp_exp_t exp; 35 mp_limb_t cy; 36 TMP_DECL; 37 38 TMP_MARK; 39 40 if (SIZ (u) == 0) 41 { 42 size = ABSIZ (v); 43 wt = TMP_ALLOC_LIMBS (size + 1); 44 MPN_COPY (wt, PTR (v), size); 45 exp = EXP (v); 46 neg = SIZ (v) < 0; 47 goto done; 48 } 49 if (SIZ (v) == 0) 50 { 51 size = ABSIZ (u); 52 wt = TMP_ALLOC_LIMBS (size + 1); 53 MPN_COPY (wt, PTR (u), size); 54 exp = EXP (u); 55 neg = SIZ (u) < 0; 56 goto done; 57 } 58 if ((SIZ (u) ^ SIZ (v)) < 0) 59 { 60 mpf_t tmp; 61 SIZ (tmp) = -SIZ (v); 62 EXP (tmp) = EXP (v); 63 PTR (tmp) = PTR (v); 64 refmpf_sub (w, u, tmp); 65 return; 66 } 67 neg = SIZ (u) < 0; 68 69 /* Compute the significance of the hi and lo end of the result. */ 70 hi = MAX (EXP (u), EXP (v)); 71 lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v)); 72 size = hi - lo; 73 ut = TMP_ALLOC_LIMBS (size + 1); 74 vt = TMP_ALLOC_LIMBS (size + 1); 75 wt = TMP_ALLOC_LIMBS (size + 1); 76 MPN_ZERO (ut, size); 77 MPN_ZERO (vt, size); 78 {int off; 79 off = size + (EXP (u) - hi) - ABSIZ (u); 80 MPN_COPY (ut + off, PTR (u), ABSIZ (u)); 81 off = size + (EXP (v) - hi) - ABSIZ (v); 82 MPN_COPY (vt + off, PTR (v), ABSIZ (v)); 83 } 84 85 cy = mpn_add_n (wt, ut, vt, size); 86 wt[size] = cy; 87 size += cy; 88 exp = hi + cy; 89 90 done: 91 if (size > PREC (w)) 92 { 93 wt += size - PREC (w); 94 size = PREC (w); 95 } 96 MPN_COPY (PTR (w), wt, size); 97 SIZ (w) = neg == 0 ? size : -size; 98 EXP (w) = exp; 99 TMP_FREE; 100 } 101 102 103 /* Add 1 "unit in last place" (ie. in the least significant limb) to f. 104 f cannot be zero, since that has no well-defined "last place". 105 106 This routine is designed for use in cases where we pay close attention to 107 the size of the data value and are using that (and the exponent) to 108 indicate the accurate part of a result, or similar. For this reason, if 109 there's a carry out we don't store 1 and adjust the exponent, we just 110 leave 100..00. We don't even adjust if there's a carry out of prec+1 111 limbs, but instead give up in that case (which we intend shouldn't arise 112 in normal circumstances). */ 113 114 void 115 refmpf_add_ulp (mpf_ptr f) 116 { 117 mp_ptr fp = PTR(f); 118 mp_size_t fsize = SIZ(f); 119 mp_size_t abs_fsize = ABSIZ(f); 120 mp_limb_t c; 121 122 if (fsize == 0) 123 { 124 printf ("Oops, refmpf_add_ulp called with f==0\n"); 125 abort (); 126 } 127 128 c = refmpn_add_1 (fp, fp, abs_fsize, CNST_LIMB(1)); 129 if (c != 0) 130 { 131 if (abs_fsize >= PREC(f) + 1) 132 { 133 printf ("Oops, refmpf_add_ulp carried out of prec+1 limbs\n"); 134 abort (); 135 } 136 137 fp[abs_fsize] = c; 138 abs_fsize++; 139 SIZ(f) = (fsize > 0 ? abs_fsize : - abs_fsize); 140 EXP(f)++; 141 } 142 } 143 144 /* Fill f with size limbs of the given value, setup as an integer. */ 145 void 146 refmpf_fill (mpf_ptr f, mp_size_t size, mp_limb_t value) 147 { 148 ASSERT (size >= 0); 149 size = MIN (PREC(f) + 1, size); 150 SIZ(f) = size; 151 EXP(f) = size; 152 refmpn_fill (PTR(f), size, value); 153 } 154 155 /* Strip high zero limbs from the f data, adjusting exponent accordingly. */ 156 void 157 refmpf_normalize (mpf_ptr f) 158 { 159 while (SIZ(f) != 0 && PTR(f)[ABSIZ(f)-1] == 0) 160 { 161 SIZ(f) = (SIZ(f) >= 0 ? SIZ(f)-1 : SIZ(f)+1); 162 EXP(f) --; 163 } 164 if (SIZ(f) == 0) 165 EXP(f) = 0; 166 } 167 168 /* refmpf_set_overlap sets up dst as a copy of src, but with PREC(dst) 169 unchanged, in preparation for an overlap test. 170 171 The full value of src is copied, and the space at PTR(dst) is extended as 172 necessary. The way PREC(dst) is unchanged is as per an mpf_set_prec_raw. 173 The return value is the new PTR(dst) space precision, in bits, ready for 174 a restoring mpf_set_prec_raw before mpf_clear. */ 175 176 unsigned long 177 refmpf_set_overlap (mpf_ptr dst, mpf_srcptr src) 178 { 179 mp_size_t dprec = PREC(dst); 180 mp_size_t ssize = ABSIZ(src); 181 unsigned long ret; 182 183 refmpf_set_prec_limbs (dst, (unsigned long) MAX (dprec, ssize)); 184 mpf_set (dst, src); 185 186 ret = mpf_get_prec (dst); 187 PREC(dst) = dprec; 188 return ret; 189 } 190 191 /* Like mpf_set_prec, but taking a precision in limbs. 192 PREC(f) ends up as the given "prec" value. */ 193 void 194 refmpf_set_prec_limbs (mpf_ptr f, unsigned long prec) 195 { 196 mpf_set_prec (f, __GMPF_PREC_TO_BITS (prec)); 197 } 198 199 200 void 201 refmpf_sub (mpf_ptr w, mpf_srcptr u, mpf_srcptr v) 202 { 203 mp_size_t hi, lo, size; 204 mp_ptr ut, vt, wt; 205 int neg; 206 mp_exp_t exp; 207 TMP_DECL; 208 209 TMP_MARK; 210 211 if (SIZ (u) == 0) 212 { 213 size = ABSIZ (v); 214 wt = TMP_ALLOC_LIMBS (size + 1); 215 MPN_COPY (wt, PTR (v), size); 216 exp = EXP (v); 217 neg = SIZ (v) > 0; 218 goto done; 219 } 220 if (SIZ (v) == 0) 221 { 222 size = ABSIZ (u); 223 wt = TMP_ALLOC_LIMBS (size + 1); 224 MPN_COPY (wt, PTR (u), size); 225 exp = EXP (u); 226 neg = SIZ (u) < 0; 227 goto done; 228 } 229 if ((SIZ (u) ^ SIZ (v)) < 0) 230 { 231 mpf_t tmp; 232 SIZ (tmp) = -SIZ (v); 233 EXP (tmp) = EXP (v); 234 PTR (tmp) = PTR (v); 235 refmpf_add (w, u, tmp); 236 if (SIZ (u) < 0) 237 mpf_neg (w, w); 238 return; 239 } 240 neg = SIZ (u) < 0; 241 242 /* Compute the significance of the hi and lo end of the result. */ 243 hi = MAX (EXP (u), EXP (v)); 244 lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v)); 245 size = hi - lo; 246 ut = TMP_ALLOC_LIMBS (size + 1); 247 vt = TMP_ALLOC_LIMBS (size + 1); 248 wt = TMP_ALLOC_LIMBS (size + 1); 249 MPN_ZERO (ut, size); 250 MPN_ZERO (vt, size); 251 {int off; 252 off = size + (EXP (u) - hi) - ABSIZ (u); 253 MPN_COPY (ut + off, PTR (u), ABSIZ (u)); 254 off = size + (EXP (v) - hi) - ABSIZ (v); 255 MPN_COPY (vt + off, PTR (v), ABSIZ (v)); 256 } 257 258 if (mpn_cmp (ut, vt, size) >= 0) 259 mpn_sub_n (wt, ut, vt, size); 260 else 261 { 262 mpn_sub_n (wt, vt, ut, size); 263 neg ^= 1; 264 } 265 exp = hi; 266 while (size != 0 && wt[size - 1] == 0) 267 { 268 size--; 269 exp--; 270 } 271 272 done: 273 if (size > PREC (w)) 274 { 275 wt += size - PREC (w); 276 size = PREC (w); 277 } 278 MPN_COPY (PTR (w), wt, size); 279 SIZ (w) = neg == 0 ? size : -size; 280 EXP (w) = exp; 281 TMP_FREE; 282 } 283 284 285 /* Validate got by comparing to want. Return 1 if good, 0 if bad. 286 287 The data in got is compared to that in want, up to either PREC(got) limbs 288 or the size of got, whichever is bigger. Clearly we always demand 289 PREC(got) of accuracy, but we go further and say that if got is bigger 290 then any extra must be correct too. 291 292 want needs to have enough data to allow this comparison. The size in 293 want doesn't have to be that big though, if it's smaller then further low 294 limbs are taken to be zero. 295 296 This validation approach is designed to allow some flexibility in exactly 297 how much data is generated by an mpf function, ie. either prec or prec+1 298 limbs. We don't try to make a reference function that emulates that same 299 size decision, instead the idea is for a validation function to generate 300 at least as much data as the real function, then compare. */ 301 302 int 303 refmpf_validate (const char *name, mpf_srcptr got, mpf_srcptr want) 304 { 305 int bad = 0; 306 mp_size_t gsize, wsize, cmpsize, i; 307 mp_srcptr gp, wp; 308 mp_limb_t glimb, wlimb; 309 310 MPF_CHECK_FORMAT (got); 311 312 if (EXP (got) != EXP (want)) 313 { 314 printf ("%s: wrong exponent\n", name); 315 bad = 1; 316 } 317 318 gsize = SIZ (got); 319 wsize = SIZ (want); 320 if ((gsize < 0 && wsize > 0) || (gsize > 0 && wsize < 0)) 321 { 322 printf ("%s: wrong sign\n", name); 323 bad = 1; 324 } 325 326 gsize = ABS (gsize); 327 wsize = ABS (wsize); 328 329 /* most significant limb of respective data */ 330 gp = PTR (got) + gsize - 1; 331 wp = PTR (want) + wsize - 1; 332 333 /* compare limb data */ 334 cmpsize = MAX (PREC (got), gsize); 335 for (i = 0; i < cmpsize; i++) 336 { 337 glimb = (i < gsize ? gp[-i] : 0); 338 wlimb = (i < wsize ? wp[-i] : 0); 339 340 if (glimb != wlimb) 341 { 342 printf ("%s: wrong data starting at index %ld from top\n", 343 name, (long) i); 344 bad = 1; 345 break; 346 } 347 } 348 349 if (bad) 350 { 351 printf (" prec %d\n", PREC(got)); 352 printf (" exp got %ld\n", (long) EXP(got)); 353 printf (" exp want %ld\n", (long) EXP(want)); 354 printf (" size got %d\n", SIZ(got)); 355 printf (" size want %d\n", SIZ(want)); 356 printf (" limbs (high to low)\n"); 357 printf (" got "); 358 for (i = ABSIZ(got)-1; i >= 0; i--) 359 { 360 gmp_printf ("%MX", PTR(got)[i]); 361 if (i != 0) 362 printf (","); 363 } 364 printf ("\n"); 365 printf (" want "); 366 for (i = ABSIZ(want)-1; i >= 0; i--) 367 { 368 gmp_printf ("%MX", PTR(want)[i]); 369 if (i != 0) 370 printf (","); 371 } 372 printf ("\n"); 373 return 0; 374 } 375 376 return 1; 377 } 378 379 380 int 381 refmpf_validate_division (const char *name, mpf_srcptr got, 382 mpf_srcptr n, mpf_srcptr d) 383 { 384 mp_size_t nsize, dsize, sign, prec, qsize, tsize; 385 mp_srcptr np, dp; 386 mp_ptr tp, qp, rp; 387 mpf_t want; 388 int ret; 389 390 nsize = SIZ (n); 391 dsize = SIZ (d); 392 ASSERT_ALWAYS (dsize != 0); 393 394 sign = nsize ^ dsize; 395 nsize = ABS (nsize); 396 dsize = ABS (dsize); 397 398 np = PTR (n); 399 dp = PTR (d); 400 prec = PREC (got); 401 402 EXP (want) = EXP (n) - EXP (d) + 1; 403 404 qsize = prec + 2; /* at least prec+1 limbs, after high zero */ 405 tsize = qsize + dsize - 1; /* dividend size to give desired qsize */ 406 407 /* dividend n, extended or truncated */ 408 tp = refmpn_malloc_limbs (tsize); 409 refmpn_copy_extend (tp, tsize, np, nsize); 410 411 qp = refmpn_malloc_limbs (qsize); 412 rp = refmpn_malloc_limbs (dsize); /* remainder, unused */ 413 414 ASSERT_ALWAYS (qsize == tsize - dsize + 1); 415 refmpn_tdiv_qr (qp, rp, (mp_size_t) 0, tp, tsize, dp, dsize); 416 417 PTR (want) = qp; 418 SIZ (want) = (sign >= 0 ? qsize : -qsize); 419 refmpf_normalize (want); 420 421 ret = refmpf_validate (name, got, want); 422 423 free (tp); 424 free (qp); 425 free (rp); 426 427 return ret; 428 }