* catan errors @ 2018-04-10 19:50 Rich Felker 2018-04-10 20:23 ` dgutson . 2018-04-10 23:08 ` Szabolcs Nagy 0 siblings, 2 replies; 10+ messages in thread From: Rich Felker @ 2018-04-10 19:50 UTC (permalink / raw) To: musl [-- Attachment #1: Type: text/plain, Size: 994 bytes --] The OpenBSD catan implementation we're using has a number of nonsensical "overflow" (goto ovrf) conditions that aren't errors, reported by mepholic on irc. I think the attached patch fixes them without introducing new problems, but I'm not sure if any other problems remain. Note that, of the three cases removed: 1. Is not an exceptional case at all, and made no sense to begin with. 2. Is only exceptional if x and a are both zero; atan(2x,0) is perfectly well-defined. 3. Is only possible if y==1.0 and x==0.0, which is the only real exceptional case for atan: z==I. I opted to replace the non-obvious (3) with an explicit check for z==I but this isn't necessary. As written the patch does not address raising exception flags; that should be fixed before it's committed but right now I'm just submitting this for comments. Prior to the patch, at least the following (utterly dumb -- the first is a very real input!) cases are broken: - catan(1.0) - catan(2*I) - catan(I) Rich [-- Attachment #2: catan.diff --] [-- Type: text/plain, Size: 687 bytes --] diff --git a/src/complex/catan.c b/src/complex/catan.c index 39ce6cf..9c2fccf 100644 --- a/src/complex/catan.c +++ b/src/complex/catan.c @@ -91,21 +91,17 @@ double complex catan(double complex z) x = creal(z); y = cimag(z); - if (x == 0.0 && y > 1.0) + if (x == 0.0 && y == 1.0) goto ovrf; x2 = x * x; a = 1.0 - x2 - (y * y); - if (a == 0.0) - goto ovrf; t = 0.5 * atan2(2.0 * x, a); w = _redupi(t); t = y - 1.0; a = x2 + (t * t); - if (a == 0.0) - goto ovrf; t = y + 1.0; a = (x2 + t * t)/a; @@ -113,7 +109,6 @@ double complex catan(double complex z) return w; ovrf: - // FIXME - w = MAXNUM + MAXNUM * I; + w = CMPLX(0, INFINITY); return w; } ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: catan errors 2018-04-10 19:50 catan errors Rich Felker @ 2018-04-10 20:23 ` dgutson . 2018-04-10 20:32 ` Rich Felker 2018-04-10 23:08 ` Szabolcs Nagy 1 sibling, 1 reply; 10+ messages in thread From: dgutson . @ 2018-04-10 20:23 UTC (permalink / raw) To: musl [-- Attachment #1: Type: text/plain, Size: 1714 bytes --] On Tue, Apr 10, 2018 at 4:50 PM, Rich Felker <dalias@libc.org> wrote: > The OpenBSD catan implementation we're using has a number of > nonsensical "overflow" (goto ovrf) conditions that aren't errors, > reported by mepholic on irc. I think the attached patch fixes them > without introducing new problems, but I'm not sure if any other > problems remain. > > Note that, of the three cases removed: > > 1. Is not an exceptional case at all, and made no sense to begin with. > > 2. Is only exceptional if x and a are both zero; atan(2x,0) is > perfectly well-defined. > > 3. Is only possible if y==1.0 and x==0.0, which is the only real > exceptional case for atan: z==I. > Besides the trigonometric case, are you considering de-normalized numbers, such as 4.94066e-324 as divisor? For example: double x = 1.0; double y = 5E-324; x / y is inf, and y != 0.0. Shouldn't 'a' be checked against that number or its absolute value >= minimum? > > I opted to replace the non-obvious (3) with an explicit check for z==I > but this isn't necessary. > > As written the patch does not address raising exception flags; that > should be fixed before it's committed but right now I'm just > submitting this for comments. > > Prior to the patch, at least the following (utterly dumb -- the first > is a very real input!) cases are broken: > > - catan(1.0) > - catan(2*I) > - catan(I) > > Rich > -- Who’s got the sweetest disposition? One guess, that’s who? Who’d never, ever start an argument? Who never shows a bit of temperament? Who's never wrong but always right? Who'd never dream of starting a fight? Who get stuck with all the bad luck? [-- Attachment #2: Type: text/html, Size: 2528 bytes --] ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: catan errors 2018-04-10 20:23 ` dgutson . @ 2018-04-10 20:32 ` Rich Felker 2018-04-10 20:41 ` dgutson . 0 siblings, 1 reply; 10+ messages in thread From: Rich Felker @ 2018-04-10 20:32 UTC (permalink / raw) To: musl On Tue, Apr 10, 2018 at 05:23:12PM -0300, dgutson . wrote: > On Tue, Apr 10, 2018 at 4:50 PM, Rich Felker <dalias@libc.org> wrote: > > > The OpenBSD catan implementation we're using has a number of > > nonsensical "overflow" (goto ovrf) conditions that aren't errors, > > reported by mepholic on irc. I think the attached patch fixes them > > without introducing new problems, but I'm not sure if any other > > problems remain. > > > > Note that, of the three cases removed: > > > > 1. Is not an exceptional case at all, and made no sense to begin with. > > > > 2. Is only exceptional if x and a are both zero; atan(2x,0) is > > perfectly well-defined. > > > > 3. Is only possible if y==1.0 and x==0.0, which is the only real > > exceptional case for atan: z==I. > > > > > Besides the trigonometric case, are you considering de-normalized numbers, > such as 4.94066e-324 as divisor? > For example: > double x = 1.0; > double y = 5E-324; > x / y is inf, and y != 0.0. > Shouldn't 'a' be checked against that number or its absolute value >= > minimum? Can you clarify where you think something goes wrong? Rich ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: catan errors 2018-04-10 20:32 ` Rich Felker @ 2018-04-10 20:41 ` dgutson . 2018-04-10 20:50 ` Rich Felker 0 siblings, 1 reply; 10+ messages in thread From: dgutson . @ 2018-04-10 20:41 UTC (permalink / raw) To: musl [-- Attachment #1: Type: text/plain, Size: 1902 bytes --] On Tue, Apr 10, 2018 at 5:32 PM, Rich Felker <dalias@libc.org> wrote: > On Tue, Apr 10, 2018 at 05:23:12PM -0300, dgutson . wrote: > > On Tue, Apr 10, 2018 at 4:50 PM, Rich Felker <dalias@libc.org> wrote: > > > > > The OpenBSD catan implementation we're using has a number of > > > nonsensical "overflow" (goto ovrf) conditions that aren't errors, > > > reported by mepholic on irc. I think the attached patch fixes them > > > without introducing new problems, but I'm not sure if any other > > > problems remain. > > > > > > Note that, of the three cases removed: > > > > > > 1. Is not an exceptional case at all, and made no sense to begin with. > > > > > > 2. Is only exceptional if x and a are both zero; atan(2x,0) is > > > perfectly well-defined. > > > > > > 3. Is only possible if y==1.0 and x==0.0, which is the only real > > > exceptional case for atan: z==I. > > > > > > > > > Besides the trigonometric case, are you considering de-normalized > numbers, > > such as 4.94066e-324 as divisor? > > For example: > > double x = 1.0; > > double y = 5E-324; > > x / y is inf, and y != 0.0. > > Shouldn't 'a' be checked against that number or its absolute value >= > > minimum? > > Can you clarify where you think something goes wrong? > - if (a == 0.0) - goto ovrf; t = y + 1.0; a = (x2 + t * t)/a; The check you removed does not look correct for me because what I mentioned. However, shouldn't you check, before the division, that a is not the nearest to zero (+ or -) denormalized representable double, in order to avoid ending in inf? > > Rich > -- Who’s got the sweetest disposition? One guess, that’s who? Who’d never, ever start an argument? Who never shows a bit of temperament? Who's never wrong but always right? Who'd never dream of starting a fight? Who get stuck with all the bad luck? [-- Attachment #2: Type: text/html, Size: 3354 bytes --] ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: catan errors 2018-04-10 20:41 ` dgutson . @ 2018-04-10 20:50 ` Rich Felker 2018-04-10 21:27 ` dgutson . 0 siblings, 1 reply; 10+ messages in thread From: Rich Felker @ 2018-04-10 20:50 UTC (permalink / raw) To: musl On Tue, Apr 10, 2018 at 05:41:46PM -0300, dgutson . wrote: > On Tue, Apr 10, 2018 at 5:32 PM, Rich Felker <dalias@libc.org> wrote: > > > On Tue, Apr 10, 2018 at 05:23:12PM -0300, dgutson . wrote: > > > On Tue, Apr 10, 2018 at 4:50 PM, Rich Felker <dalias@libc.org> wrote: > > > > > > > The OpenBSD catan implementation we're using has a number of > > > > nonsensical "overflow" (goto ovrf) conditions that aren't errors, > > > > reported by mepholic on irc. I think the attached patch fixes them > > > > without introducing new problems, but I'm not sure if any other > > > > problems remain. > > > > > > > > Note that, of the three cases removed: > > > > > > > > 1. Is not an exceptional case at all, and made no sense to begin with.. > > > > > > > > 2. Is only exceptional if x and a are both zero; atan(2x,0) is > > > > perfectly well-defined. > > > > > > > > 3. Is only possible if y==1.0 and x==0.0, which is the only real > > > > exceptional case for atan: z==I. > > > > > > > > > > > > > Besides the trigonometric case, are you considering de-normalized > > numbers, > > > such as 4.94066e-324 as divisor? > > > For example: > > > double x = 1.0; > > > double y = 5E-324; > > > x / y is inf, and y != 0.0. > > > Shouldn't 'a' be checked against that number or its absolute value >= > > > minimum? > > > > Can you clarify where you think something goes wrong? > > > > - if (a == 0.0) > - goto ovrf; > > t = y + 1.0; > a = (x2 + t * t)/a; > > > The check you removed does not look correct for me because what I mentioned.. > However, shouldn't you check, before the division, that a is not the > nearest to zero (+ or -) denormalized representable double, > in order to avoid ending in inf? Here a=x²+(y-1)², so unless both x==0 and y==1, the smallest a can be is DBL_EPSILON². When a is small, the numerator in the last line is also small (x²+(1+y)² < 2) so dividing by a does not overflow. Rich ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: catan errors 2018-04-10 20:50 ` Rich Felker @ 2018-04-10 21:27 ` dgutson . 2018-04-10 22:14 ` Rich Felker 0 siblings, 1 reply; 10+ messages in thread From: dgutson . @ 2018-04-10 21:27 UTC (permalink / raw) To: musl [-- Attachment #1: Type: text/plain, Size: 2877 bytes --] On Tue, Apr 10, 2018 at 5:50 PM, Rich Felker <dalias@libc.org> wrote: > On Tue, Apr 10, 2018 at 05:41:46PM -0300, dgutson . wrote: > > On Tue, Apr 10, 2018 at 5:32 PM, Rich Felker <dalias@libc.org> wrote: > > > > > On Tue, Apr 10, 2018 at 05:23:12PM -0300, dgutson . wrote: > > > > On Tue, Apr 10, 2018 at 4:50 PM, Rich Felker <dalias@libc.org> > wrote: > > > > > > > > > The OpenBSD catan implementation we're using has a number of > > > > > nonsensical "overflow" (goto ovrf) conditions that aren't errors, > > > > > reported by mepholic on irc. I think the attached patch fixes them > > > > > without introducing new problems, but I'm not sure if any other > > > > > problems remain. > > > > > > > > > > Note that, of the three cases removed: > > > > > > > > > > 1. Is not an exceptional case at all, and made no sense to begin > with.. > > > > > > > > > > 2. Is only exceptional if x and a are both zero; atan(2x,0) is > > > > > perfectly well-defined. > > > > > > > > > > 3. Is only possible if y==1.0 and x==0.0, which is the only real > > > > > exceptional case for atan: z==I. > > > > > > > > > > > > > > > > > Besides the trigonometric case, are you considering de-normalized > > > numbers, > > > > such as 4.94066e-324 as divisor? > > > > For example: > > > > double x = 1.0; > > > > double y = 5E-324; > > > > x / y is inf, and y != 0.0. > > > > Shouldn't 'a' be checked against that number or its absolute value >= > > > > minimum? > > > > > > Can you clarify where you think something goes wrong? > > > > > > > - if (a == 0.0) > > - goto ovrf; > > > > t = y + 1.0; > > a = (x2 + t * t)/a; > > > > > > The check you removed does not look correct for me because what I > mentioned.. > > However, shouldn't you check, before the division, that a is not the > > nearest to zero (+ or -) denormalized representable double, > > in order to avoid ending in inf? > > Here a=x²+(y-1)², so unless both x==0 and y==1, the smallest a can be > I was worried by the case when x is 0 and y is the next (or previous) representable value nearest to 1.0; the y == 1.0 check will fail, but the division may get big; so I did a small program and verified that the result of the division is about 3.24519e+32 when going towards negative and 8.11296e+31 when going towards positive, so everything is OK (I didn't dig in the atan2 arguments though). > is DBL_EPSILON². When a is small, the numerator in the last line is > also small (x²+(1+y)² < 2) so dividing by a does not overflow. > > Rich > -- Who’s got the sweetest disposition? One guess, that’s who? Who’d never, ever start an argument? Who never shows a bit of temperament? Who's never wrong but always right? Who'd never dream of starting a fight? Who get stuck with all the bad luck? [-- Attachment #2: Type: text/html, Size: 4214 bytes --] ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: catan errors 2018-04-10 21:27 ` dgutson . @ 2018-04-10 22:14 ` Rich Felker 0 siblings, 0 replies; 10+ messages in thread From: Rich Felker @ 2018-04-10 22:14 UTC (permalink / raw) To: musl On Tue, Apr 10, 2018 at 06:27:24PM -0300, dgutson . wrote: > On Tue, Apr 10, 2018 at 5:50 PM, Rich Felker <dalias@libc.org> wrote: > > > On Tue, Apr 10, 2018 at 05:41:46PM -0300, dgutson . wrote: > > > On Tue, Apr 10, 2018 at 5:32 PM, Rich Felker <dalias@libc.org> wrote: > > > > > > > On Tue, Apr 10, 2018 at 05:23:12PM -0300, dgutson . wrote: > > > > > On Tue, Apr 10, 2018 at 4:50 PM, Rich Felker <dalias@libc.org> > > wrote: > > > > > > > > > > > The OpenBSD catan implementation we're using has a number of > > > > > > nonsensical "overflow" (goto ovrf) conditions that aren't errors, > > > > > > reported by mepholic on irc. I think the attached patch fixes them > > > > > > without introducing new problems, but I'm not sure if any other > > > > > > problems remain. > > > > > > > > > > > > Note that, of the three cases removed: > > > > > > > > > > > > 1. Is not an exceptional case at all, and made no sense to begin > > with.. > > > > > > > > > > > > 2. Is only exceptional if x and a are both zero; atan(2x,0) is > > > > > > perfectly well-defined. > > > > > > > > > > > > 3. Is only possible if y==1.0 and x==0.0, which is the only real > > > > > > exceptional case for atan: z==I. > > > > > > > > > > > > > > > > > > > > > Besides the trigonometric case, are you considering de-normalized > > > > numbers, > > > > > such as 4.94066e-324 as divisor? > > > > > For example: > > > > > double x = 1.0; > > > > > double y = 5E-324; > > > > > x / y is inf, and y != 0.0. > > > > > Shouldn't 'a' be checked against that number or its absolute value >= > > > > > minimum? > > > > > > > > Can you clarify where you think something goes wrong? > > > > > > > > > > - if (a == 0.0) > > > - goto ovrf; > > > > > > t = y + 1.0; > > > a = (x2 + t * t)/a; > > > > > > > > > The check you removed does not look correct for me because what I > > mentioned.. > > > However, shouldn't you check, before the division, that a is not the > > > nearest to zero (+ or -) denormalized representable double, > > > in order to avoid ending in inf? > > > > Here a=x²+(y-1)², so unless both x==0 and y==1, the smallest a can be > > > > I was worried by the case when x is 0 and y is the next (or previous) > representable value nearest to 1.0; the y == 1.0 check will fail, but the > division may get big; so I did a small program and verified that the result > of the division is about 3.24519e+32 when going towards negative and > 8.11296e+31 when going towards positive, so everything is OK (I didn't dig > in the atan2 arguments though). That's basically what I said here, without working out any actual numbers: > > is DBL_EPSILON². When a is small, the numerator in the last line is > > also small (x²+(1+y)² < 2) so dividing by a does not overflow. The result of a subtraction cannot be smaller than the order of DBL_EPSILON times the larger-magnitude operand. In particular numbers like (y-1) are never small in any absolute sense. Rich ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: catan errors 2018-04-10 19:50 catan errors Rich Felker 2018-04-10 20:23 ` dgutson . @ 2018-04-10 23:08 ` Szabolcs Nagy 2018-04-10 23:23 ` Rich Felker 1 sibling, 1 reply; 10+ messages in thread From: Szabolcs Nagy @ 2018-04-10 23:08 UTC (permalink / raw) To: musl * Rich Felker <dalias@libc.org> [2018-04-10 15:50:07 -0400]: > > 3. Is only possible if y==1.0 and x==0.0, which is the only real > exceptional case for atan: z==I. > > I opted to replace the non-obvious (3) with an explicit check for z==I > but this isn't necessary. > ... > t = y - 1.0; > a = x2 + (t * t); > - if (a == 0.0) > - goto ovrf; why does a==0 imply x==0? if |x| < sqrt(2)*0x1p-538, x2 underflows to 0 in nearest rounding mode. to handle this correctly extra work would need to be done, so i think either way is fine (leaving the goto there or not are both wrong, but we dont guarantee correct complex functions yet) ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: catan errors 2018-04-10 23:08 ` Szabolcs Nagy @ 2018-04-10 23:23 ` Rich Felker 2018-04-11 0:52 ` Szabolcs Nagy 0 siblings, 1 reply; 10+ messages in thread From: Rich Felker @ 2018-04-10 23:23 UTC (permalink / raw) To: musl [-- Attachment #1: Type: text/plain, Size: 1035 bytes --] On Wed, Apr 11, 2018 at 01:08:50AM +0200, Szabolcs Nagy wrote: > * Rich Felker <dalias@libc.org> [2018-04-10 15:50:07 -0400]: > > > > 3. Is only possible if y==1.0 and x==0.0, which is the only real > > exceptional case for atan: z==I. > > > > I opted to replace the non-obvious (3) with an explicit check for z==I > > but this isn't necessary. > > > .... > > t = y - 1.0; > > a = x2 + (t * t); > > - if (a == 0.0) > > - goto ovrf; > > why does a==0 imply x==0? OK, at least by my reasoning it's only algebraically true not necessarily as doubles. > if |x| < sqrt(2)*0x1p-538, x2 underflows to 0 in nearest rounding mode. > > to handle this correctly extra work would need to be done, so i think > either way is fine (leaving the goto there or not are both wrong, but > we dont guarantee correct complex functions yet) I'm fine with just moving the check back here. But the special case (maybe not near-special cases with internal over/underflow though) works fine just replacing the *I with CMPLX. See attached. Rich [-- Attachment #2: catan2.diff --] [-- Type: text/plain, Size: 636 bytes --] diff --git a/src/complex/catan.c b/src/complex/catan.c index 39ce6cf..7dc2afe 100644 --- a/src/complex/catan.c +++ b/src/complex/catan.c @@ -91,29 +91,17 @@ double complex catan(double complex z) x = creal(z); y = cimag(z); - if (x == 0.0 && y > 1.0) - goto ovrf; - x2 = x * x; a = 1.0 - x2 - (y * y); - if (a == 0.0) - goto ovrf; t = 0.5 * atan2(2.0 * x, a); w = _redupi(t); t = y - 1.0; a = x2 + (t * t); - if (a == 0.0) - goto ovrf; t = y + 1.0; a = (x2 + t * t)/a; - w = w + (0.25 * log(a)) * I; - return w; - -ovrf: - // FIXME - w = MAXNUM + MAXNUM * I; + w = CMPLX(w, 0.25 * log(a)); return w; } ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: catan errors 2018-04-10 23:23 ` Rich Felker @ 2018-04-11 0:52 ` Szabolcs Nagy 0 siblings, 0 replies; 10+ messages in thread From: Szabolcs Nagy @ 2018-04-11 0:52 UTC (permalink / raw) To: musl * Rich Felker <dalias@libc.org> [2018-04-10 19:23:28 -0400]: > On Wed, Apr 11, 2018 at 01:08:50AM +0200, Szabolcs Nagy wrote: > > * Rich Felker <dalias@libc.org> [2018-04-10 15:50:07 -0400]: > > > t = y - 1.0; > > > a = x2 + (t * t); > > > - if (a == 0.0) > > > - goto ovrf; > > > > why does a==0 imply x==0? > > OK, at least by my reasoning it's only algebraically true not > necessarily as doubles. > > > if |x| < sqrt(2)*0x1p-538, x2 underflows to 0 in nearest rounding mode. > > > > to handle this correctly extra work would need to be done, so i think > > either way is fine (leaving the goto there or not are both wrong, but > > we dont guarantee correct complex functions yet) > > I'm fine with just moving the check back here. But the special case > (maybe not near-special cases with internal over/underflow though) > works fine just replacing the *I with CMPLX. See attached. > i think this makes sense for now. > Rich > diff --git a/src/complex/catan.c b/src/complex/catan.c > index 39ce6cf..7dc2afe 100644 > --- a/src/complex/catan.c > +++ b/src/complex/catan.c > @@ -91,29 +91,17 @@ double complex catan(double complex z) > x = creal(z); > y = cimag(z); > > - if (x == 0.0 && y > 1.0) > - goto ovrf; > - > x2 = x * x; > a = 1.0 - x2 - (y * y); > - if (a == 0.0) > - goto ovrf; > > t = 0.5 * atan2(2.0 * x, a); > w = _redupi(t); > > t = y - 1.0; > a = x2 + (t * t); > - if (a == 0.0) > - goto ovrf; > > t = y + 1.0; > a = (x2 + t * t)/a; > - w = w + (0.25 * log(a)) * I; > - return w; > - > -ovrf: > - // FIXME > - w = MAXNUM + MAXNUM * I; > + w = CMPLX(w, 0.25 * log(a)); > return w; > } ^ permalink raw reply [flat|nested] 10+ messages in thread
end of thread, other threads:[~2018-04-11 0:52 UTC | newest] Thread overview: 10+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2018-04-10 19:50 catan errors Rich Felker 2018-04-10 20:23 ` dgutson . 2018-04-10 20:32 ` Rich Felker 2018-04-10 20:41 ` dgutson . 2018-04-10 20:50 ` Rich Felker 2018-04-10 21:27 ` dgutson . 2018-04-10 22:14 ` Rich Felker 2018-04-10 23:08 ` Szabolcs Nagy 2018-04-10 23:23 ` Rich Felker 2018-04-11 0:52 ` Szabolcs Nagy
Code repositories for project(s) associated with this public inbox https://git.vuxu.org/mirror/musl/ This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for NNTP newsgroup(s).