PURE INTEGER FUNCTION ISQRT(I) RESULT(K)
INTEGER, INTENT(IN) :: I
J = I
DO
   K = (J + I / J) / 2
   IF (J - K <= 1) EXIT
   J = K
END DO
END FUNCTION ISQRT
PURE RECURSIVE INTEGER FUNCTION IGCD(I, J) RESULT(K)
INTEGER, INTENT(IN) :: I, J
IF      (I == J) THEN
   K = I
ELSE IF (MOD(I, 2) == 0) THEN
   IF (MOD(J, 2) /= 0) THEN
      K = IGCD(I / 2, J)          ! I even, J odd
   ELSE
      K = 2 * IGCD(I / 2, J / 2)  ! Both even
   END IF
ELSE IF (MOD(J, 2) == 0) THEN
   K = IGCD(I, J / 2)             ! I odd, J even
ELSE IF (I > J) THEN
   K = IGCD((I - J) / 2, J)       ! Both odd, difference even
ELSE
   K = IGCD(I, (J - I) / 2)       ! Both odd, difference even
ENDIF
END FUNCTION IGCD
PROGRAM PYTHAGORAS
CHARACTER(LEN=10) :: ARG, FMT
CALL GET_COMMAND_ARGUMENT(1, ARG)
READ(ARG, *) N
WRITE(FMT, '(A,I2.2,A)') '(3(I',LEN_TRIM(ARG),',A))'
DO I = 1, N
   DO J = 1, I - 1
      L = ISQRT(I**2 - J**2)
      M = MIN(J - 1, L + 1)
      DO K = L, M
         IF (J**2 + K**2 == I**2 .AND. IGCD(J, K) == 1) &
            PRINT FMT, J,'^2 +',K,'^2 =',I,'^2'
      END DO
   END DO
END DO
END PROGRAM PYTHAGORAS