! f(a,b) = 21.b.b - 2.a.a +55.b.b.b.b -10.a.a.b.b +a/2b
! Calculates the Rump polynomial #2 (Rump 2010):
! f(a,b) = 21b^2 - 2a^2 +55b^4 -10a^2b^2 +a/2b
! Evaluate at: a = 77617 and y = 33096.
! The true sum of the first 4 terms is -2 !

! Correct value: -0.82739 60599 46821 36814 11650 95479 81629 ... (with ~40 decimal places precision)
!
! Compile with: Optimization = NONE
! Single, Double, Extended, and Quadruple Precision
!
! Pentru compilarea cu CVF:
! ------------------------
! Faceti comentarii instructiunile referitoare al tipul real(16)
!
use dflib
implicit real(8) (d)
implicit real(16) (q)
logical pressed


f(a,b) = 21.*b*b - 2.*a*a +55.*b*b*b*b -10.*a*a*b*b +a/(2*b)

f1(b) =21.*b*b
f2(a,y) =-2.*a*a;
f3(b) =55.*b*b*b*b
f4(a,b) =-10.*a*a*b*b

df(da,db) = 21.*db*db - 2.*da*da +55.*db*db*db*db -10.*da*da*db*db +da/(2*db)


df1(b) =21.*db*db
df2(a,b) =-2.*da*da;
df3(b) =55.*db*db*db*db
df4(a,b) =-10.*da*da*db*db

qf(qa,qb) = 21.*qb*qb - 2.*qa*qa +55.*qb*qb*qb*qb -10.*qa*qa*qb*qb +qa/(2*qb)


qf1(b) =21.*qb*qb
qf2(a,b) =-2.*qa*qa;
qf3(b) =55.*qb*qb*qb*qb
qf4(a,b) =-10.*qa*qa*qb*qb


open(1, file ='rump II.rez', status ='unknown', form ='formatted')
write(1,*) 'Rump''s Polynomial #2 (2010)'
write(1,'(a)') ' f(a,b) = 21.*b*b - 2.*a*a +55.*b*b*b*b -10.*a*a*b*b +a/(2*b)'
write(1,'(a)') ' Evaluated at: a = 77617.; b = 33096.'
write(1,'(a)') ' The sum of the first 4 terms should be -2 (true value)'

write(1, '(/,a)') ' * Single Precision (24):'

a = 77617_4; b = 33096_4;

!rp =f(a,b)
write(1,'(a, 1pg15.6)') ' rp(a,b) =', f(a,b)

rt1 =f1(b)
rt2 =f2(a,b)
rt3 =f3(b)
rt4 =f4(a,b)

write(1,*)
write(1,'(a, 1pg15.6)') ' rt1:', rt1
write(1,'(a, 1pg15.6)') ' rt2:', rt2
write(1,'(a, 1pg15.6)') ' rt3:', rt3
write(1,'(a, 1pg15.6)') ' rt4:', rt4

write(1, '(/, a, 1pG15.6)') ' rt1 + rt2 + rt3 + rt4 =', rt1 + rt2 + rt3 + rt4

write(1, '(/, a, 1pG15.6)') ' rt4 + rt3+ rt2 + rt1 =', rt4 +rt3 +rt2 +rt1 


write(1, '(//,a)') ' * Double Precision (53):'

da = 77617_8; db = 33096_8;

!rp =df(a,b)
write(1,'(a, 1pg24.15)') ' rp(a,b) =', df(da,db)

drt1 =df1(db)
drt2 =df2(da,db)
drt3 =df3(db)
drt4 =df4(da,db)

write(1,*)
write(1,'(a, 1pg24.15)') ' rt1:', drt1
write(1,'(a, 1pg24.15)') ' rt2:', drt2
write(1,'(a, 1pg24.15)') ' rt3:', drt3
write(1,'(a, 1pg24.15)') ' rt4:', drt4

write(1,*)
write(1,'(a, 1pg24.15)') ' rt1 + rt2 + rt3 + rt4 =', drt1 +drt2 +drt3 +drt4

write(1,*)
write(1,'(a, 1pg24.15)') ' rt4+ rt3 + rt2 + rt1 =', drt4 +drt3 +drt2 +drt1

call FPprecis(FPCW$64)
write(1, '(//,a)') ' * Extended Precision (64):'

a = 77617.; b = 33096.;

!rp =f(a,b)
write(1,'(a, 1pg30.19)') ' rp(a,b) =', f(a,b)

rt1 =f1(b)
rt2 =f2(a,b)
rt3 =f3(b)
rt4 =f4(a,b)

write(1,*)
write(1,'(a, 1pg30.19)') ' rt1:', rt1
write(1,'(a, 1pg30.19)') ' rt2:', rt2
write(1,'(a, 1pg30.19)') ' rt3:', rt3
write(1,'(a, 1pg30.19)') ' rt4:', rt4

write(1, '(/, a, 1pG30.19)') ' rt1 + rt2 + rt3 + rt4 =', rt1 + rt2 + rt3 + rt4

write(1, '(/, a, 1pG30.19)') ' rt4 + rt3+ rt2 + rt1 =', rt4 +rt3 +rt2 +rt1 

! Quadruple:
write(1, '(//,a)') ' * Quadruple Precision (113):'

qa = 77617._16; qb = 33096._16;

!rp =qf(a,b)
write(1,'(a, 1pg30.24)') ' rp(a,b) =', qf(qa,qb)

qrt1 =qf1(qb)
qrt2 =qf2(qa,qb)
qrt3 =qf3(qb)
qrt4 =qf4(qa,qb)

write(1,*)
write(1,'(a, 1pg24.15)') ' rt1:', qrt1
write(1,'(a, 1pg24.15)') ' rt2:', qrt2
write(1,'(a, 1pg24.15)') ' rt3:', qrt3
write(1,'(a, 1pg24.15)') ' rt4:', qrt4

write(1,*)
write(1,'(a, 1pg24.15)') ' rt1 + rt2 + rt3 + rt4 =', qrt1 +qrt2 +qrt3 +qrt4

write(1,*)
write(1,'(a, 1pg24.15)') ' rt4+ rt3 + rt2 + rt1 =', qrt4 +qrt3 +qrt2 +qrt1

! Correct results:
!10 ^        -1 x -8.273960599468213681411650954798162919
write(1,'(//, a,/,a)') ' ** Correct value:',' -0.82739 60599 46821 36814 11650 95479 81629 19 ... (with 36 decimal digits precision)'

write(1, '(/, a, 1pG30.16)') 'Double (Exact): a/(2*b)-2 = ', da/(2*db)-2D0
write(1, '(/, a, 1pG36.30)') 'Quadruple (Exact): a/(2*b)-2 = ', qa/(2*qb)-2Q0

write(*,*) ' Output File: "rump II.rez"'

write(*, '(/,a)') 'Apasa orice tasta'
DO WHILE (.NOT. pressed)   
pressed = PEEKCHARQQ ( )
END DO

END
