|
From: Stavros M. <mac...@us...> - 2025-12-04 18:19:36
|
I tried using SBCL native complexes (following Barton's idea) and I don't get overflow, but the answer is off by a hundred ULPs, which is the sort of thing to expect with "naive" complex float division.
~~~
(defun $cldiv (a b)
(let* ((recta (risplit a))
(ac (complex (car recta) (cdr recta))) ; will err out if not numeric
(rectb (risplit b))
(bc (complex (car rectb) (cdr rectb)))
(q (/ ac bc)))
(add (realpart q) (mul '$%i (imagpart q)))))
bfdiv(a,b):=block([fpprec:32],float(rectform(bfloat(a)/bfloat(b))));
bfcabs(a):=block([fpprex:100],float(cabs(bfloat(a))));
clq: cldiv(1.13e-100*%i+5.43e-10,5.7e-312*%i+1.2e-311);
=> 3.691993880674702e301 - 1.7536970933199895e301 %i
bfq: bfdiv(1.13e-100*%i+5.43e-10,5.7e-312*%i+1.2e-311);
=>3.6919938806746145e301 - 1.7536970933199478e301 %i
bfdiv(clq,bfq)
=> 1.0000000000000238 - 1.0515920681324065e-17 %i
qr: bfcabs(%) => 1.0000000000000238
(qr-1)/unit_in_last_place(qr) => 107
~~~
---
**[bugs:#4645] rectform/polarform complex division overflow**
**Status:** open
**Group:** None
**Labels:** rectform
**Created:** Wed Dec 03, 2025 07:23 PM UTC by Stavros Macrakis
**Last Updated:** Wed Dec 03, 2025 07:23 PM UTC
**Owner:** nobody
~~~
fpprintprec:5$
q: (1.13e-100*%i+5.43e-10)/(5.7e-312*%i+1.2e-311)
rectform(q) => OVERFLOW
float(rectform(bfloat(q))) => 3.692e+301-1.7537e+301*%i
polarform(q) => OVERFLOW
polarform(bfloat(q)) => 4.0873b301*%e^-(4.4345b-1*%i)
~~~
---
Sent from sourceforge.net because max...@li... is subscribed to https://sourceforge.net/p/maxima/bugs/
To unsubscribe from further messages, a project admin can change settings at https://sourceforge.net/p/maxima/admin/bugs/options. Or, if this is a mailing list, you can unsubscribe from the mailing list. |