Hello,
Many thanks for your comments.
Enclosed you shall find my answers to your remarks. I'll send a
revised version of this example ASAP.
Loic.
[Jason wrote:]
> Wouldn't it be nice if the variable names in examples were meaningful
> rather than single letters owing their heritage to FORTRAN?
It is perfectly clear that "c" is a _c_oefficient meant to properly
normalize the integral over the real line of the density function,
so that it represents a probability, isn't it ?? ;-)
Actually, it's a Mathematical heritage, not a FORTRAN one.
In Mathematics, you usually use short names for "coefficient",
c being a favorite one...
Although I spend 7 years of my live doing numerical simulations on
High Performance Computers, I was (and I am still) reluctant to
learn such a thing like FORTRAN...
> Anyone who needs examples for math functions is not an expert in
> numerical methods. It might be helpful in this example to explain why
> your constant was defined so as to allow the function to be computed
> with a multiplication rather than a division. An expert might be
> insulted by such an observation, but an expert wouldn't be reading the
> text in the first place.
I apologize for this "cryptic" program. Doing correct and/or efficient
numerical computations implies most of the time re-writting the
equations in such a way that the computer "is pleased" with them...
In real applications, a function like normal_density() is usually
evaluated many, many times (e.g. integral computations, or random
variable simulations...), so it pays to pre-compute everything that
can be pre-computed (in particular the sqrt(2*M_PI) coefficient).
As you might have guessed, a multiplication is used instead of a
division, because floating-point multiplications are usually computed
faster.
[Clive wrote:]
> > {
> > static double c = 0.39894228040143267794; /* 1.0/sqrt(2*M_PI) */
>
>
> Not a good idea, since this may or may not be the right sequence of
> digits for any given implementation.
IMOH, this is *not* correct. IIRC, ISO C99 Annex F defines clearly
that a double corresponds to an IEEE 754 double. The definition
ensures that enough digits are provided for a double IEEE.
> If you can't just put:
>
> return exp(-x*x/2) / sqrt (2*M_PI);
Sorry, it was a reflex ;-)
> then do something like:
>
> static double c = 0;
>
> if (c < 0.25)
> c = 1.0 / sqrt (2 * M_PI);
I guess, you have understand the problem here. The fact is that
the coefficient is a constant (in the C99 sense). However, I can't
initialize it to 1.0 / sqrt (2 * M_PI), because the compiler will
issue something like:
"error: initializer element is not constant"
If you don't want to have the value hard-coded, you could declare the
coefficient "c" as:
const double c = (M_SQRT2 * M_2_SQRTPI) / 4; /* normalizer coefficient */
Ok, I agree that's hard to recognize the 1.0/sqrt(2*pi), so comments
are definitively needed.
--
+++ GMX - Mail, Messaging & more http://www.gmx.net +++
Jetzt ein- oder umsteigen und USB-Speicheruhr als Prämie sichern!
|