Taxicab Numbers, Fortran-95

Demos for Beginners, by Hannah

back to projects

Taxicab numbers (wiki), named for a story of the continued-fraction and number theory genius Srinivasa Ramanujan correcting a fellow mathematician's assertion that 1729 was "rather a dull [number]" for a cab. This is my first Fortran program, and it's not the most efficient method, but it's simple enough.

Example Program Runs

Code

Example Program Runs

Enter minimum Taxi number degree:2
Enter search range: 1-10000
Finding Taxicab Numbers from 1 to        10000 with minimum degree of Ta(2):
        1729, Ta(2) =      1^3 +    12^3 =      9^3 +    10^3
        4104, Ta(2) =      9^3 +    15^3 =      2^3 +    16^3
Enter minimum Taxi number degree:3
Enter search range: 1-100000000
Finding Taxicab Numbers from 1 to    100000000 with minimum degree of Ta(3):
    87539319, Ta(3) =    228^3 +   423^3 =    255^3 +   414^3 =    167^3 +   436^3

This one takes about half a minute.

Code

! taxicab.f95 
!       find taxicab numbers
!       by Hannah Leitheiser
!       2018-08-27
! compile and run: gfortran taxicab.f95 -o taxicab
!  (in Linux)      ./taxicab

PROGRAM find_taxicab_numbers

! Purpose:
!     Find positive INTEGERs that can be represented 
!	by the sums of at least two pairs of INTEGER 
!	cubes. 
!     
!     T = a^3 + b^3 = c^3 + d^3
!      where a, b, c, d are unique positive
!      integers.

! Steps:
!    1) Get user parameters.
!    2) Create array of pairs of cubes and their sum
!    3) Sort table by sums
!    4) Find duplicate sums in the table.

IMPLICIT NONE

INTERFACE

SUBROUTINE print_taxicab_numbers(a, b, sum_of_cubes, i, degree)

INTEGER, DIMENSION(:), INTENT(IN) :: a, b, sum_of_cubes
INTEGER, INTENT(IN) :: degree, i
INTEGER :: ii

END SUBROUTINE print_taxicab_numbers

END INTERFACE

! pseudo-parameters, should remain unchanged after user
! supplies them.

INTEGER :: SEARCH_MAXIMUM, MINIMUM_DEGREE

! derived pseudo-parameters

INTEGER :: MAXIMUM_CUBE, ARRAY_SIZE

! a ^ 3 + b ^ 3 = sum_of_cubes

INTEGER, DIMENSION ( : ), ALLOCATABLE  :: a, b, sum_of_cubes
INTEGER :: size = 1, i, ii, previous_sum
LOGICAL :: none_found = .TRUE.

! ### Step 1) Get User Parameters. ###

WRITE (*, 1002)
READ (*,*) MINIMUM_DEGREE
1002 FORMAT ('Enter minimum Taxi number degree:', $)

WRITE (*, 1003)
READ (*,*) SEARCH_MAXIMUM
1003 FORMAT ('Enter search range: 1-', $)

! Needn't search beyond the cube root of the search maximum.
! ARRAY_SIZE will be a triangle number with base MAXIMUM_CUBE

MAXIMUM_CUBE = &
      CEILING( REAL(SEARCH_MAXIMUM) ** (1./3.))
ARRAY_SIZE = &
      (MAXIMUM_CUBE * (MAXIMUM_CUBE + 1)) / 2

ALLOCATE( a(ARRAY_SIZE) )
ALLOCATE( b(ARRAY_SIZE) )
ALLOCATE( sum_of_cubes(ARRAY_SIZE) )

WRITE (*,1000) SEARCH_MAXIMUM, MINIMUM_DEGREE
1000 FORMAT ('Finding Taxicab Numbers from 1 to ', &
      I12, ' with minimum degree of Ta(', I1, '):' )

! ### Step 2) Create array of pairs of cubes and their sum. ###

DO i=1,MAXIMUM_CUBE
    DO ii=i,MAXIMUM_CUBE
        a(size) = i
        b(size) = ii
        sum_of_cubes(size) = i ** 3 + ii ** 3
        IF (sum_of_cubes(size) <= SEARCH_MAXIMUM) THEN
            size = size + 1
        END IF
    END DO
END DO

size = size - 1

! ### Step 3) Sort table by their sums (bubble sort) ###

DO i=1, size
    DO ii=i, size
        IF (sum_of_cubes(ii) < sum_of_cubes(i)) THEN
            CALL swap( a(i), a(ii))
            CALL swap( b(i), b(ii))
            CALL swap( sum_of_cubes(i), sum_of_cubes(ii) )
            END IF
      END DO
END DO

! ### Step 4) Find duplicate sums in the table. ###

previous_sum =  -1
ii = 0
DO i=1,size
    IF (sum_of_cubes(i) == previous_sum) THEN
        ii = ii + 1
    ELSE
        IF (ii >= MINIMUM_DEGREE - 1) THEN
            none_found = .FALSE.
            ! Special case, if we want the first number.
            IF (i == 1) THEN
                CALL print_taxicab_numbers( &
                        a, b, sum_of_cubes, 1, ii+1)
                none_found = .FALSE.
            ! Normal Case
            ELSE
                CALL print_taxicab_numbers( &
                          a, b, sum_of_cubes, i-1, ii+1)
                none_found = .FALSE.
            END IF
        END IF
        ii = 0
        previous_sum = sum_of_cubes(i)
    END IF
END DO

! Special case, if we need to print the last number.
IF (ii >= MINIMUM_DEGREE - 1) THEN
    none_found = .FALSE.
    call print_taxicab_numbers( a, b, sum_of_cubes, size, ii+1)
END IF

IF (none_found .EQV. .TRUE.) THEN
    WRITE (*,*) &
         "No taxicab numbers found for seach specifications."
END IF

END PROGRAM find_taxicab_numbers

! ---------------------- SUBROUTINE swap -----------------------
! swaps values in a and b

SUBROUTINE swap (a, b)
IMPLICIT NONE

INTEGER, INTENT(INOUT) :: a, b
INTEGER :: intermediate
intermediate = a
a = b
b = intermediate

END SUBROUTINE swap

! ----------- SUBROUTINE print_taxicab_numbers ------------------
! print the taxicab number and the cubes that form it
! Example output:
!     87539319, Ta(3) =    228^3 +   423^3 =
!(...)                     255^3 +   414^3 =    167^3 +   436^3

SUBROUTINE print_taxicab_numbers (a, b, sum_of_cubes, i, degree)

INTEGER, DIMENSION(:), INTENT(IN) :: a, b, sum_of_cubes
INTEGER, INTENT(IN) :: degree, i
INTEGER :: ii

WRITE (*,1001) sum_of_cubes(i), degree
1001 FORMAT (I12, ', Ta(', I1, ') = ', $)
DO ii = degree, 1, -1
        WRITE (*, 1004) a(i-ii+1), b(i-ii+1)
        1004 FORMAT ( I6, "^3 +", I6, "^3", $)
        if (ii > 1) THEN
            WRITE (*, 1006)
        END IF
        1006 FORMAT ( " = ", $)
END DO
write (*,*)
END SUBROUTINE print_taxicab_numbers