Elementary Mathematical Operations in Chapel

[damianmoz]

One can argue that some of these operations that today interface to the C library should be written in Chapel directly. With the massive improvements in the Chapel compiler since the days when such operations really needed interface to the C library for various reasons, it makes a lot of sense. Crafting these routines in Chapel simplifies their implementation. This issue is related to #24464 although that focused on optimization. It is not inferred that these routines being in Chapel will be any faster (although they should not have a performance penalty either).

The coercion operation of two numbers to a complex(w) number is one such example. This could be written as

inline operator :(x: (?,?), type t: complex(?w)) 
{
    var r : complex(w);

    r.re = x(0):real(w/2);  r.im = x(1):real(w/2);
    return r;
}

One could also write this non-generically.

For this case, and this case alone I believe, one still needs its C equivalent seen in runtime/intrinsics/chpltypes.h because they are referenced by cg-symbol.cpp . Somebody way smarter than I might be able to avoid that reference.

The conjugate and phase of a complex(w) number are further examples.

inline proc conj(in x : complex(?w))
{
    x.im = -x.im;
    return x;
}
inline proc phase(x : complex(?w))
{
   return atan2(x.im, x.re);
}

A further example is the more complicated (but computationally trivial) projection of a complex(w) onto the extended complex plane that is modelled by a Riemann Sphere. This particular functionality could not have been done properly in the past using Chapel because it needs the (only recently) introduced transmute method on floating point numbers.

enum CmplxProj
{
    ExtdCmplxPlane = 1
};

// return a complex projection
// the projection onto the extended complex plane is the default so
// it needs no qualifying argument and can be done as just proj(z).
// NOT 100% compatible with the C17 routine but makes more sense

inline proc proj(x : complex(?w), space = CmplxProj.ExtdCmplxPlane)
{
    param _w = w / 2;
    type Real = real(_w);
    type Uint = uint(_w);
    param infinity = inf:Real;

    // real programmers propagate NaNs so they can find their own errors
    // (if C17 compatability is required, you must change this to false)

    param propagate_NaNs = true;

    // capture the negative field of the imaginary component and
    // then (technically) copy that into +0.0 (which is a NO-OP)

    inline proc signedzero(t)
    {
        param _negve = 1:Uint << (_w - 1);

        return (t.transmute(Uint) & _negve).transmute(Real);
    }
    const re = x.re;
    const im = x.im;
    const are = abs(re);
    const aim = abs(im);

    return if are != infinity && aim != infinity then
        x
    else
        cmplx(if propagate_NaNs then are + aim else infinity, signedzero(im));
}

There may be better ways to write these routines in Chapel. Also, I have no idea why these existing C interfaces for these routines use pragmas and maybe something similar is needed (although I cannot think why).

I am very naughty. To make the above as readable as possible to new users of Chapel or those migrating to Chapel, one should really pass single complex(w) arguments as z because that is how they appear in most mathematics texts.

The preceding is meant to be neither perfect nor prescriptive. Hopefully I did not make any typos and do not need edits!

[damianmoz]

Looking back on #24464, I noticed that the square root function for a const real(w), the low level intrinsic is not called directly from AutoMath.chpl, but rather the libm routines (which are not macros). Would not calling the intrinsic directly from Chapel make the code simpler? Just curious.

[mppf]

@damianmoz - not sure what you mean? It looks to me like proc sqrt called on a real(64) variable/const will run this one, which runs the primitive, and eventually code-generates an the LLVM IR intrinsic:

https://github.com/chapel-lang/chapel/blob/3a49d119386f16ab6ec4bf6cc3505e64e139c06f/modules/standard/AutoMath.chpl#L608