Making readable modifyable code.

[jahboater]

The Pi3 when running in 64-bit mode will have available a 113-bit mantissa giving 34.016 digits.

??? Where is this 128bit floating point support?
Yes it has 128bit NEON registers but that ranges from a 16x 8bit vector to a 2x 64bit vector. AFAIK it isn't a single 128bit value for anything other than load/store.

Its in software.
If you use "long double" in C you automatically get 64-bit reals in 32-bit mode (same as double), and 128-bit reals in 64-bit mode.
Underneath its the _Float128 type. There is also _Complex128 to go with it. Apparently 64-bit hardware/OS are needed to implement it at a decent speed.

Its transparent to the user and works well, albeit a little bigger and presumably slower but I have never tried benchmarking it. So for example:-

#include <math.h>
long double x = 2.0;
long double y = sqrtl( x );
printf( "%.34Lg\n", y );

Prints:
1.414213562373095048801688724209698

which is correct to 34 places according to this
https://en.wikipedia.org/wiki/Square_root_of_2

Just use "-lm" no special libraries are needed.

[Heater]

jahboater,

Hmmm...that does not work on my amd64 Debian desktop even if long double is 128 bits:

Code: Select all

$ cat longdouble.c 
#include <stdio.h>
#include <math.h>

int main (int argc, char* argvi[]) {
    long double x = 2.0;
    long double y = sqrtl(x);
    printf("long double is %d bits.\n", sizeof(x) * 8);
    printf("Expect: 1.414213562373095048801688724209698\n");
    printf("Got:    %.34Lg\n", y);

    return 0;
}

$ gcc -o longdouble longdouble.c -lm 
$ file longdouble
longdouble: ELF 64-bit LSB executable, x86-64
$ ./longdouble 
long double is 128 bits.
Expect: 1.414213562373095048801688724209698
Got:    1.414213562373095048763788073031833
 

Sadly I don't have a Pi to hand to test on.

[[email protected]]

jahboater,

Hmmm...that does not work on my amd64 Debian desktop even if long double is 128 bits:

Code: Select all

$ cat longdouble.c 
#include <stdio.h>
#include <math.h>

int main (int argc, char* argv) {
    long double x = 2.0;
    long double y = sqrtl(x);
    printf("long double is %d bits.\n", sizeof(x) * 8);
    printf("Expect: 1.414213562373095048801688724209698\n");
    printf("Got:    %.34Lg\n", y);

    return 0;
}

$ gcc -o longdouble longdouble.c -lm 
$ file longdouble
longdouble: ELF 64-bit LSB executable, x86-64
$ ./longdouble 
long double is 128 bits.
Expect: 1.414213562373095048801688724209698
Got:    1.414213562373095048763788073031833
 

Sadly I don't have a Pi to hand to test on.

FWIW: I get the same on my 32-bit x86 desktop (Debian Wheezy) where an ld is 96 bits:

Code: Select all

long double is 96 bits.
Expect: 1.414213562373095048801688724209698
Got:    1.414213562373095048763788073031833

And on a Pi:

Code: Select all

long double is 64 bits.
Expect: 1.414213562373095048801688724209698
Got:    1.414213562373095145474621858738828

where ld is only 64-bits...

Also the same on a 128 bit Debian Jessie system where ld is 128 bits -

edit to add the 128-bit system ends the same 833 - the Pi ending is different.

-Gordon

[jahboater]

Don't use sizeof()!

Sadly "long double" on Intel/amd gets you the 80-bit x87 reals.

To get 128 bit reals on Intel 64-bit you have to use "_Float128" which is a real pain (and libquadmath https://gcc.gnu.org/onlinedocs/libquadmath/).

A good value to print is LDBL_DIG for the actual precision in decimal digits, its 18 on Intel and 33 on 64-bit arm (and 15 for DBL_DIG).

Note sizeof() tells you nothing much. On Intel 32 bit sizeof(long double) returns 12 and on 64-bit sizeof(long double) returns 16. They are both 80-bits taking 10 bytes, so there are some unused bytes.
This is because the ABI needs 80-bit reals on 16 byte boundaries for 64-bit mode.

Perhaps try:
printf("Got: %.*Lg\n", LDBL_DIG, y);
to print to the maximum reasonable precision.

[bensimmo]

Been interesting reading stuff I don't know about .

What about GPU (videocore wrt R-Pi) targeted code?

[Heater]

Grrr....can't get __float128 to produce the right answer on Intel/amd64/GCC/Debian/Jessie:

Code: Select all

$ cat float128.cpp 
#include <quadmath.h>
#include <stdio.h>
#include <math.h>
#include <iostream>     // std::cout
#include <iomanip>
#include <limits>       // std::numeric_limits


using namespace std;

int main (int argc, char* argv[]) {
    __float128 x = 2.0;
    __float128 y = sqrtl(x);
    char buf[128];
    int width = 46;
    cout << "__float128 is " << sizeof(x) * 8 << " bits." << endl;
    cout << numeric_limits<__float128>::digits10 << endl;

    cout << "Expect: +1.414213562373095048801688724209698e+00" << endl;
    quadmath_snprintf (buf, sizeof buf, "%+-*.33Qe", width, y);
    cout << "Got:    " << buf << endl;

    if (y == 1.414213562373095048801688724209698)
    {
       cout << "Pass." << endl;
    }
    else
    {
       cout << "Fail." << endl;
    }
    return 0;
}

$ g++ -Wall -o float128 ft128.cpp -lm -lquadmath
$ ./float128 
__float128 is 128 bits.
0
Expect: +1.414213562373095048801688724209698e+00
Got:    +1.414213562373095048763788073031833e+00      
Fail.

[jahboater]

And on a Pi:

Code: Select all

long double is 64 bits.
Expect: 1.414213562373095048801688724209698
Got:    1.414213562373095145474621858738828

where ld is only 64-bits...

Your Pi is running a 32-bit OS where "long double" is the same as "double".
You need to run a 64-bit OS for long double to be truly 128-bits.

[jahboater]

Grrr....can't get __float128 to produce the right answer on Intel/amd64/GCC/Debian/Jessie:

I told you it was a right PITA.
I gave up and only use the ARM "long double" implementation because its so easy and portable.

On arm its all transparent - just use long double and it all works. sqrtl() and suffix L (= 2.0L) just as normal. That all works too on Intel for 80-bits.

For _Float128 on Intel, the suffix is Q .
So you need sqrtq() and probably "= 2.0Q".
You will also need quadmath_snprintf() to print the result.
Told you it was a pain.

Is it available for C++ ? if you cant get it to work, try C.

[jahboater]

Code: Select all

   if (y == 1.414213562373095048801688724209698)
    {
       cout << "Pass." << endl;
    }
    else
    {
       cout << "Fail." << endl;
    }

I think the usual rules about not comparing floats for equality hold - especially for an irrational number at the limits of the available precision

[Heater]

jahboater,

Ah, thanks. I already had the quadmath_snprintf() in the above example. Forgot to change sqrtq() to sqrtl() though.

Code: Select all

$ ./float128 
__float128 is 128 bits.
0
Expect: +1.414213562373095048801688724209698e+00
Got:    +1.414213562373095048801688724209698e+00      
Fail.

[jahboater]

jahboater,

Ah, thanks. I already had the quadmath_snprintf() in the above example. Forgot to change sqrtq() to sqrtl() though.

Code: Select all

$ ./float128 
__float128 is 128 bits.
0
Expect: +1.414213562373095048801688724209698e+00
Got:    +1.414213562373095048801688724209698e+00      
Fail.

Well done!
Now you have to think of a practical use for that much precision...
The range by the way is 10e +/- 4932.
So factorial 1754 is possible, compared to factorial 170 which is the best 64-bit doubles can do.

[Heater]

jahboater,

I think the usual rules about not comparing floats for equality hold

Yes but I was hoping that as we know the result we want we could set up and exact equality.

Like so (Which now passes):

Code: Select all

$ cat float128.cpp       
#include <quadmath.h>
#include <stdio.h>
#include <math.h>
#include <iostream>     // std::cout
#include <iomanip>
#include <limits>       // std::numeric_limits

using namespace std;

int main (int argc, char* argv[]) {
    __float128 x = 2.0;
    __float128 y = sqrtq(x);
    char buf[128];
    int width = 46;
    cout << "__float128 is " << sizeof(x) * 8 << " bits." << endl;
    cout << numeric_limits<__float128>::digits10 << endl;

    cout << "Expect: +1.414213562373095048801688724209698176940240827882212721853306632542406594266193309294976643286645412445068359375" << endl;
    quadmath_snprintf (buf, sizeof buf, "%+-*.200Qe", width, y);
    cout << "Got:    " << buf << endl;

    if (y == 1.414213562373095048801688724209698176940240827882212721853306632542406594266193309294976643286645412445068359375Q) 
    {
       cout << "Pass." << endl;
    }
    else
    {
       cout << "Fail." << endl;
    }
    return 0;
}

$ g++ -Wall -o float128 float128.cpp -lquadmath
$ ./float128 
__float128 is 128 bits.
0
Expect: +1.414213562373095048801688724209698176940240827882212721853306632542406594266193309294976643286645412445068359375
Got:    +1.4142135623730950488016887242096981769402408278822127218533066325424065942661933092949766432866454124450683593750000000000000
Pass.

[jahboater]

That's impressive!

jahboater,

I think the usual rules about not comparing floats for equality hold

Yes but I was hoping that as we know the result we want we could set up an exact equality.

The actual value to 65 places is

Code: Select all

1.41421356237309504880168872420969807856967187537694807317667973799Q

But for the test anything after about 34 places will be ignored.

I wonder if this helps (the exact bit pattern from %La)

Code: Select all

if( y == 0x1.6a09e667f3bcc908b2fb1366ea95p+0Q ) 

I can cope with comparing a float to say 128.0 which will be representable in any binary format, but an irrational? Gives me the shivers ...

[tufty]

Of course your example is the equivalent of:

Code: Select all

> (/ 1 3)
1/3

It shows how one can side step the problem of number representation by not actually returning a number as a result!

No, 1/3 most certainly is a number. It's an exact rational number, to be precise.

Scheme provides for (at least) several numeric types : Complex numbers, reals, rationals, and integers. There are no theoretical limits to these; integers can have an arbitrary number of bits, as can reals (floating or fixed point numbers). There's also the concept of exactness, which is orthogonal to numeric type.

As it stands, there's no defined numeric type in the standard for "symbolic" numbers - pi, e, √2, etc, although an implementation could implement something like that. So, as it stands, the result of

Code: Select all

(let ((x (sqrt 2))) (* 2 2))

is, at least on the scheme implementations I have easy access to here, an inexact real, the result of multiplying the floating point result of (sqrt 2) by itself. However, an implementation might define (sqrt 2) as an exact "symbolic" number, and thus the result might well be an exact integer, as it should be. It rather depends on how mathematically correct you want to be.

My main gripe with "machine numbers" as handled by most languages isn't mathematical completeness, however - it's that they usually have no indication of exactness, and usually very little handling of the oddball edge cases.

No, the copy of Chez I was using doesn't run on the Pi (although it's been open sourced, https://github.com/cisco/ChezScheme, and has ARM support, so it should be relatively trivial to build it, in as much as building a scheme system is ever trivial), but most of the other schemes build for ARM too. The only one I have on this machine that treats 1/3 as 0.333333… (which isn't incorrect, just - unexpected) is chicken.

[jahboater]

My main gripe with "machine numbers" as handled by most languages isn't mathematical completeness, however - it's that they usually have no indication of exactness, and usually very little handling of the oddball edge cases.

I'm not quite sure what you mean, but in C yes exactness is hard to deal with.

1.0/3.0 will of course raise FE_INEXACT which is fine.

But so will things like log10(100), cbrt(27), log2(2), pow(4.0,0.5) because internally they might raise FE_INEXACT even though the result is definitely exact.

On the other hand sqrt(4.0) will not because it will use a single hardware instruction which will always set FE_INEXACT correctly.

And then sqrt(2) * 0.0 will raise FE_INEXACT because the flag is sticky.

Another difficulty: things like strtod( "0.3333333333333333333333333333333333333333" ) are happy even though some digits have been discarded.

Yes, very hard to deal with.

[[email protected]]

Yes, very hard to deal with.

But... computers are always right!



-Gordon

[jahboater]

Yes, very hard to deal with.

But... computers are always right!



-Gordon

I wish mine was ...

[tufty]

I'm not quite sure what you mean

Well, exactness is one thing, it's hairy and nasty and not helped by floating point. I read a very good paper on hardware accelerated ULP based representations, which allow one to keep a handle on the current (potential) inexactitude in a computation. As per usual, I've lost the link.

The main thing, though is edge cases. And there's loads and loads of them hiding in the most trivial program - collatz, for example, will overflow a 32 bit counter for results well under 32 bits, and what happens when numbers overflow is often undefined behaviour in certain cases for pretty much anything based on C. For example...

Code: Select all

uint32_t count_1; /* overflow behaviour is defined */
int32_t count_2; /* Overflow behaviour is nose demons */

It's horribly inconsistent. And so we end up with programs that break when compiled on different platforms, or using a different compiler, or without some special incantation of compiler flags. Given that computers are supposed to be deterministic calculating machines, it seems insane to me that we leave much of calculating down to little more than luck.

[jahboater]

The main thing, though is edge cases. And there's loads and loads of them hiding in the most trivial program - collatz, for example, will overflow a 32 bit counter for results well under 32 bits, and what happens when numbers overflow is often undefined behaviour in certain cases for pretty much anything based on C. For example...

Code: Select all

uint32_t count_1; /* overflow behaviour is defined */
int32_t count_2; /* Overflow behaviour is nose demons */

It's horribly inconsistent. And so we end up with programs that break when compiled on different platforms, or using a different compiler, or without some special incantation of compiler flags. Given that computers are supposed to be deterministic calculating machines, it seems insane to me that we leave much of calculating down to little more than luck.

Yes, indeed.
You can just do your best, make sure there is absolutely no UB in your program - but how many programmers fully understand aliasing, overflow rules etc etc. Not many I suspect. Even worse is the behaviour of compilers that assume overflow cannot happen in a correct program and do strange optimizations based on that.

-fwrapv makes signed integer overflow defined by the way.

[Heater]

Wow, lots of interesting debate here. Which I will have to think about after some sleep.

For many things I am happy with the crude data types of C.

However, it's very annoying that so much of that is "implementation dependent". That means you can never be sure your program works on all compilers and architectures. WTF?

I also find it insane that C will happily overflow an integer or whatever, produce an undefined result and then happily continue producing garbage.

They tried to fix that kind of problem with Ada. But it seems programmers roundly rejected the idea of such a "fussy" system.

[scruss]

My main gripe with "machine numbers" as handled by most languages isn't mathematical completeness, however - it's that they usually have no indication of exactness, and usually very little handling of the oddball edge cases.

I'd disagree about the edge cases. For the intended audience — general engineering computation — IEEE 754 does a remarkably good job of hiding the limitations inherent in trying to wedge the representation of a real number into a limited number of binary bits. It does so fairly quickly, and without needing acres of silicon. While there is a real science to knowing when you shouldn't be using it, 754 has brought reliable FP with well understood, reliable limitations to the masses. If you even know what Scheme is, you're pretty much not the intended user.

If you want to see bad floating point, look to the implementations in the BASICs of many 1980s microcomputers; edge cases a go-go!

I really must do something with the stick of AMD 4-bit ALUs I picked up surplus. That's some real adding up right there!

[tufty]

A scheme version of the collatz problem above would be great. Hint, hint...

A quick implementation. This is all tested on a mac, not on a Pi, under Chez scheme 9.4

Here's the trivial collatz step. Should be pretty readable, the only thing that might throw non-schemers is the use of `let name ....` which defines a named recursion point. so the `(let loop ((x x) (count 0)) …) block is actually a recursive function with two parameters, x and count. That explanation probably doesn't help, either :) This is tail recursive code, so no, it doesn't blow the stack. Performance-wise, though, it kinda blows chunks. 80 seconds on my machine for the 1-10M test (including compile time and function call overhead)

Code: Select all

(define (collatz x)
  (let loop ((x x) (count 0))
    (cond
     ((= 1 x)   count)
     ((even? x) (loop (/ x 2) (+ count 1)))
     ((odd? x)  (loop (+ (* x 3) 1) (+ count 1))))))

As we know we're dealing with integers, we can rewrite this in exactly the same form using the "fixnum" operators, as follows :

Code: Select all

(define (collatz-fx x)
  (let loop ((x x) (count 0))
    (cond
     ((fx=? 1 x)   count)
     ((fxeven? x) (loop (fxdiv x 2) (fx+ count 1)))
     ((fxodd? x)  (loop (fx+ (fx* x 3) 1) (fx+ count 1))))))

Basically looks the same, but significantly faster, 29 seconds for 10M, still including compile and function call overhead. I'm sure we can do better than that. fxdiv, for example, is number-theoretically correct. But as we're dividing an integer that is known to be positive and by a power of two, we can just bitshift it.

Code: Select all

(define (collatz-fxshift x)
  (let loop ((x x) (count 0))
    (cond
     ((fx=? 1 x)   count)
     ((fxeven? x) (loop (fxarithmetic-shift-right x 1) (fx+ count 1)))
     ((fxodd? x)  (loop (fx+ (fx* x 3) 1) (fx+ count 1))))))

And we're down to 8 seconds. The calling code, FWIW, looks like this, and uses fixnum operations throughout. It's ugly and I could do a lot better if I could be bothered.

Code: Select all

(define (collatz-test fn)
  (let ((start-time (current-time)))
    (let loop ((x 1) (max-x 0) (max-count 0))
      (if (fx>? x 10000000)
          (display (format "max : ~a, count: ~a, time ~a seconds\n" max-x max-count (time-difference (current-time) start-time)))
          (let ((count (fn x)))
            (if (fx>? count max-count)
                (loop (fx+ x 1) x count)
                (loop (fx+ x 1) max-x max-count)))))))

(display "Plain Collatz\n")
(collatz-test collatz)
(display "\nExplicitly fixnum\n")
(collatz-test collatz-fx)
(display "\nExplicitly fixnum with shifts for divide\n") 
(collatz-test collatz-fxshift)

At this point, memoising could be a good next step. We have fixnum vectors, so why not? This will, of course, push a lot of overhead onto the vector manipulation code, it might not work. Unfortunately, I have a lot of stuff to do and only one day off, so it will have to wait.

[Heater]

@jahboater

But for the test anything after about 34 places will be ignored.

Yep. It actually need 36 digits after the point to pass.

Anyway, I give up with that float128 thing.

@tufty

1/3 most certainly is a number. It's an exact rational number, to be precise.

I like to think so.

Problem is you can't map "1/3" onto reality without more work. Namely a division. For example how are you going plot 1/3 on the number line? Or how are you going to output it to some actuator?

Now, where is that scheme version of the collatz algorithm?

@scruss

I really must do something with the stick of AMD 4-bit ALUs I picked up surplus. That's some real adding up right there!

You mean those AMD 29K bit slice ALUs?

Awesome, how did you get hold of those?

[tufty]

Problem is you can't map "1/3" onto reality without more work. Namely a division.

Keeping things exact until it's absolutely necessary to approximate them is good practice if you care about the correctness of your result. Perhaps my bias, which comes from, amongst other things, a mixture of financials and airworthiness, shows here.

For example how are you going plot 1/3 on the number line?

I'd put it as close to exactly one third of the way from zero to one as I could.

Or how are you going to output it to some actuator?

That would rather depend on how the actuator works. If we were talking about, for example, a stepper motor with (grabs one from the pile) 48 steps per revolution, I'd send it 48/3 = 16 steps round. Not - ummm 0.3333333 * 48 = 15.9999984 steps. If it were the next one I grab, (400 steps per revolution) that would result in an inexact conversion, but I still wouldn't have lost anything over using the inexact representation. Plus, if I were clever, I'd be using analog drivers; 3 fits really nicely into 3 phase calculations...

Now, where is that scheme version of the collatz algorithm? :)

I think you'll find it just above the post I am replying to :)

[tufty]

Also, 2901s? If so, make yourself an Atari vector machine!