#include "log.h"

#define crumb printf

#include "lntable.frag"

#define tablen (sizeof(lntable)/sizeof(*lntable))



double
antiln (double x)
{
  double index, fract;
  int i, lower, upper;
  int antiinfloop, halfway;
  double retval;

//crumb("Antilogarithm of %f\n", x);
  if (x < 0) return (antiln(x + 1.0) / ee);
  if (x > 1.0) return (antiln(x - 1.0) * ee);
  if (x == 0) return 1;
  if (x == 1.0) return ee;

/* Find the two numbers x fits between. */
/* Get the indices. */
/* Lerp between the indices for a more precise value. */
/* Scale index value between 1 and e: */
/*  antilog = (index * (ee - 1) / (lntable.len)) + 1 */

  retval = 0;
  lower = 0; upper = tablen;
  antiinfloop = 0;
  while ((upper - lower > 1) && (antiinfloop < tablen))
    {
      halfway = lower + ((upper - lower) >> 1);
//crumb("Comparing x(%.18f) to [%d](%.18f), (%d, %d)\n", x, halfway, lntable[halfway], lower, upper);
      if (x < lntable[halfway])
        {
          upper = halfway;
        }
      else
        {
          lower = halfway;
        }
      antiinfloop++;
    }
  /* Now it should be such that: lntable[lower] <= x < lntable[upper] */

  if (antiinfloop >= tablen)  /* ugh ugh. */
      upper = lower + 1;

//crumb("For x = %.18f:\n", x);
//crumb("  lower = %d [%.18f]\n", lower, lntable[lower]);
//crumb("  upper = %d [%.18f]\n", upper, lntable[upper]);

  /* Lerp. */
  fract = (x - lntable[lower]) / (lntable[upper] - lntable[lower]);
  index = (float)lower + fract;
//  index *= 1.000978;  /* Always a little too low. */
  index *= 1.0009768;
//crumb("Calculating on index %.18f\n", index);

  /* Scale to 1..e */
  retval = (index * (ee - 1) / tablen) + 1;
//crumb("Antilog %f = %f\n", x, retval);
  return retval;
}



/*
Table-lookup + linear interpolation techinque.
Natural logarithm.

Given a value x,
multiply or divide by e until within in range 1..e.
Keep track of number of multiplys/divides (maintaining sign).

equation then is   ln (x * e^n) = L
Where x is the given value, n is the number of multiply/divides, L is the desired result in 1..e inclusive.

d = 1000 * ((x * e^n) - 1) / (e - 1);

d give a value between 0..1000 inclusive.
This is used for offset into lookup table.

i = int(d);

i is index into lookup table.


a = lntable[i];
b = lntable[i+1];

The value we look for is lerped between this two values based on the fractional part of d.


c = d - i;

c is the fractional part of d.


Lerp between a and b is:
L = a + ((b - a) * c)

minor subst:
L = lntable[i] + ((lntable[i+1] - lntable[i]) * (d - i));


ln (x * e^n) = L
ln(x) + ln(e^n) = L
ln(x) + n = L
ln(x) = L - n

Our real result.

*/

double
ln (double x)
{
  double retval;
  double d, c;
  int n, i;

  if (x == 0)  return 0;  /* FIXME: error condition. */
  n = 0;
  if (x > ee)
    {
      while (x > ee)
        {
          x /= ee;
          n--;
        }
    }
  else if (x < 1.0)
    {
      while (x < 1.0)
        {
          x *= ee;
          n++;
        }
    }

  d = ((sizeof(lntable)-1) / sizeof(*lntable)) * (x - 1.0) / (ee - 1.0);
  i = (int)d;
  c = (d - (float)i);
  retval = lntable[i] + ((lntable[i+1] - lntable[i]) * c);
  retval -= (float)n;
  return retval;
}



double
exp (double x)
{
  double retval;
  int i, n;
//  double shifting;
  int shifting = 0;

/* This function aims to be based iteration, instead of using antiln's recursion. */
  n = 0;
  if (x > 1.0)
    {
      while (x > 1.0)
        {
          shifting--;
          x -= 1.0;
        }
    }
  else if (x < 0)
    {
      while (x < 0)
        {
          shifting++;
          x += 1.0;
        }
    }
//crumb("exp: of %f after shifting %d\n", x, shifting);
  retval = antiln(x);
  if (shifting > 0)
    {
      while (shifting > 0)
        {
          shifting--;
          retval /= ee;
        }
    }
  else if (shifting < 0)
    {
      while (shifting < 0)
        {
          shifting++;
          retval *= ee;
        }
    }
/* Add or subtract 1 to/from x until  0 <= x < 1 */
/* Number of add/subtracts stored in `n'. */
/* Take antilogarithm of remaining x, store in x. */
/* Raise e to the power of n, store in shifting. */
/* Multiply the current x with shifting, store in retval. */
  return retval;
}



double
lg (double x)
{
  double retval;
  retval = ln(x) / ln(10);
  return retval;
}



double
pow (double base, double exponent)
{
/* Any exponention b^n can be rewritten as e^(n ln b). */
  double retval;
  float z;
//crumb("pow: %f to %f\n", base, exponent);
  z = ln(base);
//crumb("pow: z = %f\n", z);
  z = exponent * z;
//crumb("pow:  z = %f\n", z);
  retval = exp(z);
//crumb("pow: done %f\n", retval);
//  retval = exp(exponent * ln(base));
  return retval;
}




#ifndef Q3_VM
int
main ()
{
  printf("test ln(0.015) = %.57f\n", ln(0.015));

  printf("antilogarithm(0.224455) = %.57f\n", antiln(0.224455));
  printf("antilogarithm(0.98) = %.57f\n", antiln(0.98));
  printf("antilogarithm(1.1) = %.57f\n", antiln(1.1));
  printf("antilogarithm(100) = %.57f\n", antiln(100));
  printf("2^5 = %.57f\n", pow(2, 5));
  printf("2^10 = %.57f\n", pow(2, 10));
  printf("10^3 = %.57f\n", pow(10, 3));
}
#endif