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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

Just out of curiosity I tried it in 8th:

-1+0i 0.5 c:^ .

0.000000+1.000000i

and then:

-1+0i 0.5 c:^ dup c:* .

-1.000000+0.000000i


On 25/08/2024 11:34, ahmed wrote:

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

On Sun, 25 Aug 2024 8:34:02 +0000, ahmed wrote:

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?

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 ;

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

melahi_ahmed@yahoo.fr (ahmed) writes:

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
--
M. Anton Ertl http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
New standard: https://forth-standard.org/
EuroForth 2024: https://euro.theforth.net

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

On 8/25/24 03:34, ahmed wrote:

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

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

On 8/25/24 14:08, Krishna Myneni wrote:

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

melahi_ahmed@yahoo.fr (ahmed) writes:

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

Thanks. Fixed in the git head:

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

- anton
--
M. Anton Ertl http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
New standard: https://forth-standard.org/
EuroForth 2024: https://euro.theforth.net

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

On Wed, 28 Aug 2024 14:49:01 +0000, ahmed wrote:

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

And what does Matlab/gForth give for this?

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

-marcel

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

On Wed, 28 Aug 2024 17:39:18 +0000, mhx wrote:

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

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.

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

On Wed, 28 Aug 2024 19:49:02 +0000, ahmed wrote:

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

On Wed, 28 Aug 2024 21:42:15 +0000, mhx wrote:

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

On Wed, 28 Aug 2024 22:10:46 +0000, ahmed wrote:

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

On Wed, 28 Aug 2024 22:10:46 +0000, ahmed wrote:

On Wed, 28 Aug 2024 21:42:15 +0000, mhx wrote:

How about

[..]

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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
}

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)