On Wed, 6 Mar 2019, Damian McGuckin wrote:
> On Tue, 5 Mar 2019, Brad Chamberlain wrote:
> > and using this value seems to work for me:
> This not earth shattering. Just weird.
The value seems to work when it's available to cbrt() as a param (or
literal as in Brad's example):
fortytwo@magrathea:~/src/chapel (master)$ cat realweird.chpl
param x = 0x1.0p-1044;
var y = x;
writeln("chapel var : ", cbrt(y) == 0x1.0p-348);
writeln("chapel param : ", cbrt(x) == 0x1.0p-348);
writeln("chapel literal: ", cbrt(0x1p-1044) == 0x1.0p-348);
fortytwo@magrathea:~/src/chapel (master)$ ./realweird -nl 1
chapel var : false
chapel param : true
chapel literal: true
I get the same results when compiling with -O, or with --fast, or with
neither.
Paul
>
> 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
>
_______________________________________________
Chapel-developers mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/chapel-developers
