Imaginary part to log(2;16^3)

[Former user]

For some reason, log(2;16^3) evaluates to 12+9.12440459647775110183e-41j. This, of course, is not true; the real (pun not intended) answer is 12. I checked log(2;16^2) and log(2;16^4) which give the expected answers of 8 and 16 respectively. I also checked 2^(12+9.12440459647775110183e-41j) which correctly evaluates to 4096+2.59053785920993634877e-37j and 2^12 evaluates to 4096.

I am running v0.12 on windows.

[Helder Correia]

@Hadrien Theveneau @Pol Welter I believe you guys are in a better position to clarify this one (why the precision error when complex numbers are enabled).

[Pol Welter]

The error already occurs for the natural logarithm (ln), on which log is based.

This is the relevant code from cmath.cpp

CNumber CMath::ln(const CNumber& x)
{
    HNumber abs = CMath::abs(x).real;
    CNumber result;
    result.real = HMath::ln(abs);
    // Principal Value logarithm
    // https://en.wikipedia.org/wiki/Complex_logarithm#Definition_of_principal_value
    auto imag = HMath::arccos(x.real / abs);
    result.imag = (x.imag.isPositive() || x.imag.isZero()) ? imag : -imag;
    return result;
}

The truncation error happens either in the absolute value (here mathematically a no-op, but will actually square and then take the square root again) and the arccos, or both.

Some options to improve things here:

1. Check if x is real, and use a simplified logic in that case.
2. Use auto imag = HMath::arctan(x.imag / x.real); instead. Mathematically identical, but stays clear of arccos' singularity close to 1. This will require also treating separately the case of x.real==0, but not where it is only close.

These options are not mutually exclusive. Personally, I'd rather have rock-solid behaviour in case of real arguments, so I suggest to implement option 1 regardless of option 2.

Unfortunately I no longer have a tool chain set-up for SpeedCrunch.

[Hadrien Theveneau]

I will try to have a look on it, just for the fun.

[Hadrien Theveneau]

Pull request made.