Multiplication of integer intervals produce wrong results with GNU-Prolog

[mlliarm]

Артем Андросов (Artem Androsov) reported this issue in private communication. Attaching his printscreen to this issue:

I managed to reproduce the issue on GNU Prolog 1.5.0 (64 bits), on GNU/Linux:

| ?- interval_arithmetic::(new(33096, 33096,X), mul(X, X, X2), mul(X2, X, X3), mul(X3, X, X4), mul(X4, X, X5), mul(X5, X, X6), mul(X6, X, X7), mul(X7, X, X8), mul(X8, X, X9), mul(X9, X, X10)).

X = i(33096,33096)
X10 = i(930066664987295744,930066664987295744)
X2 = i(1095345216,1095345216)
X3 = i(36251545268736,36251545268736)
X4 = i(-1106061866999607296,-1106061866999607296)
X5 = i(-965778951611580416,-965778951611580416)
X6 = i(175611183360114688,175611183360114688)
X7 = i(-1002501741366738944,-1002501741366738944)
X8 = i(-22572697749815296,-22572697749815296)
X9 = i(27130257349804032,27130257349804032)

yes
| ?- interval_arithmetic::(new(33096.0, 33096.0,X), mul(X, X, X2), mul(X2, X, X3), mul(X3, X, X4), mul(X4, X, X5), mul(X5, X, X6), mul(X6, X, X7), mul(X7, X, X8), mul(X8, X, X9), mul(X9, X, X10)).

X = i(33096.0,33096.0)
X10 = i(1.5767218239165622e+45,1.5767218239165622e+45)
X2 = i(1095345216.0,1095345216.0)
X3 = i(36251545268736.0,36251545268736.0)
X4 = i(1.1997811422140867e+18,1.1997811422140867e+18)
X5 = i(3.970795668271741e+22,3.970795668271741e+22)
X6 = i(1.3141745343712153e+27,1.3141745343712153e+27)
X7 = i(4.3493920389549741e+31,4.3493920389549741e+31)
X8 = i(1.4394747892125382e+36,1.4394747892125382e+36)
X9 = i(4.7640857623778164e+40,4.7640857623778164e+40)

(1 ms) yes

SWI-Prolog: 8.5.5-3-gb856d332c-DIRTY seems to be producing the correct powers of the interval X in both integer and floating point calculations:

?- interval_arithmetic::(new(33096, 33096,X), mul(X, X, X2), mul(X2, X, X3), mul(X3, X, X4), mul(X4, X, X5), mul(X5, X, X6), mul(X6, X, X7), mul(X7, X, X8), mul(X8, X, X9), mul(X9, X, X10)).
X = i(33096, 33096),
X2 = i(1095345216, 1095345216),
X3 = i(36251545268736, 36251545268736),
X4 = i(1199781142214086656, 1199781142214086656),
X5 = i(39707956682717411966976, 39707956682717411966976),
X6 = i(1314174534371215466459037696, 1314174534371215466459037696),
X7 = i(43493920389549747077928311586816, 43493920389549747077928311586816),
X8 = i(1439474789212538429291115400277262336, 1439474789212538429291115400277262336),
X9 = i(47640857623778171855818755287576274272256, 47640857623778171855818755287576274272256),
X10 = i(1576721823916562375740177524997624373314584576, 1576721823916562375740177524997624373314584576).

?- interval_arithmetic::(new(33096.0, 33096.0,X), mul(X, X, X2), mul(X2, X, X3), mul(X3, X, X4), mul(X4, X, X5), mul(X5, X, X6), mul(X6, X, X7), mul(X7, X, X8), mul(X8, X, X9), mul(X9, X, X10)).
X = i(33096.0, 33096.0),
X2 = i(1095345216.0, 1095345216.0),
X3 = i(36251545268736.0, 36251545268736.0),
X4 = i(1.1997811422140867e+18, 1.1997811422140867e+18),
X5 = i(3.970795668271741e+22, 3.970795668271741e+22),
X6 = i(1.3141745343712153e+27, 1.3141745343712153e+27),
X7 = i(4.349392038954974e+31, 4.349392038954974e+31),
X8 = i(1.4394747892125382e+36, 1.4394747892125382e+36),
X9 = i(4.764085762377816e+40, 4.764085762377816e+40),
X10 = i(1.576721823916562e+45, 1.576721823916562e+45).

This limitation with GNU-Prolog should be mentioned/highlighted on the NOTES.md.

[mlliarm]

Seems like a limitation of GNU-Prolog on how it handles multiplication between large integers. Not sure if we can do much about that, besides making a disclaimer and suggestion to users not to use GNU-Prolog as a backend.

[mlliarm]

It's a limitation of the is predicate in GNU-Prolog.

Proof:

(GNU-Prolog 1.5.0):

| ?- X = 33096, X2 is X*X, X3 is X2*X, X4 is X3*X, X5 is X4*X, X6 is X5*X, X7 is X6*X, X8 is X7*X, X9 is X8*X, X10 is X9*X.

X = 33096
X10 = 930066664987295744
X2 = 1095345216
X3 = 36251545268736
X4 = -1106061866999607296
X5 = -965778951611580416
X6 = 175611183360114688
X7 = -1002501741366738944
X8 = -22572697749815296
X9 = 27130257349804032

(SWI-Prolog 8.5.5-3-gb856d332c-DIRTY):

?- X = 33096, X2 is X*X, X3 is X2*X, X4 is X3*X, X5 is X4*X, X6 is X5*X, X7 is X6*X, X8 is X7*X, X9 is X8*X, X10 is X9*X.
X = 33096,
X2 = 1095345216,
X3 = 36251545268736,
X4 = 1199781142214086656,
X5 = 39707956682717411966976,
X6 = 1314174534371215466459037696,
X7 = 43493920389549747077928311586816,
X8 = 1439474789212538429291115400277262336,
X9 = 47640857623778171855818755287576274272256,
X10 = 1576721823916562375740177524997624373314584576.

Python3.9 results:

>>> list(map(lambda x: 33096**x, list(range(1,11))))
[33096, 
1095345216, 
36251545268736, 
1199781142214086656, 
39707956682717411966976, 
1314174534371215466459037696, 
43493920389549747077928311586816, 
1439474789212538429291115400277262336, 
47640857623778171855818755287576274272256,
1576721823916562375740177524997624373314584576]

So SWI-Prolog works fine.

[pmoura]

Not all Prolog systems supports unbounded integer arithmetic. A warning in the readme file that mentions this issue should be enough. This likely only affects some uses of the library?

[mlliarm]

@pmoura

Well the power function X^n of an interval is pretty elementary (and thus important), especially if someone would like to do interval arithmetic with polynomials. I'll implement if as a separate predicate, but it will be based on the multiplication/3 predicate. So if multiplication goes bad, that will go bad too.

For example, someone should be able to use this library to calculate the roots of a polynomial say of 10th degree (O(X^10)), importing a big initial interval such as X0 = [3500, 4000], and I believe that this limitation would mess up with getting the correct results.

I think a warning in the NOTES.md or the README.md should be enough, I agree.

[mlliarm]

NOTES.md updated. Closing issue.