Compile Time Constants Using Compile Time Procedures

[damianmoz]

By qualifying a relatively simple proc with the param keyword, Chapel provides for the compile-time evaluation of that proc. This provides a highly readable way to express compile-time constants such as those in the Math module. With significant enhancement by @bradcray and @mppf and others over the past few years to Chapel's handling of what can and cannot be specified as param, it is proposed to exploit this feature to more clearly express such run-time constants and reduce namespace pollution.

Chapel's Math module currently contains various symbolic constants, all of which trace their heritage to equivalents in C's <math.h> header.

Rather than having a compile time constant

sqrt2

One could, by exploiting the power of Chapel's param concept, just say

sqrt(2.0)

This is facilitated by defining a param variant of the existing_sqrt_() function:

proc sqrt(param x : real(?w)) param where x == 2
{
    param A002193 = 1.4142135623730950488016887242096980785696718753;

    return A002193:real(w);
}

the first being evaluated as a 32-bit real, the second as a 64-bit real, and the last, with 2.0 having type real(64) by default, likewise.

That A002193 is the sequence within the On-Line Encyclopedia of Integer Sequences that comprises the decimal digits of the square root of 2.

One could exploit the long (30 years+) accepted name of the inverse square root function as defined by the C (and C++) standards, i.e. rsqrt(), to provide what Chapel currently calls

recipSqrt2

by providing the routine::

proc rsqrt(param x : real(?w)) param where x == 2
{
    return 0.5:real(w) * sqrt(x);
}

This can then be evaluated with no run-time overhead as

rsqrt(2.0:real(32))
rsqrt(2.0:real(64))
rsqrt(2.0);

the first being evaluated as a 32-bit real, the second as a 64-bit real, and the last, with 2.0 having type real(64) by default, likewise.

One could also, maybe contentiously, declare

proc exp(param x : real(?w)) param where x == 1
{
    param A001113 = 2.71828182845904523536028747135266249775724709369995;

    return A001113:real(w);
}

to allow a programmer zero-overhead run-time access to Napier's constant (or Euler's number or e) and provide future justification to drop e from the Math.chpl namespace.

Multiples of the constant pi such as occur in many basic complex(w) number routines and a host of other uses could be provided with

// returns the value of pi * x as a w-bit floating point number

proc pix(param x : real(?w)) param 
{
        param A000796 = 3.1415926535897932384626433832795028841971693993;

        return (A000796 * x):real(w);
}

// returns the value of x / pi as a w-bit floating point number

proc rpix(param x : real(?w)) param
{
        param A049541 = 0.3183098861837906715377675267450287240689192914809;

        return (A049541 *  x):real(w);
}

This allows Chapel to deprecate the constant pi along with several other constants, replacing their use with (for the default floating point type)

pix(1.0); // pi
pix(0.5); // pi / 2
pix(0.25); // pi / 4
pix(2.0); // 2 * pi
rpix(1.0); // 1 / pi
rpix(2.0); // 2 / pi

The expression

pix(0.25:real(32))

provides the value of pi/4 as a 32-bit floating point number.

The above does not handle twiceRecipSqrtPi. I would argue the limited use of this value outside the field of statistics means that it can be deprecated in Chapel. Even most of those who program statistical routines accurately will probably not want to use twiceRecipSqrtPi. A more detailed argument to justify that statement can be provided should it be needed.

Similarly log2E, log10E, ln2 and ln10 probably find little use outside of the family of log() and exp() routines. For the same reason, I would suggest these be deprecated. Here too, a more detailed argument to justify that statement can be provided should it be needed.

Please review and comment

The only numbers are Archimedes' constant, and Pythagorus' constant, and Napier's as these are fundamental is so multiple fields. Theodorus', i.e. sqrt(3), and others, are omitted. What do user's want/need and are they general enough to merit inclusion.

[damianmoz]

The eventual internals of the routines may not be exactly as I have quoted above.

[damianmoz]

The param proc approach provides compile-time evaluated constants to a Chapel user. The compile time values for inclusion in Math would be:

sqrt(2.0:real(w));  // Pythagoras' constant - sqrt(2)
sqrt(3.0:real(w));  // Theodorus' constant - sqrt(3)
rsqrt(2.0:real(w)); // the reciprocal of Pythagoras' constant - sqrt(1/2)
rsqrt(3.0:real(w)); // the reciprocal of Theodorus' constant - sqrt(1/3)
pix(1.0:real(w));   // Archimedes constant - pi (i.e. pi * 1)
pix(2.0:real(w));   // pi * 2
pix(0.5:real(w));   // pi / 2
pix(0.25:real(w));  // pi / 4
rpix(1.0:real(w));  // the reciprocal of Archimedes constant - 1 / pi
rpix(2.0:real(w));  // 2 / pi

These constants are useful across a wide range of fields including geometry, mathematics, physics, and engineering. This is in seriously stark contrast to other Chapel constants which are associated with logarithms and the inverse of the square root of pi which have a far more niche use.

Both pix() and rpix() actually take any real(w) number as argument, a multiplier of either pi or 1/pi respectively. Importantly, only when that argument is an integral power of 2 is the result of either of these routines not subject to any error due to rounding.

Both rpix() and rsqrt() follow the same naming scheme, that leading 'r' implying "reciprocal of" as used in the C and C++ standard libraries. That rsqrt() is a name with a nearly 40 year pedigree which first appeared in an unpublished paper by Ng and Kahan which found its way into published papers and libraries shortly thereafter and eventually standards.

Most technically focused programmers go through life without ever having to quote Napier's constant in their program, Chapel's existing constant e. It does not even appear inside the user routines of any of the recognized Math libraries.

The limited applicability of the other constants currently within Math precludes their inclusion in this approach. After thorough investigation of their use in Chapel programs, their deprecation is something for a lot of further back-and-forth discussion, review, haggling, and discussion and review and haggling, .....

[jabraham17]

Just noting the related issue #19062, which called for deprecating a few param variables, especially ones related to pi and sqrt.

That issue ended by marking a few of those constants unstable, specifically to leave us room in the future to remove them, with users instead using an equivalent param expression.

I interpret this issue as being of wider scope.

[lydia-duncan]

I find this idea appealing, thanks Damian for writing it up!

In terms of feedback, I'd probably do some cycling on the pi names. I might even be tempted to leave pi itself defined just because it is less typing than the function call with 1.0 would be, while following up on the others

Importantly, only when that argument is an integral power of 2 is the result of either of these routines not subject to any error due to rounding.

I find this intriguing. My understanding of reals was that there was always a possibility of rounding errors, though I'll definitely admit that my understanding of reals and param reals is far from complete. I think what you're saying is that the result of such calls would be completely reliable, is that accurate? That'd be awesome, if so.

[damianmoz]

Lydia, how about

proc pix(x : real(?w) = 1:real(64)) param

which means one only needs to say pix() to get the 64-bit constant, Or should it just say 1.0?

Note my edits

[damianmoz]

I should note that most of the param procs I am suggesting, the multiples of pi and _sqrt(2), appear in simple work with complex arithmtic. You might find them in Peter Harding's and my paper at ChapelCon where I first mentioned these param procs. So those constants I have targetted are extremely fundamental mathematically, not just niche ones like say that involving the inverse of the square root of pi.

[lydia-duncan]

Lydia, how about

proc pix(x : real(?w) = 1:real(64)) param

which means one only needs to say pix() to gt the 64-bit constant,

Note my edits

That's a good idea :) We'll still want to check if people are okay with that, but at least it's only a single character longer.

In terms of discussing and implementing this, we've recently finished an internal discussion about how to make breaking changes going forward. Once I have feedback from the larger team, I'm planning on sending that information out to our user base and getting their feedback. I would expect us to use a (potentially modified) version of that process to enable this change to move forward in some way (we're leaning towards test driving better community discussions soon with multidimensional array literals, so that'd be first, but I don't currently have a huge backlog of things to discuss so would anticipate being able to slot this in if that goes well. So stay tuned!)

[damianmoz]

I actually think that the concept of a param proc is extremely useful at least for me. It certainly enhances readability of the code. Maybe it was not something that presentations at conferences could extol (shout?) from the rooftops in the early days due to the limitations of what could and could not be (or be done with) a param. But it certainly is a very powerful feature these days. My 2c.

[damianmoz]

Lydia, you mentioned

I find this intriguing. My understanding of reals was
that there was always a possibility of rounding errors,

Not when multiplying an IEEE 754 number by an integral power of two, well as long as the result is a Normal Finite IEEE 754 floating point. If you were at university before 1985 and number which did not comply to the IEEE 754 floating point standard, you could think something different but I do not think anybody in your team is that old.

If you look at the hexadecimal bit pattern of Chapel's own pi and say 0.25 * pi, i.e. use transmute() to make them an uint(64), and then print them out as 'hex', you will see that the last (or rightmost) 52 bits, the significand or fractional bits, are always the same. No rounding error will have occurred due to that multiplication. Only the exponent field, which appear in the next highest 11 bits, or 3 hexadecimal digits, will change. Make sure you are working with 64-bit reals.

You also said

I think what you're saying is that the result of such calls would
be completely reliable, is that accurate? That'd be awesome, if so.

Basically, as long as you work with the reference value of pi, or what I now call pix(1.0), multiples of it like 2.0, 0.5, 0.25, or even 256.0, any integral power of 2.0, suffer no rounding error.

I hope my explanation is clear. Let me know if I need to improve it.

[damianmoz]

Has anybody surveyed the popularity of params like twiceReciprSqrtPi, ln2, ln10, log2e log10e,and even e itself within the Chapel user population?

[lydia-duncan]

Not recently. We mostly kept them around because they already existed. I have a vague memory of people using e, but that may be getting mixed in with all the times people accidentally used Math.e instead of a local e they'd defined themselves because of scoping weirdness

[damianmoz]

My experience is limited to myself and the clients for whom I work and none of us have never used e,, Our several hundred of years of experience and a few millions lines of code is large. but not a large enough population for statistical sampling purposes. But importantly, this constant never appears in any of the common or not-so-common elementary mathematical libraries. Same comment for those involving the logarithms, these appearing as a more accurate tuple, a high part and a low part. And the square root of pi that is the normalizing constant for the normal probability distribution is not even used in erf(), the library routine for that can be used for that distribution's evaluation.

Trying to go back to the first inclusion of those routines in C in BSD Unix, the opinion is that was not done rigorously. But unfortunately, the people involved seem to have long ago move on, or in a few cases, sadly left us quite recently. The only comments I have from those involved peripherally are that the people involved were far more interested in Unix system stuff and they probably just grabbed them from a text book. So a process with no rigor;

So, I would suggest that the dubious use of any of those 6 constants does not need to be perpetuated in Chapel. But only after checking with users first. That eliminates the need to document them. My 2c.

[damianmoz]

This could be the revised documentation for sqrt().

proc sqrt(param x: real(64)) param : real(64)

Returns the square root of the argument x.

It is an error if the x is less than zero.

A param argument of 2.0:real(64) or 3.0:real(64) has no run-time overhead.

proc sqrt(param x: real(32)) param : real(32)

Returns the square root of the argument x.

It is an error if the x is less than zero.

A param atgument of 2.0:real(32) or 3.0:real(32) has no run-time overhead.

NOTE Any reference to sqrt2 would disappear from the Math documentation.