Closed achan001 closed 2 years ago
b = xim / (a 2) may underflow to 0.0 if xim is tiny, (a 2) big (if we test signbit(b), underflow is not a problem, but we are not)
If xre < 0, it is possible yim have the wrong sign.
Free42 Decimal: -1.e5000-1.e-5000i SQRT → 0+1.e2500i // wrong sign !
Above example, b = -1e-5000 / (1e2500 * 2) = -5e-7501, underflow to 0.0 This is more of an issue for Free42 Binary, with much smaller exponent range.
static int mappable_sqrt_c(phloat xre, phloat xim, phloat *yre, phloat *yim) { ... phloat r = hypot(xre, xim); phloat a = sqrt((r + fabs(xre)) / 2); phloat b = xim / (a * 2); ... if (xre >= 0) { *yre = a; *yim = b; } else if (b >= 0) { // should test xim instead *yre = b; *yim = a; } else { *yre = -b; *yim = -a; } return ERR_NONE; }
Confirmed and fixed: 364d7dd17bc1a3e1ead50d6fc6041c64cf889e87
Thanks!
b = xim / (a 2) may underflow to 0.0 if xim is tiny, (a 2) big (if we test signbit(b), underflow is not a problem, but we are not)
If xre < 0, it is possible yim have the wrong sign.
Free42 Decimal: -1.e5000-1.e-5000i SQRT → 0+1.e2500i // wrong sign !
Above example, b = -1e-5000 / (1e2500 * 2) = -5e-7501, underflow to 0.0 This is more of an issue for Free42 Binary, with much smaller exponent range.