numeric complex sqrt

[Qian Yun]

Commit 2f9e0af1 is totally correct, and LLM suggests
2 improvements:

1. use a function to compute sqrt(c^2+d^2), to avoid
overflow, it is typically called "hypot".

2. avoids accuracy loss in the computation of xa or xb.
This improves the test case you added.

--- a/src/algebra/gaussian.spad
+++ b/src/algebra/gaussian.spad
@@ -353,10 +353,10 @@
s2 := sqrt(s)
if d < 0 then
s2 := -s2
- complex(s2*xa/d, s2)
+ complex(d/(2::R*s2), s2)
s := xa/(2::R)
s2 := sqrt(s)
- complex(s2, s2*xb/d)
+ complex(s2, d/(2::R*s2))

if R has TranscendentalFunctionCategory then
half := recip(2::R)::R

- Qian

[Waldek Hebisch]

On Tue, Jul 07, 2026 at 07:49:10PM +0800, Qian Yun wrote:
> Commit 2f9e0af1 is totally correct, and LLM suggests
> 2 improvements:
>
> 1. use a function to compute sqrt(c^2+d^2), to avoid
> overflow, it is typically called "hypot".

I was thinking about this, but ATM I decided to skip it. There
3 main users of 'sqrt': R=RealClosure(FRAC(INT)), R=Float
and R=DoubleFlat. RealClosure is exact, so no accuracy problems.
Float does not overflow, so also no problem. In case of
DoubleFlat our code in general does not make significant effort
to avoid overflow, rather assumes that users will keep numbers
in reasonable range. This could change in the future, but
ATM I prefer to keep this.

> 2. avoids accuracy loss in the computation of xa or xb.
> This improves the test case you added.

Thanks, please commit.

> --- a/src/algebra/gaussian.spad
> +++ b/src/algebra/gaussian.spad
> @@ -353,10 +353,10 @@
> s2 := sqrt(s)
> if d < 0 then
> s2 := -s2
> - complex(s2*xa/d, s2)
> + complex(d/(2::R*s2), s2)
> s := xa/(2::R)
> s2 := sqrt(s)
> - complex(s2, s2*xb/d)
> + complex(s2, d/(2::R*s2))
>
> if R has TranscendentalFunctionCategory then
> half := recip(2::R)::R
>
> - Qian
>

> --
> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.
> To view this discussion visit https://groups.google.com/d/msgid/fricas-devel/3974ccd4-cd1e-4a75-a8de-247b8d7ebaef%40gmail.com.

--
Waldek Hebisch

[Qian Yun]

How about calling Common Lisp version "sqrt" for DoubleFloat?

https://github.com/sbcl/sbcl/blob/sbcl-2.6.6/src/code/irrat.lisp#L920

I just realized, for
(complex MOST-POSITIVE-DOUBLE-FLOAT MOST-POSITIVE-DOUBLE-FLOAT)
compute the "hypot" will overflow, so some kind of scaling is needed.

- Qian

[Waldek Hebisch]

On Wed, Jul 08, 2026 at 07:42:29PM +0800, Qian Yun wrote:
> How about calling Common Lisp version "sqrt" for DoubleFloat?
>
> https://github.com/sbcl/sbcl/blob/sbcl-2.6.6/src/code/irrat.lisp#L920
>
> I just realized, for
> (complex MOST-POSITIVE-DOUBLE-FLOAT MOST-POSITIVE-DOUBLE-FLOAT)
> compute the "hypot" will overflow, so some kind of scaling is needed.

Note that problem is bigger: already computing 1/z may overflow
for the same reason. Our code is generic, that is we use the
same code for various domains. For most domains this causes
no trouble. DoubleFloat may overflow, I do not consider this
as a big problem because usually range of exponents is big
enough and if not either manual scaling or extened precision
allow the result.

Note that the problem we had was different: for resonable argument
we silently got completely wrong result, it is important to fix such
problems. Overflows usually mean computation earlier produced
some unresonable value. And overflows should be detected.

BTW: AFAICS complex multiplication may overflow in intermediate
expressions even when arguments and result are in range.

--
Waldek Hebisch

[Qian Yun]

But unlike '/' and '*', for 'sqrt', the result should never overflow
for valid complex double float input.

BTW, clisp fails to do it (overflows). sbcl/ccl/cmucl/ecl can
handle it correctly.

- Qian