/* dynamic-range.c
 *
 * Last changed in libpng 1.7.0
 *
 * COPYRIGHT: Written by John Cunningham Bowler, 2015
 * To the extent possible under law, the author has waived all copyright and
 * related or neighboring rights to this work.  This work is published from:
 * United States.
 *
 * Find the dynamic range of a given gamma encoding given a (linear) precision
 * and a maximum number of encoded values.
 */
#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include <math.h>
#include <assert.h>

double range(unsigned int steps, double factor, double gamma)
{
   return pow((steps * (pow(factor, 1/gamma) - 1)), gamma);
}

double max_range_gamma(unsigned int steps, double factor, double *max_range,
      double glo, double rlo, double gmid, double rmid, double ghi, double rhi)
{
   /* Given three values which contain a peak value (so rmid > rlo and rmid >
    * rhi) find the peak by repeated division of the range.  The algorithm is to
    * find the range for two gamma values mid-way between the two pairs
    * (glo,gmid), (ghi,gmid) then find the max; this gives us a new glo/ghi
    * which must be half the distance apart of the previous pair.
    */
   double gammas[5];
   double ranges[5];

   gammas[0] = glo;  ranges[0] = rlo;
   gammas[2] = gmid; ranges[2] = rmid;
   gammas[4] = ghi;  ranges[4] = rhi;

   for (;;)
   {
      int i, m;

      ranges[1] = range(steps, factor, gammas[1] = (gammas[0]+gammas[2])/2);
      ranges[3] = range(steps, factor, gammas[3] = (gammas[2]+gammas[4])/2);

      for (m=1, i=2; i<4; ++i)
         if (ranges[i] >= ranges[m])
            m = i;

      assert(gammas[0] < gammas[m] && gammas[m] < gammas[4]);
      assert(ranges[0] < ranges[m] && ranges[m] > ranges[4]);

      gammas[0] = gammas[m-1]; ranges[0] = ranges[m-1];
      gammas[4] = gammas[m+1]; ranges[4] = ranges[m+1];
      gammas[2] = gammas[m];   ranges[2] = ranges[m];

      if (((gammas[4] - gammas[0])/gammas[2]-1) < 3*DBL_EPSILON ||
          ((ranges[2] - ranges[0])/ranges[2]-1) < 6*DBL_EPSILON)
      {
         *max_range = ranges[2];
         return gammas[2];
      }
   }
}

double best_gamma(unsigned int values, double precision, double *best_range)
{
   /* The 'guess' gamma value is determined by the following formula, which is
    * itself derived from linear regression using values returned by this
    * program:
    */
   double gtry = values * precision / 2.736;
   double rtry;

   /* 'values' needs to be the number of steps after the first, we have to
    * reserve the first value, 0, for 0, so subtract 2 from values.  precision
    * must be adjusted to the step factor.
    */
   values -= 2U;
   precision += 1;
   rtry = range(values, precision, gtry);

   /* Now find two values either side of gtry with a lower range. */
   {
      double glo, ghi, rlo, rhi, gbest, rbest;

      glo = gtry;
      do
      {
         glo *= 0.9;
         rlo = range(values, precision, glo);
      }
      while (rlo >= rtry);

      ghi = gtry;
      do
      {
         ghi *= 1.1;
         rhi = range(values, precision, ghi);
      }
      while (rhi >= rtry);

      gbest = max_range_gamma(values, precision, &rbest,
            glo, rlo, gtry, rtry, ghi, rhi);

      *best_range = rbest / precision;
      return gbest;
   }
}

double linear_regression(double precision, double *bp)
{
   unsigned int values, count = 0;
   double g_sum = 0, g2_sum = 0, v_sum = 0, v2_sum = 0, gv_sum = 0;

   /* Perform simple linear regression to get:
    *
    *    gamma = a + b.values
    */
   for (values = 128; values < 65536; ++values, ++count)
   {
      double range;
      double gamma = best_gamma(values, precision, &range);

      g_sum += gamma;
      g2_sum += gamma * gamma;
      v_sum += values;
      v2_sum += values * (double)values;
      gv_sum += gamma * values;
      /* printf("%u %g %g\n", values, gamma, range); */
   }

   g_sum /= count;
   g2_sum /= count;
   v_sum /= count;
   v2_sum /= count;
   gv_sum /= count;

   {
      double b = (gv_sum - g_sum * v_sum) / (v2_sum - v_sum * v_sum);
      *bp = b;
      return g_sum - b * v_sum;
   }
}

int
main(int argc, const char **argv)
{
   double precision = argc == 2 ? atof(argv[1]) : 0;

   /* Perform a second linear regression here on b:
    *
    *    b = bA + bB * precision
    */
   if (precision == 0)
   {
      double b_sum = 0, b2_sum = 0, p_sum = 0, p2_sum = 0, bp_sum = 0,
             a_sum = 0, count = 0;

      for (precision = .001; precision <= 0.01; precision += .001, count += 1)
      {
         double b;
         double a = linear_regression(precision, &b);

         b_sum += b;
         b2_sum += b * b;
         p_sum += precision;
         p2_sum += precision * precision;
         bp_sum += b * precision;
         a_sum += a;
      }

      b_sum /= count;
      b2_sum /= count;
      p_sum /= count;
      p2_sum /= count;
      bp_sum /= count;
      a_sum /= count;

      {
         double bB = (bp_sum - b_sum * p_sum) / (p2_sum - p_sum * p_sum);
         double bA = b_sum - bB * p_sum;

         printf("a = %g, b = %g + precision/%g\n", a_sum, bA, 1/bB);
      }
   }

   else
   {
      unsigned int bits;
      double b;
      double a = linear_regression(precision, &b);
      printf("precision %g: gamma = %g + values*%g\n", precision, a, b);

      /* For information, given a precision: */
      for (bits=7U; bits <= 16U; ++bits)
      {
         unsigned int values = 1U<<bits;
         double gamma = values*precision/2.736;
         double r = range(values-2U, 1+precision, gamma);

         printf("bits: %u, gamma: %g, range: 1:%g\n", bits, gamma, r);
      }
   }

   return 0;
}