andreacavalli/WarpPI
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

First commit

Browse Source
This commit is contained in:
Gatecraft 2016-03-15 10:41:39 +01:00
parent b4f3bb8ba9
commit c3668828fc
42 changed files with 10401 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
.classpath Normal file
View File

@ -0,0 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<classpath>
<classpathentry kind="src" path="src"/>
<classpathentry kind="con" path="org.eclipse.jdt.launching.JRE_CONTAINER"/>
<classpathentry kind="output" path="bin"/>
</classpath>

17
.project Normal file
View File

@ -0,0 +1,17 @@
<?xml version="1.0" encoding="UTF-8"?>
<projectDescription>
<name>PICalculator</name>
<comment></comment>
<projects>
</projects>
<buildSpec>
<buildCommand>
<name>org.eclipse.jdt.core.javabuilder</name>
<arguments>
</arguments>
</buildCommand>
</buildSpec>
<natures>
<nature>org.eclipse.jdt.core.javanature</nature>
</natures>
</projectDescription>

8
.settings/org.eclipse.core.resources.prefs Normal file
View File

@ -0,0 +1,8 @@
eclipse.preferences.version=1
encoding//src/org/nevec/rjm/BigSurd.java=UTF-8
encoding//src/org/nevec/rjm/BigSurdVec.java=UTF-8
encoding//src/org/warpgate/pi/calculator/Main.java=UTF-8
encoding//src/org/warpgate/pi/calculator/Parentesi.java=UTF-8
encoding//src/org/warpgate/pi/calculator/Radice.java=UTF-8
encoding//src/org/warpgate/pi/calculator/RadiceQuadrata.java=UTF-8
encoding//src/org/warpgate/pi/calculator/Simboli.java=UTF-8

88
out.txt Normal file
View File

@ -0,0 +1,88 @@
•Analyzing expression:((+35(2√7))+(9(2√13)))/21
•Added implicit multiplications:((35*(2√7))+(9*(2√13)))/21
•Subdivision in classes:
•Analyzing expression:(35*(2√7))+(9*(2√13))
•Added implicit multiplications:(35*(2√7))+(9*(2√13))
•Subdivision in classes:
•Analyzing expression:35*(2√7)
•Added implicit multiplications:35*(2√7)
•Subdivision in classes:
•Added variable to expression:35
•Added variable to expression:*
•Analyzing expression:2√7
•Subdivision in classes:
•Added variable to expression:2
•Added variable to expression:√
•Added variable to expression:7
•Finished the subdivision in classes.
•Pushing classes...
•Correcting classes:
•Set variable to expression:√
var1=2
var2=7
(result)=2√7
•Finished correcting classes.
•Result:2√7
•Added variable to expression:Parentesi
•Finished the subdivision in classes.
•Pushing classes...
•Correcting classes:
•Set variable to expression:*
var1=35
var2=2√7
(result)=35*(2√7)
•Finished correcting classes.
•Result:35*(2√7)
•Added variable to expression:Parentesi
•Added variable to expression:+
•Analyzing expression:9*(2√13)
•Added implicit multiplications:9*(2√13)
•Subdivision in classes:
•Added variable to expression:9
•Added variable to expression:*
•Analyzing expression:2√13
•Subdivision in classes:
•Added variable to expression:2
•Added variable to expression:√
•Added variable to expression:13
•Finished the subdivision in classes.
•Pushing classes...
•Correcting classes:
•Set variable to expression:√
var1=2
var2=13
(result)=2√13
•Finished correcting classes.
•Result:2√13
•Added variable to expression:Parentesi
•Finished the subdivision in classes.
•Pushing classes...
•Correcting classes:
•Set variable to expression:*
var1=9
var2=2√13
(result)=9*(2√13)
•Finished correcting classes.
•Result:9*(2√13)
•Added variable to expression:Parentesi
•Finished the subdivision in classes.
•Pushing classes...
•Correcting classes:
•Set variable to expression:+
var1=35*(2√7)
var2=9*(2√13)
(result)=35*(2√7)+9*(2√13)
•Finished correcting classes.
•Result:35*(2√7)+9*(2√13)
•Added variable to expression:Parentesi
•Added variable to expression:/
•Added variable to expression:21
•Finished the subdivision in classes.
•Pushing classes...
•Correcting classes:
•Set variable to expression:/
var1=35*(2√7)+9*(2√13)
var2=21
(result)=(35*(2√7)+9*(2√13))/21
•Finished correcting classes.
•Result:(35*(2√7)+9*(2√13))/21

96
src/org/nevec/rjm/Bernoulli.java Normal file
View File

@ -0,0 +1,96 @@
package org.nevec.rjm ;
import java.math.BigInteger;
import java.util.Vector;
/** Bernoulli numbers.
* @since 2006-06-25
* @author Richard J. Mathar
*/
public class Bernoulli
{
/*
* The list of all Bernoulli numbers as a vector, n=0,2,4,....
*/
static Vector<Rational> a = new Vector<Rational>() ;
public Bernoulli()
{
if ( a.size() == 0 )
{
a.add(Rational.ONE) ;
a.add(new Rational(1,6)) ;
}
}
/** Set a coefficient in the internal table.
* @param n the zero-based index of the coefficient. n=0 for the constant term.
* @param value the new value of the coefficient.
*/
protected void set(final int n, final Rational value)
{
final int nindx = n /2 ;
if ( nindx < a.size())
a.set(nindx,value) ;
else
{
while ( a.size() < nindx )
a.add( Rational.ZERO ) ;
a.add(value) ;
}
}
/** The Bernoulli number at the index provided.
* @param n the index, non-negative.
* @return the B_0=1 for n=0, B_1=-1/2 for n=1, B_2=1/6 for n=2 etc
*/
public Rational at(int n)
{
if ( n == 1)
return(new Rational(-1,2)) ;
else if ( n % 2 != 0 )
return Rational.ZERO ;
else
{
final int nindx = n /2 ;
if( a.size() <= nindx )
{
for(int i= 2*a.size() ; i <= n; i+= 2)
set(i, doubleSum(i) ) ;
}
return a.elementAt(nindx) ;
}
}
/* Generate a new B_n by a standard double sum.
* @param n The index of the Bernoulli number.
* @return The Bernoulli number at n.
*/
private Rational doubleSum(int n)
{
Rational resul = Rational.ZERO ;
for(int k=0 ; k <= n ; k++)
{
Rational jsum = Rational.ZERO ;
BigInteger bin = BigInteger.ONE ;
for(int j=0 ; j <= k ; j++)
{
BigInteger jpown = (new BigInteger(""+j)).pow(n);
if ( j % 2 == 0)
jsum = jsum.add(bin.multiply(jpown)) ;
else
jsum = jsum.subtract(bin.multiply(jpown)) ;
/* update binomial(k,j) recursively
*/
bin = bin.multiply( new BigInteger(""+(k-j))). divide( new BigInteger(""+(j+1)) ) ;
}
resul = resul.add(jsum.divide(new BigInteger(""+(k+1)))) ;
}
return resul ;
}
} /* Bernoulli */

188
src/org/nevec/rjm/BigComplex.java Normal file
View File

@ -0,0 +1,188 @@
package org.nevec.rjm ;
import java.math.BigDecimal;
import java.math.MathContext;
/** Complex numbers with BigDecimal real and imaginary components
* @since 2008-10-26
* @author Richard J. Mathar
*/
public class BigComplex
{
/** real part
*/
BigDecimal re ;
/** imaginary part
*/
BigDecimal im ;
/** The constant that equals zero
*/
final static BigComplex ZERO = new BigComplex(BigDecimal.ZERO, BigDecimal.ZERO) ;
/** Default ctor equivalent to zero.
*/
public BigComplex()
{
re= BigDecimal.ZERO ;
im= BigDecimal.ZERO ;
}
/** ctor with real and imaginary parts
* @param x real part
* @param y imaginary part
*/
public BigComplex( BigDecimal x, BigDecimal y)
{
re=x ;
im=y ;
}
/** ctor with real part.
* @param x real part.
* The imaginary part is set to zero.
*/
public BigComplex( BigDecimal x )
{
re=x ;
im= BigDecimal.ZERO ;
}
/** ctor with real and imaginary parts
* @param x real part
* @param y imaginary part
*/
public BigComplex( double x, double y)
{
re= new BigDecimal(x) ;
im= new BigDecimal(y) ;
}
/** Multiply with another BigComplex
* @param oth The BigComplex which is a factor in the product
* @param mc Defining precision and rounding mode
* @return This multiplied by oth
* @since 2010-07-19 implemented with 3 multiplications and 5 additions/subtractions
*/
BigComplex multiply(final BigComplex oth, MathContext mc)
{
final BigDecimal a = re.add(im).multiply(oth.re) ;
final BigDecimal b = oth.re.add(oth.im).multiply(im) ;
final BigDecimal c = oth.im.subtract(oth.re).multiply(re) ;
final BigDecimal x = a.subtract(b,mc) ;
final BigDecimal y = a.add(c,mc) ;
return new BigComplex(x,y) ;
}
/** Add a BigDecimal
* @param oth the value to be added to the real part.
* @return this added to oth
*/
BigComplex add(final BigDecimal oth)
{
final BigDecimal x = re.add(oth) ;
return new BigComplex(x,im) ;
}
/** Subtract another BigComplex
* @param oth the value to be subtracted from this.
* @return this minus oth
*/
BigComplex subtract(final BigComplex oth)
{
final BigDecimal x = re.subtract(oth.re) ;
final BigDecimal y = im.subtract(oth.im) ;
return new BigComplex(x,y) ;
}
/** Complex-conjugation
* @return the complex conjugate of this.
*/
BigComplex conj()
{
return new BigComplex(re,im.negate()) ;
}
/** The absolute value squared.
* @return The sum of the squares of real and imaginary parts.
* This is the square of BigComplex.abs() .
*/
BigDecimal norm()
{
return re.multiply(re).add(im.multiply(im)) ;
}
/** The absolute value.
* @return the square root of the sum of the squares of real and imaginary parts.
* @since 2008-10-27
*/
BigDecimal abs(MathContext mc)
{
return BigDecimalMath.sqrt(norm(),mc) ;
}
/** The square root.
* @return the square root of the this.
* The branch is chosen such that the imaginary part of the result has the
* same sign as the imaginary part of this.
* @see Tim Ahrendt, <a href="http://dx.doi.org/10.1145/236869.236924">Fast High-precision computation of complex square roots</a>,
* ISSAC 1996 p142-149.
* @since 2008-10-27
*/
BigComplex sqrt(MathContext mc)
{
final BigDecimal half = new BigDecimal("2") ;
/* compute l=sqrt(re^2+im^2), then u=sqrt((l+re)/2)
* and v= +- sqrt((l-re)/2 as the new real and imaginary parts.
*/
final BigDecimal l = abs(mc) ;
if ( l.compareTo(BigDecimal.ZERO) == 0 )
return new BigComplex( BigDecimalMath.scalePrec(BigDecimal.ZERO,mc),
BigDecimalMath.scalePrec(BigDecimal.ZERO,mc) ) ;
final BigDecimal u = BigDecimalMath.sqrt( l.add(re).divide(half,mc), mc );
final BigDecimal v = BigDecimalMath.sqrt( l.subtract(re).divide(half,mc), mc );
if ( im.compareTo(BigDecimal.ZERO)>= 0 )
return new BigComplex(u,v) ;
else
return new BigComplex(u,v.negate()) ;
}
/** The inverse of this.
* @return 1/this
*/
BigComplex inverse(MathContext mc)
{
final BigDecimal hyp = norm() ;
/* 1/(x+iy)= (x-iy)/(x^2+y^2 */
return new BigComplex( re.divide(hyp,mc), im.divide(hyp,mc).negate() ) ;
}
/** Divide through another BigComplex number.
* @return this/oth
*/
BigComplex divide(BigComplex oth, MathContext mc)
{
/* lazy implementation: (x+iy)/(a+ib)= (x+iy)* 1/(a+ib) */
return multiply(oth.inverse(mc),mc) ;
}
/** Human-readable Fortran-type display
* @return real and imaginary part in parenthesis, divided by a comma.
*/
public String toString()
{
return "("+re.toString()+","+im.toString()+")" ;
}
/** Human-readable Fortran-type display
* @return real and imaginary part in parenthesis, divided by a comma.
*/
public String toString(MathContext mc)
{
return "("+re.round(mc).toString()+","+im.round(mc).toString()+")" ;
}
} /* BigComplex */

2758
src/org/nevec/rjm/BigDecimalMath.java Normal file
View File

File diff suppressed because it is too large Load Diff

546
src/org/nevec/rjm/BigIntegerMath.java Normal file
View File

@ -0,0 +1,546 @@
package org.nevec.rjm ;
import java.math.BigInteger;
import java.util.Vector;
/** BigInteger special functions and Number theory.
* @since 2009-08-06
* @author Richard J. Mathar
*/
public class BigIntegerMath
{
/** Evaluate binomial(n,k).
* @param n The upper index
* @param k The lower index
* @return The binomial coefficient
*/
static public BigInteger binomial(final int n, final int k)
{
if ( k == 0 )
return(BigInteger.ONE) ;
BigInteger bin = new BigInteger(""+n) ;
BigInteger n2 = bin ;
for(BigInteger i= new BigInteger(""+(k-1)) ; i.compareTo(BigInteger.ONE) >= 0 ; i = i.subtract(BigInteger.ONE) )
bin = bin.multiply(n2.subtract(i)) ;
for(BigInteger i= new BigInteger(""+k) ; i.compareTo(BigInteger.ONE) == 1 ; i = i.subtract(BigInteger.ONE) )
bin = bin.divide(i) ;
return ( bin) ;
} /* binomial */
/** Evaluate binomial(n,k).
* @param n The upper index
* @param k The lower index
* @return The binomial coefficient
* @since 2008-10-15
*/
static public BigInteger binomial(final BigInteger n, final BigInteger k)
{
/* binomial(n,0) =1
*/
if ( k.compareTo(BigInteger.ZERO) == 0 )
return(BigInteger.ONE) ;
BigInteger bin = new BigInteger(""+n) ;
/* the following version first calculates n(n-1)(n-2)..(n-k+1)
* in the first loop, and divides this product through k(k-1)(k-2)....2
* in the second loop. This is rather slow and replaced by a faster version
* below
* BigInteger n2 = bin ;
* BigInteger i= k.subtract(BigInteger.ONE) ;
* for( ; i.compareTo(BigInteger.ONE) >= 0 ; i = i.subtract(BigInteger.ONE) )
* bin = bin.multiply(n2.subtract(i)) ;
* i= new BigInteger(""+k) ;
* for( ; i.compareTo(BigInteger.ONE) == 1 ; i = i.subtract(BigInteger.ONE) )
* bin = bin.divide(i) ;
*/
/* calculate n then n(n-1)/2 then n(n-1)(n-2)(2*3) etc up to n(n-1)..(n-k+1)/(2*3*..k)
* This is roughly the best way to keep the individual intermediate products small
* and in the integer domain. First replace C(n,k) by C(n,n-k) if n-k<k.
*/
BigInteger truek = new BigInteger(k.toString()) ;
if ( n.subtract(k).compareTo(k) < 0 )
truek = n.subtract(k) ;
/* Calculate C(num,truek) where num=n and truek is the smaller of n-k and k.
* Have already initialized bin=n=C(n,1) above. Start definining the factorial
* in the denominator, named fden
*/
BigInteger i = new BigInteger("2") ;
BigInteger num = new BigInteger(n.toString()) ;
/* a for-loop (i=2;i<= truek;i++)
*/
for( ; i.compareTo(truek) <= 0 ; i = i.add(BigInteger.ONE) )
{
/* num = n-i+1 after this operation
*/
num = num.subtract(BigInteger.ONE) ;
/* multiply by (n-i+1)/i
*/
bin = (bin.multiply(num)).divide(i) ;
}
return ( bin) ;
} /* binomial */
/** Evaluate sigma_k(n).
* @param n the main argument which defines the divisors
* @param k the lower index, which defines the power
* @return The sum of the k-th powers of the positive divisors
*/
static public BigInteger sigmak(final BigInteger n, final int k)
{
return (new Ifactor(n.abs())).sigma(k).n ;
} /* sigmak */
/** Evaluate sigma(n).
* @param n the argument for which divisors will be searched.
* @return the sigma function. Sum of the positive divisors of the argument.
* @since 2006-08-14
* @author Richard J. Mathar
*/
static public BigInteger sigma(int n)
{
return (new Ifactor(Math.abs(n))).sigma().n ;
}
/** Compute the list of positive divisors.
* @param n The integer of which the divisors are to be found.
* @return The sorted list of positive divisors.
* @since 2010-08-27
* @author Richard J. Mathar
*/
static public Vector<BigInteger> divisors(final BigInteger n)
{
return (new Ifactor(n.abs())).divisors() ;
}
/** Evaluate sigma(n).
* @param n the argument for which divisors will be searched.
* @return the sigma function. Sum of the divisors of the argument.
* @since 2006-08-14
* @author Richard J. Mathar
*/
static public BigInteger sigma(final BigInteger n)
{
return (new Ifactor(n.abs())).sigma().n ;
}
/** Evaluate floor(sqrt(n)).
* @param n The non-negative argument.
* @return The integer square root. The square root rounded down.
* @since 2010-08-27
* @author Richard J. Mathar
*/
static public int isqrt(final int n)
{
if ( n < 0 )
throw new ArithmeticException("Negative argument "+ n) ;
final double resul= Math.sqrt((double)n) ;
return (int)Math.round(resul) ;
}
/** Evaluate floor(sqrt(n)).
* @param n The non-negative argument.
* Arguments less than zero throw an ArithmeticException.
* @return The integer square root, the square root rounded down.
* @since 2010-08-27
* @author Richard J. Mathar
*/
static public long isqrt(final long n)
{
if ( n < 0 )
throw new ArithmeticException("Negative argument "+ n) ;
final double resul= Math.sqrt((double)n) ;
return Math.round(resul) ;
}
/** Evaluate floor(sqrt(n)).
* @param n The non-negative argument.
* Arguments less than zero throw an ArithmeticException.
* @return The integer square root, the square root rounded down.
* @since 2011-02-12
* @author Richard J. Mathar
*/
static public BigInteger isqrt(final BigInteger n)
{
if ( n.compareTo(BigInteger.ZERO) < 0 )
throw new ArithmeticException("Negative argument "+ n.toString()) ;
/* Start with an estimate from a floating point reduction.
*/
BigInteger x ;
final int bl = n.bitLength() ;
if ( bl > 120)
x = n.shiftRight(bl/2-1) ;
else
{
final double resul= Math.sqrt(n.doubleValue()) ;
x = new BigInteger(""+Math.round(resul)) ;
}
final BigInteger two = new BigInteger("2") ;
while ( true)
{
/* check whether the result is accurate, x^2 =n
*/
BigInteger x2 = x.pow(2) ;
BigInteger xplus2 = x.add(BigInteger.ONE).pow(2) ;
if ( x2.compareTo(n) <= 0 && xplus2.compareTo(n) > 0)
return x ;
xplus2 = xplus2.subtract(x.shiftLeft(2)) ;
if ( xplus2.compareTo(n) <= 0 && x2.compareTo(n) > 0)
return x.subtract(BigInteger.ONE) ;
/* Newton algorithm. This correction is on the
* low side caused by the integer divisions. So the value required
* may end up by one unit too large by the bare algorithm, and this
* is caught above by comparing x^2, (x+-1)^2 with n.
*/
xplus2 = x2.subtract(n).divide(x).divide(two) ;
x = x.subtract(xplus2) ;
}
}
/** Evaluate core(n).
* Returns the smallest positive integer m such that n/m is a perfect square.
* @param n The non-negative argument.
* @return The square-free part of n.
* @since 2011-02-12
* @author Richard J. Mathar
*/
static public BigInteger core(final BigInteger n)
{
if ( n.compareTo(BigInteger.ZERO) < 0 )
throw new ArithmeticException("Negative argument "+ n) ;
final Ifactor i = new Ifactor(n) ;
return i.core() ;
}
/** Minor of an integer matrix.
* @param A The matrix.
* @param r The row index of the row to be removed (0-based).
* An exception is thrown if this is outside the range 0 to the upper row index of A.
* @param c The column index of the column to be removed (0-based).
* An exception is thrown if this is outside the range 0 to the upper column index of A.
* @return The depleted matrix. This is not a deep copy but contains references to the original.
* @since 2010-08-27
* @author Richard J. Mathar
*/
static public BigInteger[][] minor(final BigInteger[][] A, final int r, final int c) throws ArithmeticException
{
/* original row count */
final int rL = A.length ;
if ( rL == 0 )
throw new ArithmeticException("zero row count in matrix") ;
if ( r < 0 || r >= rL)
throw new ArithmeticException("row number "+r + " out of range 0.." + (rL-1)) ;
/* original column count */
final int cL = A[0].length ;
if ( cL == 0 )
throw new ArithmeticException("zero column count in matrix") ;
if ( c < 0 || c >= cL)
throw new ArithmeticException("column number "+c + " out of range 0.." + (cL-1)) ;
BigInteger M[][] = new BigInteger[rL-1][cL-1] ;
int imrow =0 ;
for (int row = 0 ; row < rL ; row++)
{
if ( row != r)
{
int imcol = 0 ;
for(int col = 0 ; col < cL ;col++)
{
if ( col != c )
{
M[imrow][imcol] = A[row][col] ;
imcol ++ ;
}
}
imrow++ ;
}
}
return M ;
}
/** Replace column of a matrix with a column vector.
* @param A The matrix.
* @param c The column index of the column to be substituted (0-based).
* @param v The column vector to be inserted.
* With the current implementation, it must be at least as long as the row count, and
* its elements that exceed that count are ignored.
* @return The modified matrix. This is not a deep copy but contains references to the original.
* @since 2010-08-27
* @author Richard J. Mathar
*/
@SuppressWarnings("unused")
static private BigInteger[][] colSubs(final BigInteger[][] A, final int c, final BigInteger[] v) throws ArithmeticException
{
/* original row count */
final int rL = A.length ;
if ( rL == 0 )
throw new ArithmeticException("zero row count in matrix") ;
/* original column count */
final int cL = A[0].length ;
if ( cL == 0 )
throw new ArithmeticException("zero column count in matrix") ;
if ( c < 0 || c >= cL)
throw new ArithmeticException("column number "+c + " out of range 0.." + (cL-1)) ;
BigInteger M[][] = new BigInteger[rL][cL] ;
for (int row = 0 ; row < rL ; row++)
{
for(int col = 0 ; col < cL ;col++)
{
/* currently, v may just be longer than the row count, and surplus
* elements will be ignored. Shorter v lead to an exception.
*/
if ( col != c )
M[row][col] = A[row][col] ;
else
M[row][col] = v[row] ;
}
}
return M ;
}
/** Determinant of an integer square matrix.
* @param A The square matrix.
* If column and row dimensions are unequal, an ArithmeticException is thrown.
* @return The determinant.
* @since 2010-08-27
* @author Richard J. Mathar
*/
static public BigInteger det(final BigInteger[][] A) throws ArithmeticException
{
BigInteger d = BigInteger.ZERO ;
/* row size */
final int rL = A.length ;
if ( rL == 0 )
throw new ArithmeticException("zero row count in matrix") ;
/* column size */
final int cL = A[0].length ;
if ( cL != rL )
throw new ArithmeticException("Non-square matrix dim "+rL + " by " + cL) ;
/* Compute the low-order cases directly.
*/
if ( rL == 1 )
return A[0][0] ;
else if ( rL == 2)
{
d = A[0][0].multiply(A[1][1]) ;
return d.subtract( A[0][1].multiply(A[1][0])) ;
}
else
{
/* Work arbitrarily along the first column of the matrix */
for (int r = 0 ; r < rL ; r++)
{
/* Do not consider minors that do no contribute anyway
*/
if ( A[r][0].compareTo(BigInteger.ZERO) != 0 )
{
final BigInteger M[][] = minor(A,r,0) ;
final BigInteger m = A[r][0].multiply( det(M)) ;
/* recursive call */
if ( r % 2 == 0)
d = d.add(m) ;
else
d = d.subtract(m) ;
}
}
}
return d;
}
/** Solve a linear system of equations.
* @param A The square matrix.
* If it is not of full rank, an ArithmeticException is thrown.
* @param rhs The right hand side. The length of this vector must match the matrix size;
* else an ArithmeticException is thrown.
* @return The vector of x in A*x=rhs.
* @since 2010-08-28
* @author Richard J. Mathar
*/
static public Rational[] solve(final BigInteger[][]A, final BigInteger[] rhs) throws ArithmeticException
{
final int rL = A.length ;
if ( rL == 0 )
throw new ArithmeticException("zero row count in matrix") ;
/* column size */
final int cL = A[0].length ;
if ( cL != rL )
throw new ArithmeticException("Non-square matrix dim "+rL + " by " + cL) ;
if ( rhs.length != rL )
throw new ArithmeticException("Right hand side dim "+ rhs.length + " unequal matrix dim " + rL) ;
/* Gauss elimination
*/
Rational x[] = new Rational[rL] ;
/* copy of r.h.s ito a mutable Rationalright hand side
*/
for(int c = 0 ; c < cL ; c++)
x[c] = new Rational(rhs[c]) ;
/* Create zeros downwards column c by linear combination of row c and row r.
*/
for(int c = 0 ; c < cL-1 ; c++)
{
/* zero on the diagonal? swap with a non-zero row, searched with index r */
if ( A[c][c].compareTo(BigInteger.ZERO) == 0)
{
boolean swpd = false ;
for(int r=c+1; r< rL ; r++)
{
if ( A[r][c].compareTo(BigInteger.ZERO) != 0)
{
for(int cpr =c ; cpr < cL; cpr++)
{
BigInteger tmp = A[c][cpr] ;
A[c][cpr] = A[r][cpr] ;
A[r][cpr] = tmp ;
}
Rational tmp = x[c] ;
x[c] = x[r] ;
x[r] = tmp ;
swpd = true ;
break;
}
}
/* not swapped with a non-zero row: determinant zero and no solution
*/
if ( ! swpd)
throw new ArithmeticException("Zero determinant of main matrix") ;
}
/* create zero at A[c+1..cL-1][c] */
for( int r=c+1; r < rL ; r++)
{
/* skip the cpr=c which actually sets the zero: this element is not visited again
*/
for(int cpr = c+1; cpr < cL; cpr++)
{
BigInteger tmp = A[c][c].multiply(A[r][cpr]) .subtract ( A[c][cpr].multiply(A[r][c])) ;
A[r][cpr] = tmp ;
}
Rational tmp = x[r].multiply(A[c][c]) .subtract ( x[c].multiply(A[r][c])) ;
x[r] = tmp ;
}
}
if ( A[cL-1][cL-1].compareTo(BigInteger.ZERO) == 0)
throw new ArithmeticException("Zero determinant of main matrix") ;
/* backward elimination */
for( int r = cL-1 ; r >= 0 ; r--)
{
x[r] = x[r].divide(A[r][r]) ;
for(int rpr = r-1 ; rpr >=0 ; rpr--)
x[rpr] = x[rpr].subtract( x[r].multiply(A[rpr][r]) ) ;
}
return x ;
}
/** The lowest common multiple
* @param a The first argument
* @param b The second argument
* @return lcm(|a|,|b|)
* @since 2010-08-27
* @author Richard J. Mathar
*/
static public BigInteger lcm(final BigInteger a, final BigInteger b)
{
BigInteger g = a.gcd(b) ;
return a.multiply(b).abs().divide(g) ;
}
/** Evaluate the value of an integer polynomial at some integer argument.
* @param c Represents the coefficients c[0]+c[1]*x+c[2]*x^2+.. of the polynomial
* @param x The abscissa point of the evaluation
* @return The polynomial value.
* @since 2010-08-27
* @author Richard J. Mathar
*/
static public BigInteger valueOf(final Vector<BigInteger>c, final BigInteger x)
{
if (c.size() == 0)
return BigInteger.ZERO ;
BigInteger res = c.lastElement() ;
for(int i= c.size()-2 ; i >=0 ; i--)
res = res.multiply(x).add( c.elementAt(i) ) ;
return res ;
}
/** The central factorial number t(n,k) number at the indices provided.
* @param n the first parameter, non-negative.
* @param k the second index, non-negative.
* @return t(n,k)
* @since 2009-08-06
* @author Richard J. Mathar
* @see <a href="http://dx.doi.org/10.1080/01630568908816313">P. L. Butzer et al, Num. Funct. Anal. Opt. 10 (5)( 1989) 419-488</a>
*/
static public Rational centrlFactNumt(int n,int k)
{
if ( k > n || k < 0 || ( k % 2 ) != (n % 2) )
return Rational.ZERO ;
else if ( k == n)
return Rational.ONE ;
else
{
/* Proposition 6.2.6 */
Factorial f = new Factorial() ;
Rational jsum = new Rational(0,1) ;
int kprime = n-k ;
for ( int j =0 ; j <= kprime ; j++)
{
Rational nusum = new Rational(0,1) ;
for(int nu =0 ; nu <= j ; nu++)
{
Rational t = new Rational(j-2*nu,2) ;
t = t.pow(kprime+j) ;
t = t.multiply( binomial(j,nu) ) ;
if ( nu % 2 != 0 )
nusum = nusum.subtract(t) ;
else
nusum = nusum.add(t) ;
}
nusum = nusum.divide( f.at(j) ).divide(n+j) ;
nusum = nusum.multiply( binomial(2*kprime,kprime-j) ) ;
if ( j % 2 != 0 )
jsum = jsum.subtract(nusum) ;
else
jsum = jsum.add(nusum) ;
}
return jsum.multiply(k).multiply( binomial(n+kprime,k) ) ;
}
} /* CentralFactNumt */
/** The central factorial number T(n,k) number at the indices provided.
* @param n the first parameter, non-negative.
* @param k the second index, non-negative.
* @return T(n,k)
* @since 2009-08-06
* @author Richard J. Mathar
* @see <a href="http://dx.doi.org/10.1080/01630568908816313">P. L. Butzer et al, Num. Funct. Anal. Opt. 10 (5)( 1989) 419-488</a>
*/
static public Rational centrlFactNumT(int n,int k)
{
if ( k > n || k < 0 || ( k % 2 ) != (n % 2) )
return Rational.ZERO ;
else if ( k == n)
return Rational.ONE ;
else
{
/* Proposition 2.1 */
return centrlFactNumT(n-2,k-2).add( centrlFactNumT(n-2,k).multiply(new Rational(k*k,4)) ) ;
}
} /* CentralFactNumT */
} /* BigIntegerMath */

668
src/org/nevec/rjm/BigIntegerPoly.java Normal file
View File

@ -0,0 +1,668 @@
package org.nevec.rjm ;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.util.Scanner;
import java.util.Vector;
/** Polynomial with integer coefficients.
* Alternatively to be interpreted as a sequence which has the polynomial as an (approximate)
* generating function.
* @since 2010-08-27
* @author Richard J. Mathar
*/
public class BigIntegerPoly implements Cloneable
{
/** The list of all coefficients, starting with a0, then a1, as in
* poly=a0+a1*x+a2*x^2+a3*x^3+...
*/
Vector<BigInteger> a ;
/** Default ctor.
* Creates the polynomial p(x)=0.
*/
public BigIntegerPoly()
{
a = new Vector<BigInteger>() ;
}
/** Ctor with a comma-separated list as the list of coefficients.
* @param L the string of the form a0,a1,a2,a3 with the coefficients
*/
public BigIntegerPoly(final String L) throws NumberFormatException
{
a = new Vector<BigInteger>() ;
Scanner sc = new Scanner(L) ;
sc.useDelimiter(",") ;
while ( sc.hasNextBigInteger())
a.add(sc.nextBigInteger()) ;
simplify() ;
sc.close();
} /* ctor */
/** Ctor with a list of coefficients.
* @param c The coefficients a0, a1, a2 etc in a0+a1*x+a2*x^2+...
*/
@SuppressWarnings("unchecked")
public BigIntegerPoly(final Vector<BigInteger> c)
{
a = (Vector<BigInteger>)c.clone() ;
simplify() ;
} /* ctor */
/** Ctor with a list of coefficients.
* @param c The coefficients a0, a1, a2 etc in a0+a1*x+a2*x^2+...
*/
public BigIntegerPoly(final BigInteger[] c)
{
for(int i=0 ; i < c.length; i++)
a.add( c[i].add(BigInteger.ZERO) ) ;
simplify() ;
} /* ctor */
/** Create a copy of this.
* @since 2010-08-27
*/
public BigIntegerPoly clone()
{
return new BigIntegerPoly(a) ;
} /* clone */
/** Translate into a RatPoly copy.
* @since 2012-03-02
*/
public RatPoly toRatPoly()
{
RatPoly bd = new RatPoly() ;
for(int i=0 ; i < a.size() ; i++)
bd.set(i, a.elementAt(i) ) ;
return bd;
} /* toRatPoly */
/** Retrieve a polynomial coefficient.
* @param n the zero-based index of the coefficient. n=0 for the constant term.
* @return the polynomial coefficient in front of x^n.
*/
public BigInteger at(final int n)
{
if ( n < a.size())
return( a.elementAt(n) ) ;
else
return( BigInteger.ZERO ) ;
} /* at */
/** Evaluate at some integer argument.
* @param x The abscissa point of the evaluation
* @return The polynomial value.
* @since 2010-08-27
* @author Richard J. Mathar
*/
public BigInteger valueOf(final BigInteger x)
{
if (a.size() == 0)
return BigInteger.ZERO ;
BigInteger res = a.lastElement() ;
/* Heron casted form
*/
for(int i= a.size()-2 ; i >=0 ; i--)
res = res.multiply(x).add( a.elementAt(i) ) ;
return res ;
} /* valueOf */
/** Horner scheme to find the function value at the argument x
* @param x The argument x.
* @return Value of the polynomial at x.
* @since 2008-11-13
*/
public BigInteger valueOf( int x)
{
return valueOf(new BigInteger(""+x)) ;
} /* valueOf */
/** Set a polynomial coefficient.
* @param n the zero-based index of the coefficient. n=0 for the constant term.
* If the polynomial has not yet the degree to need this coefficient,
* the intermediate coefficients are set to zero.
* @param value the new value of the coefficient.
*/
public void set(final int n, final BigInteger value)
{
if ( n < a.size())
a.set(n,value) ;
else
{
/* fill intermediate powers with coefficients of zero
*/
while ( a.size() < n )
{
a.add(BigInteger.ZERO ) ;
}
a.add(value) ;
}
} /* set */
/** Set a polynomial coefficient.
* @param n the zero-based index of the coefficient. n=0 for the constant term.
* If the polynomial has not yet the degree to need this coefficient,
* the intermediate coefficients are implicitly set to zero.
* @param value the new value of the coefficient.
*/
public void set(final int n, final int value)
{
BigInteger val2 = new BigInteger(""+value) ;
set(n,val2) ;
} /* set */
/** Count of coefficients.
* @return the number of polynomial coefficients.
* Differs from the polynomial degree by one.
*/
public int size()
{
return a.size() ;
} /* size */
/** Polynomial degree.
* @return the polynomial degree.
*/
public int degree()
{
return a.size()-1 ;
} /* degree */
/** Polynomial lower degree.
* @return power of the smallest non-zero coefficient.
* If the polynomial is identical to 0, 0 is returned.
*/
public int ldegree()
{
for(int n=0 ; n < a.size() ; n++)
if ( a.elementAt(n).compareTo(BigInteger.ZERO) != 0 )
return n;
return 0 ;
} /* ldegree */
/** Multiply by a constant factor.
* @param val the factor
* @return the product of this with the factor.
* All coefficients of this have been multiplied individually by the factor.
* @since 2010-08-27
*/
public BigIntegerPoly multiply(final BigInteger val)
{
BigIntegerPoly resul = new BigIntegerPoly() ;
if ( val.compareTo(BigInteger.ZERO) != 0 )
for(int n=0; n < a.size() ; n++)
resul.set(n,a.elementAt(n).multiply(val) ) ;
return resul ;
} /* multiply */
/** Multiply by another polynomial
* @param val the other polynomial
* @return the product of this with the other polynomial
*/
public BigIntegerPoly multiply(final BigIntegerPoly val)
{
BigIntegerPoly resul = new BigIntegerPoly() ;
/* the degree of the result is the sum of the two degrees.
*/
final int nmax = degree()+val.degree() ;
for(int n=0; n <= nmax ; n++)
{
BigInteger coef = BigInteger.ZERO ;
for(int nleft=0; nleft <= n ; nleft++)
coef = coef.add(at(nleft).multiply(val.at(n-nleft))) ;
resul.set(n,coef) ;
}
resul.simplify() ;
return resul ;
} /* multiply */
/** Raise to a positive power.
* @param n the exponent of the power
* @return the n-th power of this.
*/
public BigIntegerPoly pow(final int n) throws ArithmeticException
{
BigIntegerPoly resul = new BigIntegerPoly("1") ;
if ( n < 0 )
throw new ArithmeticException("negative polynomial power "+n) ;
else
{
for(int i=1 ; i <= n ; i++)
resul = resul.multiply(this) ;
resul.simplify() ;
return resul ;
}
} /* pow */
/** Add another polynomial
* @param val the other polynomial
* @return the sum of this with the other polynomial
* @since 2010-08-27
*/
public BigIntegerPoly add(final BigIntegerPoly val)
{
BigIntegerPoly resul = new BigIntegerPoly() ;
/* the degree of the result is the larger of the two degrees (before simplify() at least).
*/
final int nmax = (degree()>val.degree()) ? degree() : val.degree() ;
for(int n=0; n <= nmax ; n++)
{
BigInteger coef = at(n).add(val.at(n)) ;
resul.set(n,coef) ;
}
resul.simplify() ;
return resul ;
} /* add */
/** Subtract another polynomial
* @param val the other polynomial
* @return the difference between this and the other polynomial
* @since 2008-10-25
*/
public BigIntegerPoly subtract(final BigIntegerPoly val)
{
BigIntegerPoly resul = new BigIntegerPoly() ;
/* the degree of the result is the larger of the two degrees (before simplify() at least).
*/
final int nmax = (degree()>val.degree()) ? degree() : val.degree() ;
for(int n=0; n <= nmax ; n++)
{
BigInteger coef = at(n).subtract(val.at(n)) ;
resul.set(n,coef) ;
}
resul.simplify() ;
return resul ;
} /* subtract */
/** Divide by another polynomial.
* @param val the other polynomial
* @return A vector with [0] containg the polynomial of degree which is the
* difference of the degree of this and the degree of val. [1] the remainder polynomial.
* This = returnvalue[0] + returnvalue[1]/val .
* @since 2012-03-01
*/
public BigIntegerPoly[] divideAndRemainder(final BigIntegerPoly val)
{
BigIntegerPoly[] ret = new BigIntegerPoly[2] ;
/* remove any high-order zeros. note that the clone() operation calls simplify().
*/
BigIntegerPoly valSimpl = val.clone() ;
BigIntegerPoly thisSimpl = clone() ;
/* catch the case with val equal to zero
*/
if ( valSimpl.degree() == 0 && valSimpl.a.firstElement().compareTo(BigInteger.ZERO) == 0)
throw new ArithmeticException("Division through zero polynomial") ;
/* degree of this smaller than degree of val: remainder is this
*/
if ( thisSimpl.degree() < valSimpl.degree() )
{
/* leading polynomial equals zero
*/
ret[0] = new BigIntegerPoly() ;
ret[1] = thisSimpl ;
}
else
{
/* long division. Highest degree by dividing the highest degree
* of this thru val. At this point an exception is thrown if the
* polynomial division cannot be done with integer coefficients.
*/
ret[0] = new BigIntegerPoly() ;
BigInteger[] newc = thisSimpl.a.lastElement().divideAndRemainder( valSimpl.a.lastElement()) ;
if ( newc[1].compareTo(BigInteger.ZERO) != 0)
throw new ArithmeticException("Incompatible leading term in " + this + " / " + val) ;
ret[0].set( thisSimpl.degree()-valSimpl.degree(), newc[0]) ;
/* recurrences: build this - val*(1-termresult) and feed this
* into another round of division. Have intermediate ret[0]+ret[1]/val.
*/
ret[1] = thisSimpl.subtract( ret[0].multiply( valSimpl) );
/* any remainder left ?
*/
if ( ret[1].degree() < valSimpl.degree() )
;
else
{
BigIntegerPoly rem[] = ret[1].divideAndRemainder(val) ;
ret[0] = ret[0].add(rem[0]) ;
ret[1] = rem[1] ;
}
}
return ret ;
} /* divideAndRemainder */
/** Print as a comma-separated list of coefficients.
* @return the representation a0,a1,a2,a3,...
* @since 2010-08-27
*/
public String toString()
{
String str = new String();
for(int n=0; n < a.size() ; n++)
{
if ( n == 0 )
str += a.elementAt(n).toString() ;
else
str += ","+a.elementAt(n).toString() ;
}
if ( str.length() == 0 )
str = "0" ;
return str ;
} /* toString */
/** Print as a polyomial in x.
* @return The representation a0+a1*x+a2*x^2+...
* The terms with zero coefficients are not mentioned.
* @since 2008-10-26
*/
public String toPString()
{
String str = new String();
for(int n=0; n < a.size() ; n++)
{
final BigInteger num = a.elementAt(n) ;
if ( num.compareTo(BigInteger.ZERO) != 0 )
{
str += " " ;
if ( num.compareTo(BigInteger.ZERO) > 0 && n> 0)
str += "+" ;
str += a.elementAt(n).toString() ;
if ( n > 0 )
{
str += "*x" ;
if ( n > 1 )
str += "^"+n ;
}
}
}
if ( str.length() == 0 )
str = "0" ;
return str ;
} /* toPString */
/** Simplify the representation.
* Trailing values with zero coefficients (at high powers) are deleted.
*/
protected void simplify()
{
int n = a.size()-1 ;
if ( n >= 0)
while( a.elementAt(n).compareTo(BigInteger.ZERO) == 0 )
{
a.removeElementAt(n) ;
if( --n <0)
break ;
}
} /* simplify */
/** First derivative.
* @return The first derivative with respect to the indeterminate variable.
* @since 2008-10-26
*/
public BigIntegerPoly derive()
{
if ( a.size() <= 1)
{
/* derivative of the constant is just zero
*/
return new BigIntegerPoly() ;
}
else
{
BigIntegerPoly d = new BigIntegerPoly() ;
for(int i=1 ; i <= degree() ; i++)
{
final BigInteger c = a.elementAt(i).multiply(new BigInteger(""+i)) ;
d.set(i-1,c) ;
}
return d ;
}
} /* derive */
/** Truncate polynomial degree.
* @return The polynomial with all coefficients beyond deg set to zero.
* @since 2010-08-27
*/
public BigIntegerPoly trunc(int newdeg)
{
BigIntegerPoly t = new BigIntegerPoly() ;
for(int i=0; i <= newdeg; i++)
t.set(i,at(i)) ;
t.simplify() ;
return t ;
} /* trunc */
/** Inverse Binomial transform.
* @param maxdeg the maximum polynomial degree of the result
* @return the sequence of coefficients is the inverse binomial transform of the original sequence.
* @since 2010-08-29
*/
public BigIntegerPoly binomialTInv(int maxdeg)
{
BigIntegerPoly r = new BigIntegerPoly() ;
for(int i=0; i <= maxdeg; i++)
{
BigInteger c = BigInteger.ZERO ;
for(int j=0; j <= i && j < a.size(); j++)
if ( (j+i) % 2 != 0 )
c = c.subtract( a.elementAt(j).multiply(BigIntegerMath.binomial(i,j)) ) ;
else
c = c.add( a.elementAt(j).multiply(BigIntegerMath.binomial(i,j)) ) ;
r.set(i,c) ;
}
r.simplify() ;
return r ;
} /* binomialTInv */
/** Compute the order of the root r.
* @return 1 for simple roots, 2 for order 2 etc., 0 if not a root
* @since 2010-08-27
*/
public int rootDeg(final BigInteger r)
{
int o = 0 ;
BigIntegerPoly d = clone() ;
BigInteger roo = d.valueOf(r) ;
while ( roo.compareTo(BigInteger.ZERO) == 0 )
{
o++ ;
d = d.derive() ;
roo = d.valueOf(r) ;
}
return o ;
} /* rootDeg */
/** Generate the integer roots of the polynomial.
* @return The vector of integer roots, without their multiplicity.
* @since 2010-08-27
*/
public Vector<BigInteger> iroots()
{
/* The vector of the roots */
Vector<BigInteger> res =new Vector<BigInteger>() ;
/* collect the zero
*/
if ( a.firstElement().compareTo(BigInteger.ZERO) == 0 )
res.add(BigInteger.ZERO) ;
/* collect the divisors of the constant element (or the reduced polynomial) */
int l = ldegree() ;
if ( a.elementAt(l).compareTo(BigInteger.ZERO) != 0 )
{
Vector<BigInteger> cand = BigIntegerMath.divisors(a.elementAt(l).abs()) ;
/* check the divisors (both signs) */
for(int i=0 ; i < cand.size() ; i++)
{
BigInteger roo = valueOf( cand.elementAt(i) ) ;
if ( roo.compareTo(BigInteger.ZERO) == 0 )
/* found a root cand[i] */
res.add(cand.elementAt(i)) ;
roo = valueOf( cand.elementAt(i).negate() ) ;
if ( roo.compareTo(BigInteger.ZERO) == 0 )
res.add(cand.elementAt(i).negate()) ;
}
}
return res;
} /* iroots */
/** Generate the factors which are 2nd degree polynomials.
* @return A (potentially empty) vector of factors, without multiplicity.
* Only factors with non-zero absolute coefficient are generated.
* This means the factors are of the form x^2+a*x+b=0 with nonzero b.
* @since 2012-03-01
*/
protected Vector<BigIntegerPoly> i2roots()
{
/* The vector of the factors to be returned
*/
Vector<BigIntegerPoly> res =new Vector<BigIntegerPoly>() ;
if ( degree() < 2)
return res ;
BigInteger bsco = a.firstElement().abs() ;
Vector<BigInteger> b = BigIntegerMath.divisors(bsco) ;
BigInteger csco = a.lastElement().abs() ;
Vector<BigInteger> c = BigIntegerMath.divisors(csco) ;
/* Generate the floating point values of roots. To have some reasonable
* accuracy in the results, add zeros to the integer coefficients, scaled
* by the expected division with values of b (which are all <= a.firstele).
* Number of decimal digits in bsco by using a log2->log10 rough estimate
* and adding 6 safety digits
*/
RatPoly thisDec = toRatPoly() ;
Vector<BigComplex> roo = thisDec.roots(6+(int)(0.3*bsco.bitCount()) ) ;
final BigDecimal half = new BigDecimal("0.5") ;
/* for each of the roots z try to see whether c*z^2+a*z+b=0 with integer a, b and c
* where b is restricted to a signed divisor of the constant coefficient.
* Solve z*(c*z+a)=-b or c*z+a = -b/z or -b/z-c*z = some integer a.
*/
for( BigComplex z : roo)
{
for(BigInteger bco : b)
for(BigInteger cco : c)
{
/* the major reason to avoid the case b=0 is that this would
* require precaution of double counting below. Note that this
* case is already covered by using iroots().
*/
if ( bco.signum() != 0 )
{
for(int sig = -1 ; sig <=1 ; sig +=2)
{
BigInteger bcosig = (sig > 0 )? bco : bco.negate() ;
/* -a = b/z+c*z has real part b*Re(z)/|z|^2+c*Re(z) = Re z *( b/|z|^2+c)
*/
BigDecimal negA = BigDecimalMath.add(BigDecimalMath.divideRound(bcosig,z.norm()),cco) ;
negA = negA.multiply(z.re) ;
/* convert to a with round-to-nearest
*/
BigInteger a = negA.negate().add(half).toBigInteger() ;
/* test the polynomial remainder. if zero, add the term
* to the results.
*/
BigIntegerPoly dtst = new BigIntegerPoly(""+bcosig+","+a+","+cco) ;
try
{
BigIntegerPoly[] rm = divideAndRemainder(dtst) ;
if ( rm[1].isZero() )
res.add(dtst) ;
}
catch ( ArithmeticException ex)
{
}
}
}
}
}
return res;
} /* i2roots */
/** Test whether this polynomial value is zero.
* @return If this is a polynomial p(x)=0 for all x.
*/
public boolean isZero()
{
simplify() ;
return (a.size() ==0 ) ;
}
/** Factorization into integer polynomials.
* The current factorization detects only factors which are polynomials of order up to 2.
* @return The vector of factors. Factors with higher multiplicity are represented by repetition.
* @since 2012-03-01
*/
public Vector<BigIntegerPoly> ifactor()
{
/* this ought be entirely rewritten in terms of the LLL algorithm
*/
Vector<BigIntegerPoly> fac = new Vector<BigIntegerPoly>() ;
/* collect integer roots (polynomial factors of degree 1) */
Vector<BigInteger> r = iroots() ;
BigIntegerPoly[] res = new BigIntegerPoly[2] ;
res[0] = this ;
for( BigInteger i : r)
{
int deg = rootDeg(i) ;
/* construct the factor x-i */
BigIntegerPoly f = new BigIntegerPoly(""+i.negate()+",1") ;
for(int mu =0 ; mu < deg ; mu++)
{
fac.add(f) ;
res = res[0].divideAndRemainder(f) ;
}
}
/* collect factors which are polynomials of degree 2
*/
Vector<BigIntegerPoly> pol2 = i2roots() ;
for( BigIntegerPoly i : pol2)
{
/* the internal loop catches cases with higher
* powers of individual polynomials (of actual degree 2 or 4...)
*/
while ( res[0].degree() >= 2)
{
try
{
BigIntegerPoly[] dtst = res[0].divideAndRemainder(i) ;
if ( dtst[1].isZero() )
{
fac.add(i) ;
res = dtst ;
}
else
break ;
}
catch(ArithmeticException ex)
{
break ;
}
}
}
/* add remaining factor, if not equal to 1
*/
if ( res[0].degree() >0 || res[0].a.firstElement().compareTo (BigInteger.ONE) != 0 )
fac.add(res[0]) ;
return fac ;
} /* ifactor */
} /* BigIntegerPoly */

478
src/org/nevec/rjm/BigSurd.java Normal file
View File

@ -0,0 +1,478 @@
package org.nevec.rjm ;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.security.ProviderException;
import org.warpgate.pi.calculator.Utils;
/** Square roots on the real line.
* These represent numbers which are a product of a (signed) fraction by
* a square root of a non-negative fraction.
* This might be extended to values on the imaginary axis by allowing negative
* values underneath the square root, but this is not yet implemented.
* @since 2011-02-12
* @author Richard J. Mathar
*/
public class BigSurd implements Cloneable, Comparable<BigSurd>
{
/** The value of zero.
*/
static public BigSurd ZERO = new BigSurd() ;
/** The value of one.
*/
static public BigSurd ONE = new BigSurd(Rational.ONE,Rational.ONE) ;
/** Prefactor
*/
Rational pref ;
/** The number underneath the square root, always non-negative.
* The mathematical object has the value pref*sqrt(disc).
*/
Rational disc ;
/** Default ctor, which represents the zero.
* @since 2011-02-12
*/
public BigSurd()
{
pref = Rational.ZERO ;
disc = Rational.ZERO ;
}
/** ctor given the prefactor and the basis of the root.
* This creates an object of value a*sqrt(b).
* @param a the prefactor.
* @param b the discriminant.
* @since 2011-02-12
*/
public BigSurd(Rational a, Rational b)
{
this.pref = a ;
/* reject attempts to use a negative b
*/
if ( b.signum() < 0 )
throw new ProviderException("Not implemented: imaginary surds") ;
this.disc = b ;
normalize() ;
normalizeG() ;
}
/** ctor given the numerator and denominator of the root.
* This creates an object of value sqrt(a/b).
* @param a the numerator
* @param b the denominator.
* @since 2011-02-12
*/
public BigSurd(int a, int b)
{
this( Rational.ONE, new Rational(a,b) ) ;
}
/** ctor given the value under the root.
* This creates an object of value sqrt(a).
* @param a the discriminant.
* @since 2011-02-12
*/
public BigSurd(BigInteger a)
{
this( Rational.ONE, new Rational(a,BigInteger.ONE) ) ;
}
public BigSurd(Rational a)
{
this( Rational.ONE, a) ;
}
/** Create a deep copy.
* @since 2011-02-12
*/
public BigSurd clone()
{
Rational fclon = pref.clone() ;
Rational dclon = disc.clone() ;
/* the main intent here is to bypass any attempt to reduce the discriminant
* by figuring out the square-free part in normalize(), which has already done
* in the current copy of the number.
*/
BigSurd cl = new BigSurd() ;
cl.pref = fclon ;
cl.disc = dclon ;
return cl ;
} /* BigSurd.clone */
/** Add two surds of compatible discriminant.
* @param val The value to be added to this.
*/
public BigSurdVec add(final BigSurd val)
{
//zero plus somethings yields something
if ( signum() == 0 )
return new BigSurdVec(val) ;
else if (val.signum() == 0 )
return new BigSurdVec(this) ;
else
// let the ctor of BigSurdVec to the work
return new BigSurdVec(this,val) ;
} /* BigSurd.add */
/** Multiply by another square root.
* @param val a second number of this type.
* @return the product of this with the val.
* @since 2011-02-12
*/
public BigSurd multiply(final BigSurd val)
{
return new BigSurd( pref.multiply(val.pref), disc.multiply(val.disc) ) ;
} /* BigSurd.multiply */
/** Multiply by a rational number.
* @param val the factor.
* @return the product of this with the val.
* @since 2011-02-15
*/
public BigSurd multiply(final Rational val)
{
return new BigSurd( pref.multiply(val), disc) ;
} /* BigSurd.multiply */
/** Multiply by a BigInteger.
* @param val a second number.
* @return the product of this with the value.
* @since 2011-02-12
*/
public BigSurd multiply(final BigInteger val)
{
return new BigSurd(pref.multiply(val), disc) ;
} /* BigSurd.multiply */
/** Multiply by an integer.
* @param val a second number.
* @return the product of this with the value.
* @since 2011-02-12
*/
public BigSurd multiply(final int val)
{
BigInteger tmp = new BigInteger(""+val) ;
return multiply(tmp) ;
} /* BigSurd.multiply */
/** Compute the square.
* @return this value squared.
* @since 2011-02-12
*/
public Rational sqr()
{
Rational res = pref.pow(2) ;
res = res.multiply(disc) ;
return res;
} /* BigSurd.sqr */
/** Divide by another square root.
* @param val A second number of this type.
* @return The value of this/val
* @since 2011-02-12
*/
public BigSurd divide(final BigSurd val)
{
if( val.signum() == 0 )
throw new ArithmeticException("Dividing "+ toFancyString() + " through zero.") ;
return new BigSurd( pref.divide(val.pref), disc.divide(val.disc) ) ;
} /* BigSurd.divide */
private String toFancyString() {
BigSurd bs = this;
BigInteger denominator = pref.b;
String s = "";
if (denominator.compareTo(BigInteger.ONE) != 0) {
s += "(";
}
if (bs.isBigInteger()) {
s += bs.BigDecimalValue(new MathContext(Utils.scale, Utils.scaleMode2)).toBigInteger().toString();
} else if (bs.isRational()) {
s += bs.toRational().toString();
} else {
BigInteger numerator = bs.pref.a;
if (numerator.compareTo(BigInteger.ONE) != 0) {
s += numerator.toString();
s += "*";
s += "(";
}
s += "2√";
if (bs.disc.isInteger()) {
s += bs.disc.toString();
} else {
s += "("+bs.disc.toString()+")";
}
if (numerator.compareTo(BigInteger.ONE) != 0) {
s += ")";
}
}
return s;
}
/** Divide by an integer.
* @param val a second number.
* @return the value of this/val
* @since 2011-02-12
*/
public BigSurd divide(final BigInteger val)
{
if( val.signum() == 0 )
throw new ArithmeticException("Dividing "+ toFancyString() + " through zero.") ;
return new BigSurd( pref.divide(val), disc ) ;
} /* BigSurd.divide */
/** Divide by an integer.
* @param val A second number.
* @return The value of this/val
* @since 2011-02-12
*/
public BigSurd divide(int val)
{
if( val == 0 )
throw new ArithmeticException("Dividing "+ toFancyString() + " through zero.") ;
return new BigSurd( pref.divide(val), disc ) ;
} /* BigSurd.divide */
/** Compute the negative.
* @return -this.
* @since 2011-02-12
*/
public BigSurd negate()
{
/* This is trying to be quick, avoiding normalize(), by toggling
* the sign in a clone()
*/
BigSurd n = clone() ;
n.pref = n.pref.negate() ;
return n ;
} /* BigSurd.negate */
/** Absolute value.
* @return The absolute (non-negative) value of this.
* @since 2011-02-12
*/
public BigSurd abs()
{
return new BigSurd(pref.abs(),disc) ;
}
/** Compares the value of this with another constant.
* @param val the other constant to compare with
* @return -1, 0 or 1 if this number is numerically less than, equal to,
* or greater than val.
* @since 2011-02-12
*/
public int compareTo(final BigSurd val)
{
/* Since we keep the discriminant positive, the rough estimate
* comes from comparing the signs of the prefactors.
*/
final int sig = signum() ;
final int sigv = val.signum() ;
if ( sig < 0 && sigv >= 0 )
return -1 ;
if ( sig > 0 && sigv <= 0 )
return 1 ;
if ( sig == 0 && sigv == 0 )
return 0 ;
if ( sig == 0 && sigv > 0 )
return -1 ;
if ( sig == 0 && sigv < 0 )
return 1 ;
/* Work out the cases of equal sign. Compare absolute values by comparison
* of the squares which is forwarded to the comparison of the Rational class.
*/
final Rational this2 = sqr() ;
final Rational val2 = val.sqr() ;
final int c = this2.compareTo(val2) ;
if ( c == 0 )
return 0 ;
/* If both values have negative sign, the one with the smaller square is the larger number.
*/
else if ( sig >0 && c >0 || sig <0 && c <0 )
return 1;
else
return -1 ;
} /* BigSurd.compareTo */
/** Return a string in the format (number/denom)*()^(1/2).
* If the discriminant equals 1, print just the prefactor.
* @return the human-readable version in base 10
* @since 2011-02-12
*/
public String toString()
{
if ( disc.compareTo(Rational.ONE) != 0 && disc.compareTo(Rational.ZERO) != 0)
return( "("+pref.toString()+")*("+disc.toString()+")^(1/2)" ) ;
else
return pref.toString() ;
} /* BigSurd.toString */
/** Return a double value representation.
* @return The value with double precision.
* @since 2011-02-12
*/
public double doubleValue()
{
/* First compute the square to prevent overflows if the two pieces of
* the prefactor and the discriminant are of very different magnitude.
*/
Rational p2 = pref.pow(2).multiply(disc) ;
System.out.println("dv sq " + p2.toString()) ;
double res = p2.doubleValue() ;
System.out.println("dv sq " + res) ;
return (pref.signum() >= 0) ? Math.sqrt(res) : -Math.sqrt(res) ;
} /* BigSurd.doubleValue */
/** Return a float value representation.
* @return The value with single precision.
* @since 2011-02-12
*/
public float floatValue()
{
return (float)(doubleValue()) ;
} /* BigSurd.floatValue */
/** True if the value is integer.
* Equivalent to the indication whether a conversion to an integer
* can be exact.
* @since 2011-02-12
*/
public boolean isBigInteger()
{
return pref.isBigInteger() && ( disc.signum() ==0 || disc.compareTo(Rational.ONE) == 0 ) ;
} /* BigSurd.isBigInteger */
/** True if the value is rational.
* Equivalent to the indication whether a conversion to a Rational can be exact.
* @since 2011-02-12
*/
public boolean isRational()
{
return ( disc.signum() ==0 || disc.compareTo(Rational.ONE) == 0 ) ;
} /* BigSurd.isRational */
/** Convert to a rational value if possible
* @since 2012-02-15
*/
public Rational toRational()
{
if ( isRational() )
return pref ;
else
throw new ArithmeticException("Undefined conversion "+ toFancyString() + " to Rational.") ;
} /* BigSurd.toRational */
/** The sign: 1 if the number is >0, 0 if ==0, -1 if <0
* @return the signum of the value.
* @since 2011-02-12
*/
public int signum()
{
/* Since the disc is kept positive, this is the same
* as the sign of the prefactor. This works because a zero discriminant
* is always copied over to the prefactor, not hidden.
*/
return pref.signum() ;
} /* BigSurd.signum */
/** Normalize to squarefree discriminant.
* @since 2011-02-12
*/
protected void normalize()
{
/* Move squares out of the numerator and denominator of the discriminant
*/
if ( disc.signum() != 0 )
{
/* square-free part of the numerator: numer = numC*some^2
*/
BigInteger numC = BigIntegerMath.core(disc.numer()) ;
/* extract the perfect square of the numerator
*/
BigInteger sq = disc.numer().divide(numC) ;
/* extract the associated square root
*/
BigInteger sqf = BigIntegerMath.isqrt(sq) ;
/* move sqf over to the pre-factor
*/
pref = pref.multiply(sqf) ;
BigInteger denC = BigIntegerMath.core(disc.denom()) ;
sq = disc.denom().divide(denC) ;
sqf = BigIntegerMath.isqrt(sq) ;
pref = pref.divide(sqf) ;
disc = new Rational(numC,denC) ;
}
else
pref = Rational.ZERO ;
} /* BigSurd.normalize */
/** Normalize to coprime numerator and denominator in prefactor and discriminant
* @since 2011-02-12
*/
protected void normalizeG()
{
/* Is there a common factor between the numerator of the prefactor
* and the denominator of the discriminant ?
*/
BigInteger d = pref.numer().abs().gcd( disc.denom()) ;
if ( d.compareTo(BigInteger.ONE) > 0 )
{
pref = pref.divide(d) ;
/* instead of multiplying with the square of d, using two steps
* offers a change to recognize the common factor..
*/
disc = disc.multiply(d) ;
disc = disc.multiply(d) ;
}
/* Is there a common factor between the denominator of the prefactor
* and the numerator of the discriminant ?
*/
d = pref.denom().gcd( disc.numer()) ;
if ( d.compareTo(BigInteger.ONE) > 0 )
{
pref = pref.multiply(d) ;
/* instead of dividing through the square of d, using two steps
* offers a change to recognize the common factor..
*/
disc = disc.divide(d) ;
disc = disc.divide(d) ;
}
} /* BigSurd.normalizeG */
/** Return the approximate floating point representation.
* @param mc Description of the accuracy needed.
* @return A representation with digits valid as described by mc
* @since 2012-02-15
*/
public BigDecimal BigDecimalValue(MathContext mc)
{
/* the relative error of the result equals the relative error of the
* prefactor plus half of the relative error of the discriminant.
* So adding 3 digits temporarily is sufficient.
*/
final MathContext locmc = new MathContext(mc.getPrecision()+3,mc.getRoundingMode()) ;
/* first the square root of the discriminant
*/
BigDecimal sqrdis = BigDecimalMath.sqrt(disc.BigDecimalValue(locmc),locmc ) ;
/* Then multiply by the prefactor. If sqrdis is a terminating decimal fraction,
* we prevent early truncation of the result by truncating later.
*/
BigDecimal res = sqrdis.multiply(pref.BigDecimalValue(mc)) ;
return BigDecimalMath.scalePrec(res,mc) ;
} /* BigDecimalValue */
} /* BigSurd */

601
src/org/nevec/rjm/BigSurdVec.java Normal file
View File

@ -0,0 +1,601 @@
package org.nevec.rjm;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.util.Vector;
import org.warpgate.pi.calculator.Utils;
/**
* A BigSurdVec represents an algebraic sum or differences of values which each
* term an instance of BigSurd. This mainly means that sums or differences of
* two BigSurd (or two BigSurdVec) can be represented (exactly) as a BigSurdVec.
*
* @since 2012-02-15
* @author Richard J. Mathar
*/
public class BigSurdVec implements Comparable<BigSurdVec> {
/**
* The value of zero.
*/
static public BigSurdVec ZERO = new BigSurdVec();
/**
* The value of one.
*/
static public BigSurdVec ONE = new BigSurdVec(BigSurd.ONE);
/**
* Internal representation: Each term as a single BigSurd. The value zero is
* represented by an empty vector.
*/
Vector<BigSurd> terms;
/**
* Default ctor, which represents the zero.
*
* @since 2012-02-15
*/
public BigSurdVec() {
terms = new Vector<BigSurd>();
} /* ctor */
/**
* ctor given the value of a BigSurd.
*
* @param a
* The value to be represented by this vector.
* @since 2012-02-15
*/
public BigSurdVec(BigSurd a) {
terms = new Vector<BigSurd>(1);
terms.add(a);
} /* ctor */
/**
* ctor given two values, which (when added) represent this number a+b.
*
* @param a
* The value to be represented by the first term of the vector.
* @param b
* The value to be represented by the second term of the vector.
* @since 2012-02-15
*/
public BigSurdVec(BigSurd a, BigSurd b) {
terms = new Vector<BigSurd>(2);
terms.add(a);
terms.add(b);
normalize();
} /* ctor */
/**
* Combine terms that can be written as a single surd. This unites for
* example the terms sqrt(90) and sqrt(10) to 4*sqrt(10).
*
* @since 2012-02-15
*/
protected void normalize() {
/*
* nothing to be done if at most one term
*/
if (terms.size() <= 1)
return;
Vector<BigSurd> newter = new Vector<BigSurd>();
newter.add(terms.firstElement());
/*
* add j-th element to the existing vector and combine were possible
*/
for (int j = 1; j < terms.size(); j++) {
BigSurd todo = terms.elementAt(j);
boolean merged = false;
for (int ex = 0; ex < newter.size(); ex++) {
BigSurd v = newter.elementAt(ex);
/*
* try to merge terms[j] and newter[ex]. todo = r * v with r a
* rational number is needed. Replaces v with v+todo = v*(1+r)
* if this reduction works.
*/
BigSurd r = todo.divide(v);
if (r.isRational()) {
/* compute r+1 */
Rational newpref = r.toRational().add(1);
/*
* eliminate accidental zeros; overwrite with v*(1+r).
*/
if (newpref.compareTo(Rational.ZERO) == 0)
newter.removeElementAt(ex);
else {
v = v.multiply(newpref);
newter.setElementAt(v, ex);
}
merged = true;
break;
}
}
/*
* append if none of the existing elements matched
*/
if (!merged)
newter.add(todo);
}
/* overwrite old version */
terms = newter;
} /* normalize */
/**
* Compare algebraic value with oth. Returns -1, 0 or +1 depending on
* whether this is smaller, equal to or larger than oth.
*
* @param oth
* The value with which this is to be compared.
* @return 0 or +-1.
* @since 2012-02-15
*/
public int compareTo(BigSurdVec oth) {
final BigSurdVec diff = this.subtract(oth);
return diff.signum();
} /* compareTo */
/**
* Sign function. Returns -1, 0 or +1 depending on whether this is smaller,
* equal to or larger than zero.
*
* @return 0 or +-1.
* @since 2012-02-15
*/
public int signum() {
/*
* the case of zero is unique, because no (reduced) vector of surds
* other than the one element 0 itself can add/subtract to zero.
*/
if (terms.size() == 0)
return 0;
/*
* if there is one term: forward to the signum function of BigSurd
*/
if (terms.size() == 1)
return terms.firstElement().signum();
/*
* if all terms have a common sign: take that one offsig is the index of
* the first "offending" term in the sense that its sign doese not agree
* with the term[0].
*/
int sig0 = terms.elementAt(0).signum();
int offsig = 1;
for (; offsig < terms.size(); offsig++)
if (terms.elementAt(offsig).signum() != sig0)
break;
if (offsig >= terms.size())
return sig0;
/*
* if there are two terms (now known to have different sign): forward to
* the comparison of the two elements as BigSurds
*/
if (terms.size() == 2)
return terms.elementAt(0).compareTo(terms.elementAt(1).negate());
/*
* if there are three terms, move the one with the offending sign to the
* other side and square both sides (which looses the sign) to remove
* all but one surds. The difference of the squared sides contains at
* most two terms, which reduces to the case above. t(0)+t(offbar) <>
* -t(offs)
*/
if (terms.size() == 3) {
BigSurdVec lhs;
if (offsig == 2)
lhs = new BigSurdVec(terms.elementAt(0), terms.elementAt(1));
else
lhs = new BigSurdVec(terms.elementAt(0), terms.elementAt(2));
lhs = lhs.sqr();
/*
* Strange line: this line isn't used, but it's present in this code!
*
*
*
BigSurd rhs = new BigSurd(terms.elementAt(offsig).sqr(), Rational.ONE);
*
*
*
*/
if (lhs.compareTo(lhs) > 0)
/*
* dominating sign was t(0)+t(offbar)
*/
return terms.elementAt(0).signum();
else
return terms.elementAt(offsig).signum();
}
/*
* for a larger number of terms: take a floating point representation
* with a small but correct number of digits, and resume with the sign
* of that one.
*/
return (floatValue() > 0.) ? 1 : -1;
} /* signum */
/**
* Construct an approximate floating point representation
*
* @param mc
* The intended accuracy of the result.
* @return A truncated version with the precision described by mc
*/
public BigDecimal BigDecimalValue(MathContext mc) {
/*
* simple cases with one term forwarded to the BigSurd class
*/
if (terms.size() == 0)
return BigDecimal.ZERO;
else if (terms.size() == 1) {
return terms.firstElement().BigDecimalValue(mc);
}
/*
* To reduce cancellation errors, loop over increasing local precision
* until we are stable to the required result. Keep the old (less
* precise) estimate in res[0], and the newer, more precise in res[1].
*/
BigDecimal[] res = new BigDecimal[2];
res[0] = BigDecimal.ZERO;
for (int addpr = 1;; addpr += 3) {
MathContext locmc = new MathContext(mc.getPrecision() + addpr, mc.getRoundingMode());
res[1] = BigDecimal.ZERO;
for (BigSurd j : terms)
res[1] = BigDecimalMath.addRound(res[1], j.BigDecimalValue(locmc));
if (addpr > 1) {
BigDecimal err = res[1].subtract(res[0]).abs();
int prec = BigDecimalMath.err2prec(res[1], err);
if (prec > mc.getPrecision())
break;
}
res[0] = res[1];
}
return BigDecimalMath.scalePrec(res[1], mc);
} /* BigDecimalValue */
/**
* Construct an approximate floating point representation
*
* @return A truncated version with the precision described by mc
*/
public double doubleValue() {
BigDecimal bd = BigDecimalValue(MathContext.DECIMAL128);
return bd.doubleValue();
} /* doubleValue */
/**
* Construct an approximate floating point representation
*
* @return A truncated version with the precision described by mc
*/
public double floatValue() {
BigDecimal bd = BigDecimalValue(MathContext.DECIMAL64);
return bd.floatValue();
} /* floatValue */
/**
* Add two vectors algebraically.
*
* @param val
* The value to be added to this.
* @return The new value representing this+val.
*/
public BigSurdVec add(final BigSurdVec val) {
BigSurdVec sum = new BigSurdVec();
/*
* concatenate the vectors and eliminate common overlaps
*/
for (BigSurd term : terms) {
if (term.compareTo(BigSurd.ZERO) != 0) {
sum.terms.add(term);
}
}
for (BigSurd term : val.terms) {
if (term.compareTo(BigSurd.ZERO) != 0) {
sum.terms.add(term);
}
}
sum.normalize();
return sum;
} /* add */
/**
* Add two vectors algebraically.
*
* @param val
* The value to be added to this.
* @return The new value representing this+val.
*/
public BigSurdVec add(final BigSurd val) {
BigSurdVec sum = new BigSurdVec();
/*
* concatenate the vectors and eliminate common overlaps
*/
sum.terms.addAll(terms);
sum.terms.add(val);
sum.normalize();
return sum;
} /* add */
/**
* Subtract another number.
*
* @param val
* The value to be subtracted from this.
* @return The new value representing this-val.
*/
public BigSurdVec subtract(final BigSurdVec val) {
BigSurdVec sum = new BigSurdVec();
/*
* concatenate the vectors and eliminate common overlaps
*/
sum.terms.addAll(terms);
for (BigSurd s : val.terms)
sum.terms.add(s.negate());
sum.normalize();
return sum;
} /* subtract */
/**
* Subtract another number.
*
* @param val
* The value to be subtracted from this.
* @return The new value representing this-val.
*/
public BigSurdVec subtract(final BigSurd val) {
BigSurdVec sum = new BigSurdVec();
/*
* concatenate the vectors and eliminate common overlaps
*/
sum.terms.addAll(terms);
sum.terms.add(val.negate());
sum.normalize();
return sum;
} /* subtract */
/**
* Compute the negative.
*
* @return -this.
* @since 2012-02-15
*/
public BigSurdVec negate() {
/*
* accumulate the negated elements of term one by one
*/
BigSurdVec resul = new BigSurdVec();
for (BigSurd s : terms)
resul.terms.add(s.negate());
/*
* no normalization step here, because the negation of all terms does
* not introduce new common factors
*/
return resul;
} /* negate */
/**
* Compute the square.
*
* @return this value squared.
* @since 2012-02-15
*/
public BigSurdVec sqr() {
/*
* Binomial expansion. First the sum of the terms squared, then 2 times
* the mixed products.
*/
BigSurdVec resul = new BigSurdVec();
for (int i = 0; i < terms.size(); i++)
resul.terms.add(new BigSurd(terms.elementAt(i).sqr(), Rational.ONE));
for (int i = 0; i < terms.size() - 1; i++)
for (int j = i + 1; j < terms.size(); j++)
resul.terms.add(terms.elementAt(i).multiply(terms.elementAt(j)).multiply(2));
resul.normalize();
return resul;
} /* sqr */
/**
* Multiply by another square root.
*
* @param val
* a second number of this type.
* @return the product of this with the val.
* @since 2011-02-12
*/
public BigSurdVec multiply(final BigSurd val) {
BigSurdVec resul = new BigSurdVec();
for (BigSurd s : terms)
resul.terms.add(s.multiply(val));
resul.normalize();
return resul;
} /* multiply */
public BigSurdVec multiply(final BigSurdVec val) {
BigSurdVec resul = new BigSurdVec();
for (BigSurd s : terms) {
resul.terms.add(s);
}
for (BigSurd s : val.terms) {
resul = resul.multiply(s);
}
return resul;
} /* multiply */
public BigSurdVec divide(final BigSurd val) {
BigSurdVec resul = new BigSurdVec();
for (BigSurd s : terms)
resul.terms.add(s.divide(val));
resul.normalize();
return resul;
} /* multiply */
public BigSurdVec divide(final BigSurdVec val) {
BigSurdVec resul = new BigSurdVec();
resul.terms = terms;
for (BigSurd s : val.terms) {
resul = resul.divide(s);
}
return resul;
} /* divide */
/**
* True if the value is rational. Equivalent to the indication whether a
* conversion to a Rational can be exact.
*
* @since 2011-02-12
*/
public boolean isRational() {
boolean val = false;
for (BigSurd s : terms) {
val = s.isRational();
if (val == false) {
break;
}
}
return val;
} /* BigSurdVec.isRational */
/**
* True if the value is BigInteger. Equivalent to the indication whether a
* conversion to a BigInteger can be exact.
*
* @since 2011-02-12
*/
public boolean isBigInteger() {
boolean val = false;
for (BigSurd s : terms) {
val = s.isBigInteger();
if (val == false) {
break;
}
}
return val;
} /* BigSurdVec.isRational */
/**
* Convert to a rational value if possible
*
* @since 2012-02-15
*/
public Rational toRational() {
Rational rat = Rational.ZERO;
if (isRational() == false)
throw new ArithmeticException("Undefined conversion " + toString() + " to Rational.");
for (BigSurd s : terms) {
rat = rat.add(s.pref);
}
return rat;
} /* BigSurd.toRational */
/**
* Convert to a BigInteger value if possible
*
* @since 2012-02-15
*/
public BigInteger toBigInteger() {
BigDecimal tmp = BigDecimal.ZERO.setScale(Utils.scale, Utils.scaleMode);
if (isBigInteger() == false)
throw new ArithmeticException("Undefined conversion " + toString() + " to Rational.");
for (BigSurd s : terms) {
tmp = BigDecimalMath.addRound(tmp, s.pref.BigDecimalValue(new MathContext(Utils.scale, Utils.scaleMode2)));
}
return tmp.toBigInteger();
} /* BigSurd.toRational */
/**
* Convert to a BigDecimal value if possible
*
* @since 2012-02-15
*/
public BigDecimal toBigDecimal() {
BigDecimal tmp = BigDecimal.ZERO.setScale(Utils.scale, Utils.scaleMode);
for (BigSurd s : terms) {
tmp = BigDecimalMath.addRound(tmp, s.BigDecimalValue(new MathContext(Utils.scale, Utils.scaleMode2)));
}
return tmp;
} /* BigSurd.toBigDecimal */
/**
* Return a string in the format (number/denom)*()^(1/2). If the
* discriminant equals 1, print just the prefactor.
*
* @return the human-readable version in base 10
* @since 2012-02-16
*/
public String toString() {
/*
* simple cases with one term forwarded to the BigSurd class
*/
if (terms.size() == 0)
return new String("0");
else {
String s = new String();
for (int t = 0; t < terms.size(); t++) {
BigSurd bs = terms.elementAt(t);
if (bs.signum() > 0)
s += "+";
s += bs.toString();
}
return s;
}
} /* toString */
public String toFancyString() {
if (terms.size() == 0)
return new String("0");
else {
BigInteger denominator = BigInteger.ONE;
for (int i = 0; i < terms.size(); i++) {
denominator = denominator.multiply(terms.elementAt(i).pref.b);
}
String s = "";
if (denominator.compareTo(BigInteger.ONE) != 0) {
s += "(";
}
for (int t = 0; t < terms.size(); t++) {
BigSurd bs = terms.elementAt(t);
if (bs.signum() > 0 && t > 0)
s += "+";
if (bs.isBigInteger()) {
s += bs.BigDecimalValue(new MathContext(Utils.scale, Utils.scaleMode2)).toBigInteger().toString();
} else if (bs.isRational()) {
s += bs.toRational().toString();
} else {
BigInteger numerator = bs.pref.multiply(denominator).numer();
if (numerator.compareTo(BigInteger.ONE) != 0) {
s += numerator.toString();
s += "*";
s += "(";
}
s += "2√";
if (bs.disc.isInteger()) {
s += bs.disc.toString();
} else {
s += "("+bs.disc.toString()+")";
}
if (numerator.compareTo(BigInteger.ONE) != 0) {
s += ")";
}
}
}
if (denominator.compareTo(BigInteger.ONE) != 0) {
s += ")";
s += "/";
s += denominator;
}
return s;
}
}
} /* BigSurdVec */

69
src/org/nevec/rjm/Euler.java Normal file
View File

@ -0,0 +1,69 @@
package org.nevec.rjm ;
import java.math.BigInteger;
import java.util.Vector;
/** Euler numbers
* @see <a href="http://oeis.org/A000364">A000364</a> in the OEIS.
* @since 2008-10-30
* @author Richard J. Mathar
*/
public class Euler
{
/*
* The list of all Euler numbers as a vector, n=0,2,4,....
*/
static protected Vector<BigInteger> a = new Vector<BigInteger>() ;
/** Ctor(). Fill the hash list initially with E_0 to E_3.
*/
public Euler()
{
if ( a.size() == 0 )
{
a.add(BigInteger.ONE) ;
a.add(BigInteger.ONE) ;
a.add(new BigInteger("5")) ;
a.add(new BigInteger("61")) ;
}
}
/** Compute a coefficient in the internal table.
* @param n the zero-based index of the coefficient. n=0 for the E_0 term.
*/
protected void set(final int n)
{
while ( n >= a.size())
{
BigInteger val = BigInteger.ZERO ;
boolean sigPos = true;
int thisn = a.size() ;
for(int i= thisn-1 ; i > 0 ; i--)
{
BigInteger f = new BigInteger(""+ a.elementAt(i).toString() ) ;
f = f.multiply( BigIntegerMath.binomial(2*thisn,2*i) );
if ( sigPos )
val = val.add(f) ;
else
val = val.subtract(f) ;
sigPos = ! sigPos ;
}
if ( thisn % 2 ==0 )
val = val.subtract(BigInteger.ONE) ;
else
val = val.add(BigInteger.ONE) ;
a.add(val) ;
}
}
/** The Euler number at the index provided.
* @param n the index, non-negative.
* @return the E_0=E_1=1 , E_2=5, E_3=61 etc
*/
public BigInteger at(int n)
{
set(n) ;
return(a.elementAt(n)) ;
}
} /* Euler */

60
src/org/nevec/rjm/EulerPhi.java Normal file
View File

@ -0,0 +1,60 @@
package org.nevec.rjm ;
import java.math.BigInteger;
/** Euler totient function.
* @see <a href="http://oeis.org/A000010">A000010</a> in the OEIS.
* @since 2008-10-14
* @since 2012-03-04 Adapted to new Ifactor representation.
* @author Richard J. Mathar
*/
public class EulerPhi
{
/** Default constructor.
* Does nothing().
*/
public EulerPhi()
{
}
/** Compute phi(n).
* @param n The positive argument of the function.
* @return phi(n)
*/
public BigInteger at(int n)
{
return at(new BigInteger(""+n) ) ;
} /* at */
/** Compute phi(n).
* @param n The positive argument of the function.
* @return phi(n)
*/
public BigInteger at(BigInteger n)
{
if ( n.compareTo(BigInteger.ZERO) <= 0 )
throw new ArithmeticException("negative argument "+n+ " of EulerPhi") ;
Ifactor prFact = new Ifactor(n) ;
BigInteger phi = n ;
if ( n.compareTo(BigInteger.ONE) > 0 )
for(int i=0 ; i < prFact.primeexp.size() ; i += 2)
{
BigInteger p = new BigInteger(prFact.primeexp.elementAt(i).toString()) ;
BigInteger p_1 = p.subtract(BigInteger.ONE) ;
phi = phi.multiply(p_1).divide(p) ;
}
return phi ;
} /* at */
/** Test program.
* It takes one argument n and prints the value phi(n).<br>
* java -cp . org.nevec.rjm.EulerPhi n<br>
* @since 2006-08-14
*/
public static void main(String[] args) throws ArithmeticException
{
EulerPhi a = new EulerPhi() ;
int n = (new Integer(args[0])).intValue() ;
System.out.println("phi("+ n + ") = " + a.at(n)) ;
}
} /* EulerPhi */

70
src/org/nevec/rjm/Factorial.java Normal file
View File

@ -0,0 +1,70 @@
package org.nevec.rjm ;
import java.math.BigInteger;
import java.util.Vector;
/** Factorials.
* @since 2006-06-25
* @since 2012-02-15 Storage of the values based on Ifactor, not BigInteger.
* @author Richard J. Mathar
*/
public class Factorial
{
/** The list of all factorials as a vector.
*/
static Vector<Ifactor> a = new Vector<Ifactor>() ;
/** ctor().
* Initialize the vector of the factorials with 0!=1 and 1!=1.
*/
public Factorial()
{
if ( a.size() == 0 )
{
a.add(Ifactor.ONE) ;
a.add(Ifactor.ONE) ;
}
} /* ctor */
/** Compute the factorial of the non-negative integer.
* @param n the argument to the factorial, non-negative.
* @return the factorial of n.
*/
public BigInteger at(int n)
{
/* extend the internal list if needed.
*/
growto(n) ;
return a.elementAt(n).n ;
} /* at */
/** Compute the factorial of the non-negative integer.
* @param n the argument to the factorial, non-negative.
* @return the factorial of n.
*/
public Ifactor toIfactor(int n)
{
/* extend the internal list if needed.
*/
growto(n) ;
return a.elementAt(n) ;
} /* at */
/** Extend the internal table to cover up to n!
* @param n The maximum factorial to be supported.
* @since 2012-02-15
*/
private void growto(int n)
{
/* extend the internal list if needed. Size to be 2 for n<=1, 3 for n<=2 etc.
*/
while ( a.size() <=n )
{
final int lastn = a.size()-1 ;
final Ifactor nextn = new Ifactor(lastn+1) ;
a.add(a.elementAt(lastn).multiply(nextn) ) ;
}
} /* growto */
} /* Factorial */

39
src/org/nevec/rjm/Harmonic.java Normal file
View File

@ -0,0 +1,39 @@
package org.nevec.rjm ;
/** Harmonic numbers.
* H(n) is the sum of the inverses of the integers from 1 to n.
* @since 2008-10-19
* @author Richard J. Mathar
*/
public class Harmonic
{
/** ctor()
* Does nothing.
*/
public Harmonic()
{
}
/** The Harmonic number at the index specified
* @param n the index, non-negative.
* @return the H_1=1 for n=1, H_2=3/2 for n=2 etc.
* For values of n less than 1, zero is returned.
*/
public Rational at(int n)
{
if ( n < 1)
return(new Rational(0,1)) ;
else
{
/* start with 1 as the result
*/
Rational a = new Rational(1,1) ;
/* add 1/i for i=2..n
*/
for( int i=2 ; i <=n ; i++)
a = a.add(new Rational(1,i)) ;
return a ;
}
}
} /* Harmonic */

747
src/org/nevec/rjm/Ifactor.java Normal file
View File

@ -0,0 +1,747 @@
package org.nevec.rjm ;
import java.math.BigInteger;
import java.util.Collections;
import java.util.Vector;
/** Factored integers.
* This class contains a non-negative integer with the prime factor decomposition attached.
* @since 2006-08-14
* @since 2012-02-14 The internal representation contains the bases, and becomes sparser if few
* prime factors are present.
* @author Richard J. Mathar
*/
public class Ifactor implements Cloneable, Comparable<Ifactor>
{
/**
* The standard representation of the number
*/
public BigInteger n ;
/*
* The bases and powers of the prime factorization.
* representation n = primeexp[0]^primeexp[1]*primeexp[2]^primeexp[3]*...
* The value 0 is represented by an empty vector, the value 1 by a vector of length 1
* with a single power of 0.
*/
public Vector<Integer> primeexp ;
final public static Ifactor ONE = new Ifactor(1) ;
final public static Ifactor ZERO = new Ifactor(0) ;
/** Constructor given an integer.
* constructor with an ordinary integer
* @param number the standard representation of the integer
*/
public Ifactor(int number)
{
n = new BigInteger(""+number) ;
primeexp = new Vector<Integer>() ;
if( number > 1 )
{
int primindx = 0 ;
Prime primes = new Prime() ;
/* Test division against all primes.
*/
while(number > 1)
{
int ex=0 ;
/* primindx=0 refers to 2, =1 to 3, =2 to 5, =3 to 7 etc
*/
int p = primes.at(primindx).intValue() ;
while( number % p == 0 )
{
ex++ ;
number /= p ;
if ( number == 1 )
break ;
}
if ( ex > 0 )
{
primeexp.add(new Integer(p)) ;
primeexp.add(new Integer(ex)) ;
}
primindx++ ;
}
}
else if ( number == 1)
{
primeexp.add(new Integer(1)) ;
primeexp.add(new Integer(0)) ;
}
} /* Ifactor */
/** Constructor given a BigInteger .
* Constructor with an ordinary integer, calling a prime factor decomposition.
* @param number the BigInteger representation of the integer
*/
public Ifactor(BigInteger number)
{
n = number ;
primeexp = new Vector<Integer>() ;
if ( number.compareTo(BigInteger.ONE) == 0 )
{
primeexp.add(new Integer(1)) ;
primeexp.add(new Integer(0)) ;
}
else
{
int primindx = 0 ;
Prime primes = new Prime() ;
/* Test for division against all primes.
*/
while(number.compareTo(BigInteger.ONE) == 1)
{
int ex=0 ;
BigInteger p = primes.at(primindx) ;
while( number.remainder(p).compareTo(BigInteger.ZERO) == 0 )
{
ex++ ;
number = number.divide(p) ;
if ( number.compareTo(BigInteger.ONE) == 0 )
break ;
}
if ( ex > 0 )
{
primeexp.add(new Integer(p.intValue()) ) ;
primeexp.add(new Integer(ex) ) ;
}
primindx++ ;
}
}
} /* Ifactor */
/** Constructor given a list of exponents of the prime factor decomposition.
* @param pows the vector with the sorted list of exponents.
* pows[0] is the exponent of 2, pows[1] the exponent of 3, pows[2] the exponent of 5 etc.
* Note that this list does not include the primes, but assumes a continuous prime-smooth basis.
*/
public Ifactor(Vector<Integer> pows)
{
primeexp = new Vector<Integer>(2* pows.size()) ;
if ( pows.size() > 0 )
{
n = BigInteger.ONE ;
Prime primes = new Prime() ;
/* Build the full number by the product of all powers of the primes.
*/
for(int primindx=0 ; primindx < pows.size() ; primindx++)
{
int ex= pows.elementAt(primindx).intValue() ;
final BigInteger p = primes.at(primindx) ;
n = n.multiply( p.pow(ex) ) ;
primeexp.add(new Integer(p.intValue()) ) ;
primeexp.add(new Integer(ex) ) ;
}
}
else
n = BigInteger.ZERO ;
} /* Ifactor */
/** Copy constructor.
* @param oth the value to be copied
*/
public Ifactor(Ifactor oth)
{
n = oth.n ;
primeexp = oth.primeexp ;
} /* Ifactor */
/** Deep copy.
* @since 2009-08-14
*/
public Ifactor clone()
{
/*
* Line not used:
*
Vector<Integer> p = (Vector<Integer>)primeexp.clone();
*
*/
Ifactor cl = new Ifactor(0) ;
cl.n = new BigInteger(""+n) ;
return cl ;
} /* Ifactor.clone */
/** Comparison of two numbers.
* The value of this method is in allowing the Vector<>.contains() calls that use the value,
* not the reference for comparison.
* @param oth the number to compare this with.
* @return true if both are the same numbers, false otherwise.
*/
public boolean equals(final Ifactor oth)
{
return ( n.compareTo(oth.n) == 0 ) ;
} /* Ifactor.equals */
/** Multiply with another positive integer.
* @param oth the second factor.
* @return the product of both numbers.
*/
public Ifactor multiply(final BigInteger oth)
{
/* the optimization is to factorize oth _before_ multiplying
*/
return( multiply(new Ifactor(oth)) ) ;
} /* Ifactor.multiply */
/** Multiply with another positive integer.
* @param oth the second factor.
* @return the product of both numbers.
*/
public Ifactor multiply(final int oth)
{
/* the optimization is to factorize oth _before_ multiplying
*/
return( multiply(new Ifactor(oth)) ) ;
} /* Ifactor.multiply */
/** Multiply with another positive integer.
* @param oth the second factor.
* @return the product of both numbers.
*/
public Ifactor multiply(final Ifactor oth)
{
/* This might be done similar to the lcm() implementation by adding
* the powers of the components and calling the constructor with the
* list of exponents. This here is the simplest implementation, but slow because
* it calls another prime factorization of the product:
* return( new Ifactor(n.multiply(oth.n))) ;
*/
return multGcdLcm(oth,0) ;
}
/** Lowest common multiple of this with oth.
* @param oth the second parameter of lcm(this,oth)
* @return the lowest common multiple of both numbers. Returns zero
* if any of both arguments is zero.
*/
public Ifactor lcm(final Ifactor oth)
{
return multGcdLcm(oth,2) ;
}
/** Greatest common divisor of this and oth.
* @param oth the second parameter of gcd(this,oth)
* @return the lowest common multiple of both numbers. Returns zero
* if any of both arguments is zero.
*/
public Ifactor gcd(final Ifactor oth)
{
return multGcdLcm(oth,1) ;
}
/** Multiply with another positive integer.
* @param oth the second factor.
* @param type 0 to multiply, 1 for gcd, 2 for lcm
* @return the product, gcd or lcm of both numbers.
*/
protected Ifactor multGcdLcm(final Ifactor oth, int type)
{
Ifactor prod = new Ifactor(0) ;
/* skip the case where 0*something =0, falling thru to the empty representation for 0
*/
if( primeexp.size() != 0 && oth.primeexp.size() != 0)
{
/* Cases of 1 times something return something.
* Cases of lcm(1, something) return something.
* Cases of gcd(1, something) return 1.
*/
if ( primeexp.firstElement().intValue() == 1 && type == 0)
return oth ;
else if ( primeexp.firstElement().intValue() == 1 && type == 2)
return oth ;
else if ( primeexp.firstElement().intValue() == 1 && type == 1)
return this ;
else if ( oth.primeexp.firstElement().intValue() == 1 && type ==0)
return this ;
else if ( oth.primeexp.firstElement().intValue() == 1 && type ==2)
return this ;
else if ( oth.primeexp.firstElement().intValue() == 1 && type ==1)
return oth ;
else
{
int idxThis = 0 ;
int idxOth = 0 ;
switch(type)
{
case 0 :
prod.n = n.multiply(oth.n) ;
break;
case 1 :
prod.n = n.gcd(oth.n) ;
break;
case 2 :
/* the awkward way, lcm = product divided by gcd
*/
prod.n = n.multiply(oth.n).divide( n.gcd(oth.n) ) ;
break;
}
/* scan both representations left to right, increasing prime powers
*/
while( idxOth < oth.primeexp.size() || idxThis < primeexp.size() )
{
if ( idxOth >= oth.primeexp.size() )
{
/* exhausted the list in oth.primeexp; copy over the remaining 'this'
* if multiplying or lcm, discard if gcd.
*/
if ( type == 0 || type == 2)
{
prod.primeexp.add( primeexp.elementAt(idxThis) ) ;
prod.primeexp.add( primeexp.elementAt(idxThis+1) ) ;
}
idxThis += 2 ;
}
else if ( idxThis >= primeexp.size() )
{
/* exhausted the list in primeexp; copy over the remaining 'oth'
*/
if ( type == 0 || type == 2)
{
prod.primeexp.add( oth.primeexp.elementAt(idxOth) ) ;
prod.primeexp.add( oth.primeexp.elementAt(idxOth+1) ) ;
}
idxOth += 2 ;
}
else
{
Integer p ;
int ex ;
switch ( primeexp.elementAt(idxThis).compareTo(oth.primeexp.elementAt(idxOth) ) )
{
case 0 :
/* same prime bases p in both factors */
p = primeexp.elementAt(idxThis) ;
switch(type)
{
case 0 :
/* product means adding exponents */
ex = primeexp.elementAt(idxThis+1).intValue() +
oth.primeexp.elementAt(idxOth+1).intValue() ;
break;
case 1 :
/* gcd means minimum of exponents */
ex = Math.min( primeexp.elementAt(idxThis+1).intValue() ,
oth.primeexp.elementAt(idxOth+1).intValue()) ;
break;
default :
/* lcm means maximum of exponents */
ex = Math.max( primeexp.elementAt(idxThis+1).intValue() ,
oth.primeexp.elementAt(idxOth+1).intValue()) ;
break;
}
prod.primeexp.add( p ) ;
prod.primeexp.add( new Integer(ex) ) ;
idxOth += 2 ;
idxThis += 2 ;
break ;
case 1:
/* this prime base bigger than the other and taken later */
if ( type == 0 || type == 2)
{
prod.primeexp.add( oth.primeexp.elementAt(idxOth) ) ;
prod.primeexp.add( oth.primeexp.elementAt(idxOth+1) ) ;
}
idxOth += 2 ;
break ;
default:
/* this prime base smaller than the other and taken now */
if ( type == 0 || type == 2)
{
prod.primeexp.add( primeexp.elementAt(idxThis) ) ;
prod.primeexp.add( primeexp.elementAt(idxThis+1) ) ;
}
idxThis += 2 ;
}
}
}
}
}
return prod ;
} /* Ifactor.multGcdLcm */
/** Integer division through another positive integer.
* @param oth the denominator.
* @return the division of this through the oth, discarding the remainder.
*/
public Ifactor divide(final Ifactor oth)
{
/* todo: it'd probably be faster to cancel the gcd(this,oth) first in the prime power
* representation, which would avoid a more strenous factorization of the integer ratio
*/
return new Ifactor(n.divide(oth.n)) ;
} /* Ifactor.divide */
/** Summation with another positive integer
* @param oth the other term.
* @return the sum of both numbers
*/
public Ifactor add(final BigInteger oth)
{
/* avoid refactorization if oth is zero...
*/
if ( oth.compareTo(BigInteger.ZERO) != 0 )
return new Ifactor(n.add(oth)) ;
else
return this ;
} /* Ifactor.add */
/** Exponentiation with a positive integer.
* @param exponent the non-negative exponent
* @return n^exponent. If exponent=0, the result is 1.
*/
public Ifactor pow(final int exponent) throws ArithmeticException
{
/* three simple cases first
*/
if ( exponent < 0 )
throw new ArithmeticException("Cannot raise "+ toString() + " to negative " + exponent) ;
else if ( exponent == 0)
return new Ifactor(1) ;
else if ( exponent == 1)
return this ;
/* general case, the vector with the prime factor powers, which are component-wise
* exponentiation of the individual prime factor powers.
*/
Ifactor pows = new Ifactor(0) ;
for(int i=0 ; i < primeexp.size() ; i += 2)
{
Integer p = primeexp.elementAt(i) ;
int ex = primeexp.elementAt(i+1).intValue() ;
pows.primeexp.add( p ) ;
pows.primeexp.add( new Integer(ex*exponent) ) ;
}
return pows ;
} /* Ifactor.pow */
/** Pulling the r-th root.
* @param r the positive or negative (nonzero) root.
* @return n^(1/r).
* The return value falls into the Ifactor class if r is positive, but if r is negative
* a Rational type is needed.
* @since 2009-05-18
*/
public Rational root(final int r) throws ArithmeticException
{
if ( r == 0 )
throw new ArithmeticException("Cannot pull zeroth root of "+ toString()) ;
else if ( r < 0 )
{
/* a^(-1/b)= 1/(a^(1/b))
*/
final Rational invRoot = root(-r) ;
return Rational.ONE.divide(invRoot) ;
}
else
{
BigInteger pows = BigInteger.ONE ;
for(int i=0 ; i < primeexp.size() ; i += 2)
{
/* all exponents must be multiples of r to succeed (that is, to
* stay in the range of rational results).
*/
int ex = primeexp.elementAt(i+1).intValue() ;
if ( ex % r != 0 )
throw new ArithmeticException("Cannot pull "+ r+"th root of "+ toString()) ;
pows.multiply( new BigInteger(""+primeexp.elementAt(i)).pow(ex/r) ) ;
}
/* convert result to a Rational; unfortunately this will loose the prime factorization */
return new Rational(pows) ;
}
} /* Ifactor.root */
/** The set of positive divisors.
* @return the vector of divisors of the absolute value, sorted.
* @since 2010-08-27
*/
public Vector<BigInteger> divisors()
{
/* Recursive approach: the divisors of p1^e1*p2^e2*..*py^ey*pz^ez are
* the divisors that don't contain the factor pz, and the
* the divisors that contain any power of pz between 1 and up to ez multiplied
* by 1 or by a product that contains the factors p1..py.
*/
Vector<BigInteger> d=new Vector<BigInteger>() ;
if ( n.compareTo(BigInteger.ZERO) == 0 )
return d ;
d.add(BigInteger.ONE) ;
if ( n.compareTo(BigInteger.ONE) > 0 )
{
/* Computes sigmaIncopml(p1^e*p2^e2...*py^ey) */
Ifactor dp = dropPrime() ;
/* get ez */
final int ez = primeexp.lastElement().intValue() ;
Vector<BigInteger> partd = dp.divisors() ;
/* obtain pz by lookup in the prime list */
final BigInteger pz = new BigInteger( primeexp.elementAt(primeexp.size()-2).toString()) ;
/* the output contains all products of the form partd[]*pz^ez, ez>0,
* and with the exception of the 1, all these are appended.
*/
for(int i =1 ; i < partd.size() ; i++)
d.add( partd.elementAt(i) ) ;
for(int e =1 ; e <= ez ; e++)
{
final BigInteger pzez = pz.pow(e) ;
for(int i =0 ; i < partd.size() ; i++)
d.add( partd.elementAt(i).multiply(pzez) ) ;
}
}
Collections.sort(d) ;
return d ;
} /* Ifactor.divisors */
/** Sum of the divisors of the number.
* @return the sum of all divisors of the number, 1+....+n.
*/
public Ifactor sigma()
{
return sigma(1) ;
} /* Ifactor.sigma */
/** Sum of the k-th powers of divisors of the number.
* @return the sum of all divisors of the number, 1^k+....+n^k.
*/
public Ifactor sigma(int k)
{
/* the question is whether keeping a factorization is worth the effort
* or whether one should simply multiply these to return a BigInteger...
*/
if( n.compareTo(BigInteger.ONE) == 0 )
return ONE ;
else if( n.compareTo(BigInteger.ZERO) == 0 )
return ZERO ;
else
{
/* multiplicative: sigma_k(p^e) = [p^(k*(e+1))-1]/[p^k-1]
* sigma_0(p^e) = e+1.
*/
Ifactor resul = Ifactor.ONE ;
for(int i=0 ; i < primeexp.size() ; i += 2)
{
int ex = primeexp.elementAt(i+1).intValue() ;
if ( k == 0 )
resul = resul.multiply(ex+1) ;
else
{
Integer p = primeexp.elementAt(i) ;
BigInteger num = (new BigInteger(p.toString())).pow(k*(ex+1)).subtract(BigInteger.ONE) ;
BigInteger deno = (new BigInteger(p.toString())).pow(k).subtract(BigInteger.ONE) ;
/* This division is of course exact, no remainder
* The costly prime factorization is hidden here.
*/
Ifactor f = new Ifactor(num.divide(deno)) ;
resul = resul.multiply(f) ;
}
}
return resul ;
}
} /* Ifactor.sigma */
/** Divide through the highest possible power of the highest prime.
* If the current number is the prime factor product p1^e1 * p2*e2* p3^e3*...*py^ey * pz^ez,
* the value returned has the final factor pz^ez eliminated, which gives
* p1^e1 * p2*e2* p3^e3*...*py^ey.
* @return the new integer obtained by removing the highest prime power.
* If this here represents 0 or 1, it is returned without change.
* @since 2006-08-20
*/
public Ifactor dropPrime()
{
/* the cases n==1 or n ==0
*/
if ( n.compareTo(BigInteger.ONE) <= 0 )
return this ;
/* The cases n>1
* Start empty. Copy all but the last factor over to the result
* the vector with the new prime factor powers, which contain the
* old prime factor powers up to but not including the last one.
*/
Ifactor pows=new Ifactor(0) ;
pows.n = BigInteger.ONE ;
for(int i = 0 ; i < primeexp.size()-2 ; i += 2)
{
pows.primeexp.add( primeexp.elementAt(i)) ;
pows.primeexp.add( primeexp.elementAt(i+1)) ;
BigInteger p = new BigInteger( primeexp.elementAt(i).toString() ) ;
int ex = primeexp.elementAt(i+1).intValue() ;
pows.n = pows.n.multiply( p.pow(ex) ) ;
}
return pows ;
} /* Ifactor.dropPrime */
/** Test whether this is a square of an integer (perfect square).
* @return true if this is an integer squared (including 0), else false
*/
public boolean issquare()
{
/* check the exponents, located at the odd-indexed positions
*/
for(int i=1 ; i < primeexp.size() ; i += 2)
{
if ( primeexp.elementAt(i).intValue() % 2 != 0)
return false ;
}
return true ;
} /* Ifactor.issquare */
/** The sum of the prime factor exponents, with multiplicity.
* @return the sum over the primeexp numbers
*/
public int bigomega()
{
int resul= 0 ;
for(int i=1 ; i < primeexp.size() ; i += 2)
resul += primeexp.elementAt(i).intValue() ;
return(resul) ;
} /* Ifactor.bigomega */
/** The sum of the prime factor exponents, without multiplicity.
* @return the number of distinct prime factors.
* @since 2008-10-16
*/
public int omega()
{
return primeexp.size()/2 ;
} /* Ifactor.omega */
/** The square-free part.
* @return the minimum m such that m times this number is a square.
* @since 2008-10-16
*/
public BigInteger core()
{
BigInteger resul = BigInteger.ONE ;
for(int i=0 ; i < primeexp.size() ; i += 2)
if ( primeexp.elementAt(i+1).intValue() % 2 != 0)
resul = resul.multiply( new BigInteger(primeexp.elementAt(i).toString()) );
return resul ;
} /* Ifactor.core */
/** The Moebius function.
* 1 if n=1, else, if k is the number of distinct prime factors, return (-1)^k,
* else, if k has repeated prime factors, return 0.
* @return the moebius function.
*/
public int moebius()
{
if( n.compareTo(BigInteger.ONE) <= 0 )
return 1 ;
/* accumulate number of different primes in k */
int k=1 ;
for(int i=0 ; i < primeexp.size() ; i += 2)
{
final int e = primeexp.elementAt(i+1).intValue() ;
if ( e > 1 )
return 0 ;
else if ( e == 1)
/* accumulates (-1)^k */
k *= -1 ;
}
return( k ) ;
} /* Ifactor.moebius */
/** Maximum of two values.
* @param oth the number to compare this with.
* @return the larger of the two values.
*/
public Ifactor max(final Ifactor oth)
{
if( n.compareTo(oth.n) >= 0 )
return this ;
else
return oth ;
} /* Ifactor.max */
/** Minimum of two values.
* @param oth the number to compare this with.
* @return the smaller of the two values.
*/
public Ifactor min(final Ifactor oth)
{
if( n.compareTo(oth.n) <= 0 )
return this ;
else
return oth ;
} /* Ifactor.min */
/** Maximum of a list of values.
* @param set list of numbers.
* @return the largest in the list.
*/
public static Ifactor max(final Vector<Ifactor> set)
{
Ifactor resul = set.elementAt(0) ;
for(int i=1; i < set.size() ; i++)
resul = resul.max(set.elementAt(i)) ;
return resul ;
} /* Ifactor.max */
/** Minimum of a list of values.
* @param set list of numbers.
* @return the smallest in the list.
*/
public static Ifactor min(final Vector<Ifactor> set)
{
Ifactor resul = set.elementAt(0) ;
for(int i=1; i < set.size() ; i++)
resul = resul.min(set.elementAt(i)) ;
return resul ;
} /* Ifactor.min */
/** Compare value against another Ifactor
* @param oth The value to be compared agains.
* @return 1, 0 or -1 according to being larger, equal to or smaller than oth.
* @since 2012-02-15
*/
public int compareTo( final Ifactor oth)
{
return n.compareTo(oth.n) ;
} /* compareTo */
/** Convert to printable format
* @return a string of the form n:prime^pow*prime^pow*prime^pow...
*/
public String toString()
{
String resul = new String(n.toString()+":") ;
if ( n.compareTo(BigInteger.ONE) == 0 )
resul += "1" ;
else
{
boolean firstMul = true ;
for(int i=0 ; i < primeexp.size() ; i += 2)
{
if ( ! firstMul)
resul += "*" ;
if ( primeexp.elementAt(i+1).intValue() > 1 )
resul += primeexp.elementAt(i).toString()+"^"+primeexp.elementAt(i+1).toString() ;
else
resul += primeexp.elementAt(i).toString() ;
firstMul = false ;
}
}
return resul ;
} /* Ifactor.toString */
/** Test program.
* It takes a single argument n and prints the integer factorizaton.<br>
* java -cp . org.nevec.rjm.Ifactor n<br>
*/
public static void main(String[] args) throws Exception
{
BigInteger n = new BigInteger(args[0]) ;
System.out.println( new Ifactor(n)) ;
} /* Ifactor.main */
} /* Ifactor */

85
src/org/nevec/rjm/PartitionsP.java Normal file
View File

@ -0,0 +1,85 @@
package org.nevec.rjm ;
import java.math.BigInteger;
import java.util.Vector;
/** Number of partitions.
* @since 2008-10-15
* @author Richard J. Mathar
*/
public class PartitionsP
{
/**
* The list of all partitions as a vector.
*/
static protected Vector<BigInteger> a = new Vector<BigInteger>() ;
/**
* The maximum integer covered by the high end of the list.
*/
static protected BigInteger nMax =new BigInteger("-1") ;
/**
* Default constructor initializing a list of partitions up to 7.
*/
public PartitionsP()
{
if ( a.size() == 0 )
{
a.add(new BigInteger(""+1)) ;
a.add(new BigInteger(""+1)) ;
a.add(new BigInteger(""+2)) ;
a.add(new BigInteger(""+3)) ;
a.add(new BigInteger(""+5)) ;
a.add(new BigInteger(""+7)) ;
}
nMax = new BigInteger(""+(a.size()-1)) ;
} /* ctor */
/** return the number of partitions of i
* @param i the zero-based index into the list of partitions
* @return the ith partition number. This is 1 if i=0 or 1, 2 if i=2 and so forth.
*/
public BigInteger at(int i)
{
/* If the current list is too small, increase in intervals
* of 3 until the list has at least i elements.
*/
while ( i > nMax.intValue() )
{
growto(nMax.add(new BigInteger(""+3))) ;
}
return ( a.elementAt(i) ) ;
} /* at */
/** extend the list of known partitions up to n
* @param n the maximum integer hashed after the call.
*/
private void growto(BigInteger n)
{
while( a.size() <= n.intValue() )
{
BigInteger per = new BigInteger("0") ;
BigInteger cursiz = new BigInteger(""+a.size()) ;
for(int k=0; k < a.size() ; k++)
{
BigInteger tmp = a.elementAt(k).multiply(BigIntegerMath.sigma(a.size()-k)) ;
per = per.add(tmp) ;
}
a.add(per.divide(cursiz)) ;
}
nMax = new BigInteger(""+(a.size()-1)) ;
} /* growto */
/** Test program.
* It takes one integer argument n and prints P(n).<br>
* java -cp . org.nevec.rjm.PartitionsP n<br>
* @since 2008-10-15
*/
public static void main(String[] args) throws Exception
{
PartitionsP a = new PartitionsP() ;
int n = (new Integer(args[0])).intValue() ;
System.out.println("P("+ n +")=" + a.at(n)) ;
}
}

279
src/org/nevec/rjm/Prime.java Normal file
View File

@ -0,0 +1,279 @@
package org.nevec.rjm ;
import java.math.BigInteger;
import java.util.Vector;
/** Prime numbers.
* The implementation is a very basic computation of the set of all primes
* on demand, growing infinitely without any defined upper limit.
* The effects of such scheme are (i) the lookup-times become shorter after
* a while as more and more primes have been used and stored. The applications
* appear to become faster. (ii) Using the implementation for factorizations
* may easily require all available memory and stall finally, because indeed
* a dense list of primes with growing upper bound is kept without any hashing or lagging scheme.
* @since 2006-08-11
* @author Richard J. Mathar
*/
public class Prime
{
/** The list of all numbers as a vector.
*/
static Vector<BigInteger> a = new Vector<BigInteger>();
/** The maximum integer covered by the high end of the list.
*/
static protected BigInteger nMax = new BigInteger("-1");
/** Default constructor initializing a list of primes up to 17.
* 17 is enough to call the Miller-Rabin tests on the first 7 primes without further
* action.
*/
public Prime()
{
if ( a.size() == 0 )
{
a.add(new BigInteger(""+2)) ;
a.add(new BigInteger(""+3)) ;
a.add(new BigInteger(""+5)) ;
a.add(new BigInteger(""+7)) ;
a.add(new BigInteger(""+11)) ;
a.add(new BigInteger(""+13)) ;
a.add(new BigInteger(""+17)) ;
}
nMax = a.lastElement() ;
}
/** Test if a number is a prime.
* @param n the integer to be tested for primality
* @return true if prime, false if not
*/
public boolean contains(BigInteger n)
{
/* not documented
* return ( n.isProbablePrime() ) ;
*/
switch ( millerRabin(n) )
{
case -1:
return false ;
case 1:
return true ;
}
growto(n) ;
return( a.contains(n) ) ;
}
/** Test whether a number n is a strong pseudoprime to base a.
* @param n the integer to be tested for primality
* @param a the base
* @return true if the test is passed, so n may be a prime.
* false if the test is not passed, so n is not a prime.
* @since 2010-02-25
*/
public boolean isSPP(final BigInteger n, final BigInteger a)
{
final BigInteger two = new BigInteger(""+2) ;
/* numbers less than 2 are not prime
*/
if ( n.compareTo(two) == -1 )
return false ;
/* 2 is prime
*/
else if ( n.compareTo(two) == 0 )
return true ;
/* even numbers >2 are not prime
*/
else if ( n.remainder(two).compareTo(BigInteger.ZERO) == 0 )
return false ;
else
{
/* q= n- 1 = d *2^s with d odd
*/
final BigInteger q = n.subtract(BigInteger.ONE) ;
int s = q.getLowestSetBit() ;
BigInteger d = q.shiftRight(s) ;
/* test whether a^d = 1 (mod n)
*/
if ( a.modPow(d,n).compareTo(BigInteger.ONE) == 0 )
return true ;
/* test whether a^(d*2^r) = -1 (mod n), 0<=r<s
*/
for(int r=0; r < s ; r++)
{
if ( a.modPow(d.shiftLeft(r),n).compareTo(q) == 0 )
return true ;
}
return false ;
}
}
/** Miller-Rabin primality tests.
* @param n The prime candidate
* @return -1 if n is a composite, 1 if it is a prime, 0 if it may be a prime.
* @since 2010-02-25
*/
public int millerRabin(final BigInteger n)
{
/* list of limiting numbers which fail tests on k primes, A014233 in the OEIS
*/
final String[] mr ={"2047", "1373653", "25326001", "3215031751", "2152302898747", "3474749660383",
"341550071728321"} ;
int mrLim = 0 ;
while( mrLim < mr.length )
{
int l = n.compareTo(new BigInteger(mr[mrLim])) ;
if ( l < 0 )
break;
/* if one of the pseudo-primes: this is a composite
*/
else if ( l == 0 )
return -1 ;
mrLim++ ;
}
/* cannot test candidates larger than the last in the mr list
*/
if ( mrLim == mr.length)
return 0;
/* test the bases prime(1), prime(2) up to prime(mrLim+1)
*/
for(int p =0 ; p <= mrLim ; p++)
if ( isSPP(n, at(p)) == false )
return -1;
return 1;
}
/** return the ithprime
* @param i the zero-based index into the list of primes
* @return the ith prime. This is 2 if i=0, 3 if i=1 and so forth.
*/
public BigInteger at(int i)
{
/* If the current list is too small, increase in intervals
* of 5 until the list has at least i elements.
*/
while ( i >= a.size() )
{
growto(nMax.add(new BigInteger(""+5))) ;
}
return ( a.elementAt(i) ) ;
}
/** return the count of primes <= n
* @param n the upper limit of the scan
* @return the ith prime. This is 2 if i=0, 3 if i=1 and so forth.
*/
public BigInteger pi(BigInteger n)
{
/* If the current list is too small, increase in intervals
* of 5 until the list has at least i elements.
*/
growto(n) ;
BigInteger r = new BigInteger("0") ;
for(int i=0 ; i<a.size() ; i++)
if ( a.elementAt(i).compareTo(n) <= 0 )
r = r.add(BigInteger.ONE) ;
return r ;
}
/** return the smallest prime larger than n
* @param n lower limit of the search
* @return the next larger prime.
* @since 2008-10-16
*/
public BigInteger nextprime(BigInteger n)
{
/* if n <=1, return 2 */
if ( n.compareTo(BigInteger.ONE) <= 0)
return ( a.elementAt(0) ) ;
/* If the currently largest element in the list is too small, increase in intervals
* of 5 until the list has at least i elements.
*/
while ( a.lastElement().compareTo(n) <= 0)
{
growto(nMax.add(new BigInteger(""+5))) ;
}
for(int i=0 ; i < a.size() ; i++)
if ( a.elementAt(i).compareTo(n) == 1)
return ( a.elementAt(i) ) ;
return ( a.lastElement() ) ;
}
/** return the largest prime smaller than n
* @param n upper limit of the search
* @return the next smaller prime.
* @since 2008-10-17
*/
public BigInteger prevprime(BigInteger n)
{
/* if n <=2, return 0 */
if ( n.compareTo(BigInteger.ONE) <= 0)
return BigInteger.ZERO ;
/* If the currently largest element in the list is too small, increase in intervals
* of 5 until the list has at least i elements.
*/
while ( a.lastElement().compareTo(n) < 0)
growto(nMax.add(new BigInteger(""+5))) ;
for(int i=0 ; i < a.size() ; i++)
if ( a.elementAt(i).compareTo(n) >= 0)
return ( a.elementAt(i-1) ) ;
return ( a.lastElement() ) ;
}
/** extend the list of known primes up to n
* @param n the maximum integer known to be prime or not prime after the call.
*/
protected void growto(BigInteger n)
{
while( nMax.compareTo(n) == -1)
{
nMax = nMax.add(BigInteger.ONE) ;
boolean isp = true ;
for(int p=0; p < a.size() ; p++)
{
/*
* Test the list of known primes only up to sqrt(n)
*/
if ( a.get(p).multiply(a.get(p)).compareTo(nMax) == 1 )
break ;
/*
* The next case means that the p'th number in the list of known primes divides
* nMax and nMax cannot be a prime.
*/
if ( nMax.remainder(a.get(p)).compareTo(BigInteger.ZERO) == 0 )
{
isp = false ;
break ;
}
}
if( isp )
a.add(nMax) ;
}
}
/** Test program.
* Usage: java -cp . org.nevec.rjm.Prime n<br>
* This takes a single argument (n) and prints prime(n), the previous and next prime, and pi(n).
* @since 2006-08-14
*/
public static void main(String[] args) throws Exception
{
Prime a = new Prime() ;
int n = (new Integer(args[0])).intValue() ;
if ( n >= 1 )
{
if ( n >= 2)
System.out.println("prime("+(n-1)+") = " + a.at(n-1)) ;
System.out.println("prime("+n+") = " + a.at(n)) ;
System.out.println("prime("+(n+1)+") = " + a.at(n+1)) ;
System.out.println("pi(" + n +") = " + a.pi(new BigInteger(""+n) )) ;
}
}
} /* Prime */

902
src/org/nevec/rjm/RatPoly.java Normal file
View File

@ -0,0 +1,902 @@
package org.nevec.rjm ;
import java.math.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.Random;
import java.util.Scanner;
import java.util.Vector;
/** A one-parameter polynomial with rational coefficients.
* Alternatively to be interpreted as a sequence which has the polynomial as an (approximate)
* generating function.
* @since 2006-06-25
* @author Richard J. Mathar
*/
class RatPoly
{
/** The list of all coefficients, ascending exponents. Starting with a0, then a1, representing
* a value a0+a1*x+a2*x^2+a3*x^3+...
*/
protected Vector<Rational> a ;
/** Default ctor.
* Initializes the zero-valued polynomial x=0.
*/
public RatPoly()
{
a = new Vector<Rational>() ;
} /* ctor */
/** Constructor with an explicit list of coefficients.
* @param L the coefficients a0, a1, a2, a3,.., A deep copy of the these is created.
*/
public RatPoly(final Vector<Rational> L)
{
a = new Vector<Rational>() ;
for(int i=0 ; i < L.size() ; i++)
a.add( L.elementAt(i).clone() ) ;
simplify() ;
} /* ctor */
/** Constructor with a comma-separated list as the list of coefficients.
* @param L the string of the form a0,a1,a2,a3 with the coefficients
*/
public RatPoly(final String L) throws NumberFormatException
{
a = new Vector<Rational>() ;
Scanner sc = new Scanner(L) ;
sc.useDelimiter(",") ;
while ( sc.hasNext())
{
String tok =sc.next() ;
a.add(new Rational(tok)) ;
}
simplify() ;
sc.close();
} /* ctor */
/** Constructor from a hypergeometric series.
* @param A the list of values in the numerator of AFB
* @param B the list of values in the denominator of AFB
* @param nmax the order of the truncated polynomial representation
* @since 2008-11-13
*/
public RatPoly(final Vector<BigInteger> A, final Vector<BigInteger> B, int nmax)
{
/* To allow common initialization with the signature below,
* the main body is assembled in a separate function.
*/
init(A,B,nmax) ;
}
/** Constructor from a hypergeometric series.
* @param A the list of values in the numerator of AFB.
* At least one of these values must be a negative integer, which implicitly determines
* the order of the new polynomial.
* @param B the list of values in the denominator of AFB
* @since 2009-08-05
*/
public RatPoly(final Vector<BigInteger> A, final Vector<BigInteger> B)
{
BigInteger Nmax = BigInteger.ONE.negate() ;
for(int j=0; j < A.size() ; j++)
{
if ( A.elementAt(j).compareTo(BigInteger.ZERO) <= 0)
{
if ( Nmax.compareTo(BigInteger.ZERO) < 0 )
Nmax = A.elementAt(j).negate() ;
else
Nmax = Nmax.min( A.elementAt(j).negate() ) ;
}
}
if ( Nmax.compareTo(BigInteger.ZERO) < 0 )
throw new ArithmeticException("Infinite Number of Terms in Series "+Nmax.toString()) ;
int nmax = Nmax.intValue()-1 ;
init(A,B,nmax) ;
} /* ctor */
/** Constructor from a hypergeometric series.
* @param A the list of values in the numerator of AFB
* @param B the list of values in the denominator of AFB
* @param nmax the order of the truncated polynomial representation
* @since 2008-11-13
*/
protected void init(final Vector<BigInteger> A, final Vector<BigInteger> B, int nmax)
{
a = new Vector<Rational>() ;
Factorial f=new Factorial() ;
for( int n=0; n <= nmax ; n++)
{
Rational c = new Rational(1,1) ;
for(int j=0; j < A.size() ; j++)
{
Rational aEl = new Rational(A.elementAt(j)) ;
c = c.multiply(aEl.Pochhammer(n)) ;
}
for(int j=0; j < B.size() ; j++)
{
Rational bEl = new Rational(B.elementAt(j)) ;
c = c.divide(bEl.Pochhammer(n)) ;
}
c =c.divide(f.at(n)) ;
a.add(c) ;
}
simplify() ;
} /* init */
/** Create a copy of this.
* @since 2008-11-07
*/
@SuppressWarnings("unchecked")
public RatPoly clone()
{
RatPoly clo = new RatPoly() ;
clo.a = (Vector<Rational>)a.clone() ;
return clo ;
} /* clone */
/** Retrieve a polynomial coefficient.
* @param n the zero-based index of the coefficient. n=0 for the constant term.
* @return the polynomial coefficient in front of x^n.
*/
public Rational at(final int n)
{
if ( n < a.size())
return( a.elementAt(n) ) ;
else
return( new Rational(0,1) ) ;
} /* at */
/** Horner scheme to find the function value at the argument x
* @param x The argument of the polynomial
* @param mc The context determining the precision of the value returned.
* @since 2008-10-26
*/
public BigComplex valueOf( BigComplex x, MathContext mc)
{
/* result is initialized to zero */
BigComplex f = new BigComplex() ;
for(int i=degree() ; i >= 0 ; i--)
f = f.multiply(x,mc).add(a.elementAt(i).BigDecimalValue(mc)) ;
return f ;
} /* valueOf */
/** Horner scheme to find the function value at the argument x
* @param x The argument of the polynomial
* @since 2008-11-13
*/
public Rational valueOf( Rational x)
{
/* result is initialized to zero */
Rational f = new Rational(0,1) ;
for(int i=degree() ; i >= 0 ; i--)
f = f.multiply(x).add(a.elementAt(i)) ;
return f ;
} /* valueOf */
/** Horner scheme to find the function value at the argument x
* @param x The argument of the polynomial
* @since 2008-11-13
*/
public Rational valueOf( int x)
{
return valueOf(new Rational(x,1)) ;
} /* valueOf */
/** Horner scheme to evaluate the function at the argument x
* @param x The argument of the polynomial
* @since 2010-08-27
*/
public Rational valueOf( BigInteger x)
{
return valueOf(new Rational(x)) ;
} /* valueOf */
/* Set a polynomial coefficient.
* @param n the zero-based index of the coefficient. n=0 for the constant term.
* If the polynomial has not yet the degree to need this coefficient,
* the intermediate coefficients are implicitly set to zero.
* @param value the new value of the coefficient.
*/
public void set(final int n, final Rational value)
{
if ( n < a.size())
a.set(n,value) ;
else
{
/* fill intermediate powers with coefficients of zero
*/
while ( a.size() < n )
a.add(new Rational(0,1)) ;
a.add(value) ;
}
} /* set */
/** Set a polynomial coefficient.
* @param n the zero-based index of the coefficient. n=0 for the constant term.
* If the polynomial has not yet the degree to need this coefficient,
* the intermediate coefficients are implicitly set to zero.
* @param value the new value of the coefficient.
*/
public void set(final int n, final BigInteger value)
{
Rational val2 = new Rational(value,BigInteger.ONE) ;
set(n,val2) ;
} /* set */
/** Set a polynomial coefficient.
* @param n the zero-based index of the coefficient. n=0 for the constant term.
* If the polynomial has not yet the degree to need this coefficient,
* the intermediate coefficients are implicitly set to zero.
* @param value the new value of the coefficient.
*/
public void set(final int n, final int value)
{
Rational val2 = new Rational(value,1) ;
set(n,val2) ;
} /* set */
/* Set to the taylor series of exp(x) up to degree nmax.
* @param nmax the maximum polynomial degree
*/
public void setExp(final int nmax)
{
a.clear() ;
Factorial factorial=new Factorial() ;
for(int n=0; n <= nmax ; n++)
set(n, new Rational(BigInteger.ONE,factorial.at(n))) ;
} /* setExp */
/** Set to the taylor series representing 0+x.
*/
public void setx()
{
a.clear() ;
/* coefficient 0/1=0 */
a.add(new Rational(0,1)) ;
/* coefficient 1/1=1 */
a.add(new Rational(1,1)) ;
} /* setx */
/** Count of coefficients. One more than the degree of the polynomial.
* @return the number of polynomial coefficients.
*/
public int size()
{
return a.size() ;
} /* size */
/** Polynomial degree.
* @return the polynomial degree.
*/
public int degree()
{
return a.size()-1 ;
} /* degree */
/** Lower Polynomial degree.
* @return The smallest exponent n such that [x^n] of the polynomial is nonzero.
* If the polynmial is identical zero, the result is (still) 0.
* @since 2010-08-27
*/
public int ldegree()
{
for(int n=0 ; n < a.size() ; n++)
if ( a.elementAt(n).compareTo(BigInteger.ZERO) != 0 )
return n;
return 0 ;
} /* ldegree */
/** Multiply by a constant factor.
* @param val the factor
* @return the product of this with the factor.
* All coefficients of this have been multiplied individually by the factor.
*/
public RatPoly multiply(final Rational val)
{
RatPoly resul = new RatPoly() ;
if ( val.compareTo(BigInteger.ZERO) != 0 )
for(int n=0; n < a.size() ; n++)
resul.set(n,a.elementAt(n).multiply(val) ) ;
return resul ;
} /* multiply */
/** Multiply by a constant factor.
* @param val the factor
* @return the product of this with the factor.
* All coefficients of this have been multiplied individually by the factor.
* @since 2010-08-27
*/
public RatPoly multiply(final BigInteger val)
{
RatPoly resul = new RatPoly() ;
if ( val.compareTo(BigInteger.ZERO) != 0 )
for(int n=0; n < a.size() ; n++)
resul.set(n,a.elementAt(n).multiply(val) ) ;
return resul ;
} /* multiply */
/** Multiply by another polynomial
* @param val the other polynomial
* @return the product of this with the other polynomial
*/
public RatPoly multiply(final RatPoly val)
{
RatPoly resul = new RatPoly() ;
/* the degree of the result is the sum of the two degrees.
*/
final int nmax = degree()+val.degree() ;
for(int n=0; n <= nmax ; n++)
{
Rational coef = new Rational(0,1) ;
for(int nleft=0; nleft <= n ; nleft++)
{
coef = coef.add(at(nleft).multiply(val.at(n-nleft))) ;
}
resul.set(n,coef) ;
}
resul.simplify() ;
return resul ;
} /* multiply */
/** Raise to a positive power.
* @param n The non-negative exponent of the power
* @return The n-th power of this.
*/
public RatPoly pow(final int n) throws ArithmeticException
{
RatPoly resul = new RatPoly("1") ;
if ( n < 0 )
throw new ArithmeticException("negative polynomial power "+n) ;
else
{
/* this ought probably be done with some binary representation
* of the power and a smaller number of multiplications.
*/
for(int i=1 ; i <= n ; i++)
resul = resul.multiply(this) ;
resul.simplify() ;
return resul ;
}
} /* pow */
/** Raise to a rational power.
* The result is the taylor expansion of this, truncated at the first
* term that remains undetermined based on the current number of coefficients.
* @param r the exponent of the power
* @return This^r .
* @since 2009-05-18
*/
public RatPoly pow(final Rational r) throws ArithmeticException
{
/* split (a0+a1*x+a2*x^2+...)^r = a0^r*(1+a1/a0*r+a2/a0*r^2+..)^r
*/
Rational f = at(0) ;
f = f.pow(r) ;
/* scale the polynomial by division through the expansion coefficient of the absolute term
*/
RatPoly red = divide(a.elementAt(0)) ;
/* and remove the leading term (now equal to 1)
*/
red.set(0,0) ;
/* Binomial expansion of the rest. sum_{l=0..infinity} binomial(r,l)*red^l
*/
RatPoly resul = new RatPoly("1") ;
final int d = degree() ;
for(int l=1 ; l <= d ; l++)
{
final Rational b = Rational.binomial(r,l) ;
resul = resul.add( red.pow(l).multiply(b) ) ;
}
return resul.multiply(f) ;
} /* pow */
/** Add another polynomial
* @param val The other polynomial
* @return The sum of this and the other polynomial
* @since 2008-10-25
*/
public RatPoly add(final RatPoly val)
{
RatPoly resul = new RatPoly() ;
/* the degree of the result is the larger of the two degrees (before simplify() at least).
*/
final int nmax = (degree()>val.degree()) ? degree() : val.degree() ;
for(int n=0; n <= nmax ; n++)
{
Rational coef = at(n).add(val.at(n)) ;
resul.set(n,coef) ;
}
resul.simplify() ;
return resul ;
} /* add */
/** Subtract another polynomial
* @param val The other polynomial
* @return The difference between this and the other polynomial
* @since 2008-10-25
*/
public RatPoly subtract(final RatPoly val)
{
RatPoly resul = new RatPoly() ;
/* the degree of the result is the larger of the two degrees (before simplify() at least).
*/
final int nmax = (degree()>val.degree()) ? degree() : val.degree() ;
for(int n=0; n <= nmax ; n++)
{
Rational coef = at(n).subtract(val.at(n)) ;
resul.set(n,coef) ;
}
resul.simplify() ;
return resul ;
} /* subtract */
/** Divide by a constant.
* @param val the constant through which the coefficients will be divided.
* @return the Taylor expansion of this/val .
* @since 2009-05-18
*/
public RatPoly divide(final Rational val)
{
if ( val.compareTo(Rational.ZERO) != 0 )
{
RatPoly resul = new RatPoly() ;
for(int n=0; n < a.size() ; n++)
resul.set(n,a.elementAt(n).divide(val) ) ;
return resul ;
}
else
throw new ArithmeticException("Cannot divide " + toPString() +" through zero.") ;
} /* divide */
/** Divide by another polynomial.
* @param val the other polynomial
* @param nmax the maximum degree of the Taylor expansion of the result.
* @return the Taylor expansion of this/val up to degree nmax.
*/
public RatPoly divide(final RatPoly val,int nmax)
{
RatPoly num = this ;
RatPoly denom = val ;
/* divide by a common smallest power/degree
*/
while( num.at(0).compareTo(BigInteger.ZERO) == 0 && denom.at(0).compareTo(BigInteger.ZERO) == 0)
{
num.a.remove(0) ;
denom.a.remove(0) ;
if( num.size() <= 1 || denom.size() <= 1)
break ;
}
RatPoly resul = new RatPoly() ;
/* todo: If the polynomial division is exact, we could leave
* the loop earlier, indeed
*/
for(int n=0; n <= nmax ; n++)
{
Rational coef = num.at(n) ;
for(int nres=0; nres < n ; nres++)
{
coef = coef.subtract(resul.at(nres).multiply(denom.at(n-nres))) ;
}
coef = coef.divide(denom.at(0)) ;
resul.set(n,coef) ;
}
resul.simplify() ;
return(resul) ;
} /* divide */
/** Divide by another polynomial.
* @param val the other polynomial
* @return A vector with [0] containg the polynomial of degree which is the
* difference of thisdegree and the degree of val. [1] the remainder polynomial.
* This = returnvalue[0] + returnvalue[1]/val .
* @since 2012-03-01
*/
public RatPoly[] divideAndRemainder(final RatPoly val)
{
RatPoly[] ret = new RatPoly[2] ;
/* remove any high-order zeros
*/
RatPoly valSimpl = val.clone() ;
valSimpl.simplify() ;
RatPoly thisSimpl = clone() ;
thisSimpl.simplify() ;
/* catch the case with val equal to zero
*/
if ( valSimpl.degree() == 0 && valSimpl.a.firstElement().compareTo(Rational.ZERO) == 0)
throw new ArithmeticException("Division through zero polynomial") ;
/* degree of this smaller than degree of val: remainder is this
*/
if ( thisSimpl.degree() < valSimpl.degree() )
{
/* leading polynomial equals zero
*/
ret[0] = new RatPoly() ;
ret[1] = thisSimpl ;
}
else
{
/* long division. Highest degree by dividing the highest degree
* of this thru val.
*/
ret[0] = new RatPoly() ;
ret[0].set(thisSimpl.degree()-valSimpl.degree(),
thisSimpl.a.lastElement().divide(valSimpl.a.lastElement()) ) ;
/* recurrences: build this - val*(1-termresult) and feed this
* into another round of division. Have intermediate ret[0]+ret[1]/val.
*/
ret[1] = thisSimpl.subtract( ret[0].multiply( valSimpl) );
/* any remainder left ?
*/
if ( ret[1].degree() < valSimpl.degree() )
;
else
{
RatPoly rem[] = ret[1].divideAndRemainder(val) ;
ret[0] = ret[0].add(rem[0]) ;
ret[1] = rem[1] ;
}
}
return ret ;
} /* divideAndRemainder */
/** Print as a comma-separated list of coefficients.
* @return The representation a0,a1,a2,a3,...
* This is a sort of opposite of the ctor that takes a string as an argument.
* @since 2008-10-25
*/
public String toString()
{
String str = new String();
for(int n=0; n < a.size() ; n++)
{
if ( n == 0 )
str += a.elementAt(n).toString() ;
else
str += ","+a.elementAt(n).toString() ;
}
/* print at least a sole zero
*/
if (str.length() == 0)
str = "0" ;
return str ;
} /* toString */
/** Print as a polyomial in x.
* @return To representation a0+a1*x+a2*x^2+...
* This does not print the terms with coefficients equal to zero.
* @since 2008-10-26
*/
public String toPString()
{
String str = new String();
for(int n=0; n < a.size() ; n++)
{
final BigInteger num = a.elementAt(n).a ;
if ( num.compareTo(BigInteger.ZERO) != 0 )
{
str += " " ;
if ( num.compareTo(BigInteger.ZERO) > 0 )
str += "+" ;
str += a.elementAt(n).toString() ;
if ( n > 0 )
{
str += "*x" ;
if ( n > 1 )
str += "^"+n ;
}
}
}
/* print at least a sole zero
*/
if (str.length() == 0)
str = "0" ;
return str ;
} /* toPString */
/** Simplify the representation.
* Trailing values with zero coefficients (at high powers) are deleted.
* This modifies the polynomial on the stop (does not return another instance)
*/
private void simplify()
{
int n = a.size()-1 ;
if ( n >= 0)
while( a.elementAt(n).compareTo(BigInteger.ZERO) == 0 )
{
a.remove(n) ;
if( --n <0)
break ;
}
} /* simplify */
/** First derivative.
* @return The first derivative with respect to the indeterminate variable.
* @since 2008-10-26
*/
public RatPoly derive()
{
if ( a.size() <= 1)
/* derivative of the constant is just zero
*/
return new RatPoly() ;
else
{
RatPoly d = new RatPoly() ;
for(int i=1 ; i <= degree() ; i++)
{
final Rational c = a.elementAt(i).multiply(i) ;
d.set(i-1,c) ;
}
return d ;
}
} /* derive */
/** Scale coefficients such that the coefficient in front of the maximum degree is unity.
* @return The scaled polynomial
* @since 2008-10-26
*/
public RatPoly monic()
{
RatPoly m = new RatPoly() ;
final int d = degree() ;
for(int i=0 ; i <= d ; i++)
{
final Rational c = a.elementAt(i).divide(a.elementAt(d) ) ;
m.set(i,c) ;
}
return m ;
} /* monic */
/** Mobius transform.
* @param maxdeg the maximum polynomial degree of the result
* @return the sequence of coefficients is the Mobius transform of the original sequence.
* @since 2008-12-02
*/
public RatPoly mobiusT(int maxdeg)
{
/* Start with the polynomial 0
*/
RatPoly r = new RatPoly() ;
for(int i=1; i <= maxdeg; i++)
{
Rational c = new Rational() ;
for(int d=1; d <= i && d < a.size(); d++)
{
if (i % d == 0)
{
final Ifactor m = new Ifactor(i/d) ;
c = c.add( a.elementAt(d).multiply( m.moebius() ) ) ;
}
}
r.set(i,c) ;
}
r.simplify() ;
return r ;
} /* mobiusT */
/** Inverse Mobius transform.
* @param maxdeg the maximum polynomial degree of the result
* @return the sequence of coefficients is the inverse Mobius transform of the original sequence.
* @since 2008-12-02
*/
public RatPoly mobiusTInv(int maxdeg)
{
/* Start with the polynomial 0
*/
RatPoly r = new RatPoly() ;
for(int i=1; i <= maxdeg; i++)
{
Rational c = new Rational() ;
for(int d=1; d <= i && d < a.size(); d++)
{
if (i % d == 0)
c = c.add( a.elementAt(d) ) ;
}
r.set(i,c) ;
}
r.simplify() ;
return r ;
} /* mobiusTInv */
/** Binomial transform.
* @param maxdeg the maximum polynomial degree of the result
* @return the sequence of coefficients is the binomial transform of the original sequence.
* @since 2008-10-26
*/
public RatPoly binomialT(int maxdeg)
{
RatPoly r = new RatPoly() ;
for(int i=0; i <= maxdeg; i++)
{
Rational c = new Rational(0,1) ;
for(int j=0; j <= i && j < a.size(); j++)
c = c.add( a.elementAt(j).multiply(BigIntegerMath.binomial(i,j)) ) ;
r.set(i,c) ;
}
r.simplify() ;
return r ;
} /* binomialT */
/** Inverse Binomial transform.
* @param maxdeg the maximum polynomial degree of the result
* @return the sequence of coefficients is the inverse binomial transform of the original sequence.
* @since 2008-10-26
*/
public RatPoly binomialTInv(int maxdeg)
{
RatPoly r = new RatPoly() ;
for(int i=0; i <= maxdeg; i++)
{
Rational c = new Rational(0,1) ;
for(int j=0; j <= i && j < a.size(); j++)
if ( (j+i) % 2 != 0 )
c = c.subtract( a.elementAt(j).multiply(BigIntegerMath.binomial(i,j)) ) ;
else
c = c.add( a.elementAt(j).multiply(BigIntegerMath.binomial(i,j)) ) ;
r.set(i,c) ;
}
r.simplify() ;
return r ;
} /* binomialTInv */
/** Truncate polynomial degree.
* @param newdeg The maximum degree of the result.
* @return The polynomial with all coefficients beyond deg set to zero.
* If newdeg =3, for example the polynomial returned has at most degree 3.
* If newdeg is larger than the degree of this, zeros (at the higher orders of x)
* are appended. That polynomial would have formal degree larger than this.
* @since 2008-10-26
*/
public RatPoly trunc(int newdeg)
{
RatPoly t = new RatPoly() ;
for(int i=0; i <= newdeg; i++)
t.set(i,at(i)) ;
t.simplify() ;
return t ;
} /* trunc */
/** Generate the roots of the polynomial in floating point arithmetic.
* @see <a href="http://en.wikipedia.org/wiki/Durand-Kerner_method">Durand Kerner method</a>
* @param the number of floating point digits
* @since 2008-10-26
*/
public Vector<BigComplex> roots(int digits)
{
RatPoly mon = monic() ;
Random rand = new Random() ;
MathContext mc = new MathContext(digits+3,RoundingMode.DOWN) ;
Vector<BigComplex> res =new Vector<BigComplex>() ;
final int d = mon.degree() ;
double randRad =0. ;
for(int i=0 ; i <= d ; i++)
{
/* scale coefficient at maximum degree */
double absi = Math.abs( mon.at(i).doubleValue() ) ;
if ( absi > randRad)
randRad = absi ;
}
randRad += 1.0 ;
/* initial values randomly in radius 1+randRad
*/
for(int i=0 ; i < d ; i++)
{
double rad = randRad*rand.nextDouble() ;
double phi = 2.0*3.14159*rand.nextDouble() ;
res.add(i, new BigComplex(rad*Math.cos(phi),rad*Math.sin(phi)) ) ;
}
/* iterate until convr indicates that all values changed by less than the digits
* precision indicates.
*/
boolean convr = false ;
for(;! convr;)//ORIGINAL LINE: for(int itr =0 ; ! convr ; itr++)
{
convr = true ;
Vector<BigComplex> resPlus =new Vector<BigComplex>() ;
for(int v=0 ; v < d; v++)
{
/* evaluate f(x)/(x-root1)/(x-root2)/... (x-rootdegr), Newton method
*/
BigComplex thisx = res.elementAt(v) ;
BigComplex nv = mon.valueOf(thisx,mc) ;
for(int j=0; j < d ; j++)
{
if ( j != v )
nv = nv.divide(thisx.subtract(res.elementAt(j)),mc) ;
}
/* is this value converged ? */
if ( nv.abs(mc).doubleValue() > thisx.abs(mc).doubleValue()*Math.pow(10.0,-digits) )
convr =false;
thisx = thisx.subtract(nv) ;
/* If unstable, start over */
if ( thisx.abs(MathContext.DECIMAL32).doubleValue() > randRad )
return roots(digits) ;
resPlus.add(thisx) ;
}
res = resPlus ;
}
return res;
} /* roots */
/** Generate the integer roots of the polynomial.
* @return The vector of integer roots, with multiplicity.
* The shows alternatingly first a root then its multiplicty, then another root and multiplicty etc.
* @since 2008-10-26
*/
public Vector<BigInteger> iroots()
{
/* The vector of the roots */
Vector<BigInteger> res =new Vector<BigInteger>() ;
int lowd = ldegree() ;
if( lowd == 0 && a.elementAt(0).compareTo(BigInteger.ZERO) == 0)
{
/* Case of polynomial identical to zero:
* reported as a simple root of value 0.
*/
res.add(BigInteger.ZERO) ;
res.add(BigInteger.ONE) ;
return res ;
}
/* multiply all coefs with the lcm() to get an integer polynomial
* start with denominator of first non-zero coefficient.
*/
BigInteger lcmDeno = a.elementAt(lowd).b ;
for(int i=lowd+1; i < degree() ; i++)
lcmDeno = BigIntegerMath.lcm(lcmDeno, a.elementAt(i).b ) ;
/* and eventually get the integer polynomial by ignoring the denominators
*/
Vector<BigInteger> ipo = new Vector<BigInteger>() ;
for(int i=0 ; i < a.size() ; i++)
{
BigInteger d = a.elementAt(i).a.multiply( lcmDeno).divide( a.elementAt(i).b) ;
ipo.add(d) ;
}
BigIntegerPoly p = new BigIntegerPoly(ipo) ;
/* collect the integer roots (multiple roots only once). Since we
* removed the zero already above, cand does not contain zeros.
*/
Vector<BigInteger> cand = p.iroots() ;
for( int i =0 ; i < cand.size() ; i++)
{
final BigInteger r = cand.elementAt(i) ;
int deg = p.rootDeg( r) ;
res.add(r) ;
res.add(new BigInteger(""+deg)) ;
}
return res;
} /* iroots */
} /* RatPoly */

747
src/org/nevec/rjm/Rational.java Normal file
View File

@ -0,0 +1,747 @@
package org.nevec.rjm ;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;
/** Fractions (rational numbers).
* They are divisions of two BigInteger numbers, reduced to coprime
* numerator and denominator.
* @since 2006-06-25
* @author Richard J. Mathar
*/
public class Rational implements Cloneable, Comparable<Rational>
{
/** numerator
*/
BigInteger a ;
/** denominator, always larger than zero.
*/
BigInteger b ;
/** The maximum and minimum value of a standard Java integer, 2^31.
* @since 2009-05-18
*/
static public BigInteger MAX_INT = new BigInteger("2147483647") ;
static public BigInteger MIN_INT = new BigInteger("-2147483648") ;
/** The constant 1.
*/
public static Rational ONE = new Rational(1,1) ;
/** The constant 0.
*/
static public Rational ZERO = new Rational() ;
/** The constant 1/2
* @since 2010-05-25
*/
static public Rational HALF = new Rational(1,2) ;
/** Default ctor, which represents the zero.
* @since 2007-11-17
*/
public Rational()
{
a = BigInteger.ZERO ;
b = BigInteger.ONE ;
}
/** ctor from a numerator and denominator.
* @param a the numerator.
* @param b the denominator.
*/
public Rational(BigInteger a, BigInteger b)
{
this.a = a ;
this.b = b ;
normalize() ;
}
/** ctor from a numerator.
* @param a the BigInteger.
*/
public Rational(BigInteger a)
{
this.a = a ;
b = new BigInteger("1") ;
}
/** ctor from a numerator and denominator.
* @param a the numerator.
* @param b the denominator.
*/
public Rational(int a, int b)
{
this(new BigInteger(""+a),new BigInteger(""+b)) ;
}
/** ctor from an integer.
* @param n the integer to be represented by the new instance.
* @since 2010-07-18
*/
public Rational(int n)
{
this(n,1) ;
}
/** ctor from a string representation.
* @param str the string.
* This either has a slash in it, separating two integers, or, if there is no slash,
* is representing the numerator with implicit denominator equal to 1.
* Warning: this does not yet test for a denominator equal to zero
*/
public Rational(String str) throws NumberFormatException
{
this(str,10) ;
}
/** ctor from a string representation in a specified base.
* @param str the string.
* This either has a slash in it, separating two integers, or, if there is no slash,
* is just representing the numerator.
* @param radix the number base for numerator and denominator
* Warning: this does not yet test for a denominator equal to zero
*/
public Rational(String str, int radix) throws NumberFormatException
{
int hasslah = str.indexOf("/") ;
if ( hasslah == -1 )
{
a = new BigInteger(str,radix) ;
b = new BigInteger("1",radix) ;
/* no normalization necessary here */
}
else
{
/* create numerator and denominator separately
*/
a = new BigInteger(str.substring(0,hasslah),radix) ;
b = new BigInteger(str.substring(hasslah+1),radix) ;
normalize() ;
}
}
/** Create a copy.
* @since 2008-11-07
*/
public Rational clone()
{
/* protected access means this does not work
* return new Rational(a.clone(), b.clone()) ;
*/
BigInteger aclon = new BigInteger(""+a) ;
BigInteger bclon = new BigInteger(""+b) ;
return new Rational(aclon,bclon) ;
} /* Rational.clone */
/** Multiply by another fraction.
* @param val a second rational number.
* @return the product of this with the val.
*/
public Rational multiply(final Rational val)
{
BigInteger num = a.multiply(val.a) ;
BigInteger deno = b.multiply(val.b) ;
/* Normalization to an coprime format will be done inside
* the ctor() and is not duplicated here.
*/
return ( new Rational(num,deno) ) ;
} /* Rational.multiply */
/** Multiply by a BigInteger.
* @param val a second number.
* @return the product of this with the value.
*/
public Rational multiply(final BigInteger val)
{
Rational val2 = new Rational(val,BigInteger.ONE) ;
return ( multiply(val2) ) ;
} /* Rational.multiply */
/** Multiply by an integer.
* @param val a second number.
* @return the product of this with the value.
*/
public Rational multiply(final int val)
{
BigInteger tmp = new BigInteger(""+val) ;
return multiply(tmp) ;
} /* Rational.multiply */
/** Power to an integer.
* @param exponent the exponent.
* @return this value raised to the power given by the exponent.
* If the exponent is 0, the value 1 is returned.
*/
public Rational pow(int exponent)
{
if ( exponent == 0 )
return new Rational(1,1) ;
BigInteger num = a.pow(Math.abs(exponent)) ;
BigInteger deno = b.pow(Math.abs(exponent)) ;
if ( exponent > 0 )
return ( new Rational(num,deno) ) ;
else
return ( new Rational(deno,num) ) ;
} /* Rational.pow */
/** Power to an integer.
* @param exponent the exponent.
* @return this value raised to the power given by the exponent.
* If the exponent is 0, the value 1 is returned.
* @since 2009-05-18
*/
public Rational pow(BigInteger exponent) throws NumberFormatException
{
/* test for overflow */
if ( exponent.compareTo(MAX_INT) == 1 )
throw new NumberFormatException("Exponent "+exponent.toString()+" too large.") ;
if ( exponent.compareTo(MIN_INT) == -1 )
throw new NumberFormatException("Exponent "+exponent.toString()+" too small.") ;
/* promote to the simpler interface above */
return pow( exponent.intValue() ) ;
} /* Rational.pow */
/** r-th root.
* @param r the inverse of the exponent.
* 2 for the square root, 3 for the third root etc
* @return this value raised to the inverse power given by the root argument, this^(1/r).
* @since 2009-05-18
*/
public Rational root(BigInteger r) throws NumberFormatException
{
/* test for overflow */
if ( r.compareTo(MAX_INT) == 1 )
throw new NumberFormatException("Root "+r.toString()+" too large.") ;
if ( r.compareTo(MIN_INT) == -1 )
throw new NumberFormatException("Root "+r.toString()+" too small.") ;
int rthroot = r.intValue() ;
/* cannot pull root of a negative value with even-valued root */
if ( compareTo(ZERO) == -1 && (rthroot % 2) ==0 )
throw new NumberFormatException("Negative basis "+ toString()+" with odd root "+r.toString()) ;
/* extract a sign such that we calculate |n|^(1/r), still r carrying any sign
*/
final boolean flipsign = ( compareTo(ZERO) == -1 && (rthroot % 2) != 0) ? true : false ;
/* delegate the main work to ifactor#root()
*/
Ifactor num = new Ifactor(a.abs()) ;
Ifactor deno = new Ifactor(b) ;
final Rational resul = num.root(rthroot).divide( deno.root(rthroot) ) ;
if ( flipsign)
return resul.negate() ;
else
return resul ;
} /* Rational.root */
/** Raise to a rational power.
* @param exponent The exponent.
* @return This value raised to the power given by the exponent.
* If the exponent is 0, the value 1 is returned.
* @since 2009-05-18
*/
public Rational pow(Rational exponent) throws NumberFormatException
{
if ( exponent.a.compareTo(BigInteger.ZERO) == 0 )
return new Rational(1,1) ;
/* calculate (a/b)^(exponent.a/exponent.b) as ((a/b)^exponent.a)^(1/exponent.b)
* = tmp^(1/exponent.b)
*/
Rational tmp = pow(exponent.a) ;
return tmp.root(exponent.b) ;
} /* Rational.pow */
/** Divide by another fraction.
* @param val A second rational number.
* @return The value of this/val
*/
public Rational divide(final Rational val)
{
if( val.compareTo(Rational.ZERO) == 0 )
throw new ArithmeticException("Dividing "+ toString() + " through zero.") ;
BigInteger num = a.multiply(val.b) ;
BigInteger deno = b.multiply(val.a) ;
/* Reduction to a coprime format is done inside the ctor,
* and not repeated here.
*/
return ( new Rational(num,deno) ) ;
} /* Rational.divide */
/** Divide by an integer.
* @param val a second number.
* @return the value of this/val
*/
public Rational divide(BigInteger val)
{
if( val.compareTo(BigInteger.ZERO) == 0 )
throw new ArithmeticException("Dividing "+ toString() + " through zero.") ;
Rational val2 = new Rational(val,BigInteger.ONE) ;
return ( divide(val2)) ;
} /* Rational.divide */
/** Divide by an integer.
* @param val A second number.
* @return The value of this/val
*/
public Rational divide(int val)
{
if( val == 0 )
throw new ArithmeticException("Dividing "+ toString() + " through zero.") ;
Rational val2 = new Rational(val,1) ;
return ( divide(val2)) ;
} /* Rational.divide */
/** Add another fraction.
* @param val The number to be added
* @return this+val.
*/
public Rational add(Rational val)
{
BigInteger num = a.multiply(val.b).add(b.multiply(val.a)) ;
BigInteger deno = b.multiply(val.b) ;
return ( new Rational(num,deno) ) ;
} /* Rational.add */
/** Add another integer.
* @param val The number to be added
* @return this+val.
*/
public Rational add(BigInteger val)
{
Rational val2 = new Rational(val,BigInteger.ONE) ;
return ( add(val2) ) ;
} /* Rational.add */
/** Add another integer.
* @param val The number to be added
* @return this+val.
* @since May 26 2010
*/
public Rational add(int val)
{
BigInteger val2 = a.add(b.multiply(new BigInteger(""+val))) ;
return new Rational(val2,b) ;
} /* Rational.add */
/** Compute the negative.
* @return -this.
*/
public Rational negate()
{
return ( new Rational(a.negate(),b) ) ;
} /* Rational.negate */
/** Subtract another fraction.
* @param val the number to be subtracted from this
* @return this - val.
*/
public Rational subtract(Rational val)
{
Rational val2 = val.negate() ;
return ( add(val2) ) ;
} /* Rational.subtract */
/** Subtract an integer.
* @param val the number to be subtracted from this
* @return this - val.
*/
public Rational subtract(BigInteger val)
{
Rational val2 = new Rational(val,BigInteger.ONE) ;
return ( subtract(val2) ) ;
} /* Rational.subtract */
/** Subtract an integer.
* @param val the number to be subtracted from this
* @return this - val.
*/
public Rational subtract(int val)
{
Rational val2 = new Rational(val,1) ;
return ( subtract(val2) ) ;
} /* Rational.subtract */
/** binomial (n choose m).
* @param n the numerator. Equals the size of the set to choose from.
* @param m the denominator. Equals the number of elements to select.
* @return the binomial coefficient.
* @since 2006-06-27
* @author Richard J. Mathar
*/
public static Rational binomial(Rational n, BigInteger m)
{
if ( m.compareTo(BigInteger.ZERO) == 0 )
return Rational.ONE ;
Rational bin = n ;
for(BigInteger i=new BigInteger("2") ; i.compareTo(m) != 1 ; i = i.add(BigInteger.ONE) )
{
bin = bin.multiply(n.subtract(i.subtract(BigInteger.ONE))).divide(i) ;
}
return bin ;
} /* Rational.binomial */
/** binomial (n choose m).
* @param n the numerator. Equals the size of the set to choose from.
* @param m the denominator. Equals the number of elements to select.
* @return the binomial coefficient.
* @since 2009-05-19
* @author Richard J. Mathar
*/
public static Rational binomial(Rational n, int m)
{
if ( m == 0 )
return Rational.ONE ;
Rational bin = n ;
for( int i=2 ; i <= m ; i++ )
{
bin = bin.multiply(n.subtract(i-1)).divide(i) ;
}
return bin ;
} /* Rational.binomial */
/** Hankel's symbol (n,k)
* @param n the first parameter.
* @param k the second parameter, greater or equal to 0.
* @return Gamma(n+k+1/2)/k!/GAMMA(n-k+1/2)
* @since 2010-07-18
* @author Richard J. Mathar
*/
public static Rational hankelSymb(Rational n, int k)
{
if ( k == 0 )
return Rational.ONE ;
else if ( k < 0)
throw new ArithmeticException("Negative parameter "+k) ;
Rational nkhalf = n.subtract(k).add(Rational.HALF) ;
nkhalf = nkhalf.Pochhammer(2*k) ;
Factorial f = new Factorial() ;
return nkhalf.divide(f.at(k)) ;
} /* Rational.binomial */
/** Get the numerator.
* @return The numerator of the reduced fraction.
*/
public BigInteger numer()
{
return a ;
}
/** Get the denominator.
* @return The denominator of the reduced fraction.
*/
public BigInteger denom()
{
return b ;
}
/** Absolute value.
* @return The absolute (non-negative) value of this.
*/
public Rational abs()
{
return( new Rational(a.abs(),b.abs())) ;
}
/** floor(): the nearest integer not greater than this.
* @return The integer rounded towards negative infinity.
*/
public BigInteger floor()
{
/* is already integer: return the numerator
*/
if ( b.compareTo(BigInteger.ONE) == 0 )
return a;
else if ( a.compareTo(BigInteger.ZERO) > 0 )
return a.divide(b);
else
return a.divide(b).subtract(BigInteger.ONE) ;
} /* Rational.floor */
/** ceil(): the nearest integer not smaller than this.
* @return The integer rounded towards positive infinity.
* @since 2010-05-26
*/
public BigInteger ceil()
{
/* is already integer: return the numerator
*/
if ( b.compareTo(BigInteger.ONE) == 0 )
return a;
else if ( a.compareTo(BigInteger.ZERO) > 0 )
return a.divide(b).add(BigInteger.ONE) ;
else
return a.divide(b) ;
} /* Rational.ceil */
/** Remove the fractional part.
* @return The integer rounded towards zero.
*/
public BigInteger trunc()
{
/* is already integer: return the numerator
*/
if ( b.compareTo(BigInteger.ONE) == 0 )
return a;
else
return a.divide(b);
} /* Rational.trunc */
/** Compares the value of this with another constant.
* @param val the other constant to compare with
* @return -1, 0 or 1 if this number is numerically less than, equal to,
* or greater than val.
*/
public int compareTo(final Rational val)
{
/* Since we have always kept the denominators positive,
* simple cross-multiplying works without changing the sign.
*/
final BigInteger left = a.multiply(val.b) ;
final BigInteger right = val.a.multiply(b) ;
return left.compareTo(right) ;
} /* Rational.compareTo */
/** Compares the value of this with another constant.
* @param val the other constant to compare with
* @return -1, 0 or 1 if this number is numerically less than, equal to,
* or greater than val.
*/
public int compareTo(final BigInteger val)
{
final Rational val2 = new Rational(val,BigInteger.ONE) ;
return ( compareTo(val2) ) ;
} /* Rational.compareTo */
/** Return a string in the format number/denom.
* If the denominator equals 1, print just the numerator without a slash.
* @return the human-readable version in base 10
*/
public String toString()
{
if ( b.compareTo(BigInteger.ONE) != 0)
return( a.toString()+"/"+b.toString() ) ;
else
return a.toString() ;
} /* Rational.toString */
/** Return a double value representation.
* @return The value with double precision.
* @since 2008-10-26
*/
public double doubleValue()
{
/* To meet the risk of individual overflows of the exponents of
* a separate invocation a.doubleValue() or b.doubleValue(), we divide first
* in a BigDecimal environment and convert the result.
*/
BigDecimal adivb = (new BigDecimal(a)).divide(new BigDecimal(b), MathContext.DECIMAL128) ;
return adivb.doubleValue() ;
} /* Rational.doubleValue */
/** Return a float value representation.
* @return The value with single precision.
* @since 2009-08-06
*/
public float floatValue()
{
BigDecimal adivb = (new BigDecimal(a)).divide(new BigDecimal(b), MathContext.DECIMAL128) ;
return adivb.floatValue() ;
} /* Rational.floatValue */
/** Return a representation as BigDecimal.
* @param mc the mathematical context which determines precision, rounding mode etc
* @return A representation as a BigDecimal floating point number.
* @since 2008-10-26
*/
public BigDecimal BigDecimalValue(MathContext mc)
{
/* numerator and denominator individually rephrased
*/
BigDecimal n = new BigDecimal(a) ;
BigDecimal d = new BigDecimal(b) ;
/* the problem with n.divide(d,mc) is that the apparent precision might be
* smaller than what is set by mc if the value has a precise truncated representation.
* 1/4 will appear as 0.25, independent of mc
*/
return BigDecimalMath.scalePrec(n.divide(d,mc),mc) ;
} /* Rational.BigDecimalValue */
/** Return a string in floating point format.
* @param digits The precision (number of digits)
* @return The human-readable version in base 10.
* @since 2008-10-25
*/
public String toFString(int digits)
{
if ( b.compareTo(BigInteger.ONE) != 0)
{
MathContext mc = new MathContext(digits,RoundingMode.DOWN) ;
BigDecimal f = (new BigDecimal(a)).divide(new BigDecimal(b),mc) ;
return( f.toString() ) ;
}
else
return a.toString() ;
} /* Rational.toFString */
/** Compares the value of this with another constant.
* @param val The other constant to compare with
* @return The arithmetic maximum of this and val.
* @since 2008-10-19
*/
public Rational max(final Rational val)
{
if ( compareTo(val) > 0 )
return this;
else
return val;
} /* Rational.max */
/** Compares the value of this with another constant.
* @param val The other constant to compare with
* @return The arithmetic minimum of this and val.
* @since 2008-10-19
*/
public Rational min(final Rational val)
{
if ( compareTo(val) < 0 )
return this;
else
return val;
} /* Rational.min */
/** Compute Pochhammer's symbol (this)_n.
* @param n The number of product terms in the evaluation.
* @return Gamma(this+n)/Gamma(this) = this*(this+1)*...*(this+n-1).
* @since 2008-10-25
*/
public Rational Pochhammer(final BigInteger n)
{
if ( n.compareTo(BigInteger.ZERO) < 0 )
return null;
else if ( n.compareTo(BigInteger.ZERO) == 0 )
return Rational.ONE ;
else
{
/* initialize results with the current value
*/
Rational res = new Rational(a,b) ;
BigInteger i = BigInteger.ONE ;
for( ; i.compareTo(n) < 0 ; i=i.add(BigInteger.ONE) )
res = res.multiply( add(i) ) ;
return res;
}
} /* Rational.pochhammer */
/** Compute pochhammer's symbol (this)_n.
* @param n The number of product terms in the evaluation.
* @return Gamma(this+n)/GAMMA(this).
* @since 2008-11-13
*/
public Rational Pochhammer(int n)
{
return Pochhammer(new BigInteger(""+n)) ;
} /* Rational.pochhammer */
/** True if the value is integer.
* Equivalent to the indication whether a conversion to an integer
* can be exact.
* @since 2010-05-26
*/
public boolean isBigInteger()
{
return ( b.abs().compareTo(BigInteger.ONE) == 0 ) ;
} /* Rational.isBigInteger */
/** True if the value is integer and in the range of the standard integer.
* Equivalent to the indication whether a conversion to an integer
* can be exact.
* @since 2010-05-26
*/
public boolean isInteger()
{
if ( ! isBigInteger() )
return false;
return ( a.compareTo(MAX_INT) <= 0 && a.compareTo(MIN_INT) >= 0 ) ;
} /* Rational.isInteger */
/** Conversion to an integer value, if this can be done exactly.
* @since 2011-02-13
*/
int intValue()
{
if ( ! isInteger() )
throw new NumberFormatException("cannot convert "+toString()+" to integer.") ;
return a.intValue() ;
}
/** Conversion to a BigInteger value, if this can be done exactly.
* @since 2012-03-02
*/
BigInteger BigIntegerValue()
{
if ( ! isBigInteger() )
throw new NumberFormatException("cannot convert "+toString()+" to BigInteger.") ;
return a ;
}
/** True if the value is a fraction of two integers in the range of the standard integer.
* @since 2010-05-26
*/
public boolean isIntegerFrac()
{
return ( a.compareTo(MAX_INT) <= 0 && a.compareTo(MIN_INT) >= 0
&& b.compareTo(MAX_INT) <= 0 && b.compareTo(MIN_INT) >= 0 ) ;
} /* Rational.isIntegerFrac */
/** The sign: 1 if the number is >0, 0 if ==0, -1 if <0
* @return the signum of the value.
* @since 2010-05-26
*/
public int signum()
{
return ( b.signum() * a.signum() ) ;
} /* Rational.signum */
/** Common lcm of the denominators of a set of rational values.
* @param vals The list/set of the rational values.
* @return LCM(denom of first, denom of second, ..,denom of last)
* @since 2012-03-02
*/
static public BigInteger lcmDenom(final Rational[] vals)
{
BigInteger l = BigInteger.ONE ;
for(int v= 0 ; v < vals.length ; v++)
l = BigIntegerMath.lcm(l,vals[v].b) ;
return l ;
} /* Rational.lcmDenom */
/** Normalize to coprime numerator and denominator.
* Also copy a negative sign of the denominator to the numerator.
* @since 2008-10-19
*/
protected void normalize()
{
/* compute greatest common divisor of numerator and denominator
*/
final BigInteger g = a.gcd(b) ;
if ( g.compareTo(BigInteger.ONE) > 0 )
{
a = a.divide(g) ;
b = b.divide(g);
}
if ( b.compareTo(BigInteger.ZERO) == -1 )
{
a = a.negate() ;
b = b.negate() ;
}
} /* Rational.normalize */
} /* Rational */

590
src/org/nevec/rjm/Wigner3j.java Normal file
View File

@ -0,0 +1,590 @@
package org.nevec.rjm;
import java.math.BigInteger;
import java.util.Scanner;
/**
* Exact representations of Wigner 3jm and 3nj values of half-integer arguments.
*
* @see R. J. Mathar, <a href="http://arxiv.org/abs/1102.5125">Corrigendum to
* "Universal factorzation fo 3n-j (j>2) symbols ..[J. Phys. A: Math. Gen.37 (2004) 3259]"
* </a>
* @see R. J. Mathar, <a href="http://vixra.org/abs/1202.0093">Symmetries in
* Wigner 18-j and 21-j Symbols</a>
* @since 2011-02-15
* @author Richard J. Mathar
*/
public class Wigner3j {
/**
* Test programs. This supports three types of direct evaluations:<br>
* java -cp . org.nevec.rjm.Wigner3j 3jm 2j1+1 2j2+1 2j3+1 2m1+1 2m2+1 2m3+1
* <br>
* java -cp . org.nevec.rjm.Wigner3j 6j 2j1+1 2j2+2 .. 2j6+1<br>
* java -cp . org.nevec.rjm.Wigner3j 9j 2j1+1 2j2+2 .. 2j9+1<br>
* The first command line argument is one of the three tags which determine
* whether a 3jm, a 6j or a 9j symbol will be computed. The other arguments
* are 6 or 9 integer values, which are the physical (half-integer) values
* multplied by 2 and augmented by 1. The order of the 6 or 9 values is as
* reading the corresponding standard symbol as first row, then second row
* (and for the 9j symbol) third row.
*
* @since 2011-02-15
* @author Richard J. Mathar
*/
static public void main(String args[]) {
if (args[0].compareTo("6j") == 0) {
try {
String m1 = "6";
String t1 = "1 2 -3 -1 5 6";
String t2 = "4 -5 3 -4 -2 -6";
String j = "";
for (int i = 1; i <= 6; i++)
j += args[i] + " ";
BigSurdVec w = wigner3j(m1, t1, t2, j);
System.out.println(w.toString());
} catch (Exception e) {
System.out.println(e.getMessage());
}
} else if (args[0].compareTo("9j") == 0) {
try {
String m1 = "9";
String t1 = "1 3 2 4 6 5 7 9 8";
String t2 = "2 8 5 6 3 9 7 4 1";
String j = "";
for (int i = 1; i <= 9; i++)
j += args[i] + " ";
BigSurdVec w = wigner3j(m1, t1, t2, j);
System.out.println(w.toString());
} catch (Exception e) {
System.out.println(e.getMessage());
}
} else if (args[0].compareTo("3jm") == 0) {
int j1 = (new Integer(args[1])).intValue();
int j2 = (new Integer(args[2])).intValue();
int j3 = (new Integer(args[3])).intValue();
int m1 = (new Integer(args[4])).intValue();
int m2 = (new Integer(args[5])).intValue();
int m3 = (new Integer(args[6])).intValue();
try {
BigSurd w = wigner3jm(j1, j2, j3, m1, m2, m3);
System.out.println(w.toString());
w = w.multiply(new BigSurd(j3 + 1, 1));
System.out.println("CG factor sqrt" + (j3 + 1) + "sign " + ((j2 - j2 - m3) / 2) + " " + w.toString());
} catch (Exception e) {
System.out.println(e.getMessage());
}
} else {
System.out.println("usage:");
System.out.println(args[0] + " 6j 2j1+1 2j2+1 2j3+1 2j4+1 2j5+1 2j6+1");
System.out.println(args[0] + " 9j 2j1+1 2j2+1 2j3+1 2j4+1 2j5+1 2j6+1.. 2j9+1 ");
System.out.println(args[0] + " 3jm 2j1+1 2j2+1 2j3+1 2m1+1 2m2+1 2m3+1 ");
}
} /* Wigner3j.main */
/**
* The Wigner 3jm symbol (j1,j2,j3,m1,m2,m3). All arguments of the function
* are the actual parameters multiplied by 2, so they all allow an integer
* representation.
*
* @param j1
* integer representing 2*j1
* @param j2
* integer representing 2*j2
* @param j3
* integer representing 2*j3
* @param m1
* integer representing 2*m1
* @param m2
* integer representing 2*m2
* @param m3
* integer representing 2*m3
* @return The value of the symbol. Zero if any of the triangular
* inequalities is violated or some parameters are out of range.
* @since 2011-02-13
* @author Richard J. Mathar
*/
static public BigSurd wigner3jm(int j1, int j2, int j3, int m1, int m2, int m3) {
Rational J1 = new Rational(j1, 2);
Rational J2 = new Rational(j2, 2);
Rational J3 = new Rational(j3, 2);
Rational M1 = new Rational(m1, 2);
Rational M2 = new Rational(m2, 2);
Rational M3 = new Rational(m3, 2);
return wigner3jm(J1, J2, J3, M1, M2, M3);
} /* wigner3jm */
/**
* Wigner 3jn symbol. For the 6j symbol, the input of the 3 lines is
* "1 2 3 1 5 6", "4 5 3 4 2 6" "2j1+1 2j2+1 2j3+1 2l1+1 2l2+1 2l3+1"
*
* @param m1
* The information on the number of angular momenta.
* @param t1
* The list of one half of the triads, indexing j, whitespace
* separated
* @param t2
* The list of the second half of the triads, indexing j,
* whitespace separated
* @param j
* The list of the integer values of the angular momenta. They
* are actually the doubled j-values plus 1, whitespace
* separated. Only as many as announced by the m1 parameter are
* used; trailing numbers are ignored.
* @see A. Bar-Shalom and M. Klapisch,
* <a href="http://dx.doi.org/10.1016/0010-4655(88)90192-0">NJGRAF...
* </a>, Comp. Phys Comm. 50 (3) (1988) 375
* @since 2011-02-13
* @since 2012-02-15 Upgraded return value to BigSurdVec
* @author Richard J. Mathar
*/
static public BigSurdVec wigner3j(String m1, String t1, String t2, String j) {
/*
* The first number in the line "m" is the number of angular momenta.
* The rest of the line is ignored.
*/
Scanner s = new Scanner(m1);
int m = s.nextInt();
if (m % 3 != 0) {
s.close();
throw new IllegalArgumentException("Angular momenta " + m + " not a multiple of three.");
}
/*
* Scan the numbers in the line "j". Excess numbers beyond what has been
* announced in the "m" line are ignored.
*/
int[] jvec = new int[m];
int[] tvec = new int[2 * m];
s.close();
/*
* the third row contains positive 2j+1.
*/
s = new Scanner(j);
int ji = 0;
while (s.hasNextInt() && ji < m) {
jvec[ji++] = s.nextInt();
if (jvec[ji - 1] < 1) {
s.close();
throw new IllegalArgumentException("Illegal value " + jvec[ji - 1] + " for 2j+1.");
}
}
s.close();
/*
* the first two rows contain signed values of indices into the j list
*/
s = new Scanner(t1);
int ti = 0;
while (s.hasNextInt())
tvec[ti++] = s.nextInt();
s.close();
s = new Scanner(t2);
while (s.hasNextInt())
tvec[ti++] = s.nextInt();
/*
* Basic sanity checks. All indices in the first two lines address a
* number in the third line, and each index occurs exactly twice.
*/
if (ji % 3 != 0) {
s.close();
throw new IllegalArgumentException("j-count " + ji + " not a multiple of three.");
}
if (ti != 2 * ji) {
s.close();
throw new IllegalArgumentException("triad-count " + ti + " not twice j-count " + ji);
}
int[] jfreq = new int[m];
for (ji = 0; ji < jfreq.length; ji++)
jfreq[ji] = 0;
/*
* maintain a 0-based index which shows where the j-value has its first
* and second occurrence in the flattened list of triads.
*/
int[][] jhash = new int[m][2];
for (ti = 0; ti < 2 * m; ti++) {
int t = tvec[ti];
if (t == 0 || Math.abs(t) > jvec.length) {
s.close();
throw new IllegalArgumentException("Triad index " + t + " out of bounds");
}
if (jfreq[Math.abs(t) - 1] >= 2) {
s.close();
throw new IllegalArgumentException("Node " + t + " referenced more than twice");
}
jhash[Math.abs(t) - 1][jfreq[Math.abs(t) - 1]] = ti;
jfreq[Math.abs(t) - 1]++;
}
/*
* Move on from the 2j+1 values of the input to the j-values. Subtract
* one and divide through 2.
*/
Rational[] J = new Rational[jvec.length];
for (ji = 0; ji < jvec.length; ji++) {
J[ji] = new Rational(jvec[ji] - 1, 2);
}
/*
* Convert the 1-based indices to 0-based indices, loosing the sign
* information.
*/
int[] triadidx = new int[tvec.length];
for (ti = 0; ti < tvec.length; ti++)
triadidx[ti] = Math.abs(tvec[ti]) - 1;
/*
* The M-values are all null (undetermined) at the start.
*/
Rational[] M = new Rational[J.length];
s.close();
return wigner3j(tvec, J, M, triadidx);
} /* wigner3j */
/**
* Wigner 3jn symbol. Computes sum_{mi} (-1)^(j1-m1+j2-m2+...)
* triad(triadidx[0..2])*triad(triadidx[3..5])*... where each factor is a
* Wigner-3jm symbol with each sign of m_i occurring once at the
* corresponding l-value.
*
* @param triadidx
* 0-based indices into the list of J
* @param J
* The list of J-values
* @param M
* The list of M-values associated with the J. This contains null
* where the parameter has not yet been set by an outer loop.
* @since 2011-02-13
* @since 2012-02-15 Upgraded to return BigSurdVec
* @author Richard J. Mathar
*/
static private BigSurdVec wigner3j(final int[] tvec, final Rational[] J, final Rational[] M, final int[] triadidx) {
/*
* The result of the computation. The sum over all m-combinations of the
* triads.
*/
BigSurdVec res = new BigSurdVec();
/*
* First step is to monitor the triangular conditions on the J. If at
* least one is violated, the result is zero. Loop over the triads.
*/
for (int t = 0; t < triadidx.length; t += 3) {
/* Ensure |J[t]-J[t+1]| <= J[t+2] <= J[t]+J[t+1] */
if (J[triadidx[t]].subtract(J[triadidx[t + 1]]).abs().compareTo(J[triadidx[t + 2]]) > 0)
return res;
if (J[triadidx[t]].add(J[triadidx[t + 1]]).compareTo(J[triadidx[t + 2]]) < 0)
return res;
}
/*
* the index of the preferred member of the triad list. Preference given
* to those dangling in triads where alreaday two others are fixed, then
* to members where at least one is fixed, then to smallest associated
* J-values.
*/
int freeM = -1;
int freeMrank = -1;
for (int i = 0; i < triadidx.length; i++) {
/*
* found an m-value which has not yet been summed over.
*/
if (M[triadidx[i]] == null) {
/*
* two cases: value is fixed implicitly because already two
* others values are set in the triad. or it is still to
* maintain its own explicit loop.
*/
int triadn = i / 3;
int triadr = i % 3;
/*
* the neighbors in the triad have indices triadn*3+ (tiradr+1)
* mod 3 and triadn*3+(triadr+2) mod3
*/
int nei1 = 3 * triadn + (triadr + 1) % 3;
int nei2 = 3 * triadn + (triadr + 2) % 3;
/*
* found a candidate for which the two other values are already
* set.
*/
if (M[triadidx[nei1]] != null && M[triadidx[nei2]] != null) {
freeM = i;
break;
} else {
/*
* rough work load estimator: basically (2J1+1)*(2J2+1)
*/
Rational wt = J[triadidx[i]].multiply(2).add(1);
if (M[triadidx[nei1]] == null)
wt = wt.multiply(J[triadidx[nei1]].multiply(2).add(1));
if (M[triadidx[nei2]] == null)
wt = wt.multiply(J[triadidx[nei2]].multiply(2).add(1));
int thiswt = wt.intValue();
if (freeM < 0 || thiswt < freeMrank) {
freeM = i;
freeMrank = thiswt;
}
}
}
}
if (freeM >= 0) {
/*
* found an m-value which has not yet been summed over.
*/
if (M[triadidx[freeM]] == null) {
Rational[] childM = new Rational[M.length];
for (int ji = 0; ji < M.length; ji++)
if (M[ji] != null)
childM[ji] = M[ji];
/*
* two cases: value is fixed implicitly because already two
* others values are set in the triad. or it is still to
* maintain its own explicit loop.
*/
int triadn = freeM / 3;
int triadr = freeM % 3;
/*
* the neighbors in the triad have indices triadn*3+ (triadr+1)
* mod 3 and triadn*3+(triadr+2) mod3
*/
int nei1 = 3 * triadn + (triadr + 1) % 3;
int nei2 = 3 * triadn + (triadr + 2) % 3;
if (M[triadidx[nei1]] == null || M[triadidx[nei2]] == null) {
/*
* The J-value is J[triadidx[freeM]]. Loop from -J to +J,
* the allowed range.
*/
Rational newm = J[triadidx[freeM]].negate();
while (newm.compareTo(J[triadidx[freeM]]) <= 0) {
childM[triadidx[freeM]] = tvec[freeM] > 0 ? newm : newm.negate();
res = res.add(wigner3j(tvec, J, childM, triadidx));
newm = newm.add(Rational.ONE);
}
} else {
/*
* Set its value and the value at its companion j-value. Sum
* of the three m-values in the triad is to be zero for a
* non-zero contribution.
*/
Rational m1 = M[triadidx[nei1]];
Rational m2 = M[triadidx[nei2]];
/*
* negate if these are the second occurrences of the J in
* the triads
*/
if (tvec[nei1] < 0)
m1 = m1.negate();
if (tvec[nei2] < 0)
m2 = m2.negate();
/* m3 = -(m1+m2) */
Rational newm = tvec[freeM] > 0 ? m1.add(m2).negate() : m1.add(m2);
/*
* No contribution if the m-value enforced by the other two
* entries is outside the range -|J|..|J| enforced by its
* associated J-value. One could essentially remove this
* branching and let wigner3j() decide on this, but this is
* inefficient.
*/
if (newm.abs().compareTo(J[triadidx[freeM]]) <= 0) {
childM[triadidx[freeM]] = newm;
res = res.add(wigner3j(tvec, J, childM, triadidx));
}
/*
* zero contribution if this m-value cannot be set to any
* value compatible with the triangular conditions.
*/
}
return res;
}
}
/*
* reached the bottom of the loop where all M-values are assigned. Build
* the product over all Wigner-3jm values and the associated sign.
*/
res = BigSurdVec.ONE;
for (int ji = 0; ji < triadidx.length; ji += 3) {
Rational m1 = M[triadidx[ji]];
Rational m2 = M[triadidx[ji + 1]];
Rational m3 = M[triadidx[ji + 2]];
/*
* negate if these are associated with in-flowing vectors in the
* triads
*/
if (tvec[ji] < 0)
m1 = m1.negate();
if (tvec[ji + 1] < 0)
m2 = m2.negate();
if (tvec[ji + 2] < 0)
m3 = m3.negate();
res = res.multiply(wigner3jm(J[triadidx[ji]], J[triadidx[ji + 1]], J[triadidx[ji + 2]], m1, m2, m3));
/*
* if a partial product yields zero, the total product is zero, too,
* and offers an early exit.
*/
if (res.signum() == 0)
return BigSurdVec.ZERO;
}
/*
* The overal sign is product_{J-Mpairs} (-1)^(J-M). This is an integer
* because all the J-M are integer.
*/
Rational sig = new Rational();
for (int ji = 0; ji < J.length; ji++)
sig = sig.add(J[ji]).subtract(M[ji]);
/*
* sign depends on the sum being even or odd. We assume that "sig" is
* integer and look only at the numerator
*/
if (sig.a.abs().testBit(0))
res = res.negate();
return res;
} /* wigner3j */
/**
* The Wigner 3jm symbol (j1,j2,j3,m1,m2,m3). Warning: there is no check
* that each argument is indeed half-integer.
*
* @param j1
* integer or half-integer j1
* @param j2
* integer or half-integer j2
* @param j3
* integer or half-integer j3
* @param m1
* integer or half-integer m1
* @param m2
* integer or half-integer m2
* @param m3
* integer or half-integer m3
* @return The value of the symbol. Zero if any of the triangular
* inequalities is violated or some parameters are out of range.
* @since 2011-02-13
* @author Richard J. Mathar
*/
static protected BigSurd wigner3jm(Rational j1, Rational j2, Rational j3, Rational m1, Rational m2, Rational m3) {
/*
* Check that m1+m2+m3 = 0
*/
if (m1.add(m2).add(m3).signum() != 0)
return BigSurd.ZERO;
/*
* Check that j1+j2+j3 is integer
*/
if (j1.add(j2).add(j3).isBigInteger() == false)
return BigSurd.ZERO;
/*
* Check that |j1-j2|<=j3 <= |j1+j2|
*/
Rational j1m2 = j1.subtract(j2);
if (j1m2.abs().compareTo(j3) > 0)
return BigSurd.ZERO;
Rational j1p2 = j1.add(j2);
if (j1p2.abs().compareTo(j3) < 0)
return BigSurd.ZERO;
/*
* Check that |m_i| <= j_i
*/
if (m1.abs().compareTo(j1) > 0 || m2.abs().compareTo(j2) > 0 || m3.abs().compareTo(j3) > 0)
return BigSurd.ZERO;
/*
* Check that m_i-j_i are integer.
*/
if (!m1.subtract(j1).isBigInteger() || !m2.subtract(j2).isBigInteger() || !m3.subtract(j3).isBigInteger())
return BigSurd.ZERO;
/*
* (-)^(j1-j2-m3)*delta(-m3,m1+m2)*sqrt[ (j3+j1-j2)! (j3-j1+j2)!
* (j1+j2-j3)! /(j1+j2+j3+1)!
* *(j3-m)!*(j3+m)!(j1-m1)!*(j1+m1)!*(j2-m2)!*(j2+m2)!] *sum_k
* (-1)^k/[k!(j1+j2-j3-k)!(j1-m1-k)!(j2+m2-k)!(j3-j2+m1+k)!)*(j3-j1-m2+k
* )!]
*/
/*
* It is tacitly assumed that all the major j_i, m_i values are in the
* integer range. This is implicitly plausible since otherwise the
* execution times of the following loop over the k-values would be
* immense.
*/
int j1j2jk = j1p2.subtract(j3).intValue();
int j1m1k = j1.subtract(m1).intValue();
int j2m2k = j2.add(m2).intValue();
int jj2m1k = j3.subtract(j2).add(m1).intValue();
int jj1m2k = j3.subtract(j1).subtract(m2).intValue();
int k = Math.max(0, -jj2m1k);
k = Math.max(k, -jj1m2k);
if (k > 0) {
j1j2jk -= k;
j1m1k -= k;
j2m2k -= k;
jj2m1k += k;
jj1m2k += k;
}
Factorial f = new Factorial();
Rational sumk = new Rational();
while (true) {
BigInteger d = f.at(k).multiply(f.at(j1j2jk)).multiply(f.at(j1m1k)).multiply(f.at(j2m2k))
.multiply(f.at(jj2m1k)).multiply(f.at(jj1m2k));
if (k % 2 == 0)
sumk = sumk.add(new Rational(BigInteger.ONE, d));
else
sumk = sumk.subtract(new Rational(BigInteger.ONE, d));
j1j2jk--;
j1m1k--;
j2m2k--;
jj2m1k++;
jj1m2k++;
if (j1j2jk < 0 || j1m1k < 0 || j2m2k < 0)
break;
k++;
}
/*
* sign factor (-1)^(j1-j2-m3)
*/
if (j1m2.subtract(m3).intValue() % 2 != 0)
sumk = sumk.negate();
k = j1m2.add(j3).intValue();
BigInteger s = f.at(k);
k = j3.subtract(j1m2).intValue();
s = s.multiply(f.at(k));
k = j1p2.subtract(j3).intValue();
s = s.multiply(f.at(k));
k = j3.add(m3).intValue();
s = s.multiply(f.at(k));
k = j3.subtract(m3).intValue();
s = s.multiply(f.at(k));
k = j1.add(m1).intValue();
s = s.multiply(f.at(k));
k = j1.subtract(m1).intValue();
s = s.multiply(f.at(k));
k = j2.add(m2).intValue();
s = s.multiply(f.at(k));
k = j2.subtract(m2).intValue();
s = s.multiply(f.at(k));
k = j1p2.add(j3).intValue();
k++;
Rational disc = new Rational(s, f.at(k));
return new BigSurd(sumk, disc);
} /* wigner3jm */
} /* Wigner3j */

329
src/org/nevec/rjm/Wigner3jGUI.java Normal file
View File

@ -0,0 +1,329 @@
package org.nevec.rjm;
import java.awt.Color;
import java.awt.Font;
import java.awt.GridBagConstraints;
import java.awt.GridBagLayout;
import java.awt.Label;
import java.awt.TextArea;
import java.awt.event.ActionEvent;
import java.awt.event.ActionListener;
import java.util.Scanner;
import javax.swing.JButton;
import javax.swing.JFrame;
import javax.swing.JList;
import javax.swing.event.ListSelectionEvent;
import javax.swing.event.ListSelectionListener;
/**
* An interactive interface to the Wigner3j class. The GUI allows to preselect
* one of the symbols if the number of j-terms is small (6j up to 15j), or to
* enter any other connectivity for the triads of j-values. The actual j-values
* are entered as integers (2j+1) and the computation of one value (in exact
* square root representation) is started manually.
*
* @since 2011-02-15
*/
public class Wigner3jGUI implements ActionListener, ListSelectionListener {
/**
* The master window of the session
*/
JFrame fram;
/*
* global labels
*/
Label Lbl0;
Label Lbl1;
JButton sear;
JList<?> searJ;
String[] searOpt = { "6j", "9j", "12j 1st", "12j 2nd (not symm)", "15j 1st", "15j 2nd", "15j 3rd", "15j 4th",
"15j 5th" };
/**
* Field with the triads inputs
*/
TextArea inpGtria;
/**
* Field with the J-value inputs
*/
TextArea inpGjval;
/**
* Field of the outputs.
*/
TextArea outG;
GridBagLayout gridbag;
GridBagConstraints gridconstr;
/**
* @since 2011-02-15
*/
public void init() {
fram = new JFrame("Wigner3jGUI");
Lbl0 = new Label("Input: (Triads upper area, values 2J+1 second area");
Lbl1 = new Label("Output:");
sear = new JButton("Compute");
sear.setActionCommand("compute");
sear.addActionListener(this);
sear.setToolTipText("Compute a general 3jn value");
searJ = new JList<Object>(searOpt);
searJ.setLayoutOrientation(JList.HORIZONTAL_WRAP);
searJ.addListSelectionListener(this);
Font defFont = new Font("Monospaced", Font.PLAIN, 11);
fram.setBackground(new Color(250, 250, 250));
fram.setForeground(new Color(0, 0, 0));
Color fg = new Color(0, 200, 0);
Color bg = new Color(10, 10, 10);
gridbag = new GridBagLayout();
fram.setLayout(gridbag);
gridconstr = new GridBagConstraints();
gridconstr.gridx = 0;
gridconstr.gridy = GridBagConstraints.RELATIVE;
inpGtria = new TextArea("", 4, 80);
inpGtria.setFont(defFont);
inpGtria.setForeground(fg);
inpGtria.setBackground(bg);
inpGjval = new TextArea("", 10, 80);
inpGjval.setFont(defFont);
inpGjval.setForeground(fg);
inpGjval.setBackground(bg);
outG = new TextArea("", 12, 80);
outG.setEditable(false);
outG.setFont(defFont);
outG.setForeground(fg);
outG.setBackground(bg);
fram.add(Lbl0);
gridbag.setConstraints(Lbl0, gridconstr);
fram.add(inpGtria);
gridbag.setConstraints(inpGtria, gridconstr);
fram.add(inpGjval);
gridbag.setConstraints(inpGjval, gridconstr);
fram.add(sear);
gridbag.setConstraints(sear, gridconstr);
fram.add(searJ);
gridbag.setConstraints(searJ, gridconstr);
fram.add(Lbl1);
gridbag.setConstraints(Lbl1, gridconstr);
fram.add(outG);
gridbag.setConstraints(outG, gridconstr);
fram.pack();
fram.setVisible(true);
} /* init */
/**
* @since 2010-08-27
*/
public void compute() {
String tr = inpGtria.getText();
String[] trias = new String[4];
/*
* Read the trias configuration from inpGtria into trias[0..2], skipping
* lines that start with a hash mark.
*/
Scanner s = new Scanner(tr);
for (int l = 0; l < 3;) {
try {
trias[l] = s.nextLine().trim();
if (!trias[l].startsWith("#"))
l++;
} catch (Exception e) {
s.close();
outG.setText("ERROR: less than 3 lines in the triad definition");
return;
}
}
s.close();
/*
* Read the J values from inpGjval into trias[3] in a loop
*/
String j = inpGjval.getText();
s = new Scanner(j);
while (true) {
try {
trias[3] = s.nextLine().trim();
} catch (Exception e) {
s.close();
return;
}
if (!trias[3].startsWith("#")) {
try {
BigSurdVec w = Wigner3j.wigner3j(trias[0], trias[1], trias[2], trias[3]);
outG.append(w.toString() + " = " + w.doubleValue());
} catch (Exception e) {
outG.append(e.toString());
e.printStackTrace();
}
outG.append(" # J = ");
Scanner num = new Scanner(trias[3]);
while (num.hasNextInt()) {
int twoj1 = num.nextInt();
Rational jfrac = new Rational(twoj1 - 1, 2);
outG.append(jfrac.toString() + " ");
}
outG.append("\n");
num.close();
}
}
} /* compute */
/**
* Interpreter parser loop.
*
* @param e
* the information on which button had been pressed in the GUI
* @since 2011-02-15
*/
public void actionPerformed(ActionEvent e) {
String lin = e.getActionCommand();
/*
* debugging System.out.println("Ac"+e.paramString()) ;
* System.out.println(lin) ;
*/
if (lin == "compute") {
outG.setText("");
compute();
}
} /* actionPerformed */
/**
* Interpreter parser loop.
*
* @param e
* the information on which of the 3jn templates had been
* selected in the Menu
* @since 2011-02-18
*/
public void valueChanged(ListSelectionEvent e) {
switch (searJ.getMinSelectionIndex()) {
case 0:
inpGtria.setText("6\n");
inpGtria.append("1 2 -3 -1 5 6\n");
inpGtria.append("4 -5 3 -4 -2 -6");
outG.setText("");
break;
case 1:
/*
* Yutsis Figure 18.1 index map. j1=1, j2=2, j3=3 k1=4, k2=5, k3=6
* l1=7, l2=8, l3=9
*/
inpGtria.setText("9\n");
inpGtria.append("1 3 2 4 6 5 7 9 8 # (j1 j3 j2) (k1 k3 k2) (l1 l3 l2)\n");
inpGtria.append("-2 -8 -5 -6 -3 -9 -7 -4 -1 # (j2 l2 k2) (k3 j3 l3) (l1 k1 j1)");
outG.setText("");
break;
case 2:
/*
* Yutsis Figure 19.1 and 19.2, index map, including the sign
* reveral of the l. Assume input order j1..j4, l1..l4, k1..k4.
* j1=1, j2=2, j3=3, j4=4 l1=5, l2=6, l3=7, l4=8 k1=9, k2=10, k3=11,
* k4=12
*/
inpGtria.setText("12\n");
inpGtria.append("1 12 -8 -1 5 -2 2 6 -3 3 7 -4 # (j1 k4 l4) (j1 l1 j2) (j2 l2 j3) (j3 l3 j4)\n");
inpGtria.append("4 8 -9 9 -5 -10 10 -6 -11 11 -7 -12 # (j4 l4 k1) (k1 l1 k2) (k2 l2 k3) (k3 l3 k4)");
outG.setText("");
break;
case 3:
inpGtria.setText("12\n");
inpGtria.append("1 5 9 -9 -2 -7 2 11 8 -8 -12 -4 # (j1 l1 k1) (k1 j2 l3 ) (j2 k3 l4) (l4 k4 j4)\n");
inpGtria.append("4 7 10 -10 -3 -5 3 6 12 -6 -11 -1 # (j4 l3 k2) (k2 j3 l1) (j3 l2 k4) (l2 k3 j1)");
outG.setText("");
break;
case 4:
/*
* Yutsis Figure 20.2 to 20.3, index map. j1=1, j2=2, j3=3, j4=4,
* j5=5 l1=6, l2=7, l3=8, l4=9, l5=10 k1=11, k2=12, k3=13, k4=14,
* k5=15
*/
inpGtria.setText("15\n");
inpGtria.append(
"1 -6 2 -2 -7 3 -3 -8 4 -4 -9 5 -5 -10 11 # (j1 l1 j2)(j2 l2 j3)(j3 l3 j4)(j4 l4 j5)(j5 l5 k1)\n");
inpGtria.append(
"-11 6 12 -12 7 13 -13 8 14 -14 9 15 -15 10 -1 # (k1 l1 k2)(k2 l2 k3)(k3 l3 k4)(k4 l4 k5)(k5 l5 j1)");
outG.setText("");
break;
case 5:
inpGtria.setText("15\n");
inpGtria.append(
"-1 -6 2 -2 -7 3 -3 -8 4 -4 -9 5 1 -5 -10 # (j1 l1 j2)(j2 l2 j3)(j3 l3 j4)(j4 l4 j5)(j1 j5 l5)\n");
inpGtria.append(
"11 -15 10 9 15 -14 8 14 -13 7 13 -12 6 12 -11 # (k1 k5 l5)(l4 k5 k4)(l3 k4 k3)(l2 k3 k2)(l1 k2 k1)");
outG.setText("");
break;
case 6:
/*
* Yutsis Figure 20.4a, index map. k1=1, k1'=2, k=3, k'=4, k2=5,
* k2'=6 p1=7, p=8, p2=9, j1=10, j1'=11 j=12 j'=13 j2=14 j2'=15
*/
inpGtria.setText("15\n");
inpGtria.append(
"-13 -12 -8 12 14 10 -10 -1 7 -7 -11 -2 2 4 6 # (j' j p)(j j2 j1)(j1 k1 p1)(p1 j1' k1')(k1' k' k2')\n");
inpGtria.append(
"-4 -3 8 1 3 5 -14 -5 9 -15 -6 -9 15 11 13 # (k' k p)(k1 k k2)(j2 k2 p2)(j2' k2' p2)(j2' j1' j')");
outG.setText("");
break;
case 7:
/*
* Yutsis Figure 20.5a, index map. j1=1, k1=2 s1=3 k1'=4 j1'=5 p=6
* l=7 s=8 l'=9 p'=10 j2=11 k2=12 s2=13 k2'=14 j2'=15
*/
inpGtria.setText("15\n");
inpGtria.append(
"-14 -12 -8 12 11 -10 -11 13 -7 7 -1 3 2 1 6 # (k2' k2 s)(k2 j2 p')(j2 s2 l)(l j1 s1)(k1 j1 p)\n");
inpGtria.append(
"-4 -2 8 10 4 5 9 -5 -3 -13 -9 -15 15 -6 14 # (k1' k1 s)(p' k1' j1')(l' j1' s1)(s2 l' j2')(j2' p k2')");
outG.setText("");
break;
case 8:
/*
* Yutsis Figure 20.6, index map. k1=1 k1'=2 j1=3 l1=4 l1'=5 k2=6
* k2'=7 j2=8 l2=9 l2'=10 k3=11 k3'=12 j3=13 l3=14 l3'=15
*/
inpGtria.setText("15\n");
inpGtria.append(
"-15 1 -7 -4 -11 7 5 4 -3 -12 -5 6 12 -9 -1 # (l3' k1 k2')(l1 k3 k2')(l1' l1 j1)(k3' l1' k2)(k3' l2 k1)\n");
inpGtria.append(
"9 -8 10 -10 11 -2 -14 -6 2 14 -13 15 3 8 13 # (l2 j2 l2')(l2' k3 k1')(l3 k2 k1')(l3 j3 l3')(j1 j2 j3)");
outG.setText("");
break;
}
} /* valueChanged */
/**
* Main entry point. not taking any command line options:<br>
* java -jar Wigner3jGUI.jar<br>
*
* @since 2012-02-16
* @author Richard J. Mathar
*/
public static void main(String[] args) {
Wigner3jGUI g = new Wigner3jGUI();
g.init();
} /* main */
} /* Wigner3jGUI */

36
src/org/warpgate/pi/calculator/Calculator.java Normal file
View File

@ -0,0 +1,36 @@
package org.warpgate.pi.calculator;
import org.nevec.rjm.BigSurdVec;
public class Calculator {
public Calculator(boolean b) {
Utils.debugOn = b;
}
public Termine calcolarisultato(String string) throws Errore {
System.out.println("INPUT: " + string);
Parentesi espressione = new Parentesi(string, "");
return espressione.calcola();
}
public RisultatoEquazione calcolaequazione(String string) throws Errore {
if (string.split("=").length == 0) {
return new RisultatoEquazione(new Termine("0"), true);
}
if (string.split("=").length <= 2) {
if (string.split("=").length == 1) {
string = string + "=0";
}
Termine res1 = calcolarisultato(string.split("=")[0]);
Termine res2 = calcolarisultato(string.split("=")[1]);
Termine res = res1.add(res2.multiply(new Termine("-1")));
if (res.calcola().getTerm().toString().equals("0")) {
return new RisultatoEquazione(res.calcola(), true);
}
return new RisultatoEquazione(res.calcola(), false);
}
return new RisultatoEquazione(null, false);
}
}

23
src/org/warpgate/pi/calculator/Divisione.java Normal file
View File

@ -0,0 +1,23 @@
package org.warpgate.pi.calculator;
import org.nevec.rjm.BigSurdVec;
public class Divisione extends FunzioneDueValori {
public Divisione(Funzione value1, Funzione value2) {
super(value1, value2);
}
@Override
public String simbolo() {
return Simboli.DIVISION;
}
@Override
public Termine calcola() throws Errore {
if (getVariable2().calcola().getTerm().compareTo(BigSurdVec.ZERO) == 0) {
throw new Errore(Errori.DIVISION_BY_ZERO);
}
return getVariable1().calcola().divide(getVariable2().calcola());
}
}

15
src/org/warpgate/pi/calculator/Errore.java Normal file
View File

@ -0,0 +1,15 @@
package org.warpgate.pi.calculator;
public class Errore extends java.lang.Throwable {
/**
*
*/
private static final long serialVersionUID = -1014947815755694651L;
public Errore(Errori IDErrore) {
id = IDErrore;
}
public Errori id = Errori.ERROR;
}

8
src/org/warpgate/pi/calculator/Errori.java Normal file
View File

@ -0,0 +1,8 @@
package org.warpgate.pi.calculator;
public enum Errori {
ERROR,
DIVISION_BY_ZERO,
UNBALANCED_BRACKETS,
NOT_IMPLEMENTED
}

6
src/org/warpgate/pi/calculator/Funzione.java Normal file
View File

@ -0,0 +1,6 @@
package org.warpgate.pi.calculator;
public interface Funzione {
public String simbolo();
public Termine calcola() throws Errore;
}

28
src/org/warpgate/pi/calculator/FunzioneDueValori.java Normal file
View File

@ -0,0 +1,28 @@
package org.warpgate.pi.calculator;
import org.nevec.rjm.Rational;
public abstract class FunzioneDueValori implements Funzione {
public FunzioneDueValori(Funzione value1, Funzione value2) {
setVariable1(value1);
setVariable2(value2);
}
protected Funzione variable1 = new Termine(Rational.ZERO);
public Funzione getVariable1() {
return variable1;
}
public void setVariable1(Funzione value) {
variable1 = value;
}
protected Funzione variable2 = new Termine(Rational.ZERO);
public Funzione getVariable2() {
return variable2;
}
public void setVariable2(Funzione value) {
variable2 = value;
}
@Override
public abstract String simbolo();
@Override
public abstract Termine calcola() throws Errore;
}

54
src/org/warpgate/pi/calculator/FunzioneMultipla.java Normal file
View File

@ -0,0 +1,54 @@
package org.warpgate.pi.calculator;
import java.util.Arrays;
import java.util.List;
public abstract class FunzioneMultipla implements Funzione {
public FunzioneMultipla() {
setVariables(new Funzione[]{});
}
public FunzioneMultipla(Funzione[] values) {
setVariables(values);
}
protected Funzione[] variables;
public Funzione[] getVariables() {
return variables;
}
public void setVariables(Funzione[] value) {
variables = value;
}
public void setVariables(final List<Funzione> value) {
int vsize = value.size();
Funzione[] tmp = new Funzione[vsize];
for (int i = 0; i < vsize; i++) {
tmp[i] = value.get(i);
}
variables = tmp;
}
public Funzione getVariable(int index) {
return variables[index];
}
public void setVariable(int index, Funzione value) {
variables[index] = value;
}
public void addVariableToEnd(Funzione value) {
int index = variables.length;
setVariablesLength(index+1);
variables[index] = value;
}
public int getVariablesLength() {
return variables.length;
}
public void setVariablesLength(int length) {
variables = Arrays.copyOf(variables, length);
}
@Override
public abstract String simbolo();
@Override
public abstract Termine calcola() throws Errore;
}

20
src/org/warpgate/pi/calculator/FunzioneSingola.java Normal file
View File

@ -0,0 +1,20 @@
package org.warpgate.pi.calculator;
import org.nevec.rjm.BigSurdVec;
public abstract class FunzioneSingola implements Funzione {
public FunzioneSingola(Funzione value) {
setVariable(value);
}
protected Funzione variable = new Termine(BigSurdVec.ZERO);
public Funzione getVariable() {
return variable;
}
public void setVariable(Funzione value) {
variable = value;
}
@Override
public abstract String simbolo();
@Override
public abstract Termine calcola() throws Errore;
}

27
src/org/warpgate/pi/calculator/Main.java Normal file
View File

@ -0,0 +1,27 @@
package org.warpgate.pi.calculator;
public class Main {
public static void main(String[] args) {
Calculator c = new Calculator(false);
try {
long start = System.nanoTime();
Termine result = c.calcolarisultato("((5^2+3√(100/0.1))*Ⓐ7+9/15*2√(26/2))/21");
long end = System.nanoTime();
long timeElapsed = end-start;
System.out.println("RESULT: " + result);
System.out.println("DECIMAl RESULT: " + result.getTerm().toBigDecimal());
System.out.println("Time elapsed: " + (double) timeElapsed / 1000000000 + "\n");
start = System.nanoTime();
RisultatoEquazione eresult = c.calcolaequazione("((5^2+3√(100/0.1))*Ⓐ7+9/15*2√(26/2))/21=(175*(2√7)+3*(2√13))/105");
end = System.nanoTime();
timeElapsed = end-start;
System.out.println("Is an equation: " + eresult.isAnEquation);
System.out.println("L-R: " + eresult.LR);
System.out.println("Time elapsed: " + (((double) timeElapsed / 1000000000)) + "\n");
} catch (Errore e) {
System.err.println(e.id);
}
}
}

18
src/org/warpgate/pi/calculator/Moltiplicazione.java Normal file
View File

@ -0,0 +1,18 @@
package org.warpgate.pi.calculator;
public class Moltiplicazione extends FunzioneDueValori {
public Moltiplicazione(Funzione value1, Funzione value2) {
super(value1, value2);
}
@Override
public String simbolo() {
return Simboli.MULTIPLICATION;
}
@Override
public Termine calcola() throws Errore {
return getVariable1().calcola().multiply(getVariable2().calcola());
}
}

347
src/org/warpgate/pi/calculator/Parentesi.java Normal file
View File

@ -0,0 +1,347 @@
package org.warpgate.pi.calculator;
import static org.warpgate.pi.calculator.Utils.ArrayToRegex;
import static org.warpgate.pi.calculator.Utils.concat;
import java.util.ArrayList;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import org.nevec.rjm.BigSurdVec;
public class Parentesi extends FunzioneMultipla {
public Parentesi() {
super();
}
public Parentesi(Funzione[] values) {
super(values);
}
public Parentesi(String string, String debugSpaces) throws Errore {
super();
try{
//Se l'espressione è già un numero:
setVariables(new Funzione[]{new Termine(string)});
Utils.debug.println(debugSpaces+"•Result:"+string);
} catch (NumberFormatException ex) {
String processExpression = string;
Utils.debug.println(debugSpaces+"•Analyzing expression:"+processExpression);
debugSpaces+= " ";
//Se l'espressione non è già un numero:
//Correggi i segni ++ e -- in eccesso
Pattern pattern = Pattern.compile("\\+\\++?|\\-\\-+?");
Matcher matcher = pattern.matcher(processExpression);
boolean cambiati = false;
while (matcher.find()) {
cambiati = true;
String correzione = "+";
processExpression = processExpression.substring(0, matcher.start(0)+1)+correzione+processExpression.substring(matcher.start(0)+matcher.group(0).length(), processExpression.length());
matcher = pattern.matcher(processExpression);
}
//Correggi i segni +- e -+ in eccesso
pattern = Pattern.compile("\\+\\-|\\-\\+");
matcher = pattern.matcher(processExpression);
while (matcher.find()) {
cambiati = true;
String correzione = "-";
processExpression = processExpression.substring(0, matcher.start(0)+1)+correzione+processExpression.substring(matcher.start(0)+matcher.group(0).length(), processExpression.length());
matcher = pattern.matcher(processExpression);
}
//Rimuovi i segni appena dopo le parentesi
if (processExpression.contains("(+")) {
cambiati = true;
processExpression = processExpression.replace("(+", "(");
}
//Rimuovi i segni appena dopo l'inizio
if (processExpression.startsWith("+")) {
cambiati = true;
processExpression = processExpression.substring(1, processExpression.length());
}
//Rimuovi i + in eccesso
pattern = Pattern.compile("["+
ArrayToRegex(Utils.add(concat(Simboli.segni(true), Simboli.funzioni()), "("))
+"]\\+[^"+
ArrayToRegex(concat(concat(Simboli.segni(true), Simboli.funzioni()), new String[]{"(", ")"}))
+"]+?["+
ArrayToRegex(concat(Simboli.segni(true), Simboli.funzioni()))
+"]|["+
ArrayToRegex(concat(Simboli.segni(true), Simboli.funzioni()))
+"]+?\\+[^"+
ArrayToRegex(concat(concat(Simboli.segni(true), Simboli.funzioni()), new String[]{"(", ")"}))
+"]");
matcher = pattern.matcher(processExpression);
cambiati = false;
while (matcher.find()) {
cambiati = true;
String correzione = matcher.group(0).replaceFirst(Matcher.quoteReplacement("+"), "");
processExpression = processExpression.substring(0, matcher.start(0)+1)+correzione+processExpression.substring(matcher.start(0)+matcher.group(0).length(), processExpression.length());
matcher = pattern.matcher(processExpression);
}
//Correggi i segni - in +-
pattern = Pattern.compile("[^"+Utils.ArrayToRegex(concat(concat(Simboli.funzioni(),Simboli.parentesi()), Simboli.segni(true)))+"]-");
matcher = pattern.matcher(processExpression);
while (matcher.find()) {
cambiati = true;
String correzione = "+-";
processExpression = processExpression.substring(0, matcher.start(0)+1)+correzione+processExpression.substring(matcher.start(0)+matcher.group(0).length(), processExpression.length());
matcher = pattern.matcher(processExpression);
}
if (cambiati) {
Utils.debug.println(debugSpaces+"•Resolved signs:"+processExpression);
}
//Aggiungi i segni * accanto alle parentesi
pattern = Pattern.compile("\\([^\\(]+?\\)");
matcher = pattern.matcher(processExpression);
cambiati = false;
while (matcher.find()) {
cambiati = true;
//sistema i segni * impliciti prima e dopo l'espressione.
String beforeexp = processExpression.substring(0, matcher.start(0));
String newexp = matcher.group(0).substring(1, matcher.group(0).length()-1);
String afterexp = processExpression.substring(matcher.start(0)+matcher.group(0).length(), processExpression.length());
if (Pattern.compile("[^"+Utils.ArrayToRegex(Utils.add(concat(Simboli.funzioni(), Simboli.segni(true)), "("))+"]$").matcher(beforeexp).find()) {
//Se la stringa precedente finisce con un numero
beforeexp += "*";
}
if (Pattern.compile("^[^"+Utils.ArrayToRegex(Utils.add(concat(Simboli.funzioni(), Simboli.segni(true)), ")"))+"]").matcher(afterexp).find()) {
//Se la stringa successiva inizia con un numero
afterexp = "*"+afterexp;
}
processExpression = beforeexp+"⑴"+newexp+"⑵"+afterexp;
matcher = pattern.matcher(processExpression);
}
processExpression = processExpression.replace("⑴", "(").replace("⑵", ")");
if (cambiati) {
Utils.debug.println(debugSpaces+"•Added implicit multiplications:"+processExpression);
}
Utils.debug.println(debugSpaces+"•Subdivision in classes:");
debugSpaces += " ";
//Suddividi tutto
Parentesi parentesiNonSuddivisaCorrettamente = new Parentesi();
parentesiNonSuddivisaCorrettamente.setVariables(new Funzione[]{});
String tmp = "";
final String[] funzioni = concat(concat(Simboli.funzioni(), Simboli.parentesi()), Simboli.segni(true));
for (int i = 0; i < processExpression.length(); i++) {
//Per ogni carattere cerca se è un numero o una funzione:
String charI = processExpression.charAt(i)+"";
if (Utils.isInArray(charI, funzioni)) {
//Cerca il tipo di funzione tra le esistenti
Funzione f = null;
switch (charI) {
case Simboli.SUM:
f = new Somma(null, null);
break;
case Simboli.MULTIPLICATION:
f = new Moltiplicazione(null, null);
break;
case Simboli.DIVISION:
f = new Divisione(null, null);
break;
case Simboli.NTH_ROOT:
f = new Radice(null, null);
break;
case Simboli.SQUARE_ROOT:
f = new RadiceQuadrata(null);
break;
case Simboli.POTENZA:
f = new Potenza(null, null);
break;
case Simboli.PARENTHESIS_OPEN:
//cerca l'ultima parentesi chiusa
int startIndex = i;
int endIndex = -1;
int jumps = -1;
for (int i2 = startIndex; i2 < processExpression.length(); i2++) {
if ((processExpression.charAt(i2)+"").equals(Simboli.PARENTHESIS_CLOSE)) {
if (jumps == 0) {
endIndex = i2;
break;
} else if (jumps > 0) {
jumps -= 1;
} else if (jumps < 0) {
throw new Errore(Errori.UNBALANCED_BRACKETS);
}
} else if ((processExpression.charAt(i2)+"").equals(Simboli.PARENTHESIS_OPEN)) {
jumps +=1;
}
}
if (endIndex == -1 || endIndex < startIndex) {
throw new Errore(Errori.UNBALANCED_BRACKETS);
}
startIndex += 1;
i = startIndex;
String tmpExpr = "";
while (i < endIndex) {
tmpExpr += processExpression.charAt(i);
i++;
}
f = new Parentesi(tmpExpr, debugSpaces);
break;
default:
throw new java.lang.RuntimeException("Il carattere "+charI+" non è tra le funzioni designate!\nAggiungerlo ad esse o rimuovere il carattere dall'espressione!");
}
if (f instanceof Parentesi) {
tmp = "";
} else {
if (tmp.length() != 0) {
parentesiNonSuddivisaCorrettamente.addVariableToEnd(new Termine(tmp));
Utils.debug.println(debugSpaces+"•Added variable to expression:"+tmp);
}
}
parentesiNonSuddivisaCorrettamente.addVariableToEnd(f);
Utils.debug.println(debugSpaces+"•Added variable to expression:"+f.simbolo());
tmp = "";
} else {
try{
if (!(charI.equals("-") || charI.equals("."))) {
Double.parseDouble(tmp + charI);
}
//Se il carattere è un numero intero, un segno negativo, o un punto
tmp += charI;
} catch (NumberFormatException exc) {
throw new java.lang.RuntimeException("Il carattere "+tmp+charI+" non è nè un numero nè un espressione presente nella lista completa!\nAggiungerlo ad essa o rimuovere il carattere dall'espressione!");
}
}
}
if (tmp.length() > 0) {
Utils.debug.println(debugSpaces+"•Added variable to expression:"+tmp);
parentesiNonSuddivisaCorrettamente.addVariableToEnd(new Termine(tmp));
tmp = "";
}
int dsl = debugSpaces.length();debugSpaces = "";for (int i = 0; i < dsl-2; i++) {debugSpaces += " ";}
Utils.debug.println(debugSpaces+"•Finished the subdivision in classes.");
//Fine suddivisione di insieme
//Inizia l'affinazione dell'espressione
Utils.debug.println(debugSpaces+"•Pushing classes...");
Funzione[] funzioniOLDArray = parentesiNonSuddivisaCorrettamente.getVariables();
ArrayList<Funzione> funzioniOLD = new ArrayList<Funzione>();
for (int i = 0; i < funzioniOLDArray.length; i++) {
if (funzioniOLDArray[i] != null) {
funzioniOLD.add(funzioniOLDArray[i]);
}
}
Utils.debug.println(debugSpaces+"•Correcting classes:");
debugSpaces += " ";
int before = 0;
String fase = "funzioniSN";
int n = 0;
do {
before = funzioniOLD.size();
int i = 0;
boolean change = false;
if (Utils.ciSonoSoloNumeriESomme(funzioniOLD)) {
fase = "somme"; //QUARTA FASE
} else if (Utils.ciSonoFunzioniSNnonImpostate(funzioniOLD)) {
fase = "funzioniSN"; // PRIMA FASE
} else if (Utils.ciSonoAltreFunzioni(funzioniOLD)) {
fase = "funzioniNSN"; //SECONDA FASE
} else {
fase = "moltiplicazioni"; //TERZA FASE
}
while (i < funzioniOLD.size() && change == false) {
Funzione funzioneTMP = funzioniOLD.get(i);
if (funzioneTMP instanceof FunzioneDueValori) {
if (fase != "funzioniSN") {
if (
fase == "somme" && (funzioneTMP instanceof Somma) == true
||
(fase == "moltiplicazioni" && (funzioneTMP instanceof Somma) == false)
||
(fase == "funzioniNSN" && (funzioneTMP instanceof Somma) == false && (funzioneTMP instanceof Moltiplicazione) == false && (funzioneTMP instanceof Divisione) == false)
) {
change = true;
if (i+1 < funzioniOLD.size() && i-1 >= 0 ) {
((FunzioneDueValori) funzioneTMP).setVariable1(funzioniOLD.get(i-1));
((FunzioneDueValori) funzioneTMP).setVariable2(funzioniOLD.get(i+1));
funzioniOLD.set(i, funzioneTMP.calcola());
//è importante togliere prima gli elementi in fondo e poi quelli davanti, perché gli indici scalano da destra a sinistra.
funzioniOLD.remove(i+1);
funzioniOLD.remove(i-1);
Utils.debug.println(debugSpaces+"•Set variable to expression:"+funzioneTMP.simbolo());
Utils.debug.println(debugSpaces+" "+"var1="+((FunzioneDueValori) funzioneTMP).getVariable1().calcola());
Utils.debug.println(debugSpaces+" "+"var2="+((FunzioneDueValori) funzioneTMP).getVariable2().calcola());
Utils.debug.println(debugSpaces+" "+"(result)="+((FunzioneDueValori) funzioneTMP).calcola());
} else {
throw new java.lang.RuntimeException("Argomenti mancanti! Sistemare l'equazione!");
}
}
}
} else if (funzioneTMP instanceof FunzioneSingola) {
if (fase == "funzioniSN") {
change = true;
if (i+1 < funzioniOLD.size()) {
((FunzioneSingola) funzioneTMP).setVariable(funzioniOLD.get(i+1));
funzioniOLD.set(i, funzioneTMP);
//è importante togliere prima gli elementi in fondo e poi quelli davanti, perché gli indici scalano da destra a sinistra.
funzioniOLD.remove(i+1);
Utils.debug.println(debugSpaces+"•Set variable to expression:"+funzioneTMP.simbolo());
Utils.debug.println(debugSpaces+" "+"var="+((FunzioneSingola) funzioneTMP).getVariable().calcola());
} else {
throw new java.lang.RuntimeException("Argomenti mancanti! Sistemare l'equazione!");
}
}
} else if (funzioneTMP instanceof Termine || funzioneTMP instanceof Parentesi) {
if (n < 300) {
//Utils.debug.println(debugSpaces+"•Set variable to number:"+funzioneTMP.calcola());
}
} else {
throw new java.lang.RuntimeException("Tipo sconosciuto");
}
i++;
n++;
}
} while (funzioniOLD.size() != before || fase != "somme");
setVariables(funzioniOLD);
dsl = debugSpaces.length();debugSpaces = "";for (int i = 0; i < dsl-2; i++) {debugSpaces += " ";}
Utils.debug.println(debugSpaces+"•Finished correcting classes.");
Utils.debug.println(debugSpaces+"•Result:"+calcola());
}
}
@Override
public String simbolo() {
return "Parentesi";
}
@Override
public Termine calcola() throws Errore {
Termine result = new Termine(BigSurdVec.ZERO);
for (Funzione f : variables) {
result = result.add(f.calcola());
}
return result;
}
}

18
src/org/warpgate/pi/calculator/Potenza.java Normal file
View File

@ -0,0 +1,18 @@
package org.warpgate.pi.calculator;
public class Potenza extends FunzioneDueValori {
public Potenza(Funzione value1, Funzione value2) {
super(value1, value2);
}
@Override
public String simbolo() {
return Simboli.POTENZA;
}
@Override
public Termine calcola() throws NumberFormatException, Errore {
return getVariable1().calcola().pow(getVariable2().calcola());
}
}

22
src/org/warpgate/pi/calculator/Radice.java Normal file
View File

@ -0,0 +1,22 @@
package org.warpgate.pi.calculator;
import org.nevec.rjm.BigSurdVec;
public class Radice extends FunzioneDueValori {
public Radice(Funzione value1, Funzione value2) {
super(value1, value2);
}
@Override
public String simbolo() {
return Simboli.NTH_ROOT;
}
@Override
public Termine calcola() throws NumberFormatException, Errore {
Termine exponent = new Termine(BigSurdVec.ONE);
exponent = exponent.divide(getVariable1().calcola());
return getVariable2().calcola().pow(exponent);
}
}

20
src/org/warpgate/pi/calculator/RadiceQuadrata.java Normal file
View File

@ -0,0 +1,20 @@
package org.warpgate.pi.calculator;
import org.nevec.rjm.Rational;
public class RadiceQuadrata extends FunzioneSingola {
public RadiceQuadrata(Funzione value) {
super(value);
}
@Override
public String simbolo() {
return Simboli.SQUARE_ROOT;
}
@Override
public Termine calcola() throws NumberFormatException, Errore {
return getVariable().calcola().pow(new Termine(Rational.ONE.divide(new Rational(2,1))));
}
}

11
src/org/warpgate/pi/calculator/RisultatoEquazione.java Normal file
View File

@ -0,0 +1,11 @@
package org.warpgate.pi.calculator;
public class RisultatoEquazione {
public boolean isAnEquation = false;
public Termine LR = new Termine("0");
public RisultatoEquazione(Termine LR, boolean isAnEquation) {
this.LR = LR;
this.isAnEquation = isAnEquation;
}
}

34
src/org/warpgate/pi/calculator/Simboli.java Normal file
View File

@ -0,0 +1,34 @@
package org.warpgate.pi.calculator;
import static org.warpgate.pi.calculator.Utils.concat;
public class Simboli {
public static final String SUM = "+";
public static final String MULTIPLICATION = "*";
public static final String DIVISION = "/";
public static final String NTH_ROOT = "√";
public static final String SQUARE_ROOT = "Ⓐ";
public static final String PARENTHESIS_OPEN = "(";
public static final String PARENTHESIS_CLOSE = ")";
public static final String POTENZA = "^";
public static final String[] funzioni() {
return concat(funzioniNSN(), funzioniSN());
}
public static final String[] funzioniNSN() {
return new String[]{NTH_ROOT, POTENZA};
}
public static final String[] funzioniSN() {
return new String[]{SQUARE_ROOT};
}
public static final String[] segni(boolean withMultiplication) {
String[] ret = new String[]{SUM, DIVISION};
if (withMultiplication) {
ret = Utils.add(ret, MULTIPLICATION);
}
return ret;
}
public static final String[] parentesi() {
return new String[]{PARENTHESIS_OPEN, PARENTHESIS_CLOSE};
}
}

19
src/org/warpgate/pi/calculator/Somma.java Normal file
View File

@ -0,0 +1,19 @@
package org.warpgate.pi.calculator;
public class Somma extends FunzioneDueValori {
public Somma(Funzione value1, Funzione value2) {
super(value1, value2);
}
@Override
public String simbolo() {
return Simboli.SUM;
}
@Override
public Termine calcola() throws Errore {
return getVariable1().calcola().add(getVariable2().calcola());
}
}

143
src/org/warpgate/pi/calculator/Termine.java Normal file
View File

@ -0,0 +1,143 @@
package org.warpgate.pi.calculator;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.util.ArrayList;
import org.nevec.rjm.BigDecimalMath;
import org.nevec.rjm.BigSurd;
import org.nevec.rjm.BigSurdVec;
import org.nevec.rjm.Rational;
public class Termine implements Funzione {
protected BigSurdVec term = BigSurdVec.ZERO;
protected ArrayList<VariabileEdEsponente> variables = new ArrayList<VariabileEdEsponente>();
public Termine(BigSurdVec val) {
term = val;
}
public Termine(String s) {
term = new BigSurdVec(new BigSurd(Utils.getRational(s), Rational.ONE));
}
public Termine(Rational r) {
term = new BigSurdVec(new BigSurd(r, Rational.ONE));
}
public Termine(BigInteger r) {
term = new BigSurdVec(new BigSurd(new Rational(r, BigInteger.ONE), Rational.ONE));
}
public Termine(BigDecimal r) {
term = new BigSurdVec(new BigSurd(Utils.getRational(r), Rational.ONE));
}
public BigSurdVec getTerm() {
return term;
}
public void setTerm(BigSurdVec val) {
term = val;
}
public ArrayList<VariabileEdEsponente> getVariables() {
return variables;
}
public void setVariables(ArrayList<VariabileEdEsponente> variables) {
this.variables = variables;
}
@Override
public Termine calcola() {
return this;
}
@Override
public String simbolo() {
return null;
}
public Termine add(Termine f) throws Errore {
Termine ret = new Termine(getTerm().add(f.getTerm()));
ret.setVariables(getVariables());
if (f.getVariables().size() > 0) {
System.err.println("Moltiplicazioni con variabili non implementato");
throw new Errore(Errori.NOT_IMPLEMENTED);
}
return ret;
}
public Termine multiply(Termine f) throws Errore {
Termine ret = new Termine(getTerm().multiply(f.getTerm()));
ret.setVariables(getVariables());
if (f.getVariables().size() > 0) {
System.err.println("Moltiplicazioni con variabili non implementato");
throw new Errore(Errori.NOT_IMPLEMENTED);
}
return ret;
}
public Termine divide(Termine f) throws Errore {
Termine ret = new Termine(getTerm().divide(f.getTerm()));
ret.setVariables(getVariables());
if (f.getVariables().size() > 0) {
System.err.println("Divisioni con variabili non implementato");
throw new Errore(Errori.NOT_IMPLEMENTED);
}
return ret;
}
public Termine pow(Termine f) throws Errore {
Termine ret = new Termine(BigSurdVec.ONE);
if (f.getTerm().isBigInteger()) {
for (BigInteger i = BigInteger.ZERO; i.compareTo(f.getTerm().toBigInteger()) < 0; i = i.add(BigInteger.ONE)) {
ret = ret.multiply(new Termine(getTerm()));
ret.setVariables(getVariables());
if (f.getVariables().size() > 0) {
System.err.println("Potenze con variabili non implementato");
throw new Errore(Errori.NOT_IMPLEMENTED);
}
}
} else if (getTerm().isRational() && f.getTerm().isRational() && f.getTerm().toRational().denom().compareTo(BigInteger.ONE) == 0) {
ret = new Termine(getTerm().toRational().pow(f.getTerm().toRational()));
ret.setVariables(getVariables());
if (f.getVariables().size() > 0) {
System.err.println("Potenze con variabili non implementato");
throw new Errore(Errori.NOT_IMPLEMENTED);
}
} else if (getTerm().isRational() && f.getTerm().isRational() && f.getTerm().toRational().compareTo(Rational.HALF) == 0) {
//Rational originalExponent = f.getTerm().toRational();
//Rational rootExponent = new Rational(originalExponent.denom(), originalExponent.numer());
Rational numberToRoot = getTerm().toRational();
ret = new Termine(new BigSurdVec(new BigSurd(Rational.ONE, numberToRoot)));
ret.setVariables(getVariables());
if (f.getVariables().size() > 0) {
System.err.println("Potenze con variabili non implementato");
throw new Errore(Errori.NOT_IMPLEMENTED);
}
} else {
ret = new Termine(BigDecimalMath.pow(getTerm().BigDecimalValue(new MathContext(Utils.scale, Utils.scaleMode2)), f.getTerm().BigDecimalValue(new MathContext(Utils.scale, Utils.scaleMode2))));
ret.setVariables(getVariables());
if (f.getVariables().size() > 0) {
System.err.println("Potenze con variabili non implementato");
throw new Errore(Errori.NOT_IMPLEMENTED);
}
}
return ret;
}
@Override
public String toString() {
if (getTerm().isBigInteger()) {
return getTerm().toBigInteger().toString();
} else if (getTerm().isRational()) {
return getTerm().toRational().toString();
} else {
return getTerm().toFancyString();
}
}
}

156
src/org/warpgate/pi/calculator/Utils.java Normal file
View File

@ -0,0 +1,156 @@
package org.warpgate.pi.calculator;
import java.io.StringWriter;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.RoundingMode;
import java.util.ArrayList;
import org.nevec.rjm.BigDecimalMath;
import org.nevec.rjm.Rational;
public class Utils {
public static final int scale = 130;
public static final int scaleMode = BigDecimal.ROUND_HALF_UP;
public static final RoundingMode scaleMode2 = RoundingMode.HALF_UP;
public static DebugStream debug = new DebugStream();
public static boolean debugOn;
public static final class DebugStream extends StringWriter {
public void println(String str) {
if (debugOn) {
System.err.println(str);
}
}
int before = 0;
boolean due = false;
}
public static boolean isInArray(String ch, String[] a) {
boolean contains = false;
for (String c : a) {
if (c.equals(ch)) {
contains = true;
break;
}
}
return contains;
}
public static String ArrayToRegex(String[] array) {
String regex = null;
for(String symbol : array) {
if (regex != null) {
regex+="|\\"+symbol;
} else {
regex = "\\"+symbol;
}
}
return regex;
}
public static String[] concat(String[] a, String[] b) {
int aLen = a.length;
int bLen = b.length;
String[] c= new String[aLen+bLen];
System.arraycopy(a, 0, c, 0, aLen);
System.arraycopy(b, 0, c, aLen, bLen);
return c;
}
public static String[] add(String[] a, String b) {
int aLen = a.length;
String[] c= new String[aLen+1];
System.arraycopy(a, 0, c, 0, aLen);
c[aLen] = b;
return c;
}
public static boolean ciSonoSoloNumeriESomme(ArrayList<Funzione> fl) {
for (int i = 0; i < fl.size(); i++) {
if (!(fl.get(i) instanceof Termine || fl.get(i) instanceof Somma || fl.get(i) instanceof Parentesi)) {
return false;
}
}
return true;
}
public static boolean ciSonoFunzioniSNnonImpostate(ArrayList<Funzione> fl) {
for (int i = 0; i < fl.size(); i++) {
if (fl.get(i) instanceof FunzioneSingola) {
if (((FunzioneSingola)fl.get(i)).variable == null) {
return true;
}
}
}
return false;
}
public static boolean ciSonoAltreFunzioni(ArrayList<Funzione> fl) {
for (int i = 0; i < fl.size(); i++) {
if (!(fl.get(i) instanceof Termine || fl.get(i) instanceof Somma || fl.get(i) instanceof Parentesi || fl.get(i) instanceof FunzioneSingola || fl.get(i) instanceof Moltiplicazione || fl.get(i) instanceof Divisione)) {
return true;
}
}
return false;
}
public static Rational getRational(BigDecimal str) {
return getRational(str.toString());
}
public static Rational getRational(String str) {
try {
return new Rational(str);
} catch (NumberFormatException ex) {
long bits = Double.doubleToLongBits(Double.parseDouble(str));
long sign = bits >>> 63;
long exponent = ((bits >>> 52) ^ (sign << 11)) - 1023;
long fraction = bits << 12; // bits are "reversed" but that's not a problem
long a = 1L;
long b = 1L;
for (int i = 63; i >= 12; i--) {
a = a * 2 + ((fraction >>> i) & 1);
b *= 2;
}
if (exponent > 0)
a *= 1 << exponent;
else
b *= 1 << -exponent;
if (sign == 1)
a *= -1;
if (b == 0) {
a = 0;
b = 1;
}
return new Rational(new BigInteger(a+""),new BigInteger(b+""));
}
}
public static BigDecimal rationalToIrrationalString(Rational r) {
return BigDecimalMath.divideRound(new BigDecimal(r.numer()).setScale(Utils.scale, Utils.scaleMode), new BigDecimal(r.denom()).setScale(Utils.scale, Utils.scaleMode));
}
public static boolean variabiliUguali(ArrayList<VariabileEdEsponente> variables, ArrayList<VariabileEdEsponente> variables2) {
if (variables.size() != variables2.size()) {
return false;
} else {
for (VariabileEdEsponente v : variables) {
if (!variables2.contains(v)) {
return false;
}
}
return true;
}
}
}

25
src/org/warpgate/pi/calculator/VariabileEdEsponente.java Normal file
View File

@ -0,0 +1,25 @@
package org.warpgate.pi.calculator;
import java.math.BigInteger;
public class VariabileEdEsponente {
public char simbolo = 'X';
public BigInteger esponente = BigInteger.ONE;
public VariabileEdEsponente(char simbolo, BigInteger esponente) {
this.simbolo = simbolo;
this.esponente = esponente;
}
@Override
public boolean equals(Object o) {
if (o instanceof VariabileEdEsponente) {
if (this.simbolo == ((VariabileEdEsponente) o).simbolo) {
if (this.esponente == ((VariabileEdEsponente) o).esponente) {
return true;
}
}
}
return false;
}
}