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