high number exponentiation of complex number

[RA Prism]

Based on some discussion on Swissmicros forum on C47 it was seen

that (1+1i)^n for high n may give problems.

Free42 / Plus42 looks quite good here, but only up to n~200.

(1+1i)^200 gives expected result of 1.267651e30 + i 0.000000

((1+1i)^200)^2 gives expected result of 1.606938e60 + i 0.000000

(1+1i)^400 gives 1.606938e60 + i 1.860879e28

square root of (1+1i)^400 has imaginary part of 0.007340

[Vincent Weber]

Indeed, good catch. 

I could verify that the real 42S does not have this issue.

Certainly Thomas will fix that... 

Brgds

[Thomas Okken]

I don't think there's anything to fix here. Free42 (and Plus42) performs raising complex numbers to integer powers using repeated squaring for powers up to ±256. For powers outside that range, it switches to a complex implementation of exp(x*log(y)), because repeated squaring becomes less accurate for higher exponents. The exp(x*log(y)) logic may return a nonzero imaginary part while the exact result would be zero, but that's to be expected with finite precision.

Maybe the HP-42S gets the correct exact imaginary part in this case thanks to its use of extended precision internally, which Free42 does not do, or maybe it handles numbers of the form a+ai as a special case. The result returned by Free42 has 32 correct digits in the real part, and in the imaginary part as well if you consider the absolute rather than relative error in both parts. Thus, the only way to do better seems to me to be using extended precision internally, and given that Free42 already uses the highest-precision type offered by the Intel library, that's not very practical.

[Vincent Weber]

But the imaginary part of (1+i)^400 should be zero, and we get something E28. 

The exp(y*ln(x)) formula also gives something like that on the real 42S.

Why such a huge divergence? I don't get the "...and in the imaginary part as well if you consider the absolute rather than relative error in both parts" statement. 

[Vincent Weber]

Wait, I actually see what you mean. The real part is in E60, the imaginary part in E28, so there is such an order of magnitude between them, that when converted to polar, the angle is only in E-32 radians, and the ABS is almost equal to the real part, at least with 32 digits precision... 

[RA Prism]

Hallo Thomas,

thank you for your explanations on this topic. It was not clear to me that there is a change in algorithms for n>256. Indeed some calculators seem to apply some rules for special cases as you mentioned for a+ai -- notably HP Prime in non-CAS mode gives surprisingly quite clean results. Wonder if such special rules might apply also to HP42s.

Best regards

Ralf.

[RA Prism]

Hallo Didier, thanks for the info. I read that the Prime Home (non-CAS) has 'only' 'traditional HP 12-digit-mantissa BCD', and it's interesting to read also just now on on hpmuseum.org about HP42S history.

So might be an idea to have the same behaviour for Free42, if HP42s is applying some special rules.

Am 25.02.25 um 11:33 schrieb Didier Lachieze:

The non-CAS math routines in the Prime have been ported to C from the 48 Saturn assembler code by Cyrille de Brebisson and for some of them date back to the HP-71B , so I would expect that the 42 as the 28, 48 and the Prime behaves the same way for these special cases.

Didier

PS: Seems that Google groups block answering directly by email client to the group.

[Thomas Okken]

Special-case logic for (a+ai)^n would be pretty simple, especially in binary. Assuming n >= 0:

n % 4 == 0 : a^n * (-4)^floor(n/4)

n % 4 == 1 : a^n * (-4)^floor(n/4) * (1+i)

n % 4 == 2 : a^n * (-4)^floor(n/4) * (2i)

n % 4 == 3 : a^n * (-4)^floor(n/4) * (-2+2i)

I bet that's what the 42S does. I don't see how it could get all the way to (1+i)^3320 with the imaginary part being exactly zero if it used the log/exp algorithm.

PS: Seems that Google groups block answering directly by email client to the group.

 

 Yes, I configured it that way in order to make spamming a bit more difficult. But now that I think about it, that's kind of redundant when only members are allowed to post anyway... So, I just enabled email posting. If that turns out to have undesirable side effects, I can always turn it back off later.

[Thomas Okken]

I made the change:

https://github.com/thomasokken/free42/commit/359ed5271e5bad5b58988f2f30a316bad4a4a130

The relevant part is lines 701-721.

[RA Prism]

Great.

[Thomas Okken]

Looking into this some more, I realized that the logic for (a+ai)^n should also be used for (a-ai)^n, and I added special cases for (a+0i)^n and (ai)^n as well. All this applies to both Free42 and Plus42 and will be in their respective next releases.

And further, while testing the new logic, I noticed that one of my reference test cases, (pi-pi*i)*(pi-pi*i), returned a result with a small but nonzero real part. This turned out to be caused by Apple's ARM compiler using FMA for the second part of expressions like a*b-c*d, which effectively causes the c*d term to be evaluated with more precision than the a*b term. I added a compiler flag to disable this optimization (which is independent of the -On flag), and that fixed the problem. This affects only the MacOS ARM Binary build; the iOS ARM Binary build is affected the same way, but since I don't distribute that build, I don't expect anyone was affected by this issue on iOS.

The Decimal builds are not affected since they do all their math using the Intel library, and the overloaded arithmetic operators are not affected by the FMA optimization. And the x86_64 aren't affected, either, even though those CPUs do have an FMA instruction. Apparently, the -ffp-contract option is only enabled by default when compiling for ARM, not for x86_64.

Finally, given that this optimization, or it being enabled by default, appears to be a relatively recent change in Apple's ARM compiler, I forced a rebuild of the Intel library and re-ran its test suite, but that showed no errors.

All in all, quite a few changes as a result of something which I initially thought to be a non-issue!

[Vincent Weber]

Thanks Thomas, very interesting! 

[Thomas Okken]

I just released Free42 3.2.4, which includes all the changes mentioned above.
Those changes will also be in the next Plus42 release, but I'm still working on that one.

[RA Prism]

Thanks for the update!

A short test with the current Free42 showed only expected results. Normally I use only Plus42 (mostly on Android) - so I will check the update when available.