53 lines
903 B
Fortran
53 lines
903 B
Fortran
PROGRAM Pell
|
|
IMPLICIT NONE
|
|
INTEGER n(4)
|
|
DATA n /61, 109, 181, 277/
|
|
INTEGER*16 x, y
|
|
INTEGER i
|
|
|
|
DO 10, i = 1, size(n)
|
|
CALL solve(n(i), x, y)
|
|
10 PRINT 100, n(i), x, y
|
|
|
|
100 FORMAT ('x² - ', i3, 'y² = 1; x = ', i24, ' and y = ', i24)
|
|
END PROGRAM
|
|
|
|
|
|
SUBROUTINE solve(n, a, b)
|
|
IMPLICIT INTEGER*16 (A-Z)
|
|
INTEGER n
|
|
|
|
x = INT(sqrt(REAL(n)))
|
|
IF (x**2 .ne. n) GOTO 10
|
|
|
|
a = 1 ! perfect square; only trivial solution exists
|
|
b = 0
|
|
RETURN
|
|
|
|
10 y = x
|
|
z = 1
|
|
r = 2*x
|
|
e1 = 1
|
|
e2 = 0
|
|
f1 = 0
|
|
f2 = 1
|
|
|
|
20 y = r*z - y
|
|
z = (n - y**2) / z
|
|
r = (x + y) / z
|
|
|
|
e3 = r*e2 + e1
|
|
e1 = e2
|
|
e2 = e3
|
|
|
|
f3 = r*f2 + f1
|
|
f1 = f2
|
|
f2 = f3
|
|
|
|
a = e2 + x*f2
|
|
b = f2
|
|
|
|
IF (a**2 - n*b**2 .eq. 1) RETURN
|
|
GOTO 20
|
|
|
|
END SUBROUTINE
|