forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Precision query - how to solve KIND=16

 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> General
View previous topic :: View next topic  
Author Message
arctica



Joined: 10 Sep 2006
Posts: 105
Location: United Kingdom
PostPosted: Tue Dec 05, 2023 3:04 pm    Post subject: Precision query - how to solve KIND=16 Reply with quote
Hello

Is there a way to get the real(kind=16) or a workaround for this in ftn95?

The following code works if real(kind=2) (i.e. double precision):

Code:
program feigenbaum
      implicit none

      integer i, j, k
      real(kind=2) :: x, y, a, a1, a2, d1 ! kind=2 double precision
      !real ( KIND = 16 ) x, y, a, b, a1, a2, d1

      print '(a4,a13)', 'i', 'd'

      a1 = 1.0;
      a2 = 0.0;
      d1 = 3.2;

      do i=2,20
         a = a1 + (a1 - a2) / d1;
         do j=1,10
            x = 0
            y = 0
            do k=1,2**i
                y = 1 - 2 * y * x;
                x = a - x**2;
            end do
            a = a - x / y;
         end do

         d1 = (a1 - a2) / (a - a1);
         a2 = a1;
         a1 = a;
         print '(i4,f13.10)', i, d1
     end do
     end



I cannot reproduce the results below with real(kind=2):

i d
2 3.2185114220
3 4.3856775986
4 4.6009492765
5 4.6551304954
6 4.6661119478
7 4.6685485814
8 4.6690606606
9 4.6691715554
10 4.6691951560
11 4.6692002291
12 4.6692013133
13 4.6692015458
14 4.6692015955
15 4.6692016062
16 4.6692016085
17 4.6692016090
18 4.6692016091
19 4.6692016091
20 4.6692016091

Thanks
Lester
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 7927
Location: Salford, UK
PostPosted: Tue Dec 05, 2023 3:50 pm    Post subject: Reply with quote
REAL(KIND=3) can be used for 32 bit applications.
This gives 80 bit extended precision.
Back to top
View user's profile Send private message AIM Address
arctica



Joined: 10 Sep 2006
Posts: 105
Location: United Kingdom
PostPosted: Tue Dec 05, 2023 5:22 pm    Post subject: Re: Reply with quote
Get a bit closer to the published value at iteration 14 using kind=3 and win32 option:

i d
2 3.21851142203808791280
3 4.38567759856833909170
4 4.60094927653807532040
5 4.65513049539197968780
6 4.66611194782857030690
7 4.66854858144686456730
8 4.66906066064807793730
9 4.66917155538151432430
10 4.66919515602275888650
11 4.66920022908290646170
12 4.66920131315636925260
13 4.66920154625723728170
14 4.66920159550432773820
15 4.66920148627286091600
16 4.66920254043202110470
17 4.66919736236123777780
18 4.66920562189315892730
19 4.66896400405368285560
20 4.67168279116043707170

PaulLaidler wrote:
REAL(KIND=3) can be used for 32 bit applications.
This gives 80 bit extended precision.
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 7927
Location: Salford, UK
PostPosted: Tue Dec 05, 2023 5:52 pm    Post subject: Reply with quote
At a quick glance it looks like you may be dividing one very small value by another very small value. If you are then there will be problems with round-off error (no matter what the precision).
Back to top
View user's profile Send private message AIM Address
wahorger



Joined: 13 Oct 2014
Posts: 1217
Location: Morrison, CO, USA
PostPosted: Wed Dec 06, 2023 12:56 am    Post subject: Reply with quote
In this line
Code:
y = 1 - 2 * y * x;

the default is 32-bit. Since you wish 80-bit, then define two parameters as
Code:
real (kind=3),parameter:: one=1.0_3,two=2._3

and recode this line as
Code:
y = one - two * y * x;

or
Code:
y = 1.0_3 - 2.0_3 * y * x;


This should help.

[edited]
BTW, it turns out if you put, say, 1.1 as a constant without the _3, the value of the constant will not be 1.1, rather 1.1000000238418579.
Back to top
View user's profile Send private message Visit poster's website
JohnCampbell



Joined: 16 Feb 2006
Posts: 2554
Location: Sydney
PostPosted: Wed Dec 06, 2023 7:28 am    Post subject: Reply with quote
Just a note to "clarify".
If you use a (real*4) constant, such as in (real*8) d1 = 3.2, this will set the extended bits beyond 3.2_1 to zero.
This gives a different value where the extended bits are non zero for the extended precision constant.
The most common occurrence for this is the constant 0.1

Code:
  real*10 x,y,z
!  the following are different
  x = 1.1
  y = 1.1d0
  z = 1.1_3
  write (*,*) x,y,z,'  all different values'

!  the following are the same
  x = 1.0
  y = 1.0_3
  z = 1
  write (*,*) x,y,z,'  all same values'
end


Statements like "X = 2" will not cause any problems as the real*10 value has zero bits in the extended precision bytes. You just have to check what the default kind is during the computation.
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> General All times are GMT + 1 Hour
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group