Hi Brad, Louis, Ben,
This not earth shattering. Just weird.
On Tue, 5 Mar 2019, Brad Chamberlain wrote:
Most of Chapel's math functions do simply wrap C's 'math.h' and 'libm'
(glibc or otherwise depending on what back-end compiler you're using I
suppose). The cbrt() function itself is just a pass-through to C's
cbrt() so it seems it should not behave differently.
That's what I thought. Maybe I have been looking at IEEE 754 encodings for
too long!
According to the C program I just wrote up correctly, it looks like the
result should be:
0x1.0p-348
not:
0x1.0p-340
Sorry. Fat fingers. Typo in my email. Was testing several, 0x1.0p-1020,
0x1.0p-1032 and 0x1.0p-1044, where the last two are subnormal numbers,
during some exhaustive tests of my Chapel-ized and generically programmed
cbrt.
and using this value seems to work for me:
I would expect so. But see below.
if cbrt(0x1.0p-1044) == 0x1.0p-348 then
writeln("Matches");
else
writeln("Doesn't match");
Very weird. Maybe things work strangely when you are upside down relative
to the equator! I have a GCC backend for now.
Note that in the following, everything is real(64) and
x = 0x1.0p-1044;
{
const chapel = cbrt(x);
const damian = fpCubeRoot(x);
const exact = 0x1.0p-348;
const ypsilon = 0x1.0p+52; // == 1.0/epsilon == 1/0x1.0p-52;
const brad = if chapel == exact then "Matches" else "Does Not Match";
var y : x.type;
writeln("Brad's Test ", brad);
y = ypsilon * abs(exact - damian)/exact;
writeln("Damian's cubeRoot(", x, ") rel.error = ", y);
y = ypsilon * abs(exact - chapel)/exact;
writeln("Chapel's cubeRoot(", x, ") rel.error = ", y);
}
The output is
Brad's Test Does Not Match
Damian's cubeRoot(5.30499e-315) rel.error = 0.0
Chapel's cubeRoot(5.30499e-315) rel.error = 0.5
The relative error (as a multiple of epsilon) should be zero if the
computation is exact or 0.5 if rounded correctly but still approximate.
Weird. As I said.
Note that Damian's fpCubeRoot, when written in C, agrees with the system's
own math's library. There is nothing revolutionary in my routine except
that it is generic, i.e. the 32-bit version is the same as the 64-bit one.
It is almost identical with the C routine and produces the same results
when coded in C, well at least for the 60 million test cases in my tests.
It uses similar well identical formulae so I would expect that.
Not worth pursuing for now.
Louis, Ben, saw your notes about what is included. I was just noting it
for future reference. I have my own, less LOAC (lines-of-assembler-code)
version of nextafter written in C and a generic version for any sized
real(w) in Chapel.
Writing low level routines in Chapel, with the recent update for real(w)
expressions in params, now realizes Chapel's potential for writing truly
generic code for scientific computations for floating point precision of
any type!
// end-of-advertisement (oops, where did the advertisement start??)
Work hard - Damian
Pacific Engineering Systems International, 277-279 Broadway, Glebe NSW 2037
Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
Views & opinions here are mine and not those of any past or present employer
_______________________________________________
Chapel-developers mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/chapel-developers
