1#!/bin/sh
  2#
  3# intgamma.sh
  4#
  5# Last changed in libpng 1.6.0 [February 14, 2013]
  6#
  7# COPYRIGHT: Written by John Cunningham Bowler, 2013.
  8# To the extent possible under law, the author has waived all copyright and
  9# related or neighboring rights to this work.  This work is published from:
 10# United States.
 11#
 12# Shell script to generate png.c 8-bit and 16-bit log tables (see the code in
 13# png.c for details).
 14#
 15# This script uses the "bc" arbitrary precision calculator to calculate 32-bit
 16# fixed point values of logarithms appropriate to finding the log of an 8-bit
 17# (0..255) value and a similar table for the exponent calculation.
 18#
 19# "bc" must be on the path when the script is executed, and the math library
 20# (-lm) must be available
 21#
 22# function to print out a list of numbers as integers; the function truncates
 23# the integers which must be one-per-line
 24function print(){
 25   awk 'BEGIN{
 26      str = ""
 27   }
 28   {
 29      sub("\\.[0-9]*$", "")
 30      if ($0 == "")
 31         $0 = "0"
 32
 33      if (str == "")
 34         t = "   " $0 "U"
 35      else
 36         t = str ", " $0 "U"
 37
 38      if (length(t) >= 80) {
 39         print str ","
 40         str = "   " $0 "U"
 41      } else
 42         str = t
 43   }
 44   END{
 45      print str
 46   }'
 47}
 48#
 49# The logarithm table.
 50cat <<END
 51/* 8-bit log table: png_8bit_l2[128]
 52 * This is a table of -log(value/255)/log(2) for 'value' in the range 128 to
 53 * 255, so it's the base 2 logarithm of a normalized 8-bit floating point
 54 * mantissa.  The numbers are 32-bit fractions.
 55 */
 56static const png_uint_32
 57png_8bit_l2[128] =
 58{
 59END
 60#
 61bc -lqws <<END | print
 62f=65536*65536/l(2)
 63for (i=128;i<256;++i) { .5 - l(i/255)*f; }
 64END
 65echo '};'
 66echo
 67#
 68# The exponent table.
 69cat <<END
 70/* The 'exp()' case must invert the above, taking a 20-bit fixed point
 71 * logarithmic value and returning a 16 or 8-bit number as appropriate.  In
 72 * each case only the low 16 bits are relevant - the fraction - since the
 73 * integer bits (the top 4) simply determine a shift.
 74 *
 75 * The worst case is the 16-bit distinction between 65535 and 65534; this
 76 * requires perhaps spurious accuracy in the decoding of the logarithm to
 77 * distinguish log2(65535/65534.5) - 10^-5 or 17 bits.  There is little chance
 78 * of getting this accuracy in practice.
 79 *
 80 * To deal with this the following exp() function works out the exponent of the
 81 * frational part of the logarithm by using an accurate 32-bit value from the
 82 * top four fractional bits then multiplying in the remaining bits.
 83 */
 84static const png_uint_32
 85png_32bit_exp[16] =
 86{
 87END
 88#
 89bc -lqws <<END | print
 90f=l(2)/16
 91for (i=0;i<16;++i) {
 92   x = .5 + e(-i*f)*2^32;
 93   if (x >= 2^32) x = 2^32-1;
 94   x;
 95}
 96END
 97echo '};'
 98echo
 99#
100# And the table of adjustment values.
101cat <<END
102/* Adjustment table; provided to explain the numbers in the code below. */
103#if 0
104END
105bc -lqws <<END | awk '{ printf "%5d %s\n", 12-NR, $0 }'
106for (i=11;i>=0;--i){
107   (1 - e(-(2^i)/65536*l(2))) * 2^(32-i)
108}
109END
110echo '#endif'