literal integers in double precision expressions
alexei....@gmail.com
What is the todays opinion on the use of literal integers
in double precision expressions? I always start to doubt
after some time that they are safe to use, especially if I look
at someones code. For instance, are the two funcitons,
see below, equivalent under all circumstances?
Alexei
(rk == double precision)
function arcsinh(x)
use constants, only: one
implicit none
real(rk),intent(in) :: x(:)
real(rk) :: arcsinh(size(x))
! *** end of interface ***
real(rk) :: res(size(x),3)
where( abs(x) >= 1.0E-04_rk ) ! FIXME: double precision
arcsinh = LOG( x + sqrt( one + x**2 ) )
elsewhere
res(:,1)= x !
1.000000000000000E-004
res(:,2)= - ( 1.0_rk / 6.0_rk) * x**3 !
-1.666666666666667E-013
res(:,3)= + ( 3.0_rk / 40.0_rk) * x**5 !
7.500000000000001E-022
!es(:,4)= - ( 5.0_rk / 112.0_rk) * x**7 !
-4.464285714285714E-030
!es(:,5)= + (35.0_rk / 1152.0_rk) * x**9 !
3.038194444444445E-038
arcsinh = res(:,1) + res(:,2) + res(:,3)
endwhere
end function arcsinh
function arcsinh(x)
implicit none
real(rk),intent(in) :: x(:)
real(rk) :: arcsinh(size(x))
! *** end of interface ***
where( abs(x) >= 1.0D-04 )
arcsinh = LOG( x + sqrt( 1 + x**2 ) )
elsewhere
arcsinh = x & ! 1.000000000000000E-004
- 1 * x**3 / 6 & ! -1.666666666666667E-013
+ 3 * x**5 / 40 ! 7.500000000000001E-022
! - 5 * x**7 / 112 & ! -4.464285714285714E-030
! + 35 * x**9 / 1152 & ! 3.038194444444445E-038
endwhere
end function arcsinh
alexei....@gmail.com
> see below, equivalent under all circumstances?
sorry for line wrapping, these are the same:
function arcsinh(x)
use constants, only: one
implicit none
real(rk),intent(in) :: x(:)
real(rk) :: arcsinh(size(x))
real(rk) :: res(size(x),3)
where( abs(x) >= 1.0E-04_rk )
arcsinh = LOG( x + sqrt( one + x**2 ) )
elsewhere
res(:,1)= x
res(:,2)= - ( 1.0_rk / 6.0_rk) * x**3
res(:,3)= + ( 3.0_rk / 40.0_rk) * x**5
arcsinh = res(:,1) + res(:,2) + res(:,3)
endwhere
end function arcsinh
function arcsinh(x)
implicit none
real(rk),intent(in) :: x(:)
real(rk) :: arcsinh(size(x))
where( abs(x) >= 1.0D-04 )
arcsinh = LOG( x + sqrt( 1 + x**2 ) )
elsewhere
arcsinh = x &
- 1 * x**3 / 6 &
+ 3 * x**5 / 40
endwhere
end function arcsinh
glen herrmannsfeldt
> What is the todays opinion on the use of literal integers
> in double precision expressions? I always start to doubt
> after some time that they are safe to use, especially if I look
> at someones code. For instance, are the two funcitons,
> see below, equivalent under all circumstances?
> use constants, only: one
(snip)
> arcsinh = LOG( x + sqrt( one + x**2 ) )
(snip)
> arcsinh = LOG( x + sqrt( 1 + x**2 ) )
Is one an INTEGER PARAMETER, a REAL PARAMETER, or a
DOUBLE PRECICION PARAMETER? In any case, I would say that the
result should be the same. In years past there was some question
about the conversion from INTEGER to the appropriate type being
done at run time instead of compile time. The result should be
the same in any case. (There might be possible difference for
larger integers converting to REAL with different rounding methods.)
(snip)
> elsewhere
> arcsinh = x & ! 1.000000000000000E-004
> - 1 * x**3 / 6 & ! -1.666666666666667E-013
> + 3 * x**5 / 40 ! 7.500000000000001E-022
> ! - 5 * x**7 / 112 & ! -4.464285714285714E-030
> ! + 35 * x**9 / 1152 & ! 3.038194444444445E-038
A Horner's rule expansion would be better in this case.
Better rounding, and less computation.
-- glen
kar...@comcast.net
>
> > elsewhere
> > arcsinh = x & ! 1.000000000000000E-004
> > - 1 * x**3 / 6 & ! -1.666666666666667E-013
> > + 3 * x**5 / 40 ! 7.500000000000001E-022
> > ! - 5 * x**7 / 112 & ! -4.464285714285714E-030
> > ! + 35 * x**9 / 1152 & ! 3.038194444444445E-038
>
> A Horner's rule expansion would be better in this case.
> Better rounding, and less computation.
>
There is also the possibility that an optimizing compiler
will return 'arcsinh = x' for all values of x.
--
steve
Richard Maine
> What is the todays opinion on the use of literal integers
> in double precision expressions? I always start to doubt
> after some time that they are safe to use, especially if I look
> at someones code. For instance, are the two funcitons,
> see below, equivalent under all circumstances?
[code elided]
Certainly not under all circumstances, since that would include
circumstances where your definition of the symbol "one" was other than
what is presumably implied by its name.
I'll not argue (much) the fine points of mixed-mode arithmetic as a
style. Opinions vary; they have done so for a long time and still do.
Some people find mixed-mode expressions to be easier to read than the
equivalents with explicit mode conversion; others find the opposite.
The standard does specify how mixed-mode operations behave; it does not
vary among compilers. The expressions you showed are equivalent (at
least with any of several plausible definitions of "one"). There are
also no particularly interesting special cases for mixing integer and
double precision that are any different from mixing integer and any
other kind of real.
Probably the biggest caveat of mixed-mode operations is that if you use
them too blithely, you might be prone to write things such as
y = (1/3)*x
where the 1/3 is not a mixed-mode operation, but is integer division,
which almost certainly is not what would be intended and would have very
different results (the result of 1/3 is 0).
In the other direction, one big caveat (other than readability - note as
a specific example that I can't actually tell for sure what one of your
versions does because of the question of what "one" might mean) about
going overboard with avoiding mixed-mode operations is the case of
exponentiation. The expression
x**3
is not a mixed mode expression. Changing it to x**3.0 or any variant
thereof is usually a mistake. The two do not mean the same thing. They
are most notably different when x is negative (I'm assuming here that x
is a real of some kind), in which case x**3.0 is illegal but x**3 is
fine.
My personal style has no particular heartache with mixed-mode
expressions. Other styles differ.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
Richard Maine
If you are implying that a compiler could "optimize", for example
1 * x**3 / 6
as
(1/6) * x**3
then no. That is not an "optimization" allowed by the standard. If a
compiler does that, it is just a compiler bug. In standard-speak, the
difference between those two expressions is a mathematical difference
rather than a computational difference, Fortran's integer division being
a different mathematical operation. I think there might even be a note
in the standard about that kind of case, but I'm not feeling like
looking it up.
kar...@comcast.net
> <kar...@comcast.net> wrote:
> > On May 8, 8:49 am, glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
>
> > > > elsewhere
> > > > arcsinh = x & ! 1.000000000000000E-004
> > > > - 1 * x**3 / 6 & ! -1.666666666666667E-013
> > > > + 3 * x**5 / 40 ! 7.500000000000001E-022
> > > > ! - 5 * x**7 / 112 & ! -4.464285714285714E-030
> > > > ! + 35 * x**9 / 1152 & ! 3.038194444444445E-038
>
> > > A Horner's rule expansion would be better in this case.
> > > Better rounding, and less computation.
>
> > There is also the possibility that an optimizing compiler
> > will return 'arcsinh = x' for all values of x.
>
> If you are implying that a compiler could "optimize", for example
>
> 1 * x**3 / 6
>
> as
>
> (1/6) * x**3
>
> then no. That is not an "optimization" allowed by the standard. If a
> compiler does that, it is just a compiler bug. In standard-speak, the
> difference between those two expressions is a mathematical difference
> rather than a computational difference, Fortran's integer division being
> a different mathematical operation. I think there might even be a note
> in the standard about that kind of case, but I'm not feeling like
> looking it up.
Whoops, it's in note 7.18. I missed the prohibition on X * I / J
to X * (I / J) in the forbidden list.
--
steve
robin
news:0838beeb-ea20-4d7c...@i28g2000prd.googlegroups.com...
No it can't.