开发者

generate ill-conditioned data for testing floating point summation

开发者 https://www.devze.com 2023-03-18 07:17 出处:网络
I have implemented a Kahan floating point summation algorithm in Java. I want to compare it against the built-in floating point addition in Java and infinite precision addition in Mathematica. However

I have implemented a Kahan floating point summation algorithm in Java. I want to compare it against the built-in floating point addition in Java and infinite precision addition in Mathematica. However the data set I have is not good for testing, because the numbers are close to each other. (Condition number ~开发者_如何学JAVA= 1)

Running Kahan on my data set gives all most the same result as the built-in +.

Could anyone suggest how to generate a large amount of data that can potentially cause serious rounding off error?


However the data set I have is not good for testing, because the numbers are close to each other.

It sounds like you already know what the problem is. Get to it =)

There are a few things that you will want:

  • Numbers of wildly different magnitudes, so that most of the precision of the smaller number is lost with naive summation.
  • Numbers with different signs and nearly equal (or equal) magnitudes, such that catastrophic cancellation occurs.
  • Numbers that have some low-order bits set, to increase the effects of rounding.

To get you started, you could try some simple three-term sums, which should show the effect clearly:

1.0 + 1.0e-20 - 1.0

Evaluated with simple summation, this will give 0.0; clearly incorrect. You might also look at sums of the form:

a0 + a1 + a2 + ... + an - b

Where b is the sum a0 + ... + an evaluated naively.


You want a heap of high precision numbers? Try this:

double[] nums = new double[SIZE];
for (int i = 0; i < SIZE; i++)
    nums[i] = Math.rand();


Are we talking about number pairs or sequences?

If pairs, start with 1 for both numbers, then in every iteration divide one by 3, multiply the other by 3. It's easy to calculate the theoretical sums of those pairs and you'll get a whole host of rounding errors. (Some from the division and some from the addition. If you don't want division errors, then use 2 instead of 3.)


By experiment, I found following pattern:

public static void main(String[] args) {
    System.out.println(1.0 / 3 - 0.01 / 3);

    System.out.println(1.0 / 7 - 0.01 / 7);

    System.out.println(1.0 / 9 - 0.001 / 9);
}

I've subtracted close negative powers of prime numbers (which should not have exact representation in binary form). However, there are cases then such expression evaluates correctly, for example

    System.out.println(1.0 / 9 - 0.01 / 9);

You can automate this approach by iterating power of subtrahend and stopping when multiplication by appropriate value doesn't yield integer number, for example:

    System.out.println((1.0 / 9 - 0.001 / 9) * 9000);
    if (1000 - (1.0 / 9 - 0.001 / 9) * 9000 > 1.0)
        System.out.println("Found it!");


Scalacheck might be something for you. Here is a short sample:

cat DoubleSpecification.scala
import org.scalacheck._

object DoubleSpecification extends Properties ("Doubles") {
        /* 
            (a/1000 + b/1000) = (a+b) / 1000
            (a/x    + b/x   ) = (a+b) / x
        */
        property ("distributive") = Prop.forAll { (a: Int, b: Int, c: Int) => 
            (c == 0 || a*1.0/c + b*1.0/c == (a+b) * 1.0 / c)            }
}       

object Runner { 
    def main (args: Array[String]) {
        DoubleSpecification.check
        println ("...done")
    }
}

To run it, you need scala, and the schalacheck-jar. I used version 2.8 (I don't have to say, that your c-path will vary):

 scalac -cp /opt/scala/lib/scalacheck.jar:. DoubleSpecification.scala 
 scala -cp /opt/scala/lib/scalacheck.jar:. DoubleSpecification
! Doubles.distributive: Falsified after 6 passed tests.                       
> ARG_0: 28 (orig arg: 1030341)
> ARG_1: 9 (orig arg: 2147483647)
> ARG_2: 5

Scalacheck takes some random values (orig args) and tries to simplify these, if the test fails, in order to find simple examples.

0

精彩评论

暂无评论...
验证码 换一张
取 消