Precision query - how to solve KIND=16

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):

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

REAL(KIND=3) can be used for 32 bit applications. This gives 80 bit extended precision.

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

Quoted from PaulLaidler REAL(KIND=3) can be used for 32 bit applications. This gives 80 bit extended precision.

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).

In this line

y = 1 - 2 * y * x;

the default is 32-bit. Since you wish 80-bit, then define two parameters as

real (kind=3),parameter:: one=1.0_3,two=2._3

and recode this line as

y = one - two * y * x;

or

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.

Just a note to 'clarify'. If you use a (real4) constant, such as in (real8) 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

  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.