zmath.java

Java实现的各种数学算法

// zmath.java jplewis 97
// modified
// jan02	big in eig33
// feb01	huge bugfix in inv2x2, was only working for positive det!

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Library General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.
// 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Library General Public License for more details.
// 
// You should have received a copy of the GNU Library General Public
// License along with this library; if not, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
// Boston, MA  02111-1307, USA.

package ZS;

import java.io.PrintWriter;

import VisualNumerics.math.*;	// for testing inverse 

import zlib.*;

/**
 * Misc. math routines,
 * adapted from Mmath.c.
 * routines formerly called Mxx
 *
 * @author jplewis
 */

/****************
#ifdef FPOW
float fpow(b,e)
float b,e;
{ return( exp(log(b)*e) ); }
#endif

#ifdef IMOD
int imod(a,b)
int a,b;
{ return(a - ((a/b)*b)); }
#endif 
****************/

public final class zmath 
{
  public static final  float PI = (float)Math.PI;
  public static final  float TWOPI = (float)(2.0*Math.PI);
  public static final  float DEGRAD = (float)Math.PI / 180.f;
  public static final  float RADDEG = 180.f / (float)Math.PI;

  public static final int imax(int a, int b)
  {
    if (a > b) return(a); else return(b);
  }

  public static final int imin(int a,int b)
  {
    if (a < b) return(a); else return(b);
  }

  public static final int ifloor(double f)
  {
    if ((float)(int)f == f) return((int)f);  /* otherwise floor(-1) => -2 */
    if (f >= 0.0) return((int)f);
    return( (int)f - 1 );
  }

  public static final int round(double a)
  {
    if (a >= 0.0)
	return((int)(a + 0.5));
    else
	return((int)(a - 0.5));
  }

  public static final double fmax(double a,double b)
  {
    if (a > b) return(a); else return(b);
  }

  public static final double fmin(double a, double b)
  {
    if (a < b) return(a); else return(b);
  }

  /**
   * IEEEremainder is screwy?  returns negative numbers
   */
  public static final float rem(float a, float b)
  {
    if (a >= b) {
      int n = (int)(a / b);
      a = a - n*b;
    }
    return a;
  }

  /* just guessing */
  private static float FAPX = 0.00001f;
  private static double DAPX = 0.00000001;

  public static final boolean fapprox(float a,float b)
  {
    double diff;
    diff = a - b;
    if (diff < 0.0) diff = - diff;
    return(diff < FAPX);
  }

  public static final boolean dapprox(double a,double b)
  {
    double diff;
    diff = a - b;
    if (diff < 0.0) diff = - diff;
    return(diff < DAPX);
  }


  public static final int ipower(int x, int n)
  {
    // if x==2 return 1<<n
    int p = 1;
    for( ; n>0; --n)  p *= x;
    return(p);
  }

  public static final double log10(double x)
  {
    return Math.log(x) / Math.log(10.0);
  }


  /**
   * return bounding power of 2
   */
  public static final int nextpowerof2(int len)	
  {
    int plen = 1;
    int ipow = 0;
    while(len > plen) {
	plen *= 2;
	ipow ++;
    }
    return(ipow);
  } /*nextpowerof2*/


  public static final int getpowerof2(int len)
  {
    int p = nextpowerof2(len);
    if (1<<p != len) zliberror.warning("getpowerof2");
    return(p);
  }

  /**
   */
  public static final boolean ispowerof2(int len)
  {
    return (len == 1<<nextpowerof2(len));
  }


  /**
   * round len to the next multiple of mult,
   * e.g. len=9, mult=4, result=12
   */
  public static final int nextmultipleof(int len, int mult)
  {
    return mult*((len + mult - 1) / mult);
  }


  /**
   */
  public static final int binom(int n,int k)
  {
    int b=1;
    for (int i=1; i<=k; i++) b = b*n--/i;
    return(b);
  }

  static void testBinom()
  {
    int n = 20;
    for(int k=1;k<n;k++)
      cstdio.printf("binom(%d,%d) -> %d\n",n,k,binom(n,k));
  }


  /**
   */
  public static final int fact(int n)
  {
     int f=1;
     for (; n>1; --n) f *= n;
     return(f);
  }


  /**
   * greatest common divisor
   */
  public static final int gcd(int a,int b)
  {
    return (a<b) ? gcd(b,a) : b>0 ? gcd(b,a%b) : a;
  }


  /**
   */
  public static final double hippo(double a,double b)	
  {
    return(Math.sqrt(a*a+b*b));
  }

  /**
   * fast inverse of 2x2 matrix.
   * Tested, looks good.
   */
  public static final float inv2x2(float M00, float M01, float M10, float M11,
				   float[][] MI)
  {
    float det = M00 * M11 - M10 * M01;
    if (det != 0.f) {
      MI[0][0] =  M11 / det;
      MI[0][1] = -M01 / det;
      MI[1][0] = -M10 / det;
      MI[1][1] =  M00 / det;
    }
    return det;
  } //inv2x2


  static void testInv2x2()
  {
    float[][] M = new float[2][2];
    float[][] MI = new float[2][2];
    double[][] dM = new double[2][2];

    for( int t=0; t < 100; t++ ) {

      // build random matrix
      for( int r=0; r < 2; r++ ) {
	for( int c=0; c < 2; c++ ) {
	  //M[r][c] = rnd.rndf();
	  M[r][c] = rnd.rndf11();
	  dM[r][c] = M[r][c];
	}
      }

      // symmetric case
      //dM[1][0] = M[1][0] = M[0][1];

      try {
	double[][] dMI = DoubleMatrix.inverse(dM);
	float det = inv2x2(M[0][0],M[0][1],M[1][0],M[1][1], MI);
	System.out.println("det = " + det);

	double maxdiff = 0.;
	for( int r=0; r < 2; r++ ) {
	  for( int c=0; c < 2; c++ ) {
	    double d = MI[r][c] - dMI[r][c];
	    if (d < 0.) d = - d;
	    if (d > maxdiff) maxdiff = d;
	  } //c
	} //r

	System.out.println("maxdiff = " + maxdiff);
	if (maxdiff > 0.01) {
	  matrix.print(new PrintWriter(System.out), "MI", MI);
	  matrix.print(new PrintWriter(System.out), "dMI", dMI);
	}
      }
      catch(Exception x) {
	zliberror.die(x);
      }

    } //t
  } //testInv2x2


  /**
   * Return eigenvalues of a 2x2 SYMMETRIC matrix in v[].
   * If e[] is not null, also return the leading eigenvector.
   * The second eigenvector is of course orthogonal to the first.
   * Return value is false if the matrix is near singular.
   */
  public static final boolean eigvalues22(final float[][] m,
					  float[] v, float[] e)
  {
    float a = m[0][0];
    float b = m[0][1];
    zliberror._assert(b == m[1][0]);
    float c = m[1][1];

    return eigvalues22(a,b,c, v, e);
  }

  public static final boolean eigvalues22(final float[][] m, float[] v)
  {
    return eigvalues22(m, v, null);
  }

  public static final boolean eigvalues22(float a, float b, float c, float[] v)
  {
    return eigvalues22(a, b, c, null);
  }


  /**
   * Return eigenvalues of a 2x2 SYMMETRIC matrix in v[].
   * If e[] is not null, also return the leading eigenvector.
   * The second eigenvector is of course orthogonal to the first.
   * Return value is false if the matrix is near singular.
   *
   * Tested against EigenJacobi routine for a number of random symmetric
   * positive/negative matrices, works!
   * <pre>
   *  [a b]
   *  [b c]
   */
  public static final boolean eigvalues22(float a, float b, float c,
					  float[] v, float[] e)
  {
    boolean singular = false;
    //System.out.println("\neigvalues abc = "+a+" "+b+" "+c);

    float dmin, dmax;
    if (a > c) {
      dmax = a;
      dmin = c;
    }
    else {
      dmax = c;
      dmin = a;
    }

    dmax *= 0.01f;