inaccurate mathematical constants

[ian-bruce]

<cl_platform.h> defines some standard mathematical constants, in both single- and double-precision formats:

$ ( sed -n '/stdint/,$p' | grep 'CL_M_' ) <cl_platform.h

#define  CL_M_E             2.718281828459045090796
#define  CL_M_LOG2E         1.442695040888963387005
#define  CL_M_LOG10E        0.434294481903251816668
#define  CL_M_LN2           0.693147180559945286227
#define  CL_M_LN10          2.302585092994045901094
#define  CL_M_PI            3.141592653589793115998
#define  CL_M_PI_2          1.570796326794896557999
#define  CL_M_PI_4          0.785398163397448278999
#define  CL_M_1_PI          0.318309886183790691216
#define  CL_M_2_PI          0.636619772367581382433
#define  CL_M_2_SQRTPI      1.128379167095512558561
#define  CL_M_SQRT2         1.414213562373095145475
#define  CL_M_SQRT1_2       0.707106781186547572737
#define  CL_M_E_F           2.71828174591064f
#define  CL_M_LOG2E_F       1.44269502162933f
#define  CL_M_LOG10E_F      0.43429449200630f
#define  CL_M_LN2_F         0.69314718246460f
#define  CL_M_LN10_F        2.30258512496948f
#define  CL_M_PI_F          3.14159274101257f
#define  CL_M_PI_2_F        1.57079637050629f
#define  CL_M_PI_4_F        0.78539818525314f
#define  CL_M_1_PI_F        0.31830987334251f
#define  CL_M_2_PI_F        0.63661974668503f
#define  CL_M_2_SQRTPI_F    1.12837922573090f
#define  CL_M_SQRT2_F       1.41421353816986f
#define  CL_M_SQRT1_2_F     0.70710676908493f

Bizarrely, the single- and double-precision versions DO NOT MATCH BEYOND THE SEVENTH DIGIT:

$ ( sed -n '/stdint/,$p' | grep 'CL_M_' ) <cl_platform.h |
    ( sort | sed -r -e 's/ 0\./0./' -e 's/(.{36})(.*)/\1 \2/' )

#define  CL_M_1_PI         0.3183098 86183790691216
#define  CL_M_1_PI_F       0.3183098 7334251f

#define  CL_M_2_PI         0.6366197 72367581382433
#define  CL_M_2_PI_F       0.6366197 4668503f

#define  CL_M_2_SQRTPI      1.128379 167095512558561
#define  CL_M_2_SQRTPI_F    1.128379 22573090f

#define  CL_M_E             2.718281 828459045090796
#define  CL_M_E_F           2.718281 74591064f

#define  CL_M_LN10          2.302585 092994045901094
#define  CL_M_LN10_F        2.302585 12496948f

#define  CL_M_LN2          0.6931471 80559945286227
#define  CL_M_LN2_F        0.6931471 8246460f

#define  CL_M_LOG10E       0.4342944 81903251816668
#define  CL_M_LOG10E_F     0.4342944 9200630f

#define  CL_M_LOG2E         1.442695 040888963387005
#define  CL_M_LOG2E_F       1.442695 02162933f

#define  CL_M_PI_2          1.570796 326794896557999
#define  CL_M_PI_2_F        1.570796 37050629f

#define  CL_M_PI            3.141592 653589793115998
#define  CL_M_PI_F          3.141592 74101257f    /* re-ordered */

#define  CL_M_PI_4         0.7853981 63397448278999
#define  CL_M_PI_4_F       0.7853981 8525314f

#define  CL_M_SQRT1_2      0.7071067 81186547572737
#define  CL_M_SQRT1_2_F    0.7071067 6908493f

#define  CL_M_SQRT2         1.414213 562373095145475
#define  CL_M_SQRT2_F       1.414213 53816986f

THEY CAN'T BOTH BE CORRECT!

The situation may not be as bad as it looks -- the double-precision values may be entirely correct, and since eight decimal digits would be sufficient to represent the 24-bit mantissa of the single-precision values, their inaccuracy may be very small. But it looks quite odd, and is likely to cause problems if these definitions are ever copied and reused somewhere else.

Surely this should be investigated?

[b-sumner]

There is nothing wrong. Unfortunately we can't use hex floating point format because some silly compilers don't accept it.

The closest single precision representation of PI is 3.1415927410125732421875f AKA 0x1.921fb6p+1f

The closest double precision representation of PI is 3.141592653589793115997963468544185161590576171875 AKA 0x1.921fb54442d18p+1

[keryell]

@b-sumner Your explanation is good culture information that would worth to be added as a comment to the header.

[ian-bruce]

It's still an odd way of writing it.

$ bc
bc 1.06.95
Copyright 1991-1994, 1997, 1998, 2000, 2004, 2006 Free Software Foundation, Inc.


obase=2

3141592653589793115998   # CL_M_PI
10101010010011100101101100110001100100011001010001101110100011110101\
1110

3141592650000000000000
10101010010011100101101100101110010011011100010000010110101000000000\
0000

3141592659999999999999
10101010010011100101101100110111011001100001001010001001001111111111\
1111

#23456789012345678901234^---- 25th digit



3141592741012570000000   # CL_M_PI_F
10101010010011100101101110000001000101000100011111111110010010101000\
0000

3141592740000000000000
10101010010011100101101110000000001010001000011000011110010000000000\
0000

3141592749999999999999
10101010010011100101101110001001010000001101010010010000110111111111\
1111

#23456789012345678901234^---- 25th digit



3141592653589793115997963468544185161590576171875   # "closest double precision representation of PI"
10001001100100100111100011000000100011001011110001111001011000011111\
00101011011100100110000101001110101110011000001110110010010111111011\
00010101101100101101100011

3141592650000000000000000000000000000000000000000
10001001100100100111100010111101111010011001001000111001000001011100\
01110110001111001111011111110110111111000000100000000100000000000000\
00000000000000000000000000

3141592659999999999999999999999999999999999999999
10001001100100100111100011000101010000100101110010110101011101101001\
11101000110011111111011010100101111000111101110110000111111111111111\
11111111111111111111111111

#23456789012345678901234^---- 25th digit



3141592741012573242187500000000000000000000000000   # "closest single precision representation of PI"
10001001100100100111100100000000110001110010001000111001110011110001\
01001101111111011110011001111111000000000110101100111010101110011011\
00000000000000000000000000

3141592740000000000000000000000000000000000000000
10001001100100100111100100000000000010001011000010011000111111010101\
01111101011001111110110000011101001000101000100110101000000000000000\
00000000000000000000000000

3141592749999999999999999999999999999999999999999
10001001100100100111100100000111011000010111101100010101011011100010\
11101111111110101110101011001100000010100101111100101011111111111111\
11111111111111111111111111

#23456789012345678901234^---- 25th digit

Observe that in no case do decimal digits beyond the 9th affect binary digits beyond the 24th, the size of the single-precision mantissa.

The same could be said of decimal digits beyond the 17th and binary digits beyond the 53rd, in the case of double-precision values.

Specifying more than that has no effect, and just confuses the issue.

[b-sumner]

Well, If you want to write single precision 0x1.000002p+0 accurately in decimal, you have to write it like this: 1.00000011920928955078125. There may be other decimal strings that result in the same value, but that one results in exact same double precision and extended precision value.

[ian-bruce]

The two attached small programs demonstrate that decimal digits beyond the ninth and seventeenth, respectively, have no effect whatsoever on the binary representation of single- and double-precision floating-point numbers.

The correct decimal expansion, rounded to nine or seventeen digits, produces a result just as accurate as any other possible specification. You could add digits randomly, and it would have exactly the same result, nil. This being so, specifying more digits serves no useful purpose; it is merely pointless and misleading.

$ gcc -o literal-sp literal-sp.c
$ ./literal-sp
40490FDB <--- 3.1415927410125732421875
40490FDB <--- 3.1415927400000000000000
40490FDB <--- 3.1415927499999999999999
40490FDB <--- 3.14159265

$ gcc -o literal-dp literal-dp.c
$ ./literal-dp
400921FB54442D18 <--- 3.141592653589793115997963468544185161590576171875
400921FB54442D18 <--- 3.141592653589793100000000000000000000000000000000
400921FB54442D18 <--- 3.141592653589793199999999999999999999999999999999
400921FB54442D18 <--- 3.1415926535897932

literal-sp.c
literal-dp.c

[b-sumner]

The problem is that those decimal numbers you're writing are not exactly binary floating point numbers. Don't you think that the values specified in the headers should be exactly floating point numbers?

[ian-bruce]

I don't understand the significance of "exactly binary floating-point numbers". From the two programs I provided:

#define PI_4 3.14159265    /* single-precision */

or

#define PI_4 3.1415926535897932    /* double-precision */

decimal expansion of Pi

Those are the correct decimal expansions of Pi, rounded to nine and seventeen digits. Those two programs demonstrate (try them with any compiler you like) that those definitions produce EXACTLY THE SAME BINARY VALUE, bit-for-bit, as any longer sequence. As I said before, after the ninth or seventeenth place, you could add COMPLETELY RANDOM DIGITS and get exactly the same result. It's just misleading to specify more, as that implies that it has some effect on the floating-point value, whereas it actually just doesn't, as proven above.

[keryell]

I understand your point.
On the other hand an exact representation of say 0x1.000002p+0 when you cannot write this syntax is interesting because it is just... possible :-) and is it also independent from the fact that it is stored as an OpenCL float in that case.
For example you could take the value in another context, say bc -l and you would have the same value, even if you use a double or infinite precision representation.

But at least this kind of joke for computer scientists and nerds should be explained carefully in the headers... :-)

[keryell]

I suggest we keep the values as they are but that we explain in some comments the cleverness of the approach.

But on the other hand, even if parse time and storage cost for these exact values is negligible nowadays, it seems that it breaks some compilers with some limitations #19

[Dithermaster]

IMHO, adding incorrect digits (for the constant) to match some binary representation is the wrong way to do it. There is nothing clever about it; it is just as inexact as using the correct constants (they are equally inexact). Just specify enough (correct ) decimal places of the constant to get the best binary representation. When I saw pi out to many place but with incorrect digits, the pedant in me had a fit!

[b-sumner]

I think if you are going to have a floating point number in a header, then you should have a floating point number, not a string of digits, etc. that is hopefully rounded to the desired floating point number. This wouldn't be a problem if all compilers supported hex floating point literal syntax, which permits us to write floating point numbers efficiently. Unfortunately, due to the limitations of certain compilers, we are stuck using other less useful syntax.

[hughperkins]

I think there are two ways of looking at this:

• what option(s) will give the correct results, when used?
• what will developers new to the library/headers think about the definitions, when they first view them?

But I think the best approach might be, since this is not brand-new technology, ie mathematical constants:

• pick three existing standard libraries (eg: gcc, msvc, clang)
• look at how they each define these constants
• use the same approach as at least one of these, for consistency

[ghost]

Since floating point precision is not part of the C standard, it's not good practice to make assumptions on the underlying floating point format, and to write an incorrect decimal expansion that will fit.

What I have usually seen is rather to give enough correct decimal digits, and let the compiler approximate correctly when it reads the digits.

[ian-bruce]

Since floating point precision is not part of the C standard, it's not good practice to make assumptions on the underlying floating point format, and to write an incorrect decimal expansion that will fit.

quoting from the C99 standard, Annex F, IEC 60559 floating-point arithmetic, page 444:

This annex specifies C language support for the IEC 60559 floating-point standard. The IEC 60559 floating-point standard is specifically Binary floating-point arithmetic for microprocessor systems, second edition (IEC 60559:1989), previously designated IEC 559:1989 and as IEEE Standard for Binary Floating-Point Arithmetic (ANSI/IEEE 754−1985). An implementation that defines __STDC_IEC_559__ shall conform to the specifications in this annex.

The C floating types match the IEC 60559 formats as follows:

-- The float type matches the IEC 60559 single format.

-- The double type matches the IEC 60559 double format.

http://www.open-std.org/jtc1/sc22/WG14/www/docs/n1256.pdf

However, it's obviously the OpenCL specification that actually matters:

float -- A 32-bit floating-point. The float data type must conform to the IEEE 754 single precision storage format.

double -- A 64-bit floating-point. The double data type must conform to the IEEE 754 double precision storage format.

https://www.khronos.org/registry/OpenCL/specs/opencl-2.2-cplusplus.html#builtin-scalar-data-types

As I pointed out before:

If an IEEE 754 single-precision number is converted to a decimal string with at least 9 significant digits, and then converted back to single-precision representation, the final result must match the original number.

If an IEEE 754 double-precision number is converted to a decimal string with at least 17 significant digits, and then converted back to double-precision representation, the final result must match the original number.

https://en.wikipedia.org/wiki/Single-precision_floating-point_format
https://en.wikipedia.org/wiki/Double-precision_floating-point_format
https://people.eecs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF

The ratio between the number of significant binary digits and the number of significant decimal digits is the ratio of the logarithms of 2 and 10 :

24 binary digits * (log(2) / log(10)) = 7.2247 decimal digits

53 binary digits * (log(2) / log(10)) = 15.954 decimal digits

There is nothing indeterminate about this. It is a mathematical fact. Since the OpenCL specification requires IEEE 754 single- and double-precision floating point, defining constants with more than 9 and 17 decimal digits, respectively, is completely useless. It just confuses the issue.

I don't think any other considerations are at all relevant.

[ghost]

The interesting part about this annex is the note 356 (p. 507 in the C11 standard):

"Implementations that do not define _ STDC_IEC_559 _ are not required to conform to these
specifications."

And of course, not all compilers follow completely the latest standard anyway. Otherwise, one would have preferred hex floats in the first place...

I agree that the OpenCL specification matters most.

However, I maintain that it's not good practice to do that: code tend to be copy-pasted, and this approach is error-prone. It obfuscates the code for no reason.

[jrprice]

I've opened pull request #22 which I believe implements the change that has been proposed here, to give us something concrete to review. At this stage I believe @b-sumner is the only vocal opponent for this change.

I'd appreciate if anyone in favour of the change could quickly run their eye over the PR to make sure it accurately captures what has been proposed (as well as making sure I haven't screwed up any important mathematical constants).

[jrprice]

PR #22 is now merged - closing issue.