cmd/compile: asymmetric constant arithmetic leads to hilbert test failure #13472

griesemer · 2015-12-03T22:45:12Z

hilbert.go (attached below) computes the product H*H' of a 4x4 Hilbert matrix H and its inverse H' using constant arithmetic only. The product should be the identity matrix I. The Hilbert matrix is famously ill-conditioned leading quickly to unrecoverable floating-point errors even given the large (but limited) mantissa size used for constant arithmetic by the compiler.

That said, the matrices are symmetric and thus H*H' should lead to a sequence of totally symmetric computations (but for opposite signs). However, this test fails with a non-symmetric error, indicating that the two computations for elements (3, 4) and (4, 3) are not symmetric.

$ go run hilbert.go 
+1.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000
+0.000000e+000 +1.000000e+000 +0.000000e+000 +0.000000e+000
+0.000000e+000 +0.000000e+000 +1.000000e+000 +3.818670e-152
+0.000000e+000 +0.000000e+000 +0.000000e+000 +1.000000e+000

Investigate if this is a compiler or big.Float issue.

// WARNING: GENERATED FILE - DO NOT MODIFY MANUALLY!
// (To generate, in go/types directory: go test -run=Hilbert -H=4 -out="hilbert.go")

// This program tests arbitrary precision constant arithmetic
// by generating the constant elements of a Hilbert matrix H,
// its inverse I, and the product P = H*I. The product should
// be the identity matrix.
package main

func main() {
    if !ok {
        printProduct()
        return
    }
    println("PASS")
}

// Hilbert matrix, n = 4
const (
    h0_0, h0_1, h0_2, h0_3 = 1.0/(iota + 1), 1.0/(iota + 2), 1.0/(iota + 3), 1.0/(iota + 4)
    h1_0, h1_1, h1_2, h1_3
    h2_0, h2_1, h2_2, h2_3
    h3_0, h3_1, h3_2, h3_3
)

// Inverse Hilbert matrix
const (
    i0_0 = +1 * b4_3 * b4_3 * b0_0 * b0_0
    i0_1 = -2 * b4_2 * b5_3 * b1_0 * b1_0
    i0_2 = +3 * b4_1 * b6_3 * b2_0 * b2_0
    i0_3 = -4 * b4_0 * b7_3 * b3_0 * b3_0

    i1_0 = -2 * b5_3 * b4_2 * b1_1 * b1_1
    i1_1 = +3 * b5_2 * b5_2 * b2_1 * b2_1
    i1_2 = -4 * b5_1 * b6_2 * b3_1 * b3_1
    i1_3 = +5 * b5_0 * b7_2 * b4_1 * b4_1

    i2_0 = +3 * b6_3 * b4_1 * b2_2 * b2_2
    i2_1 = -4 * b6_2 * b5_1 * b3_2 * b3_2
    i2_2 = +5 * b6_1 * b6_1 * b4_2 * b4_2
    i2_3 = -6 * b6_0 * b7_1 * b5_2 * b5_2

    i3_0 = -4 * b7_3 * b4_0 * b3_3 * b3_3
    i3_1 = +5 * b7_2 * b5_0 * b4_3 * b4_3
    i3_2 = -6 * b7_1 * b6_0 * b5_3 * b5_3
    i3_3 = +7 * b7_0 * b7_0 * b6_3 * b6_3

)

// Product matrix
const (
    p0_0 = h0_0*i0_0 + h0_1*i1_0 + h0_2*i2_0 + h0_3*i3_0
    p0_1 = h0_0*i0_1 + h0_1*i1_1 + h0_2*i2_1 + h0_3*i3_1
    p0_2 = h0_0*i0_2 + h0_1*i1_2 + h0_2*i2_2 + h0_3*i3_2
    p0_3 = h0_0*i0_3 + h0_1*i1_3 + h0_2*i2_3 + h0_3*i3_3

    p1_0 = h1_0*i0_0 + h1_1*i1_0 + h1_2*i2_0 + h1_3*i3_0
    p1_1 = h1_0*i0_1 + h1_1*i1_1 + h1_2*i2_1 + h1_3*i3_1
    p1_2 = h1_0*i0_2 + h1_1*i1_2 + h1_2*i2_2 + h1_3*i3_2
    p1_3 = h1_0*i0_3 + h1_1*i1_3 + h1_2*i2_3 + h1_3*i3_3

    p2_0 = h2_0*i0_0 + h2_1*i1_0 + h2_2*i2_0 + h2_3*i3_0
    p2_1 = h2_0*i0_1 + h2_1*i1_1 + h2_2*i2_1 + h2_3*i3_1
    p2_2 = h2_0*i0_2 + h2_1*i1_2 + h2_2*i2_2 + h2_3*i3_2
    p2_3 = h2_0*i0_3 + h2_1*i1_3 + h2_2*i2_3 + h2_3*i3_3

    p3_0 = h3_0*i0_0 + h3_1*i1_0 + h3_2*i2_0 + h3_3*i3_0
    p3_1 = h3_0*i0_1 + h3_1*i1_1 + h3_2*i2_1 + h3_3*i3_1
    p3_2 = h3_0*i0_2 + h3_1*i1_2 + h3_2*i2_2 + h3_3*i3_2
    p3_3 = h3_0*i0_3 + h3_1*i1_3 + h3_2*i2_3 + h3_3*i3_3

)

// Verify that product is the identity matrix
const ok =
    p0_0 == 1 && p0_1 == 0 && p0_2 == 0 && p0_3 == 0 &&
    p1_0 == 0 && p1_1 == 1 && p1_2 == 0 && p1_3 == 0 &&
    p2_0 == 0 && p2_1 == 0 && p2_2 == 1 && p2_3 == 0 &&
    p3_0 == 0 && p3_1 == 0 && p3_2 == 0 && p3_3 == 1 &&
    true

func printProduct() {
    println(p0_0, p0_1, p0_2, p0_3)
    println(p1_0, p1_1, p1_2, p1_3)
    println(p2_0, p2_1, p2_2, p2_3)
    println(p3_0, p3_1, p3_2, p3_3)
}

// Binomials
const (
    b0_0 = f0 / (f0*f0)

    b1_0 = f1 / (f0*f1)
    b1_1 = f1 / (f1*f0)

    b2_0 = f2 / (f0*f2)
    b2_1 = f2 / (f1*f1)
    b2_2 = f2 / (f2*f0)

    b3_0 = f3 / (f0*f3)
    b3_1 = f3 / (f1*f2)
    b3_2 = f3 / (f2*f1)
    b3_3 = f3 / (f3*f0)

    b4_0 = f4 / (f0*f4)
    b4_1 = f4 / (f1*f3)
    b4_2 = f4 / (f2*f2)
    b4_3 = f4 / (f3*f1)
    b4_4 = f4 / (f4*f0)

    b5_0 = f5 / (f0*f5)
    b5_1 = f5 / (f1*f4)
    b5_2 = f5 / (f2*f3)
    b5_3 = f5 / (f3*f2)
    b5_4 = f5 / (f4*f1)
    b5_5 = f5 / (f5*f0)

    b6_0 = f6 / (f0*f6)
    b6_1 = f6 / (f1*f5)
    b6_2 = f6 / (f2*f4)
    b6_3 = f6 / (f3*f3)
    b6_4 = f6 / (f4*f2)
    b6_5 = f6 / (f5*f1)
    b6_6 = f6 / (f6*f0)

    b7_0 = f7 / (f0*f7)
    b7_1 = f7 / (f1*f6)
    b7_2 = f7 / (f2*f5)
    b7_3 = f7 / (f3*f4)
    b7_4 = f7 / (f4*f3)
    b7_5 = f7 / (f5*f2)
    b7_6 = f7 / (f6*f1)
    b7_7 = f7 / (f7*f0)
)

// Factorials
const (
    f0 = 1
    f1 = 1
    f2 = f1 * 2
    f3 = f2 * 3
    f4 = f3 * 4
    f5 = f4 * 5
    f6 = f5 * 6
    f7 = f6 * 7
)

The text was updated successfully, but these errors were encountered:

griesemer · 2016-01-07T22:59:52Z

This is working as intended, even if unexpected. An analysis follows below.

Simplifying the constant operations to the 2 values in question, we get the two expressions:

    2800*(1.0/6) - (280.0/6 + 420)
    2800/6.0 - (280.0/6 + 420)

which are identical mathematically, but produce the observed error if computed with limited precisions:

package main

func main() {
    const a = 2800*(1.0/6) - (280.0/6 + 420)
    const b = 2800/6.0 - (280.0/6 + 420)
    println(a, b, a == b)
}

The results are:

+3.818670e-152 +0.000000e+000 false

In other words, x/3.0 is not the same as x*(1/3.0):

package main

func main() {
    const a = 1400*(1/3.0)
    const b = 1400/3.0
    println(a, b, a-b, a == b)
}

The difference is:

+4.666667e+002 +4.666667e+002 +3.818670e-152 false

Doing the constant evaluation as done by the compiler explicitly confirms the result:
http://play.golang.org/p/nR2dfat3A0

466.6666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667
466.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666665
3.8186704543745058525649742415170221247367534257154083298670173866016145331559005721273215438850190353210939408669880839441473167847553549621689163927511474e-152
1

Looking at the bits:

package main

import (
    "fmt"
    "math/big"
)

func newf(x int64) *big.Float {
    return new(big.Float).SetPrec(512).SetInt64(x)
}

func main() {
    c := newf(1400)
    one := newf(1)
    three := newf(3)

    a := newf(0).Mul(c, newf(0).Quo(one, three))
    b := newf(0).Quo(c, three)

    fmt.Println(a.Text('p', 0))
    fmt.Println(b.Text('p', 0))
    fmt.Println(newf(0).Sub(a, b).Text('p', 0))
    fmt.Println(a.Cmp(b))
}

shows an error:

0x.e9555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555556p+9
0x.e9555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555p+9
0x.8p-502
1

mdempsky · 2016-01-07T23:26:25Z

Since floating-point literals are always rational and the only constant expression operators on floating-points (+, -, *, and /) are closed over the rationals, instead of representing floating point constants as just Sign * Mantissa256 * 10^Exponent32, could we represent them as Sign * MantissaQuotient256 / MantissaDivisor256 * 10^Exponent32? (I.e., instead of 256-bit integer + floating-point radix, use 512-bit rational + floating-point radix.)

There would still be cases where values are inexactly represented, but at least expressions like 1400*(1/3.0) - 1400/3.0 would be exactly computable.

griesemer · 2016-01-07T23:42:12Z

We could but it doesn't solve the general problem, and it's yet another number representation. What we have is fine for now.

go/constants now transitions transparently from rationals to floats if need be (which is almost never the case). It would be simpler to just plug in that code instead of writing new one. That would give precise answers in probably all realistic cases.

griesemer self-assigned this Dec 3, 2015

griesemer added this to the Unplanned milestone Dec 3, 2015

griesemer closed this as completed Jan 7, 2016

golang locked and limited conversation to collaborators Jan 7, 2017

gopherbot added the FrozenDueToAge label Jan 7, 2017

rsc unassigned griesemer Jun 23, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

cmd/compile: asymmetric constant arithmetic leads to hilbert test failure #13472

cmd/compile: asymmetric constant arithmetic leads to hilbert test failure #13472

griesemer commented Dec 3, 2015

griesemer commented Jan 7, 2016

mdempsky commented Jan 7, 2016

griesemer commented Jan 7, 2016

cmd/compile: asymmetric constant arithmetic leads to hilbert test failure #13472

cmd/compile: asymmetric constant arithmetic leads to hilbert test failure #13472

Comments

griesemer commented Dec 3, 2015

griesemer commented Jan 7, 2016

mdempsky commented Jan 7, 2016

griesemer commented Jan 7, 2016