package nxm.sys.libm;

/**
  Class to hold primary functions not in Math.

  @author Jeff Schoen
*/
public class MMath {

  /** This constructor was automatically provided by Java,
   *  but it is unneeded as all the methods are static methods.
   *  @deprecated now that Javadocs are required, it is temporarily
   *              being explicitely provided but deprecated
   */
  @Deprecated
  public MMath () {}

  /** get the next power of 2 greater than or equal to the argument
   *  @param arg argument (int)
   *  @return the next power of 2 greater than or equal to the argument
   */
  public static int power2GE (int arg) {
    int value;
    if (arg<=1) value=1;
    else if (arg<32768)
      for (value=2; value<arg; value+=value);
    else if (arg<(1<<30))
      for (value=32768; value<arg; value+=value);
    else value = 0;
    return value;
  }

  /** get the next power of 2 greater than or equal to the argument
   *  @param arg argument (double)
   *  @return the next power of 2 greater than or equal to the argument
   */
  public static double power2GE (double arg) {
    double value;
    if (arg<=1) value=1;
    else if (arg<32768)
      for (value=2; value<arg; value+=value);
    else
      for (value=32768; value<arg; value+=value);
    return value;
  }

  /** perform normalized least squares fit of given order on data vector from n1 to n2
   *  @param b     input values
   *  @param order order of polynomial
   *  @param data  result of least squares
   *  @param ndata number of data points
   *  @param n1    starting n
   *  @param n2    ending n
   *  @return 0
   */
  public static int leastSquaresFit (double[] b, int order, double[] data, int ndata, int n1, int n2) {
    int ncoef=order+1; double dx = 1.0/ndata;
    double[] xstor = new double[ncoef*2], a = new double[ncoef*ncoef];
    for (int i=0; i<ncoef; i++) b[i]=0; // zero the coefficient vectors
    // create metrics
    /*
    Uncomment to use calculate faster at the cost of ndata memory usage
    Calculates vertically, rather than horizontally, so xstor holds a column of data
    and a and b are written to only once, though data is read 2*ncoef times.
    Scales much better with order, but only gives small speed improvements for very small order.
    Around order 20, there's a ~7x speedup in this section (and gaussian elimination has
    negligible cost at that order)
     */
    /*
    xstor = new double[ndata];
    int i, j, k, start = 0;
    double sum1, sum2;
    for (i = 0; i < ndata; i++) {
      xstor[i] = 1;
    }
    for (i = 0; i < 2*ncoef; i++) {
      if ( i < ncoef ) {
        sum1 = 0;
        for (j = n1, k = 0; j < n2; j++, k++) {
          sum1 += data[j]*xstor[k];
        }
        b[i] = sum1;
      }
      else { start++; }
      sum2 = 0;
      for (j = n1, k = 0; j < n2; j++, k++) {
        sum2 += xstor[k];
        xstor[k] *= j*dx;
      }
      int pos = start*(ncoef-1)+i;
      for (j = i-start; j >= start; j--) {
        a[pos] = sum2;
        pos += order; //ncoef-1, to move down one and left one--to get the diagonal
      }
    }
    /*/
    for (int j=n1; j<n2; j++) {
      double xx=1, x=j*dx;
      for (int i=0; i<2*ncoef; i++, xx*=x) xstor[i]=xx;
      for (int i=0; i<ncoef; i++) {
        b[i] += data[j]*xstor[i];
        for (int l=0; l<ncoef; l++) a[i*ncoef+l] += xstor[i+l];
      }
    }//*/
    // solve equations
    return solveByGaussianElimination(a,b,ncoef,0.0);
  }

  /** evaluate normalized least squares
   *  @param b     input values
   *  @param order order of polynomial
   *  @param data  result of least squares
   *  @param ndata number of data points
   *  @return 0
   */
  public static int leastSquaresEval (double[] b, int order, double[] data, int ndata) {
    double dx = 1.0/ndata;
    for (int j=0; j<ndata; j++) {
      double x=j*dx, xx=x, value=b[0];
      for (int i=1; i<=order; i++, xx*=x) value += xx*b[i];
      data[j]=value;
    }
    return 0;
  }

  /** denormalize coefficients
   *  @param b      output data buffer
   *  @param order  order of polynomial
   *  @param xstart X-axis start value from input DataFile's header
   *  @param xdelta X-axis delta value from input DataFile's header
   *  @param ndata  size of (number of elements in) input DataFile
   */
  public static void deNormalizeCoef (double[] b, int order, double xstart, double xdelta, int ndata) {
    int ncoef = order+1; double xdd=1;
    int[] pt = getPascalsTriangle(ncoef);
    double xd = 1.0/(xdelta*ndata), xs = -xstart*xd;
    for (int i=0; i<ncoef; i++,xdd*=xd) { double d=0;
      for (int j=order; j>i; j--) d = (d+pt[j+i*ncoef]*b[j]) * xs;
      b[i] = (d+b[i]) * xdd;  // since next iteration does not depend on d[i-1]
    }
  }

  /** Returns the Pascal's Triangle (triangular array containing values of binomial coefficient)
      @param ncoef number of coefficients (number of rows in this case)
      @return array representation of Pascal's Triangle
   */
  public static int[] getPascalsTriangle(int ncoef) {
    int[] pt = new int[ncoef*ncoef];
    for (int i=0; i<ncoef; i++) {
      pt[i]=1;
      for (int j=1; j<i; j++) pt[i+j*ncoef] = pt[i-1+j*ncoef] + pt[i-1+(j-1)*ncoef];
      pt[i+i*ncoef]=1;
    }
    return pt;
  }

  /** Solves a system of equations by gaussian elimination (also known as row reduction)
      @param a     first row of coefficients
      @param b     second row of coefficients
      @param ncoef number of coefficients in each row
      @param tolr  tolerance range of absolute value of coefficients
      @return 0 if solved; 1 if could not be solved due to tolerance
   */
  public static int solveByGaussianElimination (double[] a, double[] b, int ncoef, double tolr) {
    int i,j,k,imax=0;
    double bigA,temp;
    if (ncoef<1) return 0;
    for (j=0; j<ncoef; j++) {
      // find column maximum
      bigA=0;
      for (i=j; i<ncoef; i++) {
        double aij = a[i+j*ncoef];
        if (Math.abs(aij) > Math.abs(bigA)) { bigA=aij; imax=i; }
      }
      if (Math.abs(bigA)<tolr) return 1;
      // interchange rows and divide by maximum
      for (i=j; i<ncoef; i++) {
        int ji = j+i*ncoef;
        int imaxi = imax+i*ncoef;
        temp = a[imaxi];
        a[imaxi] = a[ji];
        a[ji] = temp/bigA;
      }
      temp = b[imax];
      b[imax] = b[j];
      b[j] = temp/bigA;
      // eliminate next variable
      for (i=j+1; i<ncoef; i++) {
        //*
        temp = a[i + j*ncoef];
        int pos1 = i + (j+1)*ncoef;
        int pos2 = j + (j+1)*ncoef;
        for (k=j+1; k<ncoef; k++) {
          a[pos1] -= temp*a[pos2];
          pos1 += ncoef;
          pos2 += ncoef;
        }
        b[i] -= temp*b[j];
        /*/
        for (k=j+1; k<ncoef; k++) a[i+k*ncoef] -= a[i+j*ncoef]*a[j+k*ncoef];
        b[i] -= a[i+j*ncoef]*b[j];
        //*/
      }
    }
    // backward solution for all other variables
    for (i=ncoef-2; i>=0; i--) {
      for (j=ncoef-1; j>i; j--) b[i] -= a[i+j*ncoef]*b[j];
    }
    return 0;
  }


  /** Floored modulo function.  The sign on the output follows the divisor.
  <pre>
    c = ModuloFloor(a,b)

    The ModuloFloor is defined as:

      c = a%b + I*b

    where
      I = if signum(a%b) = -signum(b) then 1 else 0


    Examples:

        a  |  b  |  c
        -------------
        8  |  3  |  2
        8  | -3  | -1
       -8  |  3  |  1
       -8  | -3  | -2
        ---|-----|---
        1  |  2  |  1
        1  | -2  | -1
       -1  |  2  |  1
       -1  | -2  | -1

  </pre>
    @param a The dividend
    @param b The divisor
    @return the floored modulo
  */
  public static double ModuloFloor(double a, double b) {
    //equivalent statements
    //*/
    return (a*b < 0) ? a%b + b : a%b;
    /*/
    double c = a%b;
    if ( (c>0 && b<0) || (c<0 && b>0)) c = c + b;
    return c;
    //*/
  }

  /** Euclidean modulo function.
  <pre>
    c = ModuloEuclid(a,b)

    The ModuloEuclid is defined as:

      c = a%b + I*b

    where
      I = if ((a%b)&gt;=0) then 0 else if (b&lt;0) then 1 else -1


    Examples:

        a  |  b  |  c
        -------------
        8  |  3  |  2
        8  | -3  |  2
       -8  |  3  |  1
       -8  | -3  |  1
        ---|-----|---
        1  |  2  |  1
        1  | -2  |  1
       -1  |  2  |  1
       -1  | -2  |  1

  </pre>
    @param a dividend
    @param b divisor
    @return the modulo result
  */
  public static double ModuloEuclid(double a, double b) {
    //equivalent statements
    //*/
    return (a < 0) ? a%b + Math.abs(b) : a%b;
    /*/
    double c = a%b;
    if (c < 0) {
      if (b > 0) c = c + b;
      else       c = c - b;
    }
    return c;
    //*/
  }

  //////////////////////////////////////////////////////////////////////////////////////////////////
  // Fast Log Methods
  //
  // WARNING: LOG_BITS, LOG_SHIFT, LOG_TABLE, LOG2_E must be defined in that order!
  //////////////////////////////////////////////////////////////////////////////////////////////////
  private static final int     LOG_BITS  = 14;             // Number of mantissa bits to use for fastLog
  private static final int     LOG_SHIFT = 24 - LOG_BITS;  // Bit shift to get LOG_BITS bits in mantissa
  private static final float[] LOG_TABLE = bldLogTable();  // The fastLog look-up-table
  private static final float   LOG2_E    = fastLog2((float)Math.E); // fastLog2(E) constant used by fastLog

  /** Builds the log table. */
  private static final float[] bldLogTable() {
    float[] logs = new float[0x1 << LOG_BITS];
    double  log2 = Math.log(2.0);

    for (int i = 1; i < logs.length; i++) {
      int   bits = (i << LOG_SHIFT) | 0x40000000; // 0x40000000 sets the exponent to 0
      float val  = Float.intBitsToFloat(bits);

      logs[i] = (float)((Math.log((double)val) / log2) - 1.0);
    }

    return logs;
  }

  /** Computes a fast approximation of the natural log of a value. This is computed as
      <pre>fastLog(x) = fastLog2(x) / fastLog2(Math.E)</pre>
      since the {@link #fastLog2(float)} method is fast.
      @param x The value to take the log of.
      @return An approximation of <code>log(x)</code>.
      @since NeXtMidas 2.3.0
   */
  public static final float fastLog(float x) {
    return fastLog2(x) / LOG2_E;
  }

  /** Computes a fast approximation of the log, base-2, of a value. This is computed as follows:
      <pre>
        1. Express x in 2^exp*man form - this is already encoded as part
           of the binary representation of the float.
             LG(x) = LG(2^exp * man)        |x| = 2^exp * man

        2. Use the definition of the log function to simplify.
             LG(x) = LG(2^exp) + LG(man)
             LG(x) = exp + LG(man)

        3. Since LG(man)&lt;&lt;exp, we can approximate LG(man) using a simple
           look-up-table.
             - To make things easier, the size of the look-up-table is
               2^N and thus we can use the first N bits of man as our
               index into the look-up-table.
             - For this method use N=12 so the result is accurate within
               a tolerance of +/-0.0002 for most values.
      </pre>
      since the {@link #fastLog2(float)} method is fast.
      @param x The value to take the log of.
      @return An approximation of <code>log2(x)</code>.
      @since NeXtMidas 2.3.0
   */
  public static final float fastLog2(float x) {
    if ((x <= 0.0f) || Float.isNaN(x)) { return Float.NaN; }
    int val = Float.floatToRawIntBits(x);  // The float in raw hex notation
    int man = val & 0x7FFFFF;              // The mantissa
    int exp = ((val >> 23) & 0xFF) - 127;  // The exponent (biased by 127)

    return exp + LOG_TABLE[man >>> LOG_SHIFT];
  }
}
