Fortran - Re: Random position on sphere generator -- numerics question

Author

Re: Random position on sphere generator -- numerics question

mvukovic@nycap.rr.com

2006-09-05, 10:00 pm

Ron Shepard wrote:
> In article <1157136993.969477.56210@e3g2000cwe.googlegroups.com>,
> mvukovic@nycap.rr.com wrote:
>
>
> Ok, I think I see what might be the problem. In the test inside the
> loop, the comparison is being computed as something equivalent to
>
> one-(x1**2+x2**2) >= 0.0
>
> Then, outside the loop, it is probably being computed as
>
> (one-x1**2)-x2**2
>
> In exact arithmetic, the two ways of computing that value would be
> the same, but in finite floating point, the two values may be
> different. I'm not sure, but I think the first way of computing the
> value is probably better because there is only one subtraction
> rather than two, and it is the subtractions of nearly equal values
> that causes the roundoff problems.
>
> What if you change the test inside the DO loop to
>
> sq_arg=one-(x1**2+x2**2)
> if (sq_arg .gt. 0.0) exit
>
> This way you will test on exactly the argument that is used in the
> sqrt() call. Furthermore, sq_arg is computed only once this way (it
> is effectively computed twice in your code). Also, sqrt(sq_arg)
> should be computed only once (it is computed twice in your code).
>

I was aware of the multiple sqrt evaluations, but I thought it best to
leave it to the compiler to optimize away instead of introducing
another variable.

Mirko

glen herrmannsfeldt

2006-09-06, 7:01 pm

In comp.lang.fortran Richard Maine <nospam@see.signature> wrote:
(snip regarding VOLATILE)

> For "forces" substitute some more wafflish word - perhaps "encourages".
> The details of VOLATILE are a bit ill-defined in the standard. Partly
> that is because VOLATILE is only useful in conjunction with things
> outside of the standard (such as shared memory). The standard's
> description of VOLATILE (which I don't think I'll drag out and quote
> right now, although it isn't very long) uses, as I recall, forms of the
> word "should". In standard-speak, that word means, well, pretty much
> what it does in English, except that ISO makes the distinction between
> "should" and "shall" explicit and strong, while we sometimes blur over
> it in English. In the standard, "should" indicates something that is
> recommended, but not required.

I have heard similar comments about C's volatile, claiming that
it is useless when defined that way. PL/I has ABNORMAL, though
all the compilers I know of don't do anything with it. With
multitasking and asynchronous I/O in addition to pointers, there
is at least the possibility of variables being modified at unusual
times within the standard.

There is a comment in Rich Seifert's Gigabit Ethernet book about
words in the standard. As I remember it "shall" is the important
word, while "must" doesn't mean much at all. Also, that "may" and
"may not" mean the same thing.

> Yes, you have the right idea. I'm just being nitpicky about the word
> "forces".

It seems worthwhile in this case.

-- glen

glen herrmannsfeldt

2006-09-06, 7:01 pm

In comp.lang.fortran Ron Shepard <ron-shepard@nospam.comcast.net> wrote:

(snip)

> There is also the VOLATILE keyword in f2003 (and it
> is an extension in some earlier compilers) which inhibits register
> optimization of the associated variable (i.e. it forces memory
> address reads or writes whenever it is referenced).

I think with VOLATILE it is time for someone to figure out
how to code Duff's Device ( http://en.wikipedia.org/wiki/Duff's_device )
in Fortran. Maybe this should be a challenge for Frank.

-- glen

Craig Powers

2006-09-06, 7:01 pm

glen herrmannsfeldt wrote:
> In comp.lang.fortran Ron Shepard <ron-shepard@nospam.comcast.net> wrote:
>
> (snip)
>
>
> I think with VOLATILE it is time for someone to figure out
> how to code Duff's Device ( http://en.wikipedia.org/wiki/Duff's_device )
> in Fortran. Maybe this should be a challenge for Frank.

What does VOLATILE have to do with Duff's Device? :

:

Having said that...

It isn't as "pretty" as in C, but I think this is an equivalent. I'm
using arrays source and dest as standins for the C pointers to and from.
Unlike in C, we have the advantage of the jump-table being more
explicit through the use of the computed-GO TO.

I don't remember whether it's legal to jump into a loop body. Although
Fortran does not have a C-style do-while loop, if it's OK to jump into a
loop body, you could toss in a DO / ENDDO, invert the IF condition, and
use that to EXIT.

Note -- no use of VOLATILE. :)

I believe this is standard Fortran 95, although computed GO TO is an
obsolescent feature so it may not be standard Fortran 03.

! C uses switch with 0 as a possible target; computed GO TO requires
! that our range begin with 1, so add 1 to the result.
i = 1
GO TO (100 170 160 150 140 130 120 110) MOD(size, 8) + 1

100 CONTINUE
! do {
dest(i) = source(i)
i = i + 1

170 CONTINUE
dest(i) = source(i)
i = i + 1

160 CONTINUE
dest(i) = source(i)
i = i + 1

150 CONTINUE
dest(i) = source(i)
i = i + 1

140 CONTINUE
dest(i) = source(i)
i = i + 1

130 CONTINUE
dest(i) = source(i)
i = i + 1

120 CONTINUE
dest(i) = source(i)
i = i + 1

110 CONTINUE
dest(i) = source(i)
i = i + 1

! } while ((count -= 8) > 0)
size = size - 8
IF (size > 0) GOTO 100

Craig Powers

2006-09-06, 7:01 pm

Craig Powers wrote:
> glen herrmannsfeldt wrote:
>
> What does VOLATILE have to do with Duff's Device? :

:

Oh, now I see that his destination was originally a volatile output.
That's really more of a historical detail than something intrinsic to
the device.

Anyway, to be strictly accurate, all you need to do with the code I
posted to incorporate that detail is change dest from an array to a
normal variable and assume an appropriate (VOLATILE) declaration (plus
whatever hackery is necessary to get it pointing in the right place).

Paul J Gans

2006-09-06, 10:00 pm

In sci.math.num-analysis glen herrmannsfeldt <gah@seniti.ugcs.caltech.edu> wrote:
>I think with VOLATILE it is time for someone to figure out
>how to code Duff's Device ( http://en.wikipedia.org/wiki/Duff's_device )
>in Fortran. Maybe this should be a challenge for Frank.

I think that Duff is right not to claim credit for the "device".
IIRC it was known at least a decade earlier.

---- Paul J. Gans

Dave Thompson

2006-09-20, 9:59 pm

On Wed, 06 Sep 2006 10:34:04 -0500, Ron Shepard
<ron-shepard@NOSPAM.comcast.net> wrote:

<snip: suppressing optimization>
> Also, fortran has always been required to respect parentheses in
> expressions. This eliminates the need for many of these kinds of
> directives. Many languages are not required to respect parentheses
> in arithmetic expressions. For example, C was not required to do so
> until either the 89 or the 99 revision of the language (I forget
> which).
>
C was never and still is not required to respect parentheses that are
not mathematically necessary. Parentheses that change 'routing' e.g.
(a+b)*c versus a+b*c do work of course, and those used to indicate
arguments to hence calling a function=routine impose a sequence point
-- for an actual routine; but essentially all of the standard library
routines are permitted to be implemented as macros (which might in
turn invoke compiler builtins) in which case they are not required to
have sequence points, and (fairly) simple math functions like sqrt,
exp, log etc. are the ones where an implementation is most likely to
do this. (Formally, even printf could be open-coded and not actually
be in the library, but this will usually cost so much for so little
benefit that it makes no sense.)

What C does require, I'm pretty sure newly in C89, is that floating
expressions can be computed in precision greater than that normal for
the expression's type determined by the Standard rules, but either
explicitly casting to or storing into an object of the nominal type
drops the 'extra' precision. (Although gcc for x86 with default
options, which is probably used by more programmers than any other
single compiler configuration, at least when last I looked doesn't do
to this because of the cost of 'narrowing' x87 down from 'extended'.)

- David.Thompson1 at worldnet.att.net

glen herrmannsfeldt

2006-09-20, 9:59 pm

Dave Thompson wrote:

(snip)

> What C does require, I'm pretty sure newly in C89, is that floating
> expressions can be computed in precision greater than that normal for
> the expression's type determined by the Standard rules, but either
> explicitly casting to or storing into an object of the nominal type
> drops the 'extra' precision. (Although gcc for x86 with default
> options, which is probably used by more programmers than any other
> single compiler configuration, at least when last I looked doesn't do
> to this because of the cost of 'narrowing' x87 down from 'extended'.)

Causing many posts to newsgroups when unexpected results come out.
Especially when it is done to some expressions and not others, such that
things that should be equal aren't.

-- glen

James Giles

2006-09-20, 9:59 pm

Dave Thompson wrote:
....
> C was never and still is not required to respect parentheses that are
> not mathematically necessary. [...]

I can't find justification for that statement in the C standard. Can
you post a section number (or page/line number pair) for the part
that states this?

In fact, I've talked to C implementors that state that C's requirements
on expressions are much more restrictive than Fortran's. Not only
does C have to observe the integrity of parenthesis, but also the
specified left- or right-associativity of operators. They do *not*
have a provision that permits "mathematically equivalent" expressions
to be substituted for the ones actually present like Fortran does.
That's what they say, and that's what I find in the standard document.

The "as-if" rule applies. That is, if the compiler can tell that the
result of evaluating in some other order will be the same as if it
were evaluated in the specified order, it's free to do so. But, since
few operations are computationally associative (or any of the other
really convenient properties of "real" mathematics) such optimizations
must be rare. Given the possibility of overflow, even integer addition
cannot be reordered.

Of course, "traditional" (pre-standard) C had much looser rules.
What you say was true back then. Boy I really didn't like that.

--
J. Giles

"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare

glen herrmannsfeldt

2006-10-16, 4:01 am

James Giles wrote:

(snip on C an parenthesized expressions)

> In fact, I've talked to C implementors that state that C's requirements
> on expressions are much more restrictive than Fortran's. Not only
> does C have to observe the integrity of parenthesis, but also the
> specified left- or right-associativity of operators. They do *not*
> have a provision that permits "mathematically equivalent" expressions
> to be substituted for the ones actually present like Fortran does.
> That's what they say, and that's what I find in the standard document.

C has some strict rules on evaluating integer expressions,
especially the short circuit logical operators. I would say for
integer operations, overall C is more strict.

C has tended to not to emphasize floating point, and overall
my feeling is that C is less strict on floating point rules.

Without a weighting for each rule it would be hard to make
a more quantitative comparison.

-- glen

Ron Shepard

2006-10-16, 6:59 pm

In article < RcudnQaRi7QEgK7YnZ2dnUVZ_rudnZ2d@comcast
.com>,
glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:

> C has some strict rules on evaluating integer expressions,
> especially the short circuit logical operators. I would say for
> integer operations, overall C is more strict.

The glaring exception, of course, is integer division when one or
both arguments are negative. The fortran language has always
defined what the result should be; the C language has always left it
up to the processor (although / and % must be consistent). This is
a very common situation. I have always wondered how a C programmer
could ever tolerate that kind of indefinitness in a software tool,
and there are millions of lines of C code that are written every day
that are, according to the standard, a million bugs waiting to
happen.

> C has tended to not to emphasize floating point, and overall
> my feeling is that C is less strict on floating point rules.

This is certainly true of K&R C, it was not even required to
respect parentheses in many expressions to define the order of
operations. Instead, separate parts of a long expression had to be
broken up and evaluated in separate expressions (between sequence
points).

> Without a weighting for each rule it would be hard to make
> a more quantitative comparison.

Yes, but that would not change the qualitative difference in the
approaches taken by the two languages.

$.02 -Ron Shepard <-- IMO, of course.

Elijah Cardon

2006-10-16, 6:59 pm

"Ron Shepard" <ron-shepard@NOSPAM.comcast.net> wrote in message
news:ron-shepard-071918.10370916102006@comcast.dca.giganews.com...
> In article < RcudnQaRi7QEgK7YnZ2dnUVZ_rudnZ2d@comcast
.com>,
> glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
>
>
> The glaring exception, of course, is integer division when one or
> both arguments are negative. The fortran language has always
> defined what the result should be; the C language has always left it
> up to the processor (although / and % must be consistent). This is
> a very common situation. I have always wondered how a C programmer
> could ever tolerate that kind of indefinitness in a software tool,
> and there are millions of lines of C code that are written every day
> that are, according to the standard, a million bugs waiting to
> happen.
The standard won't address the hassle it is to get numbers in and out of
these apps.

>
> This is certainly true of K&R C, it was not even required to
> respect parentheses in many expressions to define the order of
> operations. Instead, separate parts of a long expression had to be
> broken up and evaluated in separate expressions (between sequence
> points).
>
>
> Yes, but that would not change the qualitative difference in the
> approaches taken by the two languages.
My solution is to have the front end of an app in a language that takes full
advantage of the GUI. Elsewhere I'm discussing what to do with negative
input, decimal points and thousands separator. BTW, was there a consensus
on how to get a random vector on the 2 sphere? EC

Sponsored Links