Complex square root of -1 : zsqrt(-1)

[ahmed]

Hi,
With gforth, complex.fs included.

When calculating the square root of -1, I find that:
-1e 0e 0.5e 0e z** z. 0.0000000000000000612303176911189+1.i ok
-1e 0e zsqrt z. NaN ok


What's the problem with zsqrt?

NB.
-1e 1e-16 zsqrt z. 0.0000000000000000612303176911189+1.i ok
-1e 1e-100 zsqrt z. 0.0000000000000000612303176911189+1.i ok


Ahmed

[Ron AARON]

Just out of curiosity I tried it in 8th:

ok> -1+0i 0.5 c:^ .
0.000000+1.000000i

and then:
ok> -1+0i 0.5 c:^ dup c:* .
-1.000000+0.000000i

Post by ahmed
Hi,
With gforth, complex.fs included.
-1e 0e 0.5e 0e z** z. 0.0000000000000000612303176911189+1.i  ok
-1e 0e zsqrt z. NaN ok
What's the problem with zsqrt?
NB.
-1e 1e-16 zsqrt z. 0.0000000000000000612303176911189+1.i  ok
-1e 1e-100 zsqrt z. 0.0000000000000000612303176911189+1.i  ok
Ahmed

[minforth]

Post by ahmed
Hi,
With gforth, complex.fs included.
-1e 0e 0.5e 0e z** z. 0.0000000000000000612303176911189+1.i ok
-1e 0e zsqrt z. NaN ok
What's the problem with zsqrt?

Of course you know that fp operations are nearly always just
approximations.

IOW you could improve gforth's original definition to suit your
needs:
: zsqrt ( z -- sqrt[z] )
zdup z0= 0= IF
fdup f0= IF fdrop fsqrt 0e EXIT THEN
zln z2/ zexp THEN ;

[Anton Ertl]

Post by ahmed
Hi,
With gforth, complex.fs included.
-1e 0e zsqrt z. NaN ok

Thanks for the bug report. With that bug fixed, I get

-1e 0e zsqrt z. \ output: 1.i

The others seem to be normal FP rounding.

- anton

[ahmed]

Thanks,
Another one (always with complex.fs included):

0e 0e 0e 0e z** z. NaNai ok

it must be 1 ( tested with python, julia and matlab)

and
0e 0e 1e 0e z** z. NaNai ok, it must be 0

in general, with gforth: 0e 0e (a+bi) z** gives false results,
it must give 0e when (a+bi)<>0e.

So the definition of z** also must be revised.

All my respect.

Ahmed

[Krishna Myneni]

On 8/25/24 12:38, ahmed wrote:
...

0e 0e 0e 0e z** z. NaNai  ok
it must be 1 ( tested with python, julia and matlab)
and
0e 0e 1e 0e z** z. NaNai  ok, it must be 0
in general, with gforth: 0e 0e (a+bi) z** gives false results,
it must give 0e when (a+bi)<>0e.
So the definition of z** also must be revised.

Looking at this issue in more depth, there are mathematical arguments
for the result of 0^0 to be treated as 1 or to be treated as undefined [1].


--
Krishna

1. https://en.wikipedia.org/wiki/Zero_to_the_power_of_zero

[Anton Ertl]

Post by ahmed
0e 0e 0e 0e z** z. NaNai ok
it must be 1 ( tested with python, julia and matlab)
and
0e 0e 1e 0e z** z. NaNai ok, it must be 0

Thanks. Fixed in the git head:

0e 0e 0e 0e z** z. 1. ok
0e 0e 1e 0e z** z. 0 ok

- anton

[ahmed]

With this correction,
0e 0e -1e 0e z** z. 0 ok, it is false, it must be Inf (0^(-1))

With Julia and Matlab, we have:

0^(a + i*b) gives:
- 0 for a >0 and b in R
- 1 for a =0 and b = 0
- inf for a<0 and b = 0
( in this case for julia: 0^(-1) gives Inf and O^(-1+0*im) gives
NaN + im NaN)
( but Matlab gives the same result Inf)
- NaN + i NaN for a<=0 R and b<>0

Ahmed

[ahmed]

I think this version gives the same results as Matlab


: z** ( z: a+ib c+id -- e+if)
zdup z0= if zdrop zdrop 1e 0e exit then
fover f0> if zover z0= if zdrop zdrop 0e 0e exit then then
zdup f0= f0< and if zover z0= if zdrop zdrop Inf 0e exit then then
zswap zln z* zexp
;

some tests:

0e 0e 0e 0e z** z. 1. ok
0e 0e 1e 0e z** z. 0 ok
0e 0e -1e 0e z** z. inf ok
1e 1e 0e 0e z** z. 1. ok
1e 1e 0e 1e z** z. 0.428829006294368 +0.154871752464247 i ok
0e 0e 0e 1e z** z. NaN+NaNi ok
0e 0e 1e 1e z** z. 0 ok
0e 0e -1e 1e z** z. NaN+NaNi ok
-1e 0e 0.5e 0e z** z. 0.0000000000000000612303176911189 +1. i ok
-1e 0e 0e 0e z** z. 1. ok
-1e 1e 0e 0e z** z. 1. ok
-1e 1e 1e 0e z** z. -1. +1. i ok
-1e 1e 1e 1e z** z. -0.121339466446359 +0.0569501178644237 i ok

Ahmed

[mhx]

Post by ahmed
I think this version gives the same results as Matlab
: z** ( z: a+ib c+id -- e+if)
zdup z0= if zdrop zdrop 1e 0e exit then
fover f0> if zover z0= if zdrop zdrop 0e 0e exit then then
zdup f0= f0< and if zover z0= if zdrop zdrop Inf 0e exit then then
zswap zln z* zexp
;
0e 0e 0e 0e z** z. 1. ok
0e 0e 1e 0e z** z. 0 ok
0e 0e -1e 0e z** z. inf ok
1e 1e 0e 0e z** z. 1. ok
1e 1e 0e 1e z** z. 0.428829006294368 +0.154871752464247 i ok
0e 0e 0e 1e z** z. NaN+NaNi ok
0e 0e 1e 1e z** z. 0 ok
0e 0e -1e 1e z** z. NaN+NaNi ok
-1e 0e 0.5e 0e z** z. 0.0000000000000000612303176911189 +1. i ok
-1e 0e 0e 0e z** z. 1. ok
-1e 1e 0e 0e z** z. 1. ok
-1e 1e 1e 0e z** z. -1. +1. i ok
-1e 1e 1e 1e z** z. -0.121339466446359 +0.0569501178644237 i ok
Ahmed

And what does Matlab/gForth give for this?

1e-309 0e -1e 1e z** z. ( 7.188026e+0307 -9.974133e+0308 ) ok

-marcel

[ahmed]

Post by mhx
And what does Matlab/gForth give for this?
1e-309 0e -1e 1e z** z. ( 7.188026e+0307 -9.974133e+0308 ) ok
-marcel

Yes, there is a difference between matlab and gforth for this case.
Matlab gives:
» 1e-309^(-1+i)

ans =

7.1880e+307 - Infi

Julia gives:
julia> 1e-309^(-1+im)
Inf - Inf*im

and gforth gives:
1e-309 0e -1e 1e z** z. NaN+NaNi ok

I looked at complex.fs file (within gforth directory)
z** depends on zln
zln depends on >polar

Post by mhx
polar gives 0e for 1e-309 0e ( 1e-309 0e >polar z. gives 0e) and that's
a problem.

it must gives 1e-309 0e for 1e-309 0e.

Post by mhx
polar depends on |z| which gives 0e for 1e-309 0e ( 1e-309 0e |z| f.
gives 0e) and that's a problem.

|z| depends on zsqabs which gives 0e for 1e-309 0e ( 1e-309 0e zsqabs f.
gives 0e ) and that's a problem

So the problem is in |z|. it must gives 1e-309 for 1e-309 0e.

I modified the definition of |z| like hereafter:

: |z| ( z: a+ib -- m) ( f: a b--m)
fover fover ( f: a b a b)
fabs fswap fabs
fmax ( f: a b mx)
frot frot ( f: mx a b)
fover fabs fover fabs fmax ( f: mx a b mx)
ftuck ( f: mx a mx b mx)
f/ ( f: mx a mx b/mx)
fdup f* frot frot ( f: mx [b/mx]^2 a mx)
f/ fdup f* f+ fsqrt f* ;

and with this definition of |z| I got:
1e-309 0e |z| fe. 1.00000000000000E-309 ok which is good.

and then for 1e-309 0e -1e 1e z** ,
gforth gives:
1e-309 0e -1e 1e z** z. inf-ini \ inf-inf*i ok

As you can see, it is the same result as Julia.

Ahmed

[ahmed]

Another definition for |z|

: |z| ( z: a+ib -- m) ( f: a b--m)
fover fabs fover fabs ( f: a b |a| |b|)
fmax ( f: a b mx)
frot frot ( f: mx a b)
fover fabs fover fabs fmax ( f: mx a b mx)
ftuck ( f: mx a mx b mx)
f/ ( f: mx a mx b/mx)
fdup f* frot frot ( f: mx [b/mx]^2 a mx)
f/ fdup f* f+ fsqrt f* ;

Ahmed

[mhx]

Post by ahmed
Another definition for |z|
: |z| ( z: a+ib -- m) ( f: a b--m)
fover fabs fover fabs ( f: a b |a| |b|)
fmax ( f: a b mx)
frot frot ( f: mx a b)
fover fabs fover fabs fmax ( f: mx a b mx)
ftuck ( f: mx a mx b mx)
f/ ( f: mx a mx b/mx)
fdup f* frot frot ( f: mx [b/mx]^2 a mx)
f/ fdup f* f+ fsqrt f* ;

How about
: xpythag ( F: a b -- c ) \ compute sqrt(a^2+b^2) without overflow
FABS FSWAP FABS FSWAP
F2DUP F> IF FOVER ( F: a b a -- ) F/ FSQR F1+ FSQRT F* EXIT ENDIF
FDUP F0= IF 0e
ELSE FTUCK ( F: b a b -- ) F/ FSQR F1+ FSQRT F*
ENDIF ;

FORTH> 1e-309 0e xpythag +e. 1.0000000000000000000e-0309 ok
FORTH> 1e-319 0e xpythag +e. 9.9999999999999999992e-0320 ok

-marcel

[ahmed]

Post by mhx
How about
: xpythag ( F: a b -- c ) \ compute sqrt(a^2+b^2) without overflow
FABS FSWAP FABS FSWAP
F2DUP F> IF FOVER ( F: a b a -- ) F/ FSQR F1+ FSQRT F* EXIT ENDIF
FDUP F0= IF 0e
ELSE FTUCK ( F: b a b -- ) F/ FSQR F1+ FSQRT F*
ENDIF ;
FORTH> 1e-309 0e xpythag +e. 1.0000000000000000000e-0309 ok
FORTH> 1e-319 0e xpythag +e. 9.9999999999999999992e-0320 ok
-marcel

In my defintion of |z| I forgot to process the case where a=0 or b=0.

Yes, your definition captures my approach (idea).
It works but for 0e 0e xpythag gives 0e and the fstack is not consumed.


We must empty the fstack in this case.
the defintion becomes

: xpythag ( F: a b -- c ) \ compute sqrt(a^2+b^2) without overflow
FABS FSWAP FABS FSWAP
F2DUP F> IF FOVER ( F: a b a -- ) F/ FSQR F1+ FSQRT F* EXIT ENDIF
FDUP F0= IF fdrop fdrop 0e
ELSE FTUCK ( F: b a b -- ) F/ FSQR F1+ FSQRT F*
ENDIF ;
and for z** we can use this defintion of xpythag for |z|


Ahmed

[ahmed]

Post by ahmed
In my defintion of |z| I forgot to process the case where a=0 or b=0.

I meant:

In my defintion of |z| I forgot to process the case where (a=0 and b=0).

Ahmed

[mhx]

Post by ahmed

Post by mhx
How about

[..]

Post by ahmed
Yes, your definition captures my approach (idea).
It works but for 0e 0e xpythag gives 0e and the fstack is not consumed.
We must empty the fstack in this case.
the defintion becomes
: xpythag ( F: a b -- c ) \ compute sqrt(a^2+b^2) without overflow
FABS FSWAP FABS FSWAP
F2DUP F> IF FOVER ( F: a b a -- ) F/ FSQR F1+ FSQRT F* EXIT ENDIF
FDUP F0= IF fdrop fdrop 0e
ELSE FTUCK ( F: b a b -- ) F/ FSQR F1+ FSQRT F*
ENDIF ;
and for z** we can use this defintion of xpythag for |z|

Thanks! This fixes a dormant bug in my XFLOAT-svdcmp...

But I will use "FDUP F0= IF FDROP " ( stack holds 0e 0e at this point
).

-marcel

[ahmed]

Hi,

Here is another definition for |z|.
It is based on the formula : z = |z| * exp(i*arg(z)).
One can get: |z| = z * exp(-i*arg(z)).

Here is the code:

: zarg ( z: a+bi -- ) ( f: -- theta)
fswap fatan2
;

: |z| ( z: a+bi -- ) ( f: -- |z|)
zdup zarg 0e 0e -1e z* ( z: -i*arg[z])
zexp z* ( z: z*exp(-i*arg[z])
fdrop ( z: --) ( f: |z|)
;

But it is slower than the defintion using pythag( about 3-6 times
(tested gforth)).

Ahmed

[minforth]

Better avoid calling external libraries for transcendental
functions. Calculation of fhypot is much faster and more accurate.

It is based on the calculation of fp square roots for which
there are very fast algorithms. E.g.
https://en.wikipedia.org/wiki/Fast_inverse_square_root

[minforth]

But I wouldn't bother really, use some assembly à la

double inline __declspec (naked) __fastcall sqrt(double n)
{
_asm fld qword ptr [esp+4]
_asm fsqrt
_asm ret 8
}

[ahmed]

another defintion of |z|

: |z| ( z: a+ib -- m) ( f: a b--m)
fabs fswap fabs fswap ( f: |a| |b|)
fover fover fmax ( f: |a| |b| mx)
frot frot ( f: mx |a| |b|)
2 fpick ( f: mx |a| |b| mx)
ftuck ( f: mx |a| mx |b| mx)
f/ fdup f* ( f: mx |a| mx [|b|/mx]^2)
frot frot ( f: mx [|b|/mx]^2 |a| mx)
f/ fdup f* ( f: mx [|b|/mx]^2 [|a|/mx]^2)
f+ fsqrt f* ( f: mx*[[|b|/mx]^2+[|a|/mx]^2]^0.5)
;

Ahmed

[Krishna Myneni]

Post by ahmed
Hi,
With gforth, complex.fs included.
-1e 0e 0.5e 0e z** z. 0.0000000000000000612303176911189+1.i  ok
-1e 0e zsqrt z. NaN ok
What's the problem with zsqrt?
NB.
-1e 1e-16 zsqrt z. 0.0000000000000000612303176911189+1.i  ok
-1e 1e-100 zsqrt z. 0.0000000000000000612303176911189+1.i  ok

Last time I checked, Gforth did not supply the complex number arithmetic
library provided in the Forth Scientific Library -- maybe a copyright
issue? This library was developed by David N. Williams, and has
significant amount of automated test code -- see complex-test.4th in the
kForth distribution.

Under kForth, the FSL complex library of DNW gives the following for
your test cases (note Z** is Z^ in the FSL complex library):


-1e 0e 0.5e 0e z^ z.
6.12303e-17 + i1 ok \ correct to within fp double precision

-1e 0e zsqrt z.
0 + i1 ok \ correct

0e 0e 1e 0e z^ z.
nan + inan ok \ incorrect

0e 0e 1e 0e z^ z.
nan + inan ok \ incorrect

The FSL complex library's ZSQRT handles the cases above correctly, but
Z^ does not handle the corner cases for powers of Z=0 correctly.
Examining the test code, I find that there are no tests for powers of
Z=0 in complex-test.4th.

0.1414213562373095048802E1 FCONSTANT rt2

t{ 1+i0 0+i0 z^ -> 1+i0 z}t
t{ -1+i0 0+i0 z^ -> 1+i0 z}t
t{ 0+i1 0+i0 z^ -> 1+i0 z}t
t{ 0-i1 0+i0 z^ -> 1+i0 z}t
t{ rt2 rt2 0+i0 z^ -> 1+i0 z}t
t{ rt2 -rt2 0+i0 z^ -> 1+i0 z}t
t{ -rt2 rt2 0+i0 z^ -> 1+i0 z}t
t{ -rt2 -rt2 0+i0 z^ -> 1+i0 z}t


--
Krishna

[Krishna Myneni]

Post by Krishna Myneni
0e 0e 1e 0e z^ z.
 nan + inan  ok    \ incorrect
0e 0e 1e 0e z^ z.
 nan + inan  ok     \ incorrect

Oops, I copied and pasted th same test twice. Here is the other test:

0e 0e 0e 0e z^ z.
nan + inan ok

--
KM