Hi all-

I'm a newcomer to the D language, but am quite impressed with it.  I've run into an issue, though, in a little learning project.  I'm implementing a "double-double" class based on a well-known trick using standard floating-point arithmetic.  (See, for example, http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic and links).

The techniques used to get this kind of code to work correctly, rely on subtle properties of IEEE floating point numbers.  Mixing precisions or optimisations in floating point code can destroy the results.

For example, the appended code produces the following output when compiled (DMD32 D Compiler v2.063.2 under WinXP/cygwin) with no optimization:

immutable(pair)(1.1, -2.03288e-20)
pair(1, 0.1)
pair(1.1, -8.32667e-17)

and the following results when compiled with optimization (-O):

immutable(pair)(1.1, -2.03288e-20)
pair(1, 0.1)
pair(1.1, 0)

The desired result would be:

immutable(pair)(1.1, -8.32667e-17)
pair(1, 0.1)
pair(1.1, -8.32667e-17)

That is: without optimization the run-time "normalization" is correct.  With optimization it is broken.  That is pretty easy to work around by simply compiling the relevant library without optimization.  (Though it would be nice to have, for example, pragmas to mark some functions as "delicate" or "non-optimizable".)  A bigger issue is that the compile-time normalization call gives the 'wrong' answer consistently with or without optimization.  One would expect that evaluating a pure function at run-time or compile-time would give the same result...

Cheers,
--Apollo

import std.stdio;
struct pair { double hi, lo; }
pair normalize(pair q)
{
  double h = q.hi + q.lo;
  double l = q.lo + (q.hi - h);
  return pair(h,l);
}
void main()
{
  immutable static pair spn = normalize(pair(1.0,0.1));
  writeln(spn);
  writeln(pair(1.0,0.1));
  writeln(normalize(pair(1.0,0.1)));
}

Apollo Hogan:

import std.stdio;

struct Pair { double hi, lo; }

Pair normalize(Pair q) {
    immutable h = q.hi + q.lo;
    immutable l = q.lo + (q.hi - h);
    return Pair(h, l);
}
void main() {
    immutable static Pair spn = Pair(1.0, 0.1).normalize;
    spn.writeln;
    Pair(1.0, 0.1).writeln;
    Pair(1.0, 0.1).normalize.writeln;
}

...>ldmd2 -O -run temp.d
immutable(Pair)(1.1, -2.03288e-020)
Pair(1, 0.1)
Pair(1.1, -8.32667e-017)

...>ldmd2 -run temp.d
immutable(Pair)(1.1, -2.03288e-020)
Pair(1, 0.1)
Pair(1.1, -8.32667e-017)

ldc2 0.11.0

Bye,
bearophile

On 10/23/2013 8:44 AM, Apollo Hogan wrote:
> That is: without optimization the run-time "normalization" is correct.  With
> optimization it is broken.  That is pretty easy to work around by simply
> compiling the relevant library without optimization.  (Though it would be nice
> to have, for example, pragmas to mark some functions as "delicate" or
> "non-optimizable".)  A bigger issue is that the compile-time normalization call
> gives the 'wrong' answer consistently with or without optimization.  One would
> expect that evaluating a pure function at run-time or compile-time would give
> the same result...

A D compiler is allowed to compute floating point results at arbitrarily large precision - the storage size (float, double, real) only specify the minimum precision.

This behavior is fairly deeply embedded into the front end, optimizer, and various back ends.

To precisely control maximum precision, I suggest using inline assembler to use the exact sequence of instructions needed for double-double operations.

Hi Bearophile-

Interesting.  Looks like the run-time calculation in ldmd2 works fine.

The compile-time computation in both my and your examples looks like it is being done in 80-bit arithmetic.

Thanks,
--Apollo

On Wednesday, 23 October 2013 at 16:15:56 UTC, Walter Bright wrote:
> A D compiler is allowed to compute floating point results at arbitrarily large precision - the storage size (float, double, real) only specify the minimum precision.
>
> This behavior is fairly deeply embedded into the front end, optimizer, and various back ends.

I know we've had this topic before, but just for the record, I'm still not sold on the idea of allowing CTFE to yield different results than runtime execution.

David

On 10/23/2013 9:22 AM, David Nadlinger wrote:
> On Wednesday, 23 October 2013 at 16:15:56 UTC, Walter Bright wrote:
>> A D compiler is allowed to compute floating point results at arbitrarily large
>> precision - the storage size (float, double, real) only specify the minimum
>> precision.
>>
>> This behavior is fairly deeply embedded into the front end, optimizer, and
>> various back ends.
>
> I know we've had this topic before, but just for the record, I'm still not sold
> on the idea of allowing CTFE to yield different results than runtime execution.

Java initially tried to enforce a maximum precision, and it was a major disaster for them. If I have been unable to convince you, I suggest reviewing that case history.

Back when I designed and built digital electronics boards, it was beaten into my skull that chips always get faster, never slower, and the slower parts routinely became unavailable. This means that the circuits got designed with maximum propagation delays in mind, and with a minimum delay of 0. Then, when they work with a slow part, they'll still work if you swap in a faster one.

FP precision is the same concept. Swap in more precision, and your correctly designed algorithm will still work.

On Wednesday, 23 October 2013 at 16:50:52 UTC, Walter Bright wrote:
> On 10/23/2013 9:22 AM, David Nadlinger wrote:
>> On Wednesday, 23 October 2013 at 16:15:56 UTC, Walter Bright wrote:
>>> A D compiler is allowed to compute floating point results at arbitrarily large
>>> precision - the storage size (float, double, real) only specify the minimum
>>> precision.
>>>
>>> This behavior is fairly deeply embedded into the front end, optimizer, and
>>> various back ends.
>>
>> I know we've had this topic before, but just for the record, I'm still not sold
>> on the idea of allowing CTFE to yield different results than runtime execution.
>
> Java initially tried to enforce a maximum precision, and it was a major disaster for them. If I have been unable to convince you, I suggest reviewing that case history.
>
> Back when I designed and built digital electronics boards, it was beaten into my skull that chips always get faster, never slower, and the slower parts routinely became unavailable. This means that the circuits got designed with maximum propagation delays in mind, and with a minimum delay of 0. Then, when they work with a slow part, they'll still work if you swap in a faster one.
>
> FP precision is the same concept. Swap in more precision, and your correctly designed algorithm will still work.

There are a couple of points here:

- it seems that whatever the semantics of floating-point arithmetic, they should be the same at compile-time as at run-time.

- I agree that the majority of floating point code is only improved by increasing the working precision.  (If we don't worry about reproducibility across compilers/machines/etc.)  The "real" data-type seems to be designed exactly for this: use "real" in numerical code and the compiler will give you a good answer at the highest performant precision.  However there _are_ cases where it can be very useful to have precise control of the precision that one is using.  Implementing double-double or quadruple-double data types is an example here.   Viewing D as a _systems_ language, I'd like to have the ability to just have it do what I ask (and being forced to go to assembler seems unreasonable...)

Anyway, thanks for the replies.  I guess I've got to go off and design the brand new D^^2 language to conform to my whims now.

Cheers,
--Apollo

On 10/23/2013 11:39 AM, Apollo Hogan wrote:
> There are a couple of points here:
>
> - it seems that whatever the semantics of floating-point arithmetic, they should
> be the same at compile-time as at run-time.

It's not very practical, especially considering that the compile time environment may be not at all the same as the runtime one.


> - I agree that the majority of floating point code is only improved by
> increasing the working precision.  (If we don't worry about reproducibility
> across compilers/machines/etc.)

As I mentioned earlier, Java initially mandated that floating point results be exactly reproducible in multiple environments. This was a fine idea, and turned out to be completely unworkable in practice.


> The "real" data-type seems to be designed
> exactly for this: use "real" in numerical code and the compiler will give you a
> good answer at the highest performant precision.  However there _are_ cases
> where it can be very useful to have precise control of the precision that one is
> using.  Implementing double-double or quadruple-double data types is an example
> here.

It's not that bad. You can also force a reduction in precision by calling a function like this:

    double identity(double d) { return d; }

and ensuring (via separate compilation) that the compiler cannot inline calls to identity().

> Viewing D as a _systems_ language, I'd like to have the ability to just
> have it do what I ask (and being forced to go to assembler seems
> unreasonable...)

Perhaps, but I think you are treading on implementation defined behavior here for most languages, and will be hard pressed to find a language that *guarantees* the loss of precision you desire, even though it may deliver that behavior on various platforms with various compiler switches.


> Anyway, thanks for the replies.  I guess I've got to go off and design the brand
> new D^^2 language to conform to my whims now.

Join the club! :-)

On Wednesday, 23 October 2013 at 19:10:22 UTC, Walter Bright wrote:
> On 10/23/2013 11:39 AM, Apollo Hogan wrote:
>> There are a couple of points here:
>>
>> - it seems that whatever the semantics of floating-point arithmetic, they should
>> be the same at compile-time as at run-time.
>
> It's not very practical, especially considering that the compile time environment may be not at all the same as the runtime one.

Understood, but it certainly was a surprising result to me that compiling and running the program on the same platform I got different results for a static vs. non-static variable initialization... (My version of PI as 3.14159265358979311594779789241 was a bit confusing...)

> It's not that bad. You can also force a reduction in precision by calling a function like this:
>
>     double identity(double d) { return d; }
>
> and ensuring (via separate compilation) that the compiler cannot inline calls to identity().

Thanks, a useful trick.  It at least lets me confound the optimizer a bit.  (Though doesn't help with the compile vs. run headache.  Though this seems to be workaroundable by using a static constructor.  Yeah, I'm a noob.)

Thanks for the replies,
--Apollo

On 10/23/13 1:59 PM, Apollo Hogan wrote:
> Thanks, a useful trick.  It at least lets me confound the optimizer a
> bit.  (Though doesn't help with the compile vs. run headache.  Though
> this seems to be workaroundable by using a static constructor.  Yeah,
> I'm a noob.)

+1 for "workaroundable".

Andrei