Floating-point irregularities in ga68
Nelson H. F. Beebe
beebe@math.utah.edu
Sat May 9 15:20:41 GMT 2026
I'm the author of a large book on accurate numerical computation of
elementary and special functions:
The Mathematical-Function Computation Handbook: Programming
Using the {MathCW} Portable Software Library
Springer (2017)
https://doi.org/10.1007/978-3-319-64110-2
https://www.math.utah.edu/pub/tex/bib/elefunt.html#Beebe:2017:MFC
The book, and its software, addresses several historical
floating-point designs, as well as IEEE 754 arithmetic, in both binary
and decimal formats. Since 2006, I have routinely built patched
versions of gcc releases 4 to 17 to support decimal arithmetic in C,
and recently, the recent, and independently developed by a single
expert programmer, kefir C compiler has been extended with decimal
support.
My book explains why I chose to write the code in highly portable C,
and it has appendices on interfacing the MathCW C library, -lmcw (a
drop-in, but much extended, replacement for -lm, and the formatted I/O
portions of -lc), to other languages: Ada, C++, C#, Fortran, Java, and
Pascal. Since then, I have developed interfaces and virtual book
appendices for Cobol, Go, Modula-2, Oberon, and Rust. I believe that
a similar appendix for Algol 68 is certainly feasible.
When I learned recently of the inclusion of Algol68 support in recent
gcc-16 and gcc-17 versions, I immediately added that languages to the
builds, and now have ga68 on CentOS and Ubuntu systems on x86_64, and
Debian on arm64 (Raspberry Pi + some QEMU-based virtual machines).
I then set about developing sample Algol 68 code for common numerical
programming, much as I have long done for other many programming
languages. I have now successfully interfaced Algol 68 code to -lm
and -lmcw, but quickly found serious issues. Here are four examples:
------------------------------------------------------------------------
This program compiles and runs successfully with the Genie Algol68
compiler:
% cat test1.a68
BEGIN
LONG REAL x, y;
x := LONG 0.5;
y := long sqrt(x);
print(("long sqrt(x) = ", float(y, 22, 18, 0), newline));
SKIP
END
% a68g test1.a68
long sqrt(x) = +7.0710678118654752e-1
However, this (I think) equivalent code for ga68 fails to compile:
% cat test-long-sqrt.a68
begin
long real x, y;
string newline = "'n";
x := long 0.5;
y := long sqrt(x);
puts("long sqrt(x) = " + float(y, 22, 18, 0) + newline)
end
% ga68 -c test-long-sqrt.a68
test-long-sqrt.a68:5:10: error: construct beginning with %<begin%> in \
line 1 followed by a serial clause starting in line 2 and then %<:=%> \
in line 5 and a longsety and a specification and %<;%> and a \
specification in line 6 and %<end%> in line 7 is an invalid keyword
5 | y := long sqrt(x);
| ^
Changing sqrt() to other Algol 68 functions produces similar
compilation errors.
Is this perhaps a temporary compiler limitation? As far as I can tell
from Web searches for code examples, and the Algol 68 specifications,
"long sqrt(x)", and also "long long sqrt(x)", should be the correct
way to access higher precision.
------------------------------------------------------------------------
This test file gets an infinite loop in float():
% cat inf.a68
{ Bug: infinite loop from ga68 inf.a68 && ./a.out }
begin
real x, y, z;
x := 1.0;
y := 0.0;
z := x / y;
puts("If Algol 68 arithmetic follows IEEE 754, this should be Infinity: " +
float(z, 30, 20, 0) + "'n");
skip
end
The loop is also seen on an arm64 system.
------------------------------------------------------------------------
This test file core dumps at run-time:
% cat nan.a68
{ Bug: core dump from ga68 nan.a68 && ./a.out }
begin
real x, y, z;
x := 0.0;
y := 0.0;
z := x / y;
puts("If Algol 68 arithmetic follows IEEE 754, this should be NaN: " +
float(z, 30, 20, 0) + "'n");
skip
end
% ga68 -g -g3 nan.a68
% gdb a.out
(gdb) run
...
Program received signal SIGABRT, Aborted.
0x00007ffff7a43387 in raise () from /lib64/libc.so.6
(gdb) where
#0 0x00007ffff7a43387 in raise () from /lib64/libc.so.6
#1 0x00007ffff7a44a78 in abort () from /lib64/libc.so.6
#2 0x000000000040191b in _libga68_abort (
fmt=fmt@entry=0x40c5f0 "%s:%u: runtime error: bound %zd out of range [%zd:%zd]\n") at ../../../gcc-17-20260426/libga68/ga68-error.c:48
#3 0x0000000000401a2e in _libga68_bounds (
filename=filename@entry=0x40c41a "standard.a68", lineno=lineno@entry=429,
index=<optimized out>, lower_bound=lower_bound@entry=1,
upper_bound=upper_bound@entry=16)
at ../../../gcc-17-20260426/libga68/ga68-error.c:128
#4 0x0000000000402ebb in STANDARD_digchar (x=<optimized out>)
at standard.a68:429
#5 choosedig (y=<optimized out>) at standard.a68:392
#6 STANDARD_subfixed (v=..., width=width@entry=28, after=after@entry=20)
at standard.a68:398
#7 0x0000000000403f2c in STANDARD_fixed (v=..., width=29,
after=after@entry=20) at standard.a68:200
#8 0x00000000004071b0 in STANDARD_float (v=..., width=30, after=20, exp=0)
at standard.a68:268
#9 0x0000000000401502 in __algol68_main () at nan.a68:9
#10 0x000000000040122b in main (argc=<optimized out>, argv=<optimized out>)
at ../../../gcc-17-20260426/libga68/libga68.c:41
When I run the same test on arm64, it completes successsfully, but the
output is wrong:
% ./a.out
If Algol 68 arithmetic follows IEEE 754, this should be NaN: +.0000000000000000000e0
------------------------------------------------------------------------
Numeric formatting by float() is inaccurate, even when the numbers are
small and exactly representable:
% cat powtwo032.a68
begin
real x, y, z;
int k;
string newline = "'n";
x := 2.0;
z := 1.0;
for k from 0 to 127
do
y := x ^ k;
puts ("2**(" + fixed(k, 4, 0) + ") = " +
float(y, 14, 10, 0) + " vs " +
float(z, 14, 10, 0) + " exact");
if y /= z then puts(" [DIFFER]") fi;
puts(newline);
z +:= z { exact power }
od
end
% ga68 powtwo032.a68 && ./a.out
2**( +0) = +1.000000000e0 vs +1.000000000e0 exact
2**( +1) = +2.000000000e0 vs +2.000000000e0 exact
2**( +2) = +4.000000000e0 vs +4.000000000e0 exact
2**( +3) = +8.000000000e0 vs +8.000000000e0 exact
2**( +4) = +1.599999904e1 vs +1.599999904e1 exact
2**( +5) = +3.199999809e1 vs +3.199999809e1 exact
2**( +6) = +6.399999618e1 vs +6.399999618e1 exact
2**( +7) = +1.280000095e2 vs +1.280000095e2 exact
2**( +8) = +2.560000181e2 vs +2.560000181e2 exact
2**( +9) = +5.120000362e2 vs +5.120000362e2 exact
2**( +10) = +1.024000048e3 vs +1.024000048e3 exact
...
Those numbers should be 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024,
...
Similar programs for long real and long long real also show incorrect
output conversions by float(), but happily, there are no output lines
with "DIFFER" --- thus, x ^ k is computed exactly correct in all three
floating-point precisions up to just below the overflow limit.
Some years ago, gcc switched its floating-point conversions between
binary and decimal values to use the MPFR library for correctly
rounded conversions. ga68 might well do the same.
It is worth noting here that the required digit counts for correct
round-trip binary <--> decimal conversions have been known since 1968,
but most programmers are unaware of them, and almost all programming
languages get those counts wrong. Here is what my hoc (high-order
calcular) program documents:
hoc> ?matula(nbits):
matula(nbits) is the number of decimal digits needed to ensure
correct round-trip conversion between binary and decimal of
floating-point data with nbits bits in the significand.
matula(P) is that number for this implementation of hoc.
In general, for decimal-to-binary conversion of d decimal digits
to p bits, we need to ensure that
10**d < 2**(p - 1),
and for binary-to-decimal conversion, we require that
2**p < 10**(d - 1).
See David W. Matula, ``In-and-out conversions'', Comm. ACM 11(1)
47--50, January 1968. CODEN CACMA2. ISSN 0001-0782.
DOI:10.1145/362851.362887
See also help_goldberg(), help_goldbergb(), and help_matulab().
hoc> ?goldberg
goldberg(ndecdig):
goldberg(ndecdig) is the number of bits needed to ensure
correct round-trip conversion between binary and decimal of
floating-point data with ndecdig decimal digits in the
significand.
In general, for decimal-to-binary conversion of d decimal digits
to p bits, we need to ensure that
10**d < 2**(p - 1),
and for binary-to-decimal conversion, we require that
2**p < 10**(d - 1).
See I. Bennett Goldberg, ``27 Bits Are Not Enough For 8-Digit
Accuracy'', Comm. ACM, 10(2) 105--106, February 1967. CODEN CACMA2.
ISSN 0001-0782. DOI:10.1145/363067.363112
See also help_goldbergb(), help_matula(), and help_matulab().
The two functions are simple to implement:
func matula (b) { return ceil( b * log10(2) + 1) }
func goldberg (d) { return ceil( d * log2(10) + 1 ) }
and obviously, if the compiler does not already do so, the constant
logarithms can be replaced by compile-time contants, if they have the
correct precisions.
In particular, for the IEEE 754 32-, 64-, 80-, and 128-bit binary
formats, the correct counts from matula() are 9, 17, 21, and 36
decimal digits.
----------------------------------------
P.S. I mentioned support for decimal floating-point arithmetic at the
start of this long message, and I suggest that it is something that
ga68 should give serious consideration to. Decimal arithmetic was
described in an IEEE Standard in 1987, and included in IEEE 754-2009;
it has been provided in hardware in two IBM CPU families since the
1990s. Apart from support for three new numeric types (in C,
_Decimal{32,64,128}), only at the point of code generation for +, -,
*, and / does a call to a software implementation of those decimal
operations, or on Power and S/390 (z-Series) [and perhaps soon on
RISC-V], a hardware instruction, need to be issued. I suspect that
for many programming language compilers, the work needed should be
less than a day or two, and the code changes relatively small. I have
already done the far larger work for the run-time numerical and I/O
library support of decimal floating-point arithmetic.
-------------------------------------------------------------------------------
- Nelson H. F. Beebe Tel: +1 801 581 5254 -
- University of Utah -
- Department of Mathematics, 110 LCB Internet e-mail: beebe@math.utah.edu -
- 155 S 1400 E RM 233 beebe@acm.org beebe@computer.org -
- Salt Lake City, UT 84112-0090, USA URL: https://www.math.utah.edu/~beebe -
-------------------------------------------------------------------------------
More information about the Algol68
mailing list