Re: Rounding off double precision


"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'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 .

--

"That this social order with its pauperism, famines, prisons, gallows,
armies, and wars is necessary to society; that still greater disaster
would ensue if this organization were destroyed; all this is said only
by those who profit by this organization, while those who suffer from it
- and they are ten times as numerous - think and say quite the contrary."

~~ Leo Tolstoy


.

Quantcast