Loic Domaigne said:
> [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.
Firstly, annex F does not apply to all implementations - it's an optional
extra. Secondly, the requirements on conversion of constants are, um,
complicated (they took an awful lot of discussion time in WG14 meetings).
Putting a constant like that is not sufficient. At the very least, if
IEEE754 applies, you should be using a hex constant not a decimal one.
C99 6.4.4.2:
[#3] [...]
For decimal floating constants, and also for hexadecimal
floating constants when FLT_RADIX is not a power of 2, the
result is either the nearest representable value, or the
larger or smaller representable value immediately adjacent
to the nearest representable value, chosen in an
implementation-defined manner. For hexadecimal floating
constants when FLT_RADIX is a power of 2, the result is
correctly rounded.
[...]
[#5] Floating constants are converted to internal format as
if at translation-time. The conversion of a floating
constant shall not raise an exceptional condition or a
floating-point exception at execution time.
Recommended practice
[#6] The implementation should produce a diagnostic message
if a hexadecimal constant cannot be represented exactly in
its evaluation format; the implementation should then
proceed with the translation of the program.
>> 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"
Of course. So my code initializes it on the first call and then uses that
value on all subsequent calls.
> 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.
Those constants aren't in the C Standard so I can't check them. But, yes,
that's a better approach.
--
Clive D.W. Feather | Work: <yyyyy@xxxxxxxxx> | Tel: +44 20 8495 6138
Internet Expert | Home: <yyyyy@xxxxxxxxxx> | *** NOTE CHANGE ***
Demon Internet | WWW: http://www.davros.org | Fax: +44 870 051 9937
Thus plc | | Mobile: +44 7973 377646
|