1 2 /* tarith.c 3 * 4 * Copyright (c) 2011-2013 John Cunningham Bowler 5 * 6 * Last changed in libpng 1.6.0 [February 14, 2013] 7 * 8 * This code is released under the libpng license. 9 * For conditions of distribution and use, see the disclaimer 10 * and license in png.h 11 * 12 * Test internal arithmetic functions of libpng. 13 * 14 * This code must be linked against a math library (-lm), but does not require 15 * libpng or zlib to work. Because it includes the complete source of 'png.c' 16 * it tests the code with whatever compiler options are used to build it. 17 * Changing these options can substantially change the errors in the 18 * calculations that the compiler chooses! 19 */ 20 #define _POSIX_SOURCE 1 21 #define _ISOC99_SOURCE 1 22 23 /* Obtain a copy of the code to be tested (plus other things), disabling 24 * stuff that is not required. 25 */ 26 #include <math.h> 27 #include <stdlib.h> 28 #include <ctype.h> 29 #include <string.h> 30 #include <assert.h> 31 32 #include "../../pngpriv.h" 33 34 #define png_error png_warning 35 36 void png_warning(png_const_structrp png_ptr, png_const_charp msg) 37 { 38 fprintf(stderr, "validation: %s\n", msg); 39 } 40 41 #define png_fixed_error png_fixed_warning 42 43 void png_fixed_warning(png_const_structrp png_ptr, png_const_charp msg) 44 { 45 fprintf(stderr, "overflow in: %s\n", msg); 46 } 47 48 #define png_set_error_fn(pp, ep, efp, wfp) ((void)0) 49 #define png_malloc(pp, s) malloc(s) 50 #define png_malloc_warn(pp, s) malloc(s) 51 #define png_malloc_base(pp, s) malloc(s) 52 #define png_calloc(pp, s) calloc(1, (s)) 53 #define png_free(pp, s) free(s) 54 55 #define png_safecat(b, sb, pos, str) (pos) 56 #define png_format_number(start, end, format, number) (start) 57 58 #define crc32(crc, pp, s) (crc) 59 #define inflateReset(zs) Z_OK 60 61 #define png_create_struct(type) (0) 62 #define png_destroy_struct(pp) ((void)0) 63 #define png_create_struct_2(type, m, mm) (0) 64 #define png_destroy_struct_2(pp, f, mm) ((void)0) 65 66 #undef PNG_SIMPLIFIED_READ_SUPPORTED 67 #undef PNG_SIMPLIFIED_WRITE_SUPPORTED 68 #undef PNG_USER_MEM_SUPPORTED 69 70 #include "../../png.c" 71 72 /* Validate ASCII to fp routines. */ 73 static int verbose = 0; 74 75 int validation_ascii_to_fp(int count, int argc, char **argv) 76 { 77 int showall = 0; 78 double max_error=2; /* As a percentage error-in-last-digit/.5 */ 79 double max_error_abs=17; /* Used when precision is DBL_DIG */ 80 double max = 0; 81 double max_abs = 0; 82 double test = 0; /* Important to test this. */ 83 int precision = 5; 84 int nonfinite = 0; 85 int finite = 0; 86 int ok = 0; 87 int failcount = 0; 88 int minorarith = 0; 89 90 while (--argc > 0) 91 if (strcmp(*++argv, "-a") == 0) 92 showall = 1; 93 else if (strcmp(*argv, "-e") == 0 && argc > 0) 94 { 95 --argc; 96 max_error = atof(*++argv); 97 } 98 else if (strcmp(*argv, "-E") == 0 && argc > 0) 99 { 100 --argc; 101 max_error_abs = atof(*++argv); 102 } 103 else 104 { 105 fprintf(stderr, "unknown argument %s\n", *argv); 106 return 1; 107 } 108 109 do 110 { 111 png_size_t index; 112 int state, failed = 0; 113 char buffer[64]; 114 115 if (isfinite(test)) 116 ++finite; 117 else 118 ++nonfinite; 119 120 if (verbose) 121 fprintf(stderr, "%.*g %d\n", DBL_DIG, test, precision); 122 123 /* Check for overflow in the buffer by setting a marker. */ 124 memset(buffer, 71, sizeof buffer); 125 126 png_ascii_from_fp(0, buffer, precision+10, test, precision); 127 128 /* Allow for a three digit exponent, this stuff will fail if 129 * the exponent is bigger than this! 130 */ 131 if (buffer[precision+7] != 71) 132 { 133 fprintf(stderr, "%g[%d] -> '%s'[%lu] buffer overflow\n", test, 134 precision, buffer, (unsigned long)strlen(buffer)); 135 failed = 1; 136 } 137 138 /* Following are used for the number parser below and must be 139 * initialized to zero. 140 */ 141 state = 0; 142 index = 0; 143 if (!isfinite(test)) 144 { 145 /* Expect 'inf' */ 146 if (test >= 0 && strcmp(buffer, "inf") || 147 test < 0 && strcmp(buffer, "-inf")) 148 { 149 fprintf(stderr, "%g[%d] -> '%s' but expected 'inf'\n", test, 150 precision, buffer); 151 failed = 1; 152 } 153 } 154 else if (!png_check_fp_number(buffer, precision+10, &state, &index) || 155 buffer[index] != 0) 156 { 157 fprintf(stderr, "%g[%d] -> '%s' but has bad format ('%c')\n", test, 158 precision, buffer, buffer[index]); 159 failed = 1; 160 } 161 else if (PNG_FP_IS_NEGATIVE(state) && !(test < 0)) 162 { 163 fprintf(stderr, "%g[%d] -> '%s' but negative value not so reported\n", 164 test, precision, buffer); 165 failed = 1; 166 assert(!PNG_FP_IS_ZERO(state)); 167 assert(!PNG_FP_IS_POSITIVE(state)); 168 } 169 else if (PNG_FP_IS_ZERO(state) && !(test == 0)) 170 { 171 fprintf(stderr, "%g[%d] -> '%s' but zero value not so reported\n", 172 test, precision, buffer); 173 failed = 1; 174 assert(!PNG_FP_IS_NEGATIVE(state)); 175 assert(!PNG_FP_IS_POSITIVE(state)); 176 } 177 else if (PNG_FP_IS_POSITIVE(state) && !(test > 0)) 178 { 179 fprintf(stderr, "%g[%d] -> '%s' but postive value not so reported\n", 180 test, precision, buffer); 181 failed = 1; 182 assert(!PNG_FP_IS_NEGATIVE(state)); 183 assert(!PNG_FP_IS_ZERO(state)); 184 } 185 else 186 { 187 /* Check the result against the original. */ 188 double out = atof(buffer); 189 double change = fabs((out - test)/test); 190 double allow = .5/pow(10, 191 (precision >= DBL_DIG) ? DBL_DIG-1 : precision-1); 192 193 /* NOTE: if you hit this error case are you compiling with gcc 194 * and -O0? Try -O2 - the errors can accumulate if the FP 195 * code above is not optimized and may drift outside the .5 in 196 * DBL_DIG allowed. In any case a small number of errors may 197 * occur (very small ones - 1 or 2%) because of rounding in the 198 * calculations, either in the conversion API or in atof. 199 */ 200 if (change >= allow && (isfinite(out) || 201 fabs(test/DBL_MAX) <= 1-allow)) 202 { 203 double percent = (precision >= DBL_DIG) ? max_error_abs : max_error; 204 double allowp = (change-allow)*100/allow; 205 206 if (precision >= DBL_DIG) 207 { 208 if (max_abs < allowp) max_abs = allowp; 209 } 210 211 else 212 { 213 if (max < allowp) max = allowp; 214 } 215 216 if (showall || allowp >= percent) 217 { 218 fprintf(stderr, 219 "%.*g[%d] -> '%s' -> %.*g number changed (%g > %g (%d%%))\n", 220 DBL_DIG, test, precision, buffer, DBL_DIG, out, change, allow, 221 (int)round(allowp)); 222 failed = 1; 223 } 224 else 225 ++minorarith; 226 } 227 } 228 229 if (failed) 230 ++failcount; 231 else 232 ++ok; 233 234 skip: 235 /* Generate a new number and precision. */ 236 precision = rand(); 237 if (precision & 1) test = -test; 238 precision >>= 1; 239 240 /* Generate random numbers. */ 241 if (test == 0 || !isfinite(test)) 242 test = precision+1; 243 else 244 { 245 /* Derive the exponent from the previous rand() value. */ 246 int exponent = precision % (DBL_MAX_EXP - DBL_MIN_EXP) + DBL_MIN_EXP; 247 int tmp; 248 test = frexp(test * rand(), &tmp); 249 test = ldexp(test, exponent); 250 precision >>= 8; /* arbitrary */ 251 } 252 253 /* This limits the precision to 32 digits, enough for standard 254 * IEEE implementations which have at most 15 digits. 255 */ 256 precision = (precision & 0x1f) + 1; 257 } 258 while (--count); 259 260 printf("Tested %d finite values, %d non-finite, %d OK (%d failed) %d minor " 261 "arithmetic errors\n", finite, nonfinite, ok, failcount, minorarith); 262 printf(" Error with >=%d digit precision %.2f%%\n", DBL_DIG, max_abs); 263 printf(" Error with < %d digit precision %.2f%%\n", DBL_DIG, max); 264 265 return 0; 266 } 267 268 /* Observe that valid FP numbers have the forms listed in the PNG extensions 269 * specification: 270 * 271 * [+,-]{integer,integer.fraction,.fraction}[{e,E}[+,-]integer] 272 * 273 * Test each of these in turn, including invalid cases. 274 */ 275 typedef enum checkfp_state 276 { 277 start, fraction, exponent, states 278 } checkfp_state; 279 280 /* The characters (other than digits) that characterize the states: */ 281 static const char none[] = ""; 282 static const char hexdigits[16] = "0123456789ABCDEF"; 283 284 static const struct 285 { 286 const char *start; /* Characters valid at the start */ 287 const char *end; /* Valid characters that end the state */ 288 const char *tests; /* Characters to test after 2 digits seen */ 289 } 290 state_characters[states] = 291 { 292 /* start: */ { "+-.", ".eE", "+-.e*0369" }, 293 /* fraction: */ { none, "eE", "+-.E#0147" }, 294 /* exponent: */ { "+-", none, "+-.eE^0258" } 295 }; 296 297 typedef struct 298 { 299 char number[1024]; /* Buffer for number being tested */ 300 int limit; /* Command line limit */ 301 int verbose; /* Shadows global variable */ 302 int ctimes; /* Number of numbers tested */ 303 int cmillions; /* Count of millions of numbers */ 304 int cinvalid; /* Invalid strings checked */ 305 int cnoaccept; /* Characters not accepted */ 306 } 307 checkfp_command; 308 309 typedef struct 310 { 311 int cnumber; /* Index into number string */ 312 checkfp_state check_state; /* Current number state */ 313 int at_start; /* At start (first character) of state */ 314 int cdigits_in_state; /* Digits seen in that state */ 315 int limit; /* Limit on same for checking all chars */ 316 int state; /* Current parser state */ 317 int is_negative; /* Number is negative */ 318 int is_zero; /* Number is (still) zero */ 319 int number_was_valid; /* Previous character validity */ 320 } 321 checkfp_control; 322 323 static int check_all_characters(checkfp_command *co, checkfp_control c); 324 325 static int check_some_characters(checkfp_command *co, checkfp_control c, 326 const char *tests); 327 328 static int check_one_character(checkfp_command *co, checkfp_control c, int ch) 329 { 330 /* Test this character (ch) to ensure the parser does the correct thing. 331 */ 332 png_size_t index = 0; 333 const char test = (char)ch; 334 const int number_is_valid = png_check_fp_number(&test, 1, &c.state, &index); 335 const int character_accepted = (index == 1); 336 337 if (c.check_state != exponent && isdigit(ch) && ch != '0') 338 c.is_zero = 0; 339 340 if (c.check_state == start && c.at_start && ch == '-') 341 c.is_negative = 1; 342 343 if (isprint(ch)) 344 co->number[c.cnumber++] = (char)ch; 345 else 346 { 347 co->number[c.cnumber++] = '<'; 348 co->number[c.cnumber++] = hexdigits[(ch >> 4) & 0xf]; 349 co->number[c.cnumber++] = hexdigits[ch & 0xf]; 350 co->number[c.cnumber++] = '>'; 351 } 352 co->number[c.cnumber] = 0; 353 354 if (co->verbose > 1) 355 fprintf(stderr, "%s\n", co->number); 356 357 if (++(co->ctimes) == 1000000) 358 { 359 if (co->verbose == 1) 360 fputc('.', stderr); 361 co->ctimes = 0; 362 ++(co->cmillions); 363 } 364 365 if (!number_is_valid) 366 ++(co->cinvalid); 367 368 if (!character_accepted) 369 ++(co->cnoaccept); 370 371 /* This should never fail (it's a serious bug if it does): */ 372 if (index != 0 && index != 1) 373 { 374 fprintf(stderr, "%s: read beyond end of string (%lu)\n", co->number, 375 (unsigned long)index); 376 return 0; 377 } 378 379 /* Validate the new state, note that the PNG_FP_IS_ macros all return 380 * false unless the number is valid. 381 */ 382 if (PNG_FP_IS_NEGATIVE(c.state) != 383 (number_is_valid && !c.is_zero && c.is_negative)) 384 { 385 fprintf(stderr, "%s: negative when it is not\n", co->number); 386 return 0; 387 } 388 389 if (PNG_FP_IS_ZERO(c.state) != (number_is_valid && c.is_zero)) 390 { 391 fprintf(stderr, "%s: zero when it is not\n", co->number); 392 return 0; 393 } 394 395 if (PNG_FP_IS_POSITIVE(c.state) != 396 (number_is_valid && !c.is_zero && !c.is_negative)) 397 { 398 fprintf(stderr, "%s: positive when it is not\n", co->number); 399 return 0; 400 } 401 402 /* Testing a digit */ 403 if (isdigit(ch)) 404 { 405 if (!character_accepted) 406 { 407 fprintf(stderr, "%s: digit '%c' not accepted\n", co->number, ch); 408 return 0; 409 } 410 411 if (!number_is_valid) 412 { 413 fprintf(stderr, "%s: saw a digit (%c) but number not valid\n", 414 co->number, ch); 415 return 0; 416 } 417 418 ++c.cdigits_in_state; 419 c.at_start = 0; 420 c.number_was_valid = 1; 421 422 /* Continue testing characters in this state. Either test all of 423 * them or, if we have already seen one digit in this state, just test a 424 * limited set. 425 */ 426 if (c.cdigits_in_state < 1) 427 return check_all_characters(co, c); 428 429 else 430 return check_some_characters(co, c, 431 state_characters[c.check_state].tests); 432 } 433 434 /* A non-digit; is it allowed here? */ 435 else if (((ch == '+' || ch == '-') && c.check_state != fraction && 436 c.at_start) || 437 (ch == '.' && c.check_state == start) || 438 ((ch == 'e' || ch == 'E') && c.number_was_valid && 439 c.check_state != exponent)) 440 { 441 if (!character_accepted) 442 { 443 fprintf(stderr, "%s: character '%c' not accepted\n", co->number, ch); 444 return 0; 445 } 446 447 /* The number remains valid after start of fraction but nowhere else. */ 448 if (number_is_valid && (c.check_state != start || ch != '.')) 449 { 450 fprintf(stderr, "%s: saw a non-digit (%c) but number valid\n", 451 co->number, ch); 452 return 0; 453 } 454 455 c.number_was_valid = number_is_valid; 456 457 /* Check for a state change. When changing to 'fraction' if the number 458 * is valid at this point set the at_start to false to allow an exponent 459 * 'e' to come next. 460 */ 461 if (c.check_state == start && ch == '.') 462 { 463 c.check_state = fraction; 464 c.at_start = !number_is_valid; 465 c.cdigits_in_state = 0; 466 c.limit = co->limit; 467 return check_all_characters(co, c); 468 } 469 470 else if (c.check_state < exponent && (ch == 'e' || ch == 'E')) 471 { 472 c.check_state = exponent; 473 c.at_start = 1; 474 c.cdigits_in_state = 0; 475 c.limit = co->limit; 476 return check_all_characters(co, c); 477 } 478 479 /* Else it was a sign, and the state doesn't change. */ 480 else 481 { 482 if (ch != '-' && ch != '+') 483 { 484 fprintf(stderr, "checkfp: internal error (1)\n"); 485 return 0; 486 } 487 488 c.at_start = 0; 489 return check_all_characters(co, c); 490 } 491 } 492 493 /* Testing an invalid character */ 494 else 495 { 496 if (character_accepted) 497 { 498 fprintf(stderr, "%s: character '%c' [0x%.2x] accepted\n", co->number, 499 ch, ch); 500 return 0; 501 } 502 503 if (number_is_valid != c.number_was_valid) 504 { 505 fprintf(stderr, 506 "%s: character '%c' [0x%.2x] changed number validity\n", co->number, 507 ch, ch); 508 return 0; 509 } 510 511 /* Do nothing - the parser has stuck; return success and keep going with 512 * the next character. 513 */ 514 } 515 516 /* Successful return (the caller will try the next character.) */ 517 return 1; 518 } 519 520 static int check_all_characters(checkfp_command *co, checkfp_control c) 521 { 522 int ch; 523 524 if (c.cnumber+4 < sizeof co->number) for (ch=0; ch<256; ++ch) 525 { 526 if (!check_one_character(co, c, ch)) 527 return 0; 528 } 529 530 return 1; 531 } 532 533 static int check_some_characters(checkfp_command *co, checkfp_control c, 534 const char *tests) 535 { 536 int i; 537 538 --(c.limit); 539 540 if (c.cnumber+4 < sizeof co->number && c.limit >= 0) 541 { 542 if (c.limit > 0) for (i=0; tests[i]; ++i) 543 { 544 if (!check_one_character(co, c, tests[i])) 545 return 0; 546 } 547 548 /* At the end check all the characters. */ 549 else 550 return check_all_characters(co, c); 551 } 552 553 return 1; 554 } 555 556 int validation_checkfp(int count, int argc, char **argv) 557 { 558 int result; 559 checkfp_command command; 560 checkfp_control control; 561 562 command.number[0] = 0; 563 command.limit = 3; 564 command.verbose = verbose; 565 command.ctimes = 0; 566 command.cmillions = 0; 567 command.cinvalid = 0; 568 command.cnoaccept = 0; 569 570 while (--argc > 0) 571 { 572 ++argv; 573 if (argc > 1 && strcmp(*argv, "-l") == 0) 574 { 575 --argc; 576 command.limit = atoi(*++argv); 577 } 578 579 else 580 { 581 fprintf(stderr, "unknown argument %s\n", *argv); 582 return 1; 583 } 584 } 585 586 control.cnumber = 0; 587 control.check_state = start; 588 control.at_start = 1; 589 control.cdigits_in_state = 0; 590 control.limit = command.limit; 591 control.state = 0; 592 control.is_negative = 0; 593 control.is_zero = 1; 594 control.number_was_valid = 0; 595 596 result = check_all_characters(&command, control); 597 598 printf("checkfp: %s: checked %d,%.3d,%.3d,%.3d strings (%d invalid)\n", 599 result ? "pass" : "FAIL", command.cmillions / 1000, 600 command.cmillions % 1000, command.ctimes / 1000, command.ctimes % 1000, 601 command.cinvalid); 602 603 return result; 604 } 605 606 int validation_muldiv(int count, int argc, char **argv) 607 { 608 int tested = 0; 609 int overflow = 0; 610 int error = 0; 611 int error64 = 0; 612 int passed = 0; 613 int randbits = 0; 614 png_uint_32 randbuffer; 615 png_fixed_point a; 616 png_int_32 times, div; 617 618 while (--argc > 0) 619 { 620 fprintf(stderr, "unknown argument %s\n", *++argv); 621 return 1; 622 } 623 624 /* Find out about the random number generator. */ 625 randbuffer = RAND_MAX; 626 while (randbuffer != 0) ++randbits, randbuffer >>= 1; 627 printf("Using random number generator that makes %d bits\n", randbits); 628 for (div=0; div<32; div += randbits) 629 randbuffer = (randbuffer << randbits) ^ rand(); 630 631 a = 0; 632 times = div = 0; 633 do 634 { 635 png_fixed_point result; 636 /* NOTE: your mileage may vary, a type is required below that can 637 * hold 64 bits or more, if floating point is used a 64-bit or 638 * better mantissa is required. 639 */ 640 long long int fp, fpround; 641 unsigned long hi, lo; 642 int ok; 643 644 /* Check the values, png_64bit_product can only handle positive 645 * numbers, so correct for that here. 646 */ 647 { 648 long u1, u2; 649 int n = 0; 650 if (a < 0) u1 = -a, n = 1; else u1 = a; 651 if (times < 0) u2 = -times, n = !n; else u2 = times; 652 png_64bit_product(u1, u2, &hi, &lo); 653 if (n) 654 { 655 /* -x = ~x+1 */ 656 lo = ((~lo) + 1) & 0xffffffff; 657 hi = ~hi; 658 if (lo == 0) ++hi; 659 } 660 } 661 662 fp = a; 663 fp *= times; 664 if ((fp & 0xffffffff) != lo || ((fp >> 32) & 0xffffffff) != hi) 665 { 666 fprintf(stderr, "png_64bit_product %d * %d -> %lx|%.8lx not %llx\n", 667 a, times, hi, lo, fp); 668 ++error64; 669 } 670 671 if (div != 0) 672 { 673 /* Round - this is C round to zero. */ 674 if ((fp < 0) != (div < 0)) 675 fp -= div/2; 676 else 677 fp += div/2; 678 679 fp /= div; 680 fpround = fp; 681 /* Assume 2's complement here: */ 682 ok = fpround <= PNG_UINT_31_MAX && 683 fpround >= -1-(long long int)PNG_UINT_31_MAX; 684 if (!ok) ++overflow; 685 } 686 else 687 ok = 0, ++overflow, fpround = fp/*misleading*/; 688 689 if (verbose) 690 fprintf(stderr, "TEST %d * %d / %d -> %lld (%s)\n", a, times, div, 691 fp, ok ? "ok" : "overflow"); 692 693 ++tested; 694 if (png_muldiv(&result, a, times, div) != ok) 695 { 696 ++error; 697 if (ok) 698 fprintf(stderr, "%d * %d / %d -> overflow (expected %lld)\n", a, 699 times, div, fp); 700 else 701 fprintf(stderr, "%d * %d / %d -> %d (expected overflow %lld)\n", a, 702 times, div, result, fp); 703 } 704 else if (ok && result != fpround) 705 { 706 ++error; 707 fprintf(stderr, "%d * %d / %d -> %d not %lld\n", a, times, div, result, 708 fp); 709 } 710 else 711 ++passed; 712 713 /* Generate three new values, this uses rand() and rand() only returns 714 * up to RAND_MAX. 715 */ 716 /* CRUDE */ 717 a += times; 718 times += div; 719 div = randbuffer; 720 randbuffer = (randbuffer << randbits) ^ rand(); 721 } 722 while (--count > 0); 723 724 printf("%d tests including %d overflows, %d passed, %d failed (%d 64-bit " 725 "errors)\n", tested, overflow, passed, error, error64); 726 return 0; 727 } 728 729 /* When FP is on this just becomes a speed test - compile without FP to get real 730 * validation. 731 */ 732 #ifdef PNG_FLOATING_ARITHMETIC_SUPPORTED 733 #define LN2 .000010576586617430806112933839 /* log(2)/65536 */ 734 #define L2INV 94548.46219969910586572651 /* 65536/log(2) */ 735 736 /* For speed testing, need the internal functions too: */ 737 static png_uint_32 png_log8bit(unsigned x) 738 { 739 if (x > 0) 740 return (png_uint_32)floor(.5-log(x/255.)*L2INV); 741 742 return 0xffffffff; 743 } 744 745 static png_uint_32 png_log16bit(png_uint_32 x) 746 { 747 if (x > 0) 748 return (png_uint_32)floor(.5-log(x/65535.)*L2INV); 749 750 return 0xffffffff; 751 } 752 753 static png_uint_32 png_exp(png_uint_32 x) 754 { 755 return (png_uint_32)floor(.5 + exp(x * -LN2) * 0xffffffffU); 756 } 757 758 static png_byte png_exp8bit(png_uint_32 log) 759 { 760 return (png_byte)floor(.5 + exp(log * -LN2) * 255); 761 } 762 763 static png_uint_16 png_exp16bit(png_uint_32 log) 764 { 765 return (png_uint_16)floor(.5 + exp(log * -LN2) * 65535); 766 } 767 #endif /* FLOATING_ARITHMETIC */ 768 769 int validation_gamma(int argc, char **argv) 770 { 771 double gamma[9] = { 2.2, 1.8, 1.52, 1.45, 1., 1/1.45, 1/1.52, 1/1.8, 1/2.2 }; 772 double maxerr; 773 int i, silent=0, onlygamma=0; 774 775 /* Silence the output with -s, just test the gamma functions with -g: */ 776 while (--argc > 0) 777 if (strcmp(*++argv, "-s") == 0) 778 silent = 1; 779 else if (strcmp(*argv, "-g") == 0) 780 onlygamma = 1; 781 else 782 { 783 fprintf(stderr, "unknown argument %s\n", *argv); 784 return 1; 785 } 786 787 if (!onlygamma) 788 { 789 /* First validate the log functions: */ 790 maxerr = 0; 791 for (i=0; i<256; ++i) 792 { 793 double correct = -log(i/255.)/log(2.)*65536; 794 double error = png_log8bit(i) - correct; 795 796 if (i != 0 && fabs(error) > maxerr) 797 maxerr = fabs(error); 798 799 if (i == 0 && png_log8bit(i) != 0xffffffff || 800 i != 0 && png_log8bit(i) != floor(correct+.5)) 801 { 802 fprintf(stderr, "8-bit log error: %d: got %u, expected %f\n", 803 i, png_log8bit(i), correct); 804 } 805 } 806 807 if (!silent) 808 printf("maximum 8-bit log error = %f\n", maxerr); 809 810 maxerr = 0; 811 for (i=0; i<65536; ++i) 812 { 813 double correct = -log(i/65535.)/log(2.)*65536; 814 double error = png_log16bit(i) - correct; 815 816 if (i != 0 && fabs(error) > maxerr) 817 maxerr = fabs(error); 818 819 if (i == 0 && png_log16bit(i) != 0xffffffff || 820 i != 0 && png_log16bit(i) != floor(correct+.5)) 821 { 822 if (error > .68) /* By experiment error is less than .68 */ 823 { 824 fprintf(stderr, "16-bit log error: %d: got %u, expected %f" 825 " error: %f\n", i, png_log16bit(i), correct, error); 826 } 827 } 828 } 829 830 if (!silent) 831 printf("maximum 16-bit log error = %f\n", maxerr); 832 833 /* Now exponentiations. */ 834 maxerr = 0; 835 for (i=0; i<=0xfffff; ++i) 836 { 837 double correct = exp(-i/65536. * log(2.)) * (65536. * 65536); 838 double error = png_exp(i) - correct; 839 840 if (fabs(error) > maxerr) 841 maxerr = fabs(error); 842 if (fabs(error) > 1883) /* By experiment. */ 843 { 844 fprintf(stderr, "32-bit exp error: %d: got %u, expected %f" 845 " error: %f\n", i, png_exp(i), correct, error); 846 } 847 } 848 849 if (!silent) 850 printf("maximum 32-bit exp error = %f\n", maxerr); 851 852 maxerr = 0; 853 for (i=0; i<=0xfffff; ++i) 854 { 855 double correct = exp(-i/65536. * log(2.)) * 255; 856 double error = png_exp8bit(i) - correct; 857 858 if (fabs(error) > maxerr) 859 maxerr = fabs(error); 860 if (fabs(error) > .50002) /* By experiment */ 861 { 862 fprintf(stderr, "8-bit exp error: %d: got %u, expected %f" 863 " error: %f\n", i, png_exp8bit(i), correct, error); 864 } 865 } 866 867 if (!silent) 868 printf("maximum 8-bit exp error = %f\n", maxerr); 869 870 maxerr = 0; 871 for (i=0; i<=0xfffff; ++i) 872 { 873 double correct = exp(-i/65536. * log(2.)) * 65535; 874 double error = png_exp16bit(i) - correct; 875 876 if (fabs(error) > maxerr) 877 maxerr = fabs(error); 878 if (fabs(error) > .524) /* By experiment */ 879 { 880 fprintf(stderr, "16-bit exp error: %d: got %u, expected %f" 881 " error: %f\n", i, png_exp16bit(i), correct, error); 882 } 883 } 884 885 if (!silent) 886 printf("maximum 16-bit exp error = %f\n", maxerr); 887 } /* !onlygamma */ 888 889 /* Test the overall gamma correction. */ 890 for (i=0; i<9; ++i) 891 { 892 unsigned j; 893 double g = gamma[i]; 894 png_fixed_point gfp = floor(g * PNG_FP_1 + .5); 895 896 if (!silent) 897 printf("Test gamma %f\n", g); 898 899 maxerr = 0; 900 for (j=0; j<256; ++j) 901 { 902 double correct = pow(j/255., g) * 255; 903 png_byte out = png_gamma_8bit_correct(j, gfp); 904 double error = out - correct; 905 906 if (fabs(error) > maxerr) 907 maxerr = fabs(error); 908 if (out != floor(correct+.5)) 909 { 910 fprintf(stderr, "8bit %d ^ %f: got %d expected %f error %f\n", 911 j, g, out, correct, error); 912 } 913 } 914 915 if (!silent) 916 printf("gamma %f: maximum 8-bit error %f\n", g, maxerr); 917 918 maxerr = 0; 919 for (j=0; j<65536; ++j) 920 { 921 double correct = pow(j/65535., g) * 65535; 922 png_uint_16 out = png_gamma_16bit_correct(j, gfp); 923 double error = out - correct; 924 925 if (fabs(error) > maxerr) 926 maxerr = fabs(error); 927 if (fabs(error) > 1.62) 928 { 929 fprintf(stderr, "16bit %d ^ %f: got %d expected %f error %f\n", 930 j, g, out, correct, error); 931 } 932 } 933 934 if (!silent) 935 printf("gamma %f: maximum 16-bit error %f\n", g, maxerr); 936 } 937 938 return 0; 939 } 940 941 /**************************** VALIDATION TESTS ********************************/ 942 /* Various validation routines are included herein, they require some 943 * definition for png_warning and png_error, seetings of VALIDATION: 944 * 945 * 1: validates the ASCII to floating point conversions 946 * 2: validates png_muldiv 947 * 3: accuracy test of fixed point gamma tables 948 */ 949 950 /* The following COUNT (10^8) takes about 1 hour on a 1GHz Pentium IV 951 * processor. 952 */ 953 #define COUNT 1000000000 954 955 int main(int argc, char **argv) 956 { 957 int count = COUNT; 958 959 while (argc > 1) 960 { 961 if (argc > 2 && strcmp(argv[1], "-c") == 0) 962 { 963 count = atoi(argv[2]); 964 argc -= 2; 965 argv += 2; 966 } 967 968 else if (strcmp(argv[1], "-v") == 0) 969 { 970 ++verbose; 971 --argc; 972 ++argv; 973 } 974 975 else 976 break; 977 } 978 979 if (count > 0 && argc > 1) 980 { 981 if (strcmp(argv[1], "ascii") == 0) 982 return validation_ascii_to_fp(count, argc-1, argv+1); 983 else if (strcmp(argv[1], "checkfp") == 0) 984 return validation_checkfp(count, argc-1, argv+1); 985 else if (strcmp(argv[1], "muldiv") == 0) 986 return validation_muldiv(count, argc-1, argv+1); 987 else if (strcmp(argv[1], "gamma") == 0) 988 return validation_gamma(argc-1, argv+1); 989 } 990 991 /* Bad argument: */ 992 fprintf(stderr, 993 "usage: tarith [-v] [-c count] {ascii,muldiv,gamma} [args]\n"); 994 fprintf(stderr, " arguments: ascii [-a (all results)] [-e error%%]\n"); 995 fprintf(stderr, " checkfp [-l max-number-chars]\n"); 996 fprintf(stderr, " muldiv\n"); 997 fprintf(stderr, " gamma -s (silent) -g (only gamma; no log)\n"); 998 return 1; 999 } 1000
