contrib/libtests/tarith.c

   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       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 positive 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    size_t index = 0;
 333    const char test = (char)ch;
 334    int number_is_valid = png_check_fp_number(&test, 1, &c.state, &index);
 335    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 }