Imaginary Elementary Mathematical Functions - Very Low Priority

Users
1

Chapel has not defined these (except for abs and conj). These are vast simplifications of the complex(w) ones.That should be crafed at some stage.

This is a very, very low priority task but I thought I would note this while it is in my brain. Should I make it an issue? Relying on an implicit conversion of imag to complex condemns the user to significantly more computations.

Note that, given an imag(w) number x, and

z = exp(x);

then z.re == cos(x) and z.im == sin(x). The components should be computed with the same low overhead of the sincos routine seen in many C libraries.

2

For that example quoted, one might have in Math.chpl something like the untested code:

  inline proc exp :(x: imag(32)) {
    extern proc _chpl_exp_imag32(t:real(32)) : complex(64);
    return _chpl_exp_imag32(t);
  }

  inline proc exp :(x: imag(64)) {
    extern proc _chpl_exp_imag64(t:real(64)) : complex(128);
    return _chpl_exp_imag64(t);
  }

and in chplmath.h, a tweaked clone of Brad's code from chpltypes.h:

static inline _complex128 _chpl_exp_imag64(_real64 t) {
// though CMPLX works for some C++ compilers, it doesn't work for all in our
// test environments, so dodge it to be safe;  Currently, we only compile
// this header using a C++ compiler when compiling our re2 stubs, which
// don't seem to use this routine anyway.
  _real64 re, im;

  sincos(t, &im, &re);

#if defined(CMPLX) && !defined(__cplusplus)
  return CMPLX(re, im);
#else
#ifndef CHPL_DONT_USE_CMPLX_PTR_ALIASING
#define cmplx_re64(c) (((double *)&(c))[0])
#define cmplx_im64(c) (((double *)&(c))[1])
  _complex128 val;
  cmplx_re64(val) = re;
  cmplx_im64(val) = im;
  return val;
#else
  // This can generate bad values in the face of inf/nan values
  return re + im*_Complex_I;
#endif
#endif
}

static inline _complex64 _chpl_exp_imag32(_real32 t) {
// though CMPLXF works for some C++ compilers, it doesn't work for all in our
// test environments, so dodge it to be safe;  Currently, we only compile
// this header using a C++ compiler when compiling our re2 stubs, which
// don't seem to use this routine anyway.
  _real32 re, im;

  sincos(t, &im, &re);

#if defined(CMPLXF) && !defined(__cplusplus)
  return CMPLXF(re, im);
#else
#ifndef CHPL_DONT_USE_CMPLX_PTR_ALIASING
#define cmplx_re32(c) (((float *)&(c))[0])
#define cmplx_im32(c) (((float *)&(c))[1])
  _complex64 val;
  cmplx_re32(val) = re;
  cmplx_im32(val) = im;
  return val;
#else
  // This can generate bad values in the face of inf/nan values
  return re + im*_Complex_I;
#endif
#endif
}

Something a little more complicated to throw at somebody like an intern would be the square root of an imag(w).

3

I thought I would just check that.

inline proc sqrt(y : imag(?w))
{
        type Real = real(w);
        type Uint = uint(w);

        inline proc negativeBitTransfer(t : Real, s : Real)
        {
                const n = s.transmute(Uint) & (1:Uint << (w - 1));

                return (n | t.transmute(Uint)).transmute(Real);
        }
        param half = 0.5:Real;
        const realy = y:Real;
        const a = sqrt(abs(half * realy));
        const b = negativeBitTransfer(a, realy);

        return (a, b):complex(2 * w);
}

config const r = 36.0;

proc verify(const r : real(?w)) // compute complex square root and check errors
{
        const ups = (1 << (if w == 32 then 23 else 52)):real(w);
        const y = r:imag(w);
        const s = sqrt(y);

        writeln(' SQRT(y)  = ', s);
        writeln('ERROR/EPS = ', ups * (s * s - y) / (if r == 0 then 1:real(w) else r));
}

proc main
{
        verify(r);
        verify(-r);
        verify(+0.0);
        verify(-0.0);
        exit(0);
}

ALL GOOD

4

ALL GOOD

It always helps if there is no bug in the reference program.

5

Hi Damian -- It sounds like you have figured out the issue and everything is OK at this point.

I wanted to respond to just one thing

Yes, making an issue is a good way to record something for future attention. It's more typical to look at / comment on old issues than to look at old Discourse posts. As a result, I think issues are more discoverable for people looking for a bug to fix and more likely to get comments from other users who are facing the same thing.

Of course, this is moot for this specific topic.