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

Removed useless files

Browse Source
This commit is contained in:
Cavallium 2018-10-15 19:37:56 +02:00
parent 971759a6e7
commit 76985314d4
38 changed files with 124 additions and 10508 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitmodules vendored
View File

@ -1,3 +1,6 @@
[submodule "Flow"] [submodule "Flow"]
path = Flow path = Flow
url = https://github.com/Cavallium/Flow url = https://github.com/Cavallium/Flow
[submodule "bigdecimal-math"]
path = bigdecimal-math
url = https://github.com/Cavallium/bigdecimal-math

2
LICENSE
View File

@ -186,7 +186,7 @@
same "printed page" as the copyright notice for easier same "printed page" as the copyright notice for easier
identification within third-party archives. identification within third-party archives.
Copyright [yyyy] [name of copyright owner] Copyright 2018 Andrea Cavalli
Licensed under the Apache License, Version 2.0 (the "License"); Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License. you may not use this file except in compliance with the License.

33
MmapByteBuffer.java
View File

@ -1,33 +0,0 @@
package it.cavallium.warppi;
import java.nio.ByteBuffer;
public class MmapByteBuffer {
private int fd;
private int address;
private int length;
private ByteBuffer buffer;
public MmapByteBuffer(int fd, int address, int length, ByteBuffer buffer) {
this.fd = fd;
this.address = address;
this.length = length;
this.buffer = buffer;
}
public int getFd() {
return fd;
}
public int getAddress() {
return address;
}
public int getLength() {
return length;
}
public ByteBuffer getBuffer() {
return buffer;
}
}

1
bigdecimal-math Submodule

@ -0,0 +1 @@
Subproject commit 92a67122350a38af726beae58d2582fb867a4c61

14
core/pom.xml
View File

@ -6,7 +6,7 @@
<parent> <parent>
<groupId>it.cavallium</groupId> <groupId>it.cavallium</groupId>
<artifactId>warppi</artifactId> <artifactId>warppi</artifactId>
<version>0.9.0a2</version> <version>0.9.0a2</version>
</parent> </parent>
<artifactId>warppi-core</artifactId> <artifactId>warppi-core</artifactId>
@ -17,7 +17,17 @@
<dependency> <dependency>
<groupId>it.cavallium</groupId> <groupId>it.cavallium</groupId>
<artifactId>warppi-flow</artifactId> <artifactId>warppi-flow</artifactId>
<version>${project.version}</version> <version>0.9.0a2</version>
</dependency>
<dependency>
<groupId>it.cavallium</groupId>
<artifactId>warppi-util</artifactId>
<version>0.9.0a2</version>
</dependency>
<dependency>
<groupId>it.cavallium</groupId>
<artifactId>bigdecimal-math</artifactId>
<version>0.9.0a2</version>
</dependency> </dependency>
<dependency> <dependency>
<groupId>it.unimi.dsi</groupId> <groupId>it.unimi.dsi</groupId>

2
core/src/main/java/it/cavallium/warppi/extra/tetris/TetrisScreen.java
View File

@ -2,8 +2,6 @@ package it.cavallium.warppi.extra.tetris;
import java.io.IOException; import java.io.IOException;
import org.nevec.rjm.Wigner3j;
import it.cavallium.warppi.Engine; import it.cavallium.warppi.Engine;
import it.cavallium.warppi.StaticVars; import it.cavallium.warppi.StaticVars;
import it.cavallium.warppi.device.Keyboard; import it.cavallium.warppi.device.Keyboard;

102
core/src/main/java/org/nevec/rjm/Bernoulli.java
View File

@ -1,102 +0,0 @@
package org.nevec.rjm;
import java.math.BigInteger;
import java.util.Vector;
import it.cavallium.warppi.util.Error;
/**
* 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<>();
public Bernoulli() {
if (Bernoulli.a.size() == 0) {
Bernoulli.a.add(Rational.ONE);
Bernoulli.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 < Bernoulli.a.size()) {
Bernoulli.a.set(nindx, value);
} else {
while (Bernoulli.a.size() < nindx) {
Bernoulli.a.add(Rational.ZERO);
}
Bernoulli.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
* @throws Error
*/
public Rational at(final int n) throws Error {
if (n == 1) {
return new Rational(-1, 2);
} else if (n % 2 != 0) {
return Rational.ZERO;
} else {
final int nindx = n / 2;
if (Bernoulli.a.size() <= nindx) {
for (int i = 2 * Bernoulli.a.size(); i <= n; i += 2) {
set(i, doubleSum(i));
}
}
return Bernoulli.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(final int n) throws Error {
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++) {
final 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 */

220
core/src/main/java/org/nevec/rjm/BigComplex.java
View File

@ -1,220 +0,0 @@
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(final BigDecimal x, final BigDecimal y) {
re = x;
im = y;
}
/**
* ctor with real part.
*
* @param x
* real part.
* The imaginary part is set to zero.
*/
public BigComplex(final BigDecimal x) {
re = x;
im = BigDecimal.ZERO;
}
/**
* ctor with real and imaginary parts
*
* @param x
* real part
* @param y
* imaginary part
*/
public BigComplex(final double x, final 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, final 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(final 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="https://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(final 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(final 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(final BigComplex oth, final 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.
*/
@Override
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(final MathContext mc) {
return "(" + re.round(mc).toString() + "," + im.round(mc).toString() + ")";
}
} /* BigComplex */

3038
core/src/main/java/org/nevec/rjm/BigDecimalMath.java
View File

File diff suppressed because it is too large Load Diff

644
core/src/main/java/org/nevec/rjm/BigIntegerMath.java
View File

@ -1,644 +0,0 @@
package org.nevec.rjm;
import java.math.BigInteger;
import java.util.Vector;
import it.cavallium.warppi.util.Error;
/**
* 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);
final 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(final int n) {
return new Ifactor(Math.abs(n)).sigma().n;
}
/**
* Compute the list of positive divisors.<br>
* <b>Deprecated: This function is extremely slow and inefficient!</b>
*
* @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
*/
@Deprecated
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(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(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
*/
final 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));
}
final 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));
}
final 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[][] = BigIntegerMath.minor(A, r, 0);
final BigInteger m = A[r][0].multiply(BigIntegerMath.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
* @throws Error
*/
static public Rational[] solve(final BigInteger[][] A, final BigInteger[] rhs) throws ArithmeticException, Error {
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
*/
final 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++) {
final BigInteger tmp = A[c][cpr];
A[c][cpr] = A[r][cpr];
A[r][cpr] = tmp;
}
final 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++) {
final BigInteger tmp = A[c][c].multiply(A[r][cpr]).subtract(A[c][cpr].multiply(A[r][c]));
A[r][cpr] = tmp;
}
final 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) {
final 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
* @throws Error
* @see <a href="https://doi.org/10.1080/01630568908816313">P. L. Butzer
* et al, Num. Funct. Anal. Opt. 10 (5)( 1989) 419-488</a>
*/
static public Rational centrlFactNumt(final int n, final int k) throws Error {
if (k > n || k < 0 || k % 2 != n % 2) {
return Rational.ZERO;
} else if (k == n) {
return Rational.ONE;
} else {
/* Proposition 6.2.6 */
final Factorial f = new Factorial();
Rational jsum = new Rational(0, 1);
final 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(BigIntegerMath.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(BigIntegerMath.binomial(2 * kprime, kprime - j));
if (j % 2 != 0) {
jsum = jsum.subtract(nusum);
} else {
jsum = jsum.add(nusum);
}
}
return jsum.multiply(k).multiply(BigIntegerMath.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="https://doi.org/10.1080/01630568908816313">P. L. Butzer
* et al, Num. Funct. Anal. Opt. 10 (5)( 1989) 419-488</a>
*/
static public Rational centrlFactNumT(final int n, final 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 BigIntegerMath.centrlFactNumT(n - 2, k - 2).add(BigIntegerMath.centrlFactNumT(n - 2, k).multiply(new Rational(k * k, 4)));
}
} /* CentralFactNumT */
} /* BigIntegerMath */

764
core/src/main/java/org/nevec/rjm/BigIntegerPoly.java
View File

@ -1,764 +0,0 @@
package org.nevec.rjm;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.util.Scanner;
import java.util.Vector;
import it.cavallium.warppi.util.Error;
/**
* 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<>();
}
/**
* 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<>();
final 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 (final BigInteger element : c) {
a.add(element.add(BigInteger.ZERO));
}
simplify();
} /* ctor */
/**
* Create a copy of this.
*
* @since 2010-08-27
*/
@Override
public BigIntegerPoly clone() {
return new BigIntegerPoly(a);
} /* clone */
/**
* Translate into a RatPoly copy.
*
* @since 2012-03-02
*/
public RatPoly toRatPoly() {
final 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(final 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) {
final 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) {
final 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) {
final 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) {
final 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++) {
final 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) {
final 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++) {
final 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) {
final BigIntegerPoly[] ret = new BigIntegerPoly[2];
/*
* remove any high-order zeros. note that the clone() operation calls
* simplify().
*/
final BigIntegerPoly valSimpl = val.clone();
final 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();
final 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 {
final 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
*/
@Override
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 {
final 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(final int newdeg) {
final 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(final int maxdeg) {
final 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 */
final Vector<BigInteger> res = new Vector<>();
/*
* 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)
*/
final int l = ldegree();
if (a.elementAt(l).compareTo(BigInteger.ZERO) != 0) {
@SuppressWarnings("deprecation")
final 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.
* @throws Error
* @since 2012-03-01
*/
protected Vector<BigIntegerPoly> i2roots() throws Error {
/*
* The vector of the factors to be returned
*/
final Vector<BigIntegerPoly> res = new Vector<>();
if (degree() < 2) {
return res;
}
final BigInteger bsco = a.firstElement().abs();
@SuppressWarnings("deprecation")
final Vector<BigInteger> b = BigIntegerMath.divisors(bsco);
final BigInteger csco = a.lastElement().abs();
@SuppressWarnings("deprecation")
final 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
*/
final RatPoly thisDec = toRatPoly();
final 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 (final BigComplex z : roo) {
for (final BigInteger bco : b) {
for (final 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) {
final 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
*/
final BigInteger a = negA.negate().add(half).toBigInteger();
/*
* test the polynomial remainder. if zero, add the
* term
* to the results.
*/
final BigIntegerPoly dtst = new BigIntegerPoly("" + bcosig + "," + a + "," + cco);
try {
final BigIntegerPoly[] rm = divideAndRemainder(dtst);
if (rm[1].isZero()) {
res.add(dtst);
}
} catch (final 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.
* @throws Error
* @since 2012-03-01
*/
public Vector<BigIntegerPoly> ifactor() throws Error {
/*
* this ought be entirely rewritten in terms of the LLL algorithm
*/
final Vector<BigIntegerPoly> fac = new Vector<>();
/* collect integer roots (polynomial factors of degree 1) */
final Vector<BigInteger> r = iroots();
BigIntegerPoly[] res = new BigIntegerPoly[2];
res[0] = this;
for (final BigInteger i : r) {
final int deg = rootDeg(i);
/* construct the factor x-i */
final 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
*/
final Vector<BigIntegerPoly> pol2 = i2roots();
for (final 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 {
final BigIntegerPoly[] dtst = res[0].divideAndRemainder(i);
if (dtst[1].isZero()) {
fac.add(i);
res = dtst;
} else {
break;
}
} catch (final 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 */

567
core/src/main/java/org/nevec/rjm/BigSurd.java
View File

@ -1,567 +0,0 @@
package org.nevec.rjm;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.security.ProviderException;
import it.cavallium.warppi.util.Error;
import it.cavallium.warppi.util.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(final Rational a, final Rational b) {
pref = a;
/*
* reject attempts to use a negative b
*/
if (b.signum() < 0) {
throw new ProviderException("Not implemented: imaginary surds");
}
disc = b;
try {
normalize();
normalizeG();
} catch (final Error e) {
e.printStackTrace();
}
}
/**
* 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(final int a, final 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(final BigInteger a) {
this(Rational.ONE, new Rational(a, BigInteger.ONE));
}
public BigSurd(final Rational a) {
this(Rational.ONE, a);
}
/**
* Create a deep copy.
*
* @since 2011-02-12
*/
@Override
public BigSurd clone() {
final Rational fclon = pref.clone();
final 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.
*/
final 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) {
final 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
* @throws Error
* @since 2011-02-12
*/
public BigSurd divide(final BigSurd val) throws Error {
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() {
final BigSurd bs = this;
final 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 {
final 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
* @throws Error
* @since 2011-02-12
*/
public BigSurd divide(final BigInteger val) throws Error {
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
* @throws Error
* @since 2011-02-12
*/
public BigSurd divide(final int val) throws Error {
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()
*/
final 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
*/
@Override
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;
} 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
*/
@Override
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.
*/
final Rational p2 = pref.pow(2).multiply(disc);
System.out.println("dv sq " + p2.toString());
final 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.
*
* @throws Error
* @since 2011-02-12
*/
protected void normalize() throws Error {
/*
* 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
*/
final 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);
final 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
*
* @throws Error
* @since 2011-02-12
*/
protected void normalizeG() throws Error {
/*
* 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(final 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
*/
final 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.
*/
final BigDecimal res = sqrdis.multiply(pref.BigDecimalValue(mc));
return BigDecimalMath.scalePrec(res, mc);
} /* BigDecimalValue */
} /* BigSurd */

650
core/src/main/java/org/nevec/rjm/BigSurdVec.java
View File

@ -1,650 +0,0 @@
package org.nevec.rjm;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.util.Vector;
import it.cavallium.warppi.util.Error;
import it.cavallium.warppi.util.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<>();
} /* ctor */
/**
* ctor given the value of a BigSurd.
*
* @param a
* The value to be represented by this vector.
* @since 2012-02-15
*/
public BigSurdVec(final BigSurd a) {
terms = new Vector<>(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(final BigSurd a, final BigSurd b) {
terms = new Vector<>(2);
terms.add(a);
terms.add(b);
try {
normalize();
} catch (final Error e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
} /* 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).
*
* @throws Error
*
* @since 2012-02-15
*/
protected void normalize() throws Error {
/*
* nothing to be done if at most one term
*/
if (terms.size() <= 1) {
return;
}
final Vector<BigSurd> newter = new Vector<>();
newter.add(terms.firstElement());
/*
* add j-th element to the existing vector and combine were possible
*/
for (int j = 1; j < terms.size(); j++) {
final 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.
*/
final BigSurd r = todo.divide(v);
if (r.isRational()) {
/* compute r+1 */
final 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
*/
@Override
public int compareTo(final BigSurdVec oth) {
BigSurdVec diff;
try {
diff = this.subtract(oth);
return diff.signum();
} catch (final Error e) {
// TODO Auto-generated catch block
e.printStackTrace();
return 0;
}
} /* compareTo */
/**
* Sign function. Returns -1, 0 or +1 depending on whether this is smaller,
* equal to or larger than zero.
*
* @return 0 or +-1.
* @throws Error
* @since 2012-02-15
*/
public int signum() throws Error {
/*
* 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].
*/
final 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(final 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].
*/
final BigDecimal[] res = new BigDecimal[2];
res[0] = BigDecimal.ZERO;
for (int addpr = 1;; addpr += 3) {
final MathContext locmc = new MathContext(mc.getPrecision() + addpr, mc.getRoundingMode());
res[1] = BigDecimal.ZERO;
for (final BigSurd j : terms) {
res[1] = BigDecimalMath.addRound(res[1], j.BigDecimalValue(locmc));
}
if (addpr > 1) {
final BigDecimal err = res[1].subtract(res[0]).abs();
final 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() {
final 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() {
final 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.
* @throws Error
*/
public BigSurdVec add(final BigSurdVec val) throws Error {
final BigSurdVec sum = new BigSurdVec();
/*
* concatenate the vectors and eliminate common overlaps
*/
for (final BigSurd term : terms) {
if (term.compareTo(BigSurd.ZERO) != 0) {
sum.terms.add(term);
}
}
for (final 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.
* @throws Error
*/
public BigSurdVec add(final BigSurd val) throws Error {
final 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.
* @throws Error
*/
public BigSurdVec subtract(final BigSurdVec val) throws Error {
final BigSurdVec sum = new BigSurdVec();
/*
* concatenate the vectors and eliminate common overlaps
*/
sum.terms.addAll(terms);
for (final 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.
* @throws Error
*/
public BigSurdVec subtract(final BigSurd val) throws Error {
final 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
*/
final BigSurdVec resul = new BigSurdVec();
for (final 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.
* @throws Error
* @since 2012-02-15
*/
public BigSurdVec sqr() throws Error {
/*
* Binomial expansion. First the sum of the terms squared, then 2 times
* the mixed products.
*/
final 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.
* @throws Error
* @since 2011-02-12
*/
public BigSurdVec multiply(final BigSurd val) throws Error {
final BigSurdVec resul = new BigSurdVec();
for (final BigSurd s : terms) {
resul.terms.add(s.multiply(val));
}
resul.normalize();
return resul;
} /* multiply */
public BigSurdVec multiply(final BigSurdVec val) throws Error {
BigSurdVec resul = new BigSurdVec();
for (final BigSurd s : terms) {
resul.terms.add(s);
}
for (final BigSurd s : val.terms) {
resul = resul.multiply(s);
}
return resul;
} /* multiply */
public BigSurdVec divide(final BigSurd val) throws Error {
final BigSurdVec resul = new BigSurdVec();
for (final BigSurd s : terms) {
resul.terms.add(s.divide(val));
}
resul.normalize();
return resul;
} /* multiply */
public BigSurdVec divide(final BigSurdVec val) throws Error {
BigSurdVec resul = new BigSurdVec();
resul.terms = terms;
for (final 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 (final 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 (final 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 (final 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 (final 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 (final 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
*/
@Override
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++) {
final 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++) {
final 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 {
final BigInteger numerator = bs.pref.multiply(denominator).numer();
if (numerator.compareTo(BigInteger.ONE) != 0) {
s += numerator.toString();
s += "*";
// s += "("; Radice quadrata. non servono le parentesi.
}
s += "Ⓐ";
if (bs.disc.isInteger()) {
s += bs.disc.toString();
} else {
s += "(" + bs.disc.toString() + ")";
}
if (numerator.compareTo(BigInteger.ONE) != 0) {
// s += ")"; Radice quadrata. non servono le parentesi.
}
}
}
if (denominator.compareTo(BigInteger.ONE) != 0) {
s += ")";
s += "/";
s += denominator;
}
return s;
}
}
} /* BigSurdVec */

73
core/src/main/java/org/nevec/rjm/Euler.java
View File

@ -1,73 +0,0 @@
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<>();
/**
* Ctor(). Fill the hash list initially with E_0 to E_3.
*/
public Euler() {
if (Euler.a.size() == 0) {
Euler.a.add(BigInteger.ONE);
Euler.a.add(BigInteger.ONE);
Euler.a.add(new BigInteger("5"));
Euler.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 >= Euler.a.size()) {
BigInteger val = BigInteger.ZERO;
boolean sigPos = true;
final int thisn = Euler.a.size();
for (int i = thisn - 1; i > 0; i--) {
BigInteger f = new BigInteger("" + Euler.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);
}
Euler.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(final int n) {
set(n);
return Euler.a.elementAt(n);
}
} /* Euler */

66
core/src/main/java/org/nevec/rjm/EulerPhi.java
View File

@ -1,66 +0,0 @@
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(final 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(final BigInteger n) {
if (n.compareTo(BigInteger.ZERO) <= 0) {
throw new ArithmeticException("negative argument " + n + " of EulerPhi");
}
final Ifactor prFact = new Ifactor(n);
BigInteger phi = n;
if (n.compareTo(BigInteger.ONE) > 0) {
for (int i = 0; i < prFact.primeexp.size(); i += 2) {
final BigInteger p = new BigInteger(prFact.primeexp.elementAt(i).toString());
final 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(final String[] args) throws ArithmeticException {
final EulerPhi a = new EulerPhi();
final int n = new Integer(args[0]).intValue();
System.out.println("phi(" + n + ") = " + a.at(n));
}
} /* EulerPhi */

79
core/src/main/java/org/nevec/rjm/Factorial.java
View File

@ -1,79 +0,0 @@
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<>();
/**
* ctor().
* Initialize the vector of the factorials with 0!=1 and 1!=1.
*/
public Factorial() {
if (Factorial.a.size() == 0) {
Factorial.a.add(Ifactor.ONE);
Factorial.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(final int n) {
/*
* extend the internal list if needed.
*/
growto(n);
return Factorial.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(final int n) {
/*
* extend the internal list if needed.
*/
growto(n);
return Factorial.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(final int n) {
/*
* extend the internal list if needed. Size to be 2 for n<=1, 3 for n<=2
* etc.
*/
while (Factorial.a.size() <= n) {
final int lastn = Factorial.a.size() - 1;
final Ifactor nextn = new Ifactor(lastn + 1);
Factorial.a.add(Factorial.a.elementAt(lastn).multiply(nextn));
}
} /* growto */
} /* Factorial */

43
core/src/main/java/org/nevec/rjm/Harmonic.java
View File

@ -1,43 +0,0 @@
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(final 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 */

835
core/src/main/java/org/nevec/rjm/Ifactor.java
View File

@ -1,835 +0,0 @@
package org.nevec.rjm;
import java.math.BigInteger;
import java.util.Collections;
import java.util.Vector;
import it.cavallium.warppi.util.Error;
/**
* 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<>();
if (number > 1) {
int primindx = 0;
final 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
*/
final 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<>();
if (number.compareTo(BigInteger.ONE) == 0) {
primeexp.add(new Integer(1));
primeexp.add(new Integer(0));
} else {
int primindx = 0;
final Prime primes = new Prime();
/*
* Test for division against all primes.
*/
while (number.compareTo(BigInteger.ONE) == 1) {
int ex = 0;
final 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(final Vector<Integer> pows) {
primeexp = new Vector<>(2 * pows.size());
if (pows.size() > 0) {
n = BigInteger.ONE;
final 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++) {
final 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(final Ifactor oth) {
n = oth.n;
primeexp = oth.primeexp;
} /* Ifactor */
/**
* Deep copy.
*
* @since 2009-08-14
*/
@Override
public Ifactor clone() {
/*
* Line not used:
*
* Vector<Integer> p = (Vector<Integer>)primeexp.clone();
*
*/
final 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, final int type) {
final 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.
*/
final Ifactor pows = new Ifactor(0);
for (int i = 0; i < primeexp.size(); i += 2) {
final Integer p = primeexp.elementAt(i);
final 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.
* @throws Error
* @since 2009-05-18
*/
public Rational root(final int r) throws ArithmeticException, Error {
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 {
final 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).
*/
final 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.
*/
final Vector<BigInteger> d = new Vector<>();
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) */
final Ifactor dp = dropPrime();
/* get ez */
final int ez = primeexp.lastElement().intValue();
final 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(final 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 Ifactor.ONE;
} else if (n.compareTo(BigInteger.ZERO) == 0) {
return Ifactor.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) {
final int ex = primeexp.elementAt(i + 1).intValue();
if (k == 0) {
resul = resul.multiply(ex + 1);
} else {
final Integer p = primeexp.elementAt(i);
final BigInteger num = new BigInteger(p.toString()).pow(k * (ex + 1)).subtract(BigInteger.ONE);
final 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.
*/
final 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.
*/
final 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));
final BigInteger p = new BigInteger(primeexp.elementAt(i).toString());
final 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
*/
@Override
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...
*/
@Override
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(final String[] args) throws Exception {
final BigInteger n = new BigInteger(args[0]);
System.out.println(new Ifactor(n));
} /* Ifactor.main */
} /* Ifactor */

88
core/src/main/java/org/nevec/rjm/PartitionsP.java
View File

@ -1,88 +0,0 @@
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<>();
/**
* 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 (PartitionsP.a.size() == 0) {
PartitionsP.a.add(new BigInteger("" + 1));
PartitionsP.a.add(new BigInteger("" + 1));
PartitionsP.a.add(new BigInteger("" + 2));
PartitionsP.a.add(new BigInteger("" + 3));
PartitionsP.a.add(new BigInteger("" + 5));
PartitionsP.a.add(new BigInteger("" + 7));
}
PartitionsP.nMax = new BigInteger("" + (PartitionsP.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(final int i) {
/*
* If the current list is too small, increase in intervals
* of 3 until the list has at least i elements.
*/
while (i > PartitionsP.nMax.intValue()) {
growto(PartitionsP.nMax.add(new BigInteger("" + 3)));
}
return PartitionsP.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(final BigInteger n) {
while (PartitionsP.a.size() <= n.intValue()) {
BigInteger per = new BigInteger("0");
final BigInteger cursiz = new BigInteger("" + PartitionsP.a.size());
for (int k = 0; k < PartitionsP.a.size(); k++) {
final BigInteger tmp = PartitionsP.a.elementAt(k).multiply(BigIntegerMath.sigma(PartitionsP.a.size() - k));
per = per.add(tmp);
}
PartitionsP.a.add(per.divide(cursiz));
}
PartitionsP.nMax = new BigInteger("" + (PartitionsP.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(final String[] args) throws Exception {
final PartitionsP a = new PartitionsP();
final int n = new Integer(args[0]).intValue();
System.out.println("P(" + n + ")=" + a.at(n));
}
}

321
core/src/main/java/org/nevec/rjm/Prime.java
View File

@ -1,321 +0,0 @@
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<>();
/**
* 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 (Prime.a.size() == 0) {
Prime.a.add(new BigInteger("" + 2));
Prime.a.add(new BigInteger("" + 3));
Prime.a.add(new BigInteger("" + 5));
Prime.a.add(new BigInteger("" + 7));
Prime.a.add(new BigInteger("" + 11));
Prime.a.add(new BigInteger("" + 13));
Prime.a.add(new BigInteger("" + 17));
}
Prime.nMax = Prime.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(final BigInteger n) {
/*
* not documented
* return ( n.isProbablePrime() ) ;
*/
switch (millerRabin(n)) {
case -1:
return false;
case 1:
return true;
}
growto(n);
return Prime.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;
} else if (n.compareTo(two) == 0) {
return true;
} 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);
final int s = q.getLowestSetBit();
final 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) {
final int l = n.compareTo(new BigInteger(mr[mrLim]));
if (l < 0) {
break;
} 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(final int i) {
/*
* If the current list is too small, increase in intervals
* of 5 until the list has at least i elements.
*/
while (i >= Prime.a.size()) {
growto(Prime.nMax.add(new BigInteger("" + 5)));
}
return Prime.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(final 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 < Prime.a.size(); i++) {
if (Prime.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(final BigInteger n) {
/* if n <=1, return 2 */
if (n.compareTo(BigInteger.ONE) <= 0) {
return Prime.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 (Prime.a.lastElement().compareTo(n) <= 0) {
growto(Prime.nMax.add(new BigInteger("" + 5)));
}
for (int i = 0; i < Prime.a.size(); i++) {
if (Prime.a.elementAt(i).compareTo(n) == 1) {
return Prime.a.elementAt(i);
}
}
return Prime.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(final 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 (Prime.a.lastElement().compareTo(n) < 0) {
growto(Prime.nMax.add(new BigInteger("" + 5)));
}
for (int i = 0; i < Prime.a.size(); i++) {
if (Prime.a.elementAt(i).compareTo(n) >= 0) {
return Prime.a.elementAt(i - 1);
}
}
return Prime.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(final BigInteger n) {
while (Prime.nMax.compareTo(n) == -1) {
Prime.nMax = Prime.nMax.add(BigInteger.ONE);
boolean isp = true;
for (int p = 0; p < Prime.a.size(); p++) {
/*
* Test the list of known primes only up to sqrt(n)
*/
if (Prime.a.get(p).multiply(Prime.a.get(p)).compareTo(Prime.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 (Prime.nMax.remainder(Prime.a.get(p)).compareTo(BigInteger.ZERO) == 0) {
isp = false;
break;
}
}
if (isp) {
Prime.a.add(Prime.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(final String[] args) throws Exception {
final Prime a = new Prime();
final 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 */

1045
core/src/main/java/org/nevec/rjm/RatPoly.java
View File

File diff suppressed because it is too large Load Diff

900
core/src/main/java/org/nevec/rjm/Rational.java
View File

@ -1,900 +0,0 @@
package org.nevec.rjm;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;
import it.cavallium.warppi.util.Error;
import it.cavallium.warppi.util.Errors;
/**
* 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(final BigInteger a, final BigInteger b) {
this.a = a;
this.b = b;
normalize();
}
/**
* ctor from a numerator.
*
* @param a
* the BigInteger.
*/
public Rational(final 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(final int a, final 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(final 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(final String str) {
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(final String str, final int radix) {
final 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
*/
@Override
public Rational clone() {
/*
* protected access means this does not work return new
* Rational(a.clone(), b.clone()) ;
*/
final BigInteger aclon = new BigInteger("" + a);
final 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) {
final BigInteger num = a.multiply(val.a);
final 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) {
final 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) {
final 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(final int exponent) {
if (exponent == 0) {
return new Rational(1, 1);
}
final BigInteger num = a.pow(Math.abs(exponent));
final 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.
* @throws Error
* @since 2009-05-18
*/
public Rational pow(final BigInteger exponent) throws Error {
/* test for overflow */
if (exponent.compareTo(Rational.MAX_INT) == 1) {
throw new Error(Errors.NUMBER_TOO_LARGE);
}
if (exponent.compareTo(Rational.MIN_INT) == -1) {
throw new Error(Errors.NUMBER_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).
* @throws Error
* @since 2009-05-18
*/
public Rational root(final BigInteger r) throws Error {
/* test for overflow */
if (r.compareTo(Rational.MAX_INT) == 1) {
throw new Error(Errors.NUMBER_TOO_LARGE);
}
if (r.compareTo(Rational.MIN_INT) == -1) {
throw new Error(Errors.NUMBER_TOO_SMALL);
}
final int rthroot = r.intValue();
/* cannot pull root of a negative value with even-valued root */
if (compareTo(Rational.ZERO) == -1 && rthroot % 2 == 0) {
throw new Error(Errors.NEGATIVE_PARAMETER);
}
/*
* extract a sign such that we calculate |n|^(1/r), still r carrying any
* sign
*/
final boolean flipsign = compareTo(Rational.ZERO) == -1 && rthroot % 2 != 0 ? true : false;
/*
* delegate the main work to ifactor#root()
*/
final Ifactor num = new Ifactor(a.abs());
final 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.
* @throws Error
* @since 2009-05-18
*/
public Rational pow(final Rational exponent) throws Error {
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)
*/
final 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
* @throws Error
*/
public Rational divide(final Rational val) throws Error {
if (val.compareTo(Rational.ZERO) == 0) {
throw new Error(Errors.DIVISION_BY_ZERO);
}
final BigInteger num = a.multiply(val.b);
final 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
* @throws Error
*/
public Rational divide(final BigInteger val) throws Error {
if (val.compareTo(BigInteger.ZERO) == 0) {
throw new Error(Errors.DIVISION_BY_ZERO);
}
final 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
* @throws Error
*/
public Rational divide(final int val) throws Error {
if (val == 0) {
throw new Error(Errors.DIVISION_BY_ZERO);
}
final 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(final Rational val) {
final BigInteger num = a.multiply(val.b).add(b.multiply(val.a));
final 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(final BigInteger val) {
final 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(final int val) {
final 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(final Rational val) {
final 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(final BigInteger val) {
final 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(final int val) {
final 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
* @throws Error
*/
public static Rational binomial(final Rational n, final BigInteger m) throws Error {
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
* @throws Error
*/
public static Rational binomial(final Rational n, final int m) throws Error {
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
* @throws Error
*/
public static Rational hankelSymb(final Rational n, final int k) throws Error {
if (k == 0) {
return Rational.ONE;
} else if (k < 0) {
throw new Error(Errors.NEGATIVE_PARAMETER);
}
Rational nkhalf = n.subtract(k).add(Rational.HALF);
nkhalf = nkhalf.Pochhammer(2 * k);
final 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.
*/
@Override
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
*/
@Override
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.
*/
final 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() {
final 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(final MathContext mc) {
/*
* numerator and denominator individually rephrased
*/
final BigDecimal n = new BigDecimal(a);
final 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(final int digits) {
if (b.compareTo(BigInteger.ONE) != 0) {
final MathContext mc = new MathContext(digits, RoundingMode.DOWN);
final 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(final 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(Rational.MAX_INT) <= 0 && a.compareTo(Rational.MIN_INT) >= 0;
} /* Rational.isInteger */
/**
* Conversion to an integer value, if this can be done exactly.
*
* @throws Error
*
* @since 2011-02-13
*/
int intValue() throws Error {
if (!isInteger()) {
throw new Error(Errors.CONVERSION_ERROR);
}
return a.intValue();
}
/**
* Conversion to a BigInteger value, if this can be done exactly.
*
* @throws Error
*
* @since 2012-03-02
*/
BigInteger BigIntegerValue() throws Error {
if (!isBigInteger()) {
throw new Error(Errors.CONVERSION_ERROR);
}
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(Rational.MAX_INT) <= 0 && a.compareTo(Rational.MIN_INT) >= 0 && b.compareTo(Rational.MAX_INT) <= 0 && b.compareTo(Rational.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 (final Rational val : vals) {
l = BigIntegerMath.lcm(l, val.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 */

27
core/src/main/java/org/nevec/rjm/SafeMathContext.java
View File

@ -1,27 +0,0 @@
package org.nevec.rjm;
import java.math.MathContext;
import java.math.RoundingMode;
import it.cavallium.warppi.Engine;
import it.cavallium.warppi.Platform.ConsoleUtils;
public final class SafeMathContext {
public static MathContext newMathContext(int precision) {
if (precision <= 0) {
Engine.getPlatform().getConsoleUtils().out().print(ConsoleUtils.OUTPUTLEVEL_DEBUG_MIN, "Warning! MathContext precision is <= 0 (" + precision + ")");
precision = 1;
}
return new MathContext(precision);
}
public static MathContext newMathContext(int precision, final RoundingMode roundingMode) {
if (precision <= 0) {
Engine.getPlatform().getConsoleUtils().out().print(ConsoleUtils.OUTPUTLEVEL_DEBUG_MIN, "Warning! MathContext precision is <= 0 (" + precision + ")");
precision = 1;
}
return new MathContext(precision, roundingMode);
}
}

628
core/src/main/java/org/nevec/rjm/Wigner3j.java
View File

@ -1,628 +0,0 @@
package org.nevec.rjm;
import java.math.BigInteger;
import java.util.Scanner;
import it.cavallium.warppi.util.Error;
/**
* 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
* @throws Error
*/
static public void main(final String args[]) throws Error {
if (args[0].compareTo("6j") == 0) {
try {
final String m1 = "6";
final String t1 = "1 2 -3 -1 5 6";
final String t2 = "4 -5 3 -4 -2 -6";
String j = "";
for (int i = 1; i <= 6; i++) {
j += args[i] + " ";
}
final BigSurdVec w = Wigner3j.wigner3j(m1, t1, t2, j);
System.out.println(w.toString());
} catch (final Exception e) {
System.out.println(e.getMessage());
}
} else if (args[0].compareTo("9j") == 0) {
try {
final String m1 = "9";
final String t1 = "1 3 2 4 6 5 7 9 8";
final String t2 = "2 8 5 6 3 9 7 4 1";
String j = "";
for (int i = 1; i <= 9; i++) {
j += args[i] + " ";
}
final BigSurdVec w = Wigner3j.wigner3j(m1, t1, t2, j);
System.out.println(w.toString());
} catch (final Exception e) {
System.out.println(e.getMessage());
}
} else if (args[0].compareTo("3jm") == 0) {
final int j1 = new Integer(args[1]).intValue();
final int j2 = new Integer(args[2]).intValue();
final int j3 = new Integer(args[3]).intValue();
final int m1 = new Integer(args[4]).intValue();
final int m2 = new Integer(args[5]).intValue();
final int m3 = new Integer(args[6]).intValue();
try {
BigSurd w = Wigner3j.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 (final 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
* @throws Error
*/
static public BigSurd wigner3jm(final int j1, final int j2, final int j3, final int m1, final int m2, final int m3)
throws Error {
final Rational J1 = new Rational(j1, 2);
final Rational J2 = new Rational(j2, 2);
final Rational J3 = new Rational(j3, 2);
final Rational M1 = new Rational(m1, 2);
final Rational M2 = new Rational(m2, 2);
final Rational M3 = new Rational(m3, 2);
return Wigner3j.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="https://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
* @throws Error
*/
static public BigSurdVec wigner3j(final String m1, final String t1, final String t2, final String j) throws Error {
/*
* 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);
final 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.
*/
final int[] jvec = new int[m];
final 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);
}
final 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.
*/
final int[][] jhash = new int[m][2];
for (ti = 0; ti < 2 * m; ti++) {
final 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.
*/
final 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.
*/
final 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.
*/
final Rational[] M = new Rational[J.length];
s.close();
return Wigner3j.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
* @throws Error
*/
static private BigSurdVec wigner3j(final int[] tvec, final Rational[] J, final Rational[] M, final int[] triadidx)
throws Error {
/*
* 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.
*/
final int triadn = i / 3;
final int triadr = i % 3;
/*
* the neighbors in the triad have indices triadn*3+ (tiradr+1)
* mod 3 and triadn*3+(triadr+2) mod3
*/
final int nei1 = 3 * triadn + (triadr + 1) % 3;
final 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));
}
final 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) {
final 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.
*/
final int triadn = freeM / 3;
final int triadr = freeM % 3;
/*
* the neighbors in the triad have indices triadn*3+ (triadr+1)
* mod 3 and triadn*3+(triadr+2) mod3
*/
final int nei1 = 3 * triadn + (triadr + 1) % 3;
final 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.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) */
final 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.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(Wigner3j.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
* @throws Error
*/
static protected BigSurd wigner3jm(final Rational j1, final Rational j2, final Rational j3, final Rational m1,
final Rational m2, final Rational m3) throws Error {
/*
* 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|
*/
final Rational j1m2 = j1.subtract(j2);
if (j1m2.abs().compareTo(j3) > 0) {
return BigSurd.ZERO;
}
final 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;
}
final Factorial f = new Factorial();
Rational sumk = new Rational();
while (true) {
final 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++;
final Rational disc = new Rational(s, f.at(k));
return new BigSurd(sumk, disc);
} /* wigner3jm */
} /* Wigner3j */

335
core/src/main/java/org/nevec/rjm/Wigner3jGUI.java
View File

@ -1,335 +0,0 @@
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;
import it.cavallium.warppi.util.Error;
/**
* 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);
final Font defFont = new Font("Monospaced", Font.PLAIN, 11);
fram.setBackground(new Color(250, 250, 250));
fram.setForeground(new Color(0, 0, 0));
final Color fg = new Color(0, 200, 0);
final 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 */
/**
* @throws Error
* @since 2010-08-27
*/
public void compute() throws Error {
final String tr = inpGtria.getText();
final 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 (final 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
*/
final String j = inpGjval.getText();
s = new Scanner(j);
while (true) {
try {
trias[3] = s.nextLine().trim();
} catch (final Exception e) {
s.close();
return;
}
if (!trias[3].startsWith("#")) {
try {
final BigSurdVec w = Wigner3j.wigner3j(trias[0], trias[1], trias[2], trias[3]);
outG.append(w.toString() + " = " + w.doubleValue());
} catch (final Exception e) {
outG.append(e.toString());
e.printStackTrace();
}
outG.append(" # J = ");
final Scanner num = new Scanner(trias[3]);
while (num.hasNextInt()) {
final int twoj1 = num.nextInt();
final 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
*/
@Override
public void actionPerformed(final ActionEvent e) {
final String lin = e.getActionCommand();
/*
* debugging System.out.println("Ac"+e.paramString()) ;
* System.out.println(lin) ;
*/
if (lin == "compute") {
outG.setText("");
try {
compute();
} catch (final Error e1) {
// TODO Auto-generated catch block
e1.printStackTrace();
}
}
} /* actionPerformed */
/**
* Interpreter parser loop.
*
* @param e
* the information on which of the 3jn templates had been
* selected in the Menu
* @since 2011-02-18
*/
@Override
public void valueChanged(final 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(final String[] args) {
final Wigner3jGUI g = new Wigner3jGUI();
g.init();
} /* main */
} /* Wigner3jGUI */

22
desktop/pom.xml
View File

@ -10,20 +10,18 @@
</parent> </parent>
<artifactId>warppi-desktop</artifactId> <artifactId>warppi-desktop</artifactId>
<packaging>bundle</packaging>
<name>WarpPI Calculator Desktop</name> <name>WarpPI Calculator Desktop</name>
<description>WarpPI Calculator desktop project</description> <description>WarpPI Calculator desktop project</description>
<dependencies> <dependencies>
<dependency> <dependency>
<groupId>it.cavallium</groupId> <groupId>it.cavallium</groupId>
<artifactId>warppi-core</artifactId> <artifactId>warppi-core</artifactId>
<version>${project.version}</version> <version>0.9.0a2</version>
</dependency> </dependency>
<dependency> <dependency>
<groupId>it.cavallium</groupId> <groupId>it.cavallium</groupId>
<artifactId>warppi-engine-jogl</artifactId> <artifactId>warppi-engine-jogl</artifactId>
<version>${project.version}</version> <version>0.9.0a2</version>
</dependency> </dependency>
<dependency> <dependency>
<groupId>org.fusesource.jansi</groupId> <groupId>org.fusesource.jansi</groupId>
@ -55,23 +53,15 @@
</resource> </resource>
</resources> </resources>
<plugins> <plugins>
<plugin>
<groupId>org.apache.maven.plugins</groupId>
<artifactId>maven-compiler-plugin</artifactId>
</plugin>
<plugin> <plugin>
<groupId>org.apache.maven.plugins</groupId> <groupId>org.apache.maven.plugins</groupId>
<artifactId>maven-source-plugin</artifactId> <artifactId>maven-source-plugin</artifactId>
</plugin> </plugin>
<plugin>
<groupId>org.apache.felix</groupId>
<artifactId>maven-bundle-plugin</artifactId>
<extensions>true</extensions>
<configuration>
<instructions>
<Import-Package>*</Import-Package>
<Export-Package>it.cavallium.warppi.*</Export-Package>
<Bundle-SymbolicName>warppi-desktop</Bundle-SymbolicName>
</instructions>
</configuration>
</plugin>
<!-- Maven Assembly Plugin --> <!-- Maven Assembly Plugin -->
<plugin> <plugin>
<groupId>org.apache.maven.plugins</groupId> <groupId>org.apache.maven.plugins</groupId>

16
engine-jogl/pom.xml
View File

@ -10,15 +10,13 @@
</parent> </parent>
<artifactId>warppi-engine-jogl</artifactId> <artifactId>warppi-engine-jogl</artifactId>
<packaging>bundle</packaging>
<name>WarpPI Calculator JOGL Engine</name> <name>WarpPI Calculator JOGL Engine</name>
<description>WarpPI Calculator engine-jogl project</description> <description>WarpPI Calculator engine-jogl project</description>
<dependencies> <dependencies>
<dependency> <dependency>
<groupId>it.cavallium</groupId> <groupId>it.cavallium</groupId>
<artifactId>warppi-core</artifactId> <artifactId>warppi-core</artifactId>
<version>${project.version}</version> <version>0.9.0a2</version>
</dependency> </dependency>
<dependency> <dependency>
<groupId>org.jogamp.jogl</groupId> <groupId>org.jogamp.jogl</groupId>
@ -38,16 +36,8 @@
<build> <build>
<plugins> <plugins>
<plugin> <plugin>
<groupId>org.apache.felix</groupId> <groupId>org.apache.maven.plugins</groupId>
<artifactId>maven-bundle-plugin</artifactId> <artifactId>maven-compiler-plugin</artifactId>
<extensions>true</extensions>
<configuration>
<instructions>
<Import-Package>*</Import-Package>
<Export-Package>it.cavallium.warppi.*</Export-Package>
<Bundle-SymbolicName>warppi-engine-jogl</Bundle-SymbolicName>
</instructions>
</configuration>
</plugin> </plugin>
<plugin> <plugin>
<groupId>org.apache.maven.plugins</groupId> <groupId>org.apache.maven.plugins</groupId>

0
libs/raspi2fb → hardware/libs/raspi2fb
View File

4
hardware/pom.xml
View File

@ -18,12 +18,12 @@
<dependency> <dependency>
<groupId>it.cavallium</groupId> <groupId>it.cavallium</groupId>
<artifactId>warppi-core</artifactId> <artifactId>warppi-core</artifactId>
<version>${project.version}</version> <version>0.9.0a2</version>
</dependency> </dependency>
<dependency> <dependency>
<groupId>it.cavallium</groupId> <groupId>it.cavallium</groupId>
<artifactId>warppi-engine-jogl</artifactId> <artifactId>warppi-engine-jogl</artifactId>
<version>${project.version}</version> <version>0.9.0a2</version>
</dependency> </dependency>
<dependency> <dependency>
<groupId>com.pi4j</groupId> <groupId>com.pi4j</groupId>

2
pom.xml
View File

@ -39,6 +39,8 @@
<modules> <modules>
<module>core</module> <module>core</module>
<module>util</module>
<module>bigdecimal-math</module>
<module>Flow</module> <module>Flow</module>
<module>rules</module> <module>rules</module>
<!-- <module>hardware</module>--> <!-- <module>hardware</module>-->

15
rules/pom.xml
View File

@ -9,7 +9,6 @@
<version>0.9.0a2</version> <version>0.9.0a2</version>
</parent> </parent>
<artifactId>warppi-rules</artifactId> <artifactId>warppi-rules</artifactId>
<packaging>bundle</packaging>
<name>WarpPI Calculator Rules</name> <name>WarpPI Calculator Rules</name>
<description>WarpPI Calculator rules project</description> <description>WarpPI Calculator rules project</description>
@ -18,23 +17,15 @@
<dependency> <dependency>
<groupId>it.cavallium</groupId> <groupId>it.cavallium</groupId>
<artifactId>warppi-core</artifactId> <artifactId>warppi-core</artifactId>
<version>${project.version}</version> <version>0.9.0a2</version>
</dependency> </dependency>
</dependencies> </dependencies>
<build> <build>
<plugins> <plugins>
<plugin> <plugin>
<groupId>org.apache.felix</groupId> <groupId>org.apache.maven.plugins</groupId>
<artifactId>maven-bundle-plugin</artifactId> <artifactId>maven-compiler-plugin</artifactId>
<extensions>true</extensions>
<configuration>
<instructions>
<Import-Package>*</Import-Package>
<Export-Package>it.cavallium.warppi.*</Export-Package>
<Bundle-SymbolicName>warppi-rules</Bundle-SymbolicName>
</instructions>
</configuration>
</plugin> </plugin>
<plugin> <plugin>
<groupId>org.apache.maven.plugins</groupId> <groupId>org.apache.maven.plugins</groupId>

4
teavm/pom.xml
View File

@ -17,12 +17,12 @@
<dependency> <dependency>
<groupId>it.cavallium</groupId> <groupId>it.cavallium</groupId>
<artifactId>warppi-core</artifactId> <artifactId>warppi-core</artifactId>
<version>${project.version}</version> <version>0.9.0a2</version>
</dependency> </dependency>
<dependency> <dependency>
<groupId>it.cavallium</groupId> <groupId>it.cavallium</groupId>
<artifactId>warppi-rules</artifactId> <artifactId>warppi-rules</artifactId>
<version>${project.version}</version> <version>0.9.0a2</version>
</dependency> </dependency>
<dependency> <dependency>
<groupId>org.teavm</groupId> <groupId>org.teavm</groupId>

27
util/.classpath Normal file
View File

@ -0,0 +1,27 @@
<?xml version="1.0" encoding="UTF-8"?>
<classpath>
<classpathentry kind="src" output="target/classes" path="src/main/java">
<attributes>
<attribute name="optional" value="true"/>
<attribute name="maven.pomderived" value="true"/>
</attributes>
</classpathentry>
<classpathentry kind="src" output="target/test-classes" path="src/test/java">
<attributes>
<attribute name="optional" value="true"/>
<attribute name="maven.pomderived" value="true"/>
<attribute name="test" value="true"/>
</attributes>
</classpathentry>
<classpathentry kind="con" path="org.eclipse.jdt.launching.JRE_CONTAINER/org.eclipse.jdt.internal.debug.ui.launcher.StandardVMType/JavaSE-1.8">
<attributes>
<attribute name="maven.pomderived" value="true"/>
</attributes>
</classpathentry>
<classpathentry kind="con" path="org.eclipse.m2e.MAVEN2_CLASSPATH_CONTAINER">
<attributes>
<attribute name="maven.pomderived" value="true"/>
</attributes>
</classpathentry>
<classpathentry kind="output" path="target/classes"/>
</classpath>

1
util/.gitignore vendored Normal file
View File

@ -0,0 +1 @@
/target/

23
util/.project Normal file
View File

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

38
util/pom.xml Normal file
View File

@ -0,0 +1,38 @@
<?xml version="1.0" encoding="UTF-8"?>
<project xmlns="http://maven.apache.org/POM/4.0.0"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/xsd/maven-4.0.0.xsd">
<modelVersion>4.0.0</modelVersion>
<parent>
<groupId>it.cavallium</groupId>
<artifactId>warppi</artifactId>
<version>0.9.0a2</version>
</parent>
<artifactId>warppi-util</artifactId>
<name>warppi-util</name>
<description>WarpPI Util</description>
<developers>
<developer>
<id>cavallium</id>
<name>Andrea Cavalli</name>
<email>andrea@warp.ovh</email>
<url>https://cavallium.it/</url>
</developer>
</developers>
<build>
<plugins>
<plugin>
<groupId>org.apache.maven.plugins</groupId>
<artifactId>maven-compiler-plugin</artifactId>
</plugin>
<plugin>
<groupId>org.apache.maven.plugins</groupId>
<artifactId>maven-source-plugin</artifactId>
</plugin>
</plugins>
</build>
</project>

0
core/src/main/java/it/cavallium/warppi/util/Error.java → util/src/main/java/it/cavallium/warppi/util/Error.java
View File

0
core/src/main/java/it/cavallium/warppi/util/Errors.java → util/src/main/java/it/cavallium/warppi/util/Errors.java
View File