Skip to content

Commit

Permalink
Quadruple division without using BigDecimal
Browse files Browse the repository at this point in the history
  • Loading branch information
apete committed Aug 4, 2024
1 parent 23ffcce commit 35e0c31
Show file tree
Hide file tree
Showing 3 changed files with 184 additions and 10 deletions.
5 changes: 5 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,11 @@ Added / Changed / Deprecated / Fixed / Removed / Security
- The new LP solver implementation is now the default alternative.
- The two classes `LinearSolver.GeneralBuilder` and `LinearSolver.StandardBuilder` are replaced by a common `LinearSolver.Builder`. Consequently the methods `newGeneralBuilder()` and `newStandardBuilder()` are deprecated and replaced by `newBuilder()`.

#### org.ojalgo.scalar

- Re-implemented `Quadruple`'s divide method using primitives only. (Previously delegated to `BigDecimal`.)


## [54.0.0] – 2024-06-06

### Added
Expand Down
65 changes: 56 additions & 9 deletions src/main/java/org/ojalgo/scalar/Quadruple.java
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,7 @@ private static Quadruple add(final double base1, final double remainder1, final

t1 = base1 + base2;
e = t1 - base1;
t2 = ((base2 - e) + (base1 - (t1 - e))) + remainder1 + remainder2;
t2 = base2 - e + (base1 - (t1 - e)) + remainder1 + remainder2;

double base = t1 + t2;
double remainder = t2 - (base - t1);
Expand Down Expand Up @@ -206,14 +206,54 @@ private static Quadruple multiply(final double base1, final double remainder1, f

t1 = c11 + c2;
e = t1 - c11;
t2 = remainder1 * remainder2 + ((c2 - e) + (c11 - (t1 - e))) + c21;
t2 = remainder1 * remainder2 + (c2 - e + (c11 - (t1 - e))) + c21;

double base = t1 + t2;
double remainder = t2 - (base - t1);

return new Quadruple(base, remainder);
}

/**
* High-precision division.
*/
private static Quadruple divide(final double base1, final double remainder1, final double base2, final double remainder2) {

double q1, r1, r2, s1, s2, t1, t2;
double c1, c2, cona, conb, a1, a2, b1, b2;

q1 = base1 / base2;

cona = q1 * SPLIT;
a1 = cona - (cona - q1);
a2 = q1 - a1;

conb = base2 * SPLIT;
b1 = conb - (conb - base2);
b2 = base2 - b1;

// high part of q1 * base2
c1 = q1 * base2;
// rounding errors and cross terms
c2 = a1 * b1 - c1 + a1 * b2 + a2 * b1 + a2 * b2;

// first remainder
r1 = base1 - c1;
r2 = remainder1 - c2;

// next approximation for remainder1 / base2
s1 = r1 / base2;
t1 = s1 * base2;
t2 = a1 * b1 - t1 + a1 * b2 + a2 * b1 + a2 * b2;

s2 = r2 / base2;

double base = q1 + s1;
double remainder = s2 - (base - (q1 + s1));

return new Quadruple(base, remainder);
}

private final double myBase;
private transient BigDecimal myDecimal = null;
private final double myRemainder;
Expand Down Expand Up @@ -247,7 +287,7 @@ public Quadruple add(final Quadruple arg) {
}

if (Quadruple.isInfinite(this)) {
if (!Quadruple.isInfinite(arg) || (this.sign() == arg.sign())) {
if (!Quadruple.isInfinite(arg) || this.sign() == arg.sign()) {
return this;
} else {
return NaN;
Expand Down Expand Up @@ -285,8 +325,15 @@ public Quadruple divide(final Quadruple arg) {
return NaN;
}

return Quadruple.divide(this, arg);
// return Quadruple.divide(this, arg);

double base1 = this.getBase();
double remainder1 = this.getRemainder();

double base2 = arg.getBase();
double remainder2 = arg.getRemainder();

return Quadruple.divide(base1, remainder1, base2, remainder2);
}

@Override
Expand All @@ -308,8 +355,8 @@ public boolean equals(final Object obj) {
return false;
}
Quadruple other = (Quadruple) obj;
if ((Double.doubleToLongBits(myBase) != Double.doubleToLongBits(other.myBase))
|| (Double.doubleToLongBits(myRemainder) != Double.doubleToLongBits(other.myRemainder))) {
if (Double.doubleToLongBits(myBase) != Double.doubleToLongBits(other.myBase)
|| Double.doubleToLongBits(myRemainder) != Double.doubleToLongBits(other.myRemainder)) {
return false;
}
return true;
Expand All @@ -331,9 +378,9 @@ public int hashCode() {
int result = 1;
long temp;
temp = Double.doubleToLongBits(myBase);
result = prime * result + (int) (temp ^ (temp >>> 32));
result = prime * result + (int) (temp ^ temp >>> 32);
temp = Double.doubleToLongBits(myRemainder);
return prime * result + (int) (temp ^ (temp >>> 32));
return prime * result + (int) (temp ^ temp >>> 32);
}

@Override
Expand Down Expand Up @@ -432,7 +479,7 @@ public Quadruple subtract(final Quadruple arg) {
}

if (Quadruple.isInfinite(this)) {
if (!Quadruple.isInfinite(arg) || (this.sign() != arg.sign())) {
if (!Quadruple.isInfinite(arg) || this.sign() != arg.sign()) {
return this;
} else {
return NaN;
Expand Down
124 changes: 123 additions & 1 deletion src/test/java/org/ojalgo/scalar/QuadrupleTest.java
Original file line number Diff line number Diff line change
Expand Up @@ -83,9 +83,10 @@ public void testDivide() {
Quadruple denominator = Quadruple.valueOf(d);

Quadruple quotient = numerator.divide(denominator);
BigDecimal bigDecimal = quotient.toBigDecimal();

BigDecimal expected = NC_EVAL.enforce(MissingMath.divide(BigDecimal.valueOf(n), BigDecimal.valueOf(d)));
BigDecimal actual = NC_EVAL.enforce(quotient.toBigDecimal());
BigDecimal actual = NC_EVAL.enforce(bigDecimal);

TestUtils.assertEquals(expected, actual, NC_EVAL);

Expand All @@ -94,6 +95,76 @@ public void testDivide() {
}
}

@Test
public void testDivideBy3() {

Quadruple numerator = Quadruple.ONE;
Quadruple denominator = Quadruple.ONE.add(Quadruple.TWO);

Quadruple quotient = numerator.divide(denominator);

BigDecimal expected = MissingMath.divide(numerator.toBigDecimal(), denominator.toBigDecimal());
BigDecimal actual = quotient.toBigDecimal();

TestUtils.assertEquals(expected, actual, NC_EVAL);

expected = NC_EVAL.enforce(expected);
actual = NC_EVAL.enforce(actual);

TestUtils.assertTrue(actual.compareTo(expected) == 0);

// Negated denominator

numerator = Quadruple.ONE;
denominator = Quadruple.ONE.add(Quadruple.TWO).negate();

quotient = numerator.divide(denominator);

expected = MissingMath.divide(numerator.toBigDecimal(), denominator.toBigDecimal());
actual = quotient.toBigDecimal();

TestUtils.assertEquals(expected, actual, NC_EVAL);

expected = NC_EVAL.enforce(expected);
actual = NC_EVAL.enforce(actual);

TestUtils.assertTrue(actual.compareTo(expected) == 0);

// Negated numerator

numerator = Quadruple.ONE.negate();
denominator = Quadruple.ONE.add(Quadruple.TWO);

quotient = numerator.divide(denominator);

expected = MissingMath.divide(numerator.toBigDecimal(), denominator.toBigDecimal());
actual = quotient.toBigDecimal();

TestUtils.assertEquals(expected, actual, NC_EVAL);

expected = NC_EVAL.enforce(expected);
actual = NC_EVAL.enforce(actual);

TestUtils.assertTrue(actual.compareTo(expected) == 0);

// Negated numerator and denominator

numerator = Quadruple.ONE.negate();
denominator = Quadruple.ONE.add(Quadruple.TWO).negate();

quotient = numerator.divide(denominator);

expected = MissingMath.divide(numerator.toBigDecimal(), denominator.toBigDecimal());
actual = quotient.toBigDecimal();

TestUtils.assertEquals(expected, actual, NC_EVAL);

expected = NC_EVAL.enforce(expected);
actual = NC_EVAL.enforce(actual);

TestUtils.assertTrue(actual.compareTo(expected) == 0);
}

@Test
public void testInvert() {

Expand Down Expand Up @@ -130,6 +201,57 @@ public void testMultiply() {
TestUtils.assertEquals(expected, actual, NC_EVAL);

TestUtils.assertTrue(actual.compareTo(expected) == 0);

// Negated left

pi = BigMath.PI.negate();
e = BigMath.E;

left = Quadruple.valueOf(pi);
right = Quadruple.valueOf(e);

product = left.multiply(right);

expected = NC_EVAL.enforce(pi.multiply(e));
actual = NC_EVAL.enforce(product.toBigDecimal());

TestUtils.assertEquals(expected, actual, NC_EVAL);

TestUtils.assertTrue(actual.compareTo(expected) == 0);

// Negated right

pi = BigMath.PI;
e = BigMath.E.negate();

left = Quadruple.valueOf(pi);
right = Quadruple.valueOf(e);

product = left.multiply(right);

expected = NC_EVAL.enforce(pi.multiply(e));
actual = NC_EVAL.enforce(product.toBigDecimal());

TestUtils.assertEquals(expected, actual, NC_EVAL);

TestUtils.assertTrue(actual.compareTo(expected) == 0);

// Negated both

pi = BigMath.PI.negate();
e = BigMath.E.negate();

left = Quadruple.valueOf(pi);
right = Quadruple.valueOf(e);

product = left.multiply(right);

expected = NC_EVAL.enforce(pi.multiply(e));
actual = NC_EVAL.enforce(product.toBigDecimal());

TestUtils.assertEquals(expected, actual, NC_EVAL);

TestUtils.assertTrue(actual.compareTo(expected) == 0);
}

@Test
Expand Down

0 comments on commit 35e0c31

Please sign in to comment.