! Calculeaza "polinomul" Rump:
! f(x; y) = 333.75y^6 + x^2*(11x^2y^2 -y^6 -121y^4 -2) + 5.5y^8 +x/2y
!
! In: x = 77617 and y = 33096.
! Suma primilor 3 terneni ar trebui sa fie -2 (valoarea adevarata)!

! Valoarea corecta: -0.82739 60599 46821 36814 11650 95479 81629 ... (cu precizie de ~40 cifre zecimale)
!
! Compileaza cu: Optimization: Disable
! Precizie: Simpla; Dubla; Extinsa; Cvadrupla;
!
!

use dflib
implicit real(8) (d)
implicit real(16) (q)
logical pressed


f(x,y) = 333.75*y**6 + x**2*(11*(x**2)*y**2 -y**6 -121*y**4 -2) + 5.5*y**8 +x/(2*y)

f1(y) =333.75*y**6
f20(x,y) =11*(x**2)*y**2 -y**6 -121*y**4 -2; f2(x,y) =x*x*f20(x,y);
f3(y) =5.5*y**8
f4(x,y) =x/(2.*y)

df(dx,dy) = 333.75*dy**6 + dx**2*(11*(dx**2)*dy**2 -dy**6 -121*dy**4 -2) + 5.5*dy**8 +dx/(2*dy)


df1(dy) =333.75*dy**6
df20(dx,dy) =11*(dx**2)*dy**2 -dy**6 -121*dy**4 -2; df2(dx,dy) =dx*dx*df20(x,y);
df3(dy) =5.5*dy**8
df4(dx,dy) =dx/(2.*dy)

qf(x,y) = 333.75*qy**6 + qx**2*(11*(qx**2)*qy**2 -qy**6 -121*qy**4 -2) + 5.5*qy**8 +qx/(2*qy)

qf1(y) =333.75*qy**6
qf20(x,y) =11*(qx**2)*qy**2 -qy**6 -121*qy**4 -2; qf2(qx,qy) =qx*qx*qf20(qx,qy);
qf3(y) =5.5*qy**8
qf4(x,y) =qx/(2.*qy)

open(1, file ='rump.rez', status ='unknown', form ='formatted')
write(1,*) 'Polinomul Rump'
write(1,'(a)') ' rp(x,y) = 333.75y^6 + x^2*(11x^2*y^2 -y^6 -121y^4 -2) + 5.5y^8 +x/(2y)'
write(1,'(/,a)') ' Calculat in: x = 77617.; y = 33096.'
write(1,'(a)') ' Suma primilor 3 terneni ar trebui sa fie -2 (valoarea adevarata) !'

write(1, '(/,a)') ' * Precizie Simpla (24):'

x = 77617_4; y = 33096_4;

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

rt1 =f1(y)
rt2 =f2(x,y)
rt3 =f3(y)
rt4 =f4(x,y)

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)') ' * Precizie Dubla (53):'

dx = 77617_8; dy = 33096_8;

!rp =df(x,y)
write(1,'(a, 1pg24.15)') ' rp(x,y) =', df(dx,dy)

drt1 =df1(dy)
drt2 =df2(dx,dy)
drt3 =df3(dy)
drt4 =df4(dx,dy)

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 =', drt1 +drt2 +drt3

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

call FPprecis(FPCW$64)
write(1, '(//,a)') ' * Precizie Extinsa (64):'

x = 77617.; y = 33096.;

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

rt1 =f1(y)
rt2 =f2(x,y)
rt3 =f3(y)
rt4 =f4(x,y)

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 =', rt1 + rt2 + rt3

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

! Quadrupla:
write(1, '(//,a)') ' * Precizie Cvadrupla (113):'

qx = 77617._16; qy = 33096._16;

!rp =qf(x,y)
write(1,'(a, 1pg24.15)') ' rp(x,y) =', qf(qx,qy)

qrt1 =qf1(qy)
qrt2 =qf2(qx,qy)
qrt3 =qf3(qy)
qrt4 =qf4(qx,qy)

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 =', qrt1 +qrt2 +qrt3

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

! Correct results:
write(1,'(//, a,/,a)') ' ** Valoarea corecta:',' -.82739605994682136814116509547981629199903311 (43 cifre zecimale de precizie)'

!write(1, '(/, a, 1pG30.19)') ' x/(2*y)-2 =', dx/(2*dy)-2
write(1, '(/, a, 1pG22.16)') 'Dubla (Exact):    -2 + x/(2*y) = ', dx/(2*dy)-2
write(1, '(/, a, 1pG36.30)') 'Cvadrupla (Exact): -2 + x/(2*y) = ', qx/(2*qy)-2Q0


write(*, '(/,a)') ' Fisier de iesire: "rump.rez"'

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

END
