# Imaginary Elementary Mathematical Functions - Very Low Priority

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.

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).

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

ALL GOOD

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

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.