FortranCode/equal_area_triples.f90 at main · Beliavsky/FortranCode · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
! https://www.johndcook.com/blog/2025/07/30/pythagorean-triangles/
! Primitive Pythagorean triangles with the same area
! Posted on 30 July 2025 by John D. Cook
! equal_area_triples.f90
!
! search for three distinct primitive pythagorean triples
!   (a,b,c) with max(a,b,c) < 20000
! that share the same right-triangle area  area = a*b/2
!
! a primitive triple is generated by euclid's formula
!   a = m**2 - n**2 ,  b = 2*m*n ,  c = m**2 + n**2
!   with m > n > 0, gcd(m,n)=1, and opposite parity
!
! the program uses a simple array-based map
!   area -> (count, up to three triples)
! and halts once the first area reaches count = 3

program equal_area_triples
implicit none

integer, parameter :: i8 = selected_int_kind(12)   ! 64-bit integers
integer, parameter :: max_side   = 20000
integer, parameter :: max_sets   = 30000           ! plenty of room

type triple_set
   integer(kind=i8) :: area = 0_i8
   integer           :: count = 0
   integer           :: a(3) = 0, b(3) = 0, c(3) = 0
end type triple_set

type(triple_set), dimension(max_sets) :: sets
integer :: n_sets = 0

integer :: m_max, m, n
integer :: a, b, c, t, idx
integer(kind=i8) :: area
logical :: found

! -------------------------------- main search loop --------------------------------
m_max = int(sqrt(real(max_side))) + 1

do m = 2, m_max
   do n = 1, m - 1
      if (mod(m - n, 2) == 0) cycle          ! need opposite parity
      if (igcd(m, n) /= 1)  cycle            ! need gcd = 1

      a = m * m - n * n
      b = 2 * m * n
      c = m * m + n * n

      if (a <= 0 .or. b <= 0) cycle
      if (a > b) then                        ! keep a <= b
         t = a ; a = b ; b = t
      end if
      if (a >= max_side .or. b >= max_side .or. c >= max_side) cycle

      area = int(a,kind=i8) * int(b,kind=i8) / 2

      ! look for this area in the existing list
      found = .false.
      do idx = 1, n_sets
         if (sets(idx)%area == area) then
            found = .true.
            exit
         end if
      end do

      if (.not. found) then
         n_sets = n_sets + 1
         sets(n_sets)%area  = area
         sets(n_sets)%count = 0
         idx = n_sets
      end if

      if (sets(idx)%count < 3) then
         sets(idx)%count            = sets(idx)%count + 1
         sets(idx)%a(sets(idx)%count) = a
         sets(idx)%b(sets(idx)%count) = b
         sets(idx)%c(sets(idx)%count) = c
      end if

      if (sets(idx)%count == 3) then
         print *, 'found area with three primitive triples:'
         print *, 'area = ', sets(idx)%area
         do t = 1, 3
            print *, '  (', sets(idx)%a(t), ',', &
                          sets(idx)%b(t), ',', &
                          sets(idx)%c(t), ')'
         end do
         stop
      end if
   end do
end do

print *, 'no area shared by three primitive triples found.'
stop

! --------------------------------- helper routine ---------------------------------
contains

integer function igcd(x, y)
   implicit none
   integer, intent(in) :: x, y
   integer :: a, b, r
   a = abs(x)
   b = abs(y)
   do while (b /= 0)
      r = mod(a, b)
      a = b
      b = r
   end do
   igcd = a
end function igcd

end program equal_area_triples