I am analyzing the ASM out of Brian's LLVM port and ran into
an interesting piece of code I though merits some discussion.

The code is part of a CFD package originally written in the
late 1960 or early 1970. Then in 2005 it was translated from
FORTRAN (apparently 2) to C.

A snippet of the code in question follows--this code is inside
the innermost loop::

dpdx[i] = 0.0;
dpdy[i] = rho;
/*
Evaluate the residuals.
*/
ur[i] = dudt[i] - nu * ( dudxx[i] + dudyy[i] )
+ u[i] * dudx[i] + v[i] * dudy[i] + dpdx[i] / rho - f[i];

vr[i] = dvdt[i] - nu * ( dvdxx[i] + dvdyy[i] )
+ u[i] * dvdx[i] + v[i] * dvdy[i] + dpdy[i] / rho - g[i];

So, Brian's compiler happily propagates dpdx[i] = 0.0; to the
division at dpdx[i] / rho and the compiler produces:

fdiv r1,#0,r21

Then the compiler happily propagates dpdy[i] = rho; to the
division at dpdy[i] / rho and the compiler produces:

fdiv r1,r21,r21

Now, since this code was written decades before IEEE 754
the code is not anticipating NaNs and Infinities, so would
it be reasonable to replace::

dpdx[i] / rho with 0.0
and
dpdy[i] / rho with 1.0

?????

MitchAlsup <MitchAlsup@aol.com> schrieb:
> I am analyzing the ASM out of Brian's LLVM port and ran into
> an interesting piece of code I though merits some discussion.
>
> The code is part of a CFD package originally written in the
> late 1960 or early 1970. Then in 2005 it was translated from
> FORTRAN (apparently 2) to C.
>
> A snippet of the code in question follows--this code is inside
> the innermost loop::
>
> dpdx[i] = 0.0;
> dpdy[i] = rho;
> /*
> Evaluate the residuals.
> */
> ur[i] = dudt[i] - nu * ( dudxx[i] + dudyy[i] )
> + u[i] * dudx[i] + v[i] * dudy[i] + dpdx[i] / rho - f[i];
>
> vr[i] = dvdt[i] - nu * ( dvdxx[i] + dvdyy[i] )
> + u[i] * dvdx[i] + v[i] * dvdy[i] + dpdy[i] / rho - g[i];
>
> So, Brian's compiler happily propagates dpdx[i] = 0.0; to the
> division at dpdx[i] / rho and the compiler produces:
>
> fdiv r1,#0,r21
>
> Then the compiler happily propagates dpdy[i] = rho; to the
> division at dpdy[i] / rho and the compiler produces:
>
> fdiv r1,r21,r21
>
> Now, since this code was written decades before IEEE 754
> the code is not anticipating NaNs and Infinities, so would
> it be reasonable to replace::
>
> dpdx[i] / rho with 0.0
> and
> dpdy[i] / rho with 1.0

Under control of user-set flags, absolutely; otherwise not,
because rho may be NaN. -ffinite-math-only should do the trick
(I haven't checked).

In a similar vein: For your architecture, you probably do not
want to set errno on math functions with erroneous arguments. This
is permitted by the C standard, and (for example) MacOS does so, see
https://developer.apple.com/library/archive/documentation/System/Conceptual/ManPages_iPhoneOS/man3/acos.3.html#//apple_ref/doc/man/3/acos
as an example.

The advantage for My 66000 is obvious: The compiler then uses
the transcendental instructions directly instead of calling a
library function.

On Saturday, March 11, 2023 at 11:57:50 AM UTC-7, MitchAlsup wrote:
> would
> it be reasonable to replace::
>
> dpdx[i] / rho with 0.0
> and
> dpdy[i] / rho with 1.0

It certainly _seems_ like it would be reasonable to do so,
but it also seems that it would still also be reasonable
even if the floating-point numbers used followed an older
standard instead of IEEE 754.

That seems like a warning sign: what was the reason it
wasn't done that way originally?

So the question is: is this code written the way it is in
order to be self-documenting, with the programmer either
failing to optimize the inner loop, or expecting the compiler
to do so?

Or is something else strange happening that isn't apparent
from the code snippet given?

John Savard

According to MitchAlsup <MitchAlsup@aol.com>:
>Now, since this code was written decades before IEEE 754
>the code is not anticipating NaNs and Infinities, so would
>it be reasonable to replace::
>
>dpdx[i] / rho with 0.0
>and
>dpdy[i] / rho with 1.0
>
>?????

The Fortran 66 standard said that a compiler can make any
transformation that is mathematically equivalent, not just numerically
equivalent, so yes.

Not that a lot of compilers did that. There was a guy at IBM turning
Fortran H into Fortran X which did better optimization, but he set a
rule for himself that the results had to be bit identical.

--
Regards,
John Levine, johnl@taugh.com, Primary Perpetrator of "The Internet for Dummies",
Please consider the environment before reading this e-mail. https://jl.ly

On Sunday, March 12, 2023 at 3:16:50 AM UTC-5, Thomas Koenig wrote:
> MitchAlsup <Mitch...@aol.com> schrieb:
> > I am analyzing the ASM out of Brian's LLVM port and ran into
> > an interesting piece of code I though merits some discussion.
> >
> > The code is part of a CFD package originally written in the
> > late 1960 or early 1970. Then in 2005 it was translated from
> > FORTRAN (apparently 2) to C.
> >
> > A snippet of the code in question follows--this code is inside
> > the innermost loop::
> >
> > dpdx[i] = 0.0;
> > dpdy[i] = rho;
> > /*
> > Evaluate the residuals.
> > */
> > ur[i] = dudt[i] - nu * ( dudxx[i] + dudyy[i] )
> > + u[i] * dudx[i] + v[i] * dudy[i] + dpdx[i] / rho - f[i];
> >
> > vr[i] = dvdt[i] - nu * ( dvdxx[i] + dvdyy[i] )
> > + u[i] * dvdx[i] + v[i] * dvdy[i] + dpdy[i] / rho - g[i];
> >
> > So, Brian's compiler happily propagates dpdx[i] = 0.0; to the
> > division at dpdx[i] / rho and the compiler produces:
> >
> > fdiv r1,#0,r21
> >
> > Then the compiler happily propagates dpdy[i] = rho; to the
> > division at dpdy[i] / rho and the compiler produces:
> >
> > fdiv r1,r21,r21
> >
> > Now, since this code was written decades before IEEE 754
> > the code is not anticipating NaNs and Infinities, so would
> > it be reasonable to replace::
> >
> > dpdx[i] / rho with 0.0
> > and
> > dpdy[i] / rho with 1.0
> Under control of user-set flags, absolutely; otherwise not,
> because rho may be NaN. -ffinite-math-only should do the trick
> (I haven't checked).
>
> In a similar vein: For your architecture, you probably do not
> want to set errno on math functions with erroneous arguments. This
> is permitted by the C standard, and (for example) MacOS does so, see
> https://developer.apple.com/library/archive/documentation/System/Conceptual/ManPages_iPhoneOS/man3/acos.3.html#//apple_ref/doc/man/3/acos
> as an example.
>
> The advantage for My 66000 is obvious: The compiler then uses
> the transcendental instructions directly instead of calling a
> library function.
<
As I read the C standard and the IEEE 754-2019 as long as you deliver
something in accordance with 754-2019 you do not have to set errno.
Since all transcendentals deliver something in accordance, errno is not
in the picture.

On Sunday, March 12, 2023 at 3:20:29 AM UTC-5, Quadibloc wrote:
> On Saturday, March 11, 2023 at 11:57:50 AM UTC-7, MitchAlsup wrote:
> > would
> > it be reasonable to replace::
> >
> > dpdx[i] / rho with 0.0
> > and
> > dpdy[i] / rho with 1.0
> It certainly _seems_ like it would be reasonable to do so,
> but it also seems that it would still also be reasonable
> even if the floating-point numbers used followed an older
> standard instead of IEEE 754.
>
> That seems like a warning sign: what was the reason it
> wasn't done that way originally?
<
Written 2-to-1.5 decades prior to IEEE 754, and sloppily
translated into C.
>
> So the question is: is this code written the way it is in
> order to be self-documenting, with the programmer either
> failing to optimize the inner loop, or expecting the compiler
> to do so?
<
After thinking about the problem for a day, I think this is a
job for special casing the HW function unit. If FDIV runs
into either (above) is delivers a result in accordance with
754-2019 in the pipeline delay of minimal cycles rather
than the pipeline delay of FDIV--screwing up code scheduling
but being faster.
>
> Or is something else strange happening that isn't apparent
> from the code snippet given?
<
I blame sloppy translation from the original FORTRAN,
>
> John Savard

On Sunday, March 12, 2023 at 11:04:18 AM UTC-5, John Levine wrote:
> According to MitchAlsup <Mitch...@aol.com>:
> >Now, since this code was written decades before IEEE 754
> >the code is not anticipating NaNs and Infinities, so would
> >it be reasonable to replace::
> >
> >dpdx[i] / rho with 0.0
> >and
> >dpdy[i] / rho with 1.0
> >
> >?????
> The Fortran 66 standard said that a compiler can make any
> transformation that is mathematically equivalent, not just numerically
> equivalent, so yes.
<
Yes, but IEEE 754 disagrees with that compiler assertion.
>
> Not that a lot of compilers did that. There was a guy at IBM turning
> Fortran H into Fortran X which did better optimization, but he set a
> rule for himself that the results had to be bit identical.
>
Compilers are better now.
>
> --
> Regards,
> John Levine, jo...@taugh.com, Primary Perpetrator of "The Internet for Dummies",
> Please consider the environment before reading this e-mail. https://jl.ly

MitchAlsup <MitchAlsup@aol.com> schrieb:
> On Sunday, March 12, 2023 at 3:16:50 AM UTC-5, Thomas Koenig wrote:
>> MitchAlsup <Mitch...@aol.com> schrieb:
>> > I am analyzing the ASM out of Brian's LLVM port and ran into
>> > an interesting piece of code I though merits some discussion.
>> >
>> > The code is part of a CFD package originally written in the
>> > late 1960 or early 1970. Then in 2005 it was translated from
>> > FORTRAN (apparently 2) to C.
>> >
>> > A snippet of the code in question follows--this code is inside
>> > the innermost loop::
>> >
>> > dpdx[i] = 0.0;
>> > dpdy[i] = rho;
>> > /*
>> > Evaluate the residuals.
>> > */
>> > ur[i] = dudt[i] - nu * ( dudxx[i] + dudyy[i] )
>> > + u[i] * dudx[i] + v[i] * dudy[i] + dpdx[i] / rho - f[i];
>> >
>> > vr[i] = dvdt[i] - nu * ( dvdxx[i] + dvdyy[i] )
>> > + u[i] * dvdx[i] + v[i] * dvdy[i] + dpdy[i] / rho - g[i];
>> >
>> > So, Brian's compiler happily propagates dpdx[i] = 0.0; to the
>> > division at dpdx[i] / rho and the compiler produces:
>> >
>> > fdiv r1,#0,r21
>> >
>> > Then the compiler happily propagates dpdy[i] = rho; to the
>> > division at dpdy[i] / rho and the compiler produces:
>> >
>> > fdiv r1,r21,r21
>> >
>> > Now, since this code was written decades before IEEE 754
>> > the code is not anticipating NaNs and Infinities, so would
>> > it be reasonable to replace::
>> >
>> > dpdx[i] / rho with 0.0
>> > and
>> > dpdy[i] / rho with 1.0
>> Under control of user-set flags, absolutely; otherwise not,
>> because rho may be NaN. -ffinite-math-only should do the trick
>> (I haven't checked).
>>
>> In a similar vein: For your architecture, you probably do not
>> want to set errno on math functions with erroneous arguments. This
>> is permitted by the C standard, and (for example) MacOS does so, see
>> https://developer.apple.com/library/archive/documentation/System/Conceptual/ManPages_iPhoneOS/man3/acos.3.html#//apple_ref/doc/man/3/acos
>> as an example.
>>
>> The advantage for My 66000 is obvious: The compiler then uses
>> the transcendental instructions directly instead of calling a
>> library function.
><
> As I read the C standard and the IEEE 754-2019 as long as you deliver
> something in accordance with 754-2019 you do not have to set errno.
> Since all transcendentals deliver something in accordance, errno is not
> in the picture.

And yet, there are systems which set it on a domain error, for
example glibc. I would consider that a mistake, especially since
it is not prescribed by the C standard. I don't know what history
is behind that. BSD does not do so, neither does Apple.

Thomas Koenig <tkoenig@netcologne.de> writes:
>MitchAlsup <MitchAlsup@aol.com> schrieb:
>> On Sunday, March 12, 2023 at 3:16:50 AM UTC-5, Thomas Koenig wrote:
>>> MitchAlsup <Mitch...@aol.com> schrieb:
>>> > I am analyzing the ASM out of Brian's LLVM port and ran into
>>> > an interesting piece of code I though merits some discussion.
>>> >
>>> > The code is part of a CFD package originally written in the
>>> > late 1960 or early 1970. Then in 2005 it was translated from
>>> > FORTRAN (apparently 2) to C.
>>> >
>>> > A snippet of the code in question follows--this code is inside
>>> > the innermost loop::
>>> >
>>> > dpdx[i] = 0.0;
>>> > dpdy[i] = rho;
>>> > /*
>>> > Evaluate the residuals.
>>> > */
>>> > ur[i] = dudt[i] - nu * ( dudxx[i] + dudyy[i] )
>>> > + u[i] * dudx[i] + v[i] * dudy[i] + dpdx[i] / rho - f[i];
>>> >
>>> > vr[i] = dvdt[i] - nu * ( dvdxx[i] + dvdyy[i] )
>>> > + u[i] * dvdx[i] + v[i] * dvdy[i] + dpdy[i] / rho - g[i];
>>> >
>>> > So, Brian's compiler happily propagates dpdx[i] = 0.0; to the
>>> > division at dpdx[i] / rho and the compiler produces:
>>> >
>>> > fdiv r1,#0,r21
>>> >
>>> > Then the compiler happily propagates dpdy[i] = rho; to the
>>> > division at dpdy[i] / rho and the compiler produces:
>>> >
>>> > fdiv r1,r21,r21
>>> >
>>> > Now, since this code was written decades before IEEE 754
>>> > the code is not anticipating NaNs and Infinities, so would
>>> > it be reasonable to replace::
>>> >
>>> > dpdx[i] / rho with 0.0
>>> > and
>>> > dpdy[i] / rho with 1.0
>>> Under control of user-set flags, absolutely; otherwise not,
>>> because rho may be NaN. -ffinite-math-only should do the trick
>>> (I haven't checked).
>>>
>>> In a similar vein: For your architecture, you probably do not
>>> want to set errno on math functions with erroneous arguments. This
>>> is permitted by the C standard, and (for example) MacOS does so, see
>>> https://developer.apple.com/library/archive/documentation/System/Conceptual/ManPages_iPhoneOS/man3/acos.3.html#//apple_ref/doc/man/3/acos
>>> as an example.
>>>
>>> The advantage for My 66000 is obvious: The compiler then uses
>>> the transcendental instructions directly instead of calling a
>>> library function.
>><
>> As I read the C standard and the IEEE 754-2019 as long as you deliver
>> something in accordance with 754-2019 you do not have to set errno.
>> Since all transcendentals deliver something in accordance, errno is not
>> in the picture.
>
>And yet, there are systems which set it on a domain error, for
>example glibc. I would consider that a mistake, especially since
>it is not prescribed by the C standard. I don't know what history
>is behind that. BSD does not do so, neither does Apple.

This interface dates back to the SVID (System V Interface Definition),
in the trig(BA_LIB) section.

RETURN VALUE

...

If the absolute value of the argument of asin() or acos() is
greater than one or if both arguments of atan2() are zero, zero is
returned and errno is set to EDOM. In addition, a message indicating
DOMAIN error is printed on the standard error output.

...

- System V Interface Definition, Third Edition, Volume I. (c) 1983-1989.

On Sunday, March 12, 2023 at 2:03:22 PM UTC-5, Thomas Koenig wrote:
> MitchAlsup <Mitch...@aol.com> schrieb:
> > On Sunday, March 12, 2023 at 3:16:50 AM UTC-5, Thomas Koenig wrote:
> >> MitchAlsup <Mitch...@aol.com> schrieb:
> >> > I am analyzing the ASM out of Brian's LLVM port and ran into
> >> > an interesting piece of code I though merits some discussion.
> >> >
> >> > The code is part of a CFD package originally written in the
> >> > late 1960 or early 1970. Then in 2005 it was translated from
> >> > FORTRAN (apparently 2) to C.
> >> >
> >> > A snippet of the code in question follows--this code is inside
> >> > the innermost loop::
> >> >
> >> > dpdx[i] = 0.0;
> >> > dpdy[i] = rho;
> >> > /*
> >> > Evaluate the residuals.
> >> > */
> >> > ur[i] = dudt[i] - nu * ( dudxx[i] + dudyy[i] )
> >> > + u[i] * dudx[i] + v[i] * dudy[i] + dpdx[i] / rho - f[i];
> >> >
> >> > vr[i] = dvdt[i] - nu * ( dvdxx[i] + dvdyy[i] )
> >> > + u[i] * dvdx[i] + v[i] * dvdy[i] + dpdy[i] / rho - g[i];
> >> >
> >> > So, Brian's compiler happily propagates dpdx[i] = 0.0; to the
> >> > division at dpdx[i] / rho and the compiler produces:
> >> >
> >> > fdiv r1,#0,r21
> >> >
> >> > Then the compiler happily propagates dpdy[i] = rho; to the
> >> > division at dpdy[i] / rho and the compiler produces:
> >> >
> >> > fdiv r1,r21,r21
> >> >
> >> > Now, since this code was written decades before IEEE 754
> >> > the code is not anticipating NaNs and Infinities, so would
> >> > it be reasonable to replace::
> >> >
> >> > dpdx[i] / rho with 0.0
> >> > and
> >> > dpdy[i] / rho with 1.0
> >> Under control of user-set flags, absolutely; otherwise not,
> >> because rho may be NaN. -ffinite-math-only should do the trick
> >> (I haven't checked).
> >>
> >> In a similar vein: For your architecture, you probably do not
> >> want to set errno on math functions with erroneous arguments. This
> >> is permitted by the C standard, and (for example) MacOS does so, see
> >> https://developer.apple.com/library/archive/documentation/System/Conceptual/ManPages_iPhoneOS/man3/acos.3.html#//apple_ref/doc/man/3/acos
> >> as an example.
> >>
> >> The advantage for My 66000 is obvious: The compiler then uses
> >> the transcendental instructions directly instead of calling a
> >> library function.
> ><
> > As I read the C standard and the IEEE 754-2019 as long as you deliver
> > something in accordance with 754-2019 you do not have to set errno.
> > Since all transcendentals deliver something in accordance, errno is not
> > in the picture.
<
> And yet, there are systems which set it on a domain error, for
> example glibc. I would consider that a mistake, especially since
> it is not prescribed by the C standard. I don't know what history
> is behind that. BSD does not do so, neither does Apple.
<
Back when transcendental function took 150-500 cycles, setting
errno added 2 cycles, it was a reasonable decision. Now with
transcendental functions in the 20-cycle range and errno an
offset into thread local memory, setting errno might cost more
cycles than executing the transcendental !!

On Sunday, March 12, 2023 at 12:28:55 PM UTC-6, MitchAlsup wrote:
> On Sunday, March 12, 2023 at 3:20:29 AM UTC-5, Quadibloc wrote:

> > Or is something else strange happening that isn't apparent
> > from the code snippet given?
> <
> I blame sloppy translation from the original FORTRAN,

In that case, it's unfortunate you don't have a copy of the FORTRAN
original - although *if* the C code is known to deliver correct results,
I suppose you can manage from it.

Performing a wholly unnecessary arithmetic operation which is meaningless
isn't just sloppy translation - it means that something that would make it
meaningful was omitted, because someone wouldn't just substitute an expression
for 0 or 1 for no reason. So translation that's sloppy, but still correct, hardly even
seems like a possibility.

John Savard

On Monday, March 13, 2023 at 7:36:59 AM UTC-5, Quadibloc wrote:
> On Sunday, March 12, 2023 at 12:28:55 PM UTC-6, MitchAlsup wrote:
> > On Sunday, March 12, 2023 at 3:20:29 AM UTC-5, Quadibloc wrote:
>
> > > Or is something else strange happening that isn't apparent
> > > from the code snippet given?
> > <
> > I blame sloppy translation from the original FORTRAN,
<
> In that case, it's unfortunate you don't have a copy of the FORTRAN
> original - although *if* the C code is known to deliver correct results,
> I suppose you can manage from it.
<
The code reads as if this is a "canned" way of numerically performing
a spatial derivative, and then immediately applying that derivative to
another numerical calculation also of "canned" form.
>
> Performing a wholly unnecessary arithmetic operation which is meaningless
> isn't just sloppy translation - it means that something that would make it
> meaningful was omitted, because someone wouldn't just substitute an expression
> for 0 or 1 for no reason. So translation that's sloppy, but still correct, hardly even
> seems like a possibility.
>
> John Savard

Does anyone know (¿ Terje, Ivan ?) why 754-2019 specifies::

5.2.1 remainder(x, ∞) is x for finite x.
6.1 remainder(subnormal, ∞) signals underflow.
<
6.1 seems to contradict 5.2.1
<
Plus the the n of:: r = x – y×n where n = nearest integer to x / y; and whenever
| x / y | = ½, n is even; n is not bound by a 64-bit container. And ∞×0 either
raises overflow or raises Operand prior to x – ∞×0 raising underflow.

Why is underflow specified instead of overflow or operation ???

And why is the result not just x as in 5.2.1 ???

MitchAlsup wrote:
> Does anyone know (¿ Terje, Ivan ?) why 754-2019 specifies::
>
> 5.2.1 remainder(x, ∞) is x for finite x.
> 6.1 remainder(subnormal, ∞) signals underflow.

I _think_ this is just the normal rule: Subnormal result -> underflow.
> <
> 6.1 seems to contradict 5.2.1
> <
> Plus the the n of:: r = x – y×n where n = nearest integer to x / y; and whenever
> | x / y | = ½, n is even; n is not bound by a 64-bit container. And ∞×0 either
> raises overflow or raises Operand prior to x – ∞×0 raising underflow.
>
> Why is underflow specified instead of overflow or operation ???
>
> And why is the result not just x as in 5.2.1 ???
>

There are a few of these specific rules about exactly which borderline
operations deliver which signals/traps that I don't agree with, but this
is what the working group ended up with.

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

On 3/13/2023 3:09 PM, MitchAlsup wrote:
> Does anyone know (¿ Terje, Ivan ?) why 754-2019 specifies::
>
> 5.2.1 remainder(x, ∞) is x for finite x.
> 6.1 remainder(subnormal, ∞) signals underflow.
> <
> 6.1 seems to contradict 5.2.1
> <
> Plus the the n of:: r = x – y×n where n = nearest integer to x / y; and whenever
> | x / y | = ½, n is even; n is not bound by a 64-bit container. And ∞×0 either
> raises overflow or raises Operand prior to x – ∞×0 raising underflow.
>
> Why is underflow specified instead of overflow or operation ???
>
> And why is the result not just x as in 5.2.1 ???

Beats me. Terje - wanna ask the 754s?

On Tuesday, March 14, 2023 at 3:43:30 AM UTC-5, Terje Mathisen wrote:
> MitchAlsup wrote:
> > Does anyone know (¿ Terje, Ivan ?) why 754-2019 specifies::
> >
> > 5.2.1 remainder(x, ∞) is x for finite x.
> > 6.1 remainder(subnormal, ∞) signals underflow.
BTW:: 2008 has the same wording.
> I _think_ this is just the normal rule: Subnormal result -> underflow.
> > <
> > 6.1 seems to contradict 5.2.1
> > <
> > Plus the the n of:: r = x – y×n where n = nearest integer to x / y; and whenever
> > | x / y | = ½, n is even; n is not bound by a 64-bit container. And ∞×0 either
> > raises overflow or raises Operand prior to x – ∞×0 raising underflow.
> >
> > Why is underflow specified instead of overflow or operation ???
> >
> > And why is the result not just x as in 5.2.1 ???
> >
> There are a few of these specific rules about exactly which borderline
> operations deliver which signals/traps that I don't agree with, but this
> is what the working group ended up with.
>
> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"