Home > Archive > Fortran > March 2008 > Re: double precision difference between matlab and fortran 90
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
| Author |
Re: double precision difference between matlab and fortran 90
|
|
| louisrom1@gmail.com 2008-03-17, 8:18 am |
| hi
actually, in fortran
I have a variable 'alpha', when I print its value:
print*,'alpha',alpha
it gives:
0.999999999999999
and when I do :
print*,1.d0-alpha
it gives:
1.110223024625157E-015
But when I do this:
print*,1.d0-0.999999999999999d0
I actually get the same value than in Matlab which is:
9.992007221626409E-016
So I think that you are right it is matter of rounding.
The question that remains is why does fortran print 0.999999999999999
for alpha, whereas 1.d0-alpha has not the same value as
1.d0-0.999999999999999d0 ?
| |
| Kurt Kallblad 2008-03-17, 8:18 am |
|
<louisrom1@gmail.com> wrote in message
news:4fc4df4f-5053-44e0-b1ad-9ebb51dd7031@b1g2000hsg.googlegroups.com...
> hi
>
> actually, in fortran
> I have a variable 'alpha', when I print its value:
> print*,'alpha',alpha
> it gives:
> 0.999999999999999
> and when I do :
> print*,1.d0-alpha
> it gives:
> 1.110223024625157E-015
>
> But when I do this:
> print*,1.d0-0.999999999999999d0
> I actually get the same value than in Matlab which is:
> 9.992007221626409E-016
>
> So I think that you are right it is matter of rounding.
>
> The question that remains is why does fortran print 0.999999999999999
> for alpha, whereas 1.d0-alpha has not the same value as
> 1.d0-0.999999999999999d0 ?
>
With Compaq Visual Fortran 6.6 I got 9.992007221626409E-016 in all six
lines from this program.
program test
real(kind=8) :: a = 0.999999999999999d0
print *, 1 - a
print *, 1.0 - a
print *, 1d0 - a
print *, 1 - 0.999999999999999d0
print *, 1.0 - 0.999999999999999d0
print *, 1d0 - 0.999999999999999d0
end program test
The above program can't be used to test the three precision and four
rounding
modes as these only can be set to influence calculation during run time. All
calculation carried out by the CVF compiler is rounded to nearest value and
I
think this is the case for all lines above. I did another test with separate
compiled
subroutines to avoid compile time calculations but could not get any other
result
than 9.992007221626409E-016.
If you show the compleate code you use, somebody may understand what
happens.
Kurt
| |
|
| "What every computer scientist (and programmer) should know about
floating-point arithmetic"
http://docs.sun.com/source/806-3568/ncg_goldberg.html
<louisrom1@gmail.com> wrote in message
news:4fc4df4f-5053-44e0-b1ad-9ebb51dd7031@b1g2000hsg.googlegroups.com...
> hi
>
> actually, in fortran
> I have a variable 'alpha', when I print its value:
> print*,'alpha',alpha
> it gives:
> 0.999999999999999
> and when I do :
> print*,1.d0-alpha
> it gives:
> 1.110223024625157E-015
>
> But when I do this:
> print*,1.d0-0.999999999999999d0
> I actually get the same value than in Matlab which is:
> 9.992007221626409E-016
>
> So I think that you are right it is matter of rounding.
>
> The question that remains is why does fortran print 0.999999999999999
> for alpha, whereas 1.d0-alpha has not the same value as
> 1.d0-0.999999999999999d0 ?
>
>
>
| |
|
|
<louisrom1@gmail.com> wrote in message
news:6dfd3a70-0e9f-4145-81f6-6059e8bb9e6b@s12g2000prg.googlegroups.com...
> On Mar 17, 11:38 am, louisr...@gmail.com wrote:
>
>
>
> As I said above,
> I v'e come to think thanks to you all that this is a matter of
> rounding.
Not quite.
You should know that lots of floating point numbers are not representable in
binary.
They are stored as approximations.
You should also note that
0.999999999999999
is not the same as
0.999999999999999d0
> The value that 1.D0-alpha contains is rounded to
> 1.110223024625157E-015
> whereas the value of 1.D0-0.999999999999999 is rounded to
> 9.992007221626409E-016.
> The question that remains is why those two values are not rounded to
> the same value?
>
because alpha may not be exactly 0.999999999999999
it only looks like it, when printed.
You should compare the actual contents of alpha (in hex) with the hex value
of 0.999999999999999
Les
| |
| Kurt Kallblad 2008-03-19, 7:15 pm |
|
"Steven G. Kargl" <kargl@troutmask.apl.washington.edu> wrote in message
news:frrg0m$s0b$1@gnus01.u.washington.edu...
> In article
> <e63d7e90-3251-46ea-851f-ea4c6944ac1c@a1g2000hsb.googlegroups.com>,
> louisrom1@gmail.com writes:
>
> program a
> real x
> integer i
> x = asin(1.0)
> i = transfer(x,i)
> write(*,'(Z0)') i
> end program a
>
> --
> Steve
The transfer function is not needed. "Write(*,'(z0)') x" and you get the
same result as the Z-format works with the variables bits and transfer do
not changes the bit pattern.
If you want to compare the specific bits, use "Write(*,'(b0)') x".
If x is double precision you must declare i with the same size, in CVF as
integer(8) to avoid that some bits are lost.
Kurt
| |
| Richard Maine 2008-03-19, 7:15 pm |
| Kurt Kallblad <kurt.kallblad@tele2.se> wrote:
> "Steven G. Kargl" <kargl@troutmask.apl.washington.edu> wrote in message
> news:frrg0m$s0b$1@gnus01.u.washington.edu...
>
> The transfer function is not needed. "Write(*,'(z0)') x" and you get the
> same result as the Z-format works with the variables bits and transfer do
> not changes the bit pattern.
What you describe is the case for some compilers. It is *NOT* standard
Fortran. The standard specifies that Z edit descriptors are only for
integers. If you want the code to be standard and portable, you must use
something like transfer or equivalence.
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
| |
| Kurt Kallblad 2008-03-19, 7:15 pm |
|
"Richard Maine" <nospam@see.signature> wrote in message
news:1ie1r9i.88pkjodnl5wcN%nospam@see.signature...
> Kurt Kallblad <kurt.kallblad@tele2.se> wrote:
>
>
>
> What you describe is the case for some compilers. It is *NOT* standard
> Fortran. The standard specifies that Z edit descriptors are only for
> integers. If you want the code to be standard and portable, you must use
> something like transfer or equivalence.
>
> --
> Richard Maine | Good judgement comes from experience;
> email: last name at domain . net | experience comes from bad judgement.
> domain: summertriangle | -- Mark Twain
Hi Richard,
As in old time, I'm learning a lot from you.
Kurt
| |
| Richard Maine 2008-03-19, 7:15 pm |
| User-Agent: MacSOUP/2.8.2 (Mac OS X version 10.5.2 (x86))
X-Complaints-To: abuse@supernews.com
Lines: 33
Bytes: 2641
Xref: number1.nntp.dca.giganews.com comp.lang.fortran:175141
FX <coudert@alussinan.org> wrote:
>
> I thought equivalences didn't give any guarantee of the sort? Mainly,
> that if a and b are equivalenced, you can't set a and read b. Am I wrong?
You are technically right... but then the difference in this regard
between equivalence and transfer is essentially nil.
With equivalence, b is undefined, which means that referencing it is
nonstandard, which in turns means that anything can happen.
With transfer, the value of b is processor-dependent, which has pretty
much the same implications.
The reasons behind those two caveats are identical. The exact expression
of the restrictions in standard-speak is slightly different.
In theory, a compiler could detect that the equivalence trick had been
used and refuse to compile or some such thing. In practice, the number
of compilers that do this is zero. Even those compilers that try very
hard to track undefined variables do *NOT* catch this case.
I consider them to be functionally equivalent in practice. That's unlike
the case of using z edit descriptors directly on noninteger variables.
There are compilers that won't do that one.
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
| |
| glen herrmannsfeldt 2008-03-19, 7:15 pm |
| Richard Maine wrote:
> Kurt Kallblad <kurt.kallblad@tele2.se> wrote:
(snip)
[color=darkred]
> What you describe is the case for some compilers. It is *NOT* standard
> Fortran. The standard specifies that Z edit descriptors are only for
> integers. If you want the code to be standard and portable, you must use
> something like transfer or equivalence.
Most likely especially true for values with no equivalent sized
integer. While 64 bit integers are getting common, the 128 bit
integer needed to do double complex (complex(kind(1.d0))) are rare.
Equivalence to an array of integers and use the appropriate width
Z format, such as Z8.8, to force the leading zeros. It seems,
then, not completely portable as you need to know the appropriate
width for the Z, and also the number of array elements corresponding
to a given type. The the result of using Z for a real type is
already dependent on the internal representation.
(Why is there no C_SIZEOF function in C interoperability?)
-- glen
| |
| Richard Maine 2008-03-19, 7:15 pm |
| glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
> (Why is there no C_SIZEOF function in C interoperability?)
Well, there isn't even a Fortran sizeof. I proposed one (which I called
storage_size), which I thought was awfully simple, but I was too late in
the f2003 process to sucessfully sell it.
I found it  that f2003 instead introduced a "new" feature that only
supported f77 style. Seemed like a step backwards to me. In particular,
f2003 introduced named constants for the size of the numeric storage
unit and the character storage unit. So if you restricted yourself to
the types and kinds that were in f77 (those are the ones that have sizes
in terms of numeric and character storage units), you can find their
sizes. A bit of a pet peeve of mine.
At least for a while storage_size was proposed for f2008, but I haven't
kept track of all the ins and outs. I do recall that it kept getting
shot down by arguments that I thought completely bogus that kept trying
to make it out to be far more complicated than it was. I had specified
it to apply to only the already well-defined cases (the ones where the
standard already said there was a storage unit of some size, but didn't
provide a way to query the size of). Seemed like some people argued that
it would be useless unless it was extended to cover more ill-defined
cases, but then it would be too complicated. Alas. Of course, I'm sure
the real reason it got shot down is that nobody (such as myself) was at
the relevant meetings actively pushing it.
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
| |
| Tim Prince 2008-03-19, 7:15 pm |
| glen herrmannsfeldt wrote:
> Richard Maine wrote:
>
> (snip)
>
>
>
> Most likely especially true for values with no equivalent sized
> integer. While 64 bit integers are getting common, the 128 bit
> integer needed to do double complex (complex(kind(1.d0))) are rare.
>
> Equivalence to an array of integers and use the appropriate width
> Z format, such as Z8.8, to force the leading zeros. It seems,
> then, not completely portable as you need to know the appropriate
> width for the Z, and also the number of array elements corresponding
> to a given type. The the result of using Z for a real type is
> already dependent on the internal representation.
I guess this warning includes the possibility of real and integer types
having different byte storage orders, as it was on the PDP and VAX, and
some early personal computers. That may be enough reason for Fortran not
defining this usage.
| |
| glen herrmannsfeldt 2008-03-19, 7:15 pm |
| Tim Prince wrote:
(snip, I wrote)
[color=darkred]
> I guess this warning includes the possibility of real and integer types
> having different byte storage orders, as it was on the PDP and VAX, and
> some early personal computers.
If you EQUIVALENCE (or TRANSFER) between REAL and the appropriate
sized INTEGER, and the byte order are the same then the appropriate
byte order will result. Then again, one might want the Z format
output in native byte order (especially on PDP-11 and VAX systems).
VAX/VMS Fortran since the beginning supported IBM style Z constants.
The constants are hex integer values, stored in little endian format,
and then treated as VAX floating point values.
> That may be enough reason for Fortran not defining this usage.
Well, for printed output the reader should realize that they
are system dependent values. The ability to print them doesn't
need to be system dependent.
-- glen
| |
|
| > At least for a while storage_size was proposed for f2008, but I haven't
> kept track of all the ins and outs.
It is in the current draft.
--
FX
| |
| Dan Nagle 2008-03-19, 7:15 pm |
| Hello,
On 2008-03-19 16:45:27 -0400, nospam@see.signature (Richard Maine) said:
> glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
>
There is in f08. See 08-007r2 15.2.3.7 at [429].[color=darkred]
>
> Well, there isn't even a Fortran sizeof. I proposed one (which I called
> storage_size), which I thought was awfully simple, but I was too late in
> the f2003 process to sucessfully sell it.
Well, storage_size() applies to the targets of polymorphic pointers,
so some Actual Thought (tm) was required. And there are
some issues of what to do when the item is register-only
(for example, x86 80-bit real). Another complication
is the scalar versus array element issue. So it wasn't all *that* easy.
>
> I found it that f2003 instead introduced a "new" feature that only
> supported f77 style. Seemed like a step backwards to me. In particular,
> f2003 introduced named constants for the size of the numeric storage
> unit and the character storage unit. So if you restricted yourself to
> the types and kinds that were in f77 (those are the ones that have sizes
> in terms of numeric and character storage units), you can find their
> sizes. A bit of a pet peeve of mine.
But character storage units, file storage units, and numeric storage units
are all defined in the standard. Not that their sizes are specified
(that's what the constants provide), but they exist within the standard.
Other items take an undefined storage unit. That is a little tougher.
>
> At least for a while storage_size was proposed for f2008, but I haven't
> kept track of all the ins and outs. I do recall that it kept getting
> shot down by arguments that I thought completely bogus that kept trying
> to make it out to be far more complicated than it was. I had specified
> it to apply to only the already well-defined cases (the ones where the
> standard already said there was a storage unit of some size, but didn't
> provide a way to query the size of). Seemed like some people argued that
> it would be useless unless it was extended to cover more ill-defined
> cases, but then it would be too complicated. Alas. Of course, I'm sure
> the real reason it got shot down is that nobody (such as myself) was at
> the relevant meetings actively pushing it.
Alek, some others, and I pushed it. See 13.7.159 at [388]
in 08-007r2. (And Notes 13.20, 13.21!)
One of *my* pet peeves is when someone says "this is too easy". ;-) ;-)
The toughest part was to get *both* c_sizeof() and storage_size()
into the standard at the same time. Some argued they were redundant.
--
Cheers!
Dan Nagle
| |
| glen herrmannsfeldt 2008-03-19, 7:15 pm |
| Dan Nagle wrote:
(snip regarding c_sizeof)
> Alek, some others, and I pushed it. See 13.7.159 at [388]
> in 08-007r2. (And Notes 13.20, 13.21!)
> One of *my* pet peeves is when someone says "this is too easy". ;-) ;-)
> The toughest part was to get *both* c_sizeof() and storage_size()
> into the standard at the same time. Some argued they were redundant.
Once you get passed redundant, I would think it would be easier
to get them in together. The work to add them (both to the standard
and to compilers) likely overlaps. There are enough redundant
functions already that adding another one shouldn't be so bad.
-- glen
| |
| Dan Nagle 2008-03-19, 7:15 pm |
| Hello,
On 2008-03-19 19:11:49 -0400, glen herrmannsfeldt <gah@ugcs.caltech.edu> said:
> Dan Nagle wrote:
>
> Once you get passed redundant, I would think it would be easier
> to get them in together. The work to add them (both to the standard
> and to compilers) likely overlaps. There are enough redundant
> functions already that adding another one shouldn't be so bad.
The issue with c_sizeof() is that it depends upon the C compiler,
so potentially, at least, it changes. Perhaps it changes depending
upon the C compiler's options.
--
Cheers!
Dan Nagle
| |
| glen herrmannsfeldt 2008-03-19, 10:14 pm |
| Dan Nagle wrote:
(snip)
> The issue with c_sizeof() is that it depends upon the C compiler,
> so potentially, at least, it changes. Perhaps it changes depending
> upon the C compiler's options.
Yes, but so do other C interoperability features. Most important,
it depends on the number of bits in a char, with sizeof(char) by
definition always one.
-- glen
| |
| Dan Nagle 2008-03-20, 7:15 pm |
| Hello,
On 2008-03-19 22:44:44 -0400, glen herrmannsfeldt <gah@ugcs.caltech.edu> said:
>
> Yes, but so do other C interoperability features.
That doesn't make defining c_sizeof() any easier.
> Most important,
> it depends on the number of bits in a char, with sizeof(char) by
> definition always one.
And since a C compiler might support different size chars
as default according to the options used, Fortran must account
for that.
The point being, that the sizes of things have more subtlity
than you might think, especially if you start with a "size is easy" view.
And the standard must work in *all* cases, not just the easy ones.
--
Cheers!
Dan Nagle
|
|
|
|
|