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