Re: Rounding off double precision

Gerry Ford wrote:
"Herman D. Knoble" <SkipKnobleLESS@xxxxxxxxxxxxxxx> wrote in message news:oe2av3hjksjeptro4ijst0afa0ume93ug1@xxxxxxxxxx
On Thu, 03 Apr 2008 16:10:34 GMT, *** Hendrickson <***.hendrickson@xxxxxxx> wrote:

-|Gerry Ford wrote:
-|> "Dave Seaman" <dseaman@xxxxxxxxxxxx> wrote in message
-|> news:fstcm0$4eh$1@xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
-|>> On Mon, 31 Mar 2008 22:04:56 -0700, Gerry Ford wrote:
-|>>
-|>>
-|>>
-|>>> "Dave Seaman" <dseaman@xxxxxxxxxxxx> wrote in message
-|>>> news:fssah4$iv8$2@xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
-|>>>> On Mon, 31 Mar 2008 19:43:46 -0700 (PDT), Bamm wrote:
-|>>>>>> A parameter is not a literal constant and vice-versa.
-|>>>>> what's the difference? I thought parameter was fortran's way of
-|>>>>> declaring a constant.
-|>>>> A parameter is a symbolic constant, not a literal constant.
-|>>> I think illustrative source snippets go a long way:
-|>>> INTEGER, PARAMETER :: dp = kind(1d0)
-|>>> is probably the parameter I get the most mileage from.
-|>> And in your example, "dp" is a symbolic constant and "1d0" is a literal
-|>> constant.
-|> I didn't know that.
-|>
-|> INTEGER, PARAMETER :: qp = kind(???)
-|>
-|> Is there an analogous literal here? It sounds like I'm trying to buy
-|> pot.:-)
-|>
-|It depends on what you want to do. If the processor doesn't
-|support quad precision, there is no statement that will work. ;)
-|
-|But, if you need quad and want the program to fail at compile
-|time if quad isn't available, then something like
-| integer, parameter :: qp = kind (precision(1.0D0) + 5)
-| real (qp) :: xxx
-|will give you a higher precision on machines that support more than
-|single and double.
-|
-|If you want to be able to use quad if it is there, and use double
-|otherwise, something like
-|
-| integer, parameter :: qp_help = kind (precision(1.0D0) + 5)
-| integer, parameter :: qp = max(dp, qp_help)
-|
-|will work on any processor that uses increasing kind values for
-|increasing precisions. James van Buskirk (I think) has a much better
-|way that works in the general case; but it's too complicated
-|for me to remember :( .
-|
-|*** Hendrickson

! This code allows the compiler to fall back on double
! precision if quadruple precision is not available.
! Public domain 2004 by James Van Buskirk
!
module mykinds
implicit none
integer, parameter :: dp = selected_real_kind(15,300)
integer, parameter :: qp_preferred = selected_real_kind(30,1000)
integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1-sign(1,qp_preferred))/2*dp
end module mykinds

Skip Knoble

Thanks, Skip. I've used this source to put my own machine through its paces:

program list_kinds

real (kind=selected_real_kind(1, 3)):: var1
real*4 var2
integer :: number

integer, parameter :: sp = selected_real_kind(3,3)
integer, parameter :: dp = selected_real_kind(15,30)
!$$$$$$ I've changed James' (15, 300) to (15,30)
!$$$$$$ I think it's wrong to have the second value be
!$$$$$$ a binary order of magnitude greater than the first.
integer, parameter :: qp_preferred = selected_real_kind(30,1000)
integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1-sign(1,qp_preferred))/2*dp

real (kind=dp):: double1
integer :: buskirk1, buskirk2

!$$$$$$ do i = 1,50
!$$$$$$ print *, i, selected_real_kind(i, 100), selected_real_kind(i), &
!$$$$$$ & selected_real_kind(i,5)
!$$$$$$ enddo

var1 = .33333333333333_sp
var2 = var1-(1.0)/(3.0)
print *, "var 1 is", var1
number = selected_real_kind(1, 3)
print *, "number is ", number
! type is real*4
! number is one
print *, "var 2 is", var2
print *, "sp is of kind", sp

number = dp
print *, "dp reals are of kind number", number
double1 = .33333333333333
print *, double1 - var1

buskirk1 = 1+sign(1,qp_preferred)
print *, "(1+sign(1,qp_preferred) equals ", buskirk1

buskirk1 = ( 2*qp_preferred )+ 1-sign(1,qp_preferred)
print *, "expression equals ", buskirk1

print *, "qp is of kind", qp

buskirk1 = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1-sign(1,qp_preferred))/2*dp

print *, "expression equals ", buskirk1

buskirk2 = 42
! replace negative one with 42:
buskirk1 = (1+sign(1,buskirk2))/2*buskirk2+ &
(1-sign(1,buskirk2))/2*dp

print *, "expression equals ", buskirk1

! goocher is irrelevant in IDE
! to redirect, use dos

end program
! end source, begin output:
var 1 is 0.333333
number is 1
var 2 is 0.00000
sp is of kind 1
dp reals are of kind number 2
0.00000000000
(1+sign(1,qp_preferred) equals 0
expression equals 0
qp is of kind 2
expression equals 2
expression equals 42
! end output, begin comment

There's a lot to learn and comment on here. If I understand Richard correctly, the second number in the selected real kind invocation is to accomodate a real of size 10** your input. By making these numbers intentionally different, I think you could confuse a good compiler. Heck, it confuses me.

I used the debugger to figure out what type my machine was giving me. Instead of having to track this information down in a debugger, can you make an implementation tell you the type of any given variable? It looks like the real * power of 2 type declaration works here.

I don't think you mean the "type". That's obvious from the
declaration of the variable. If you want to know the properties of
a number, look at the precision and range intrinsics (and maybe the
kind).
print *, range(1), range(1.0), range(1.0d0)
etc, with kind and precision (there's no integer precision)

For real numbers precision and range are two different and unrelated
concepts. The SELECTED_REAL_KIND(precision,range) intrinsic
lets you specify either or both of what you need. The problem with
understanding them is that existing machines only support a couple
of versions of real numbers. So, there isn't much practical
variation in what different values of precision or range you
can specify. But, there is no contradiction in asking for a
precision of 20 decimal digits and an exponent range of 20
(meaning zero and at least 10**(-20) through 10**20).

Try the program from before with something like
do I = 1,1000,25
print *, I, selected_int_kind(range=i)
enddo
and you'll see that the kind value changes slowly and compare
that with the results of the prints from 15 lines above.

*** Hendrickson

I'm mystified that the ultimate difference here is not the higher bits in double precision:
number = dp
print *, "dp reals are of kind number", number
double1 = .33333333333333
print *, double1 - var1

What a talent is James van Buskirk! You hand him negative one, or he'll hand you the number you wanted. Since I don't know the sign syntax I'm still a little fuzzy with how he does this with arithmetic relationships. Richard Maine has written of misgivings of having these notions be integers. The upside is seen here. I always like a picture:
http://zaxfuuq.net/fortran196913.jpg .

.