How do I do big integers in Fortran?

I need to generate some big integers. See example below.

Input   Result    
 40 165580141    
 80 37889062373143906    
120 8670007398507948658051921    
160 1983924214061919432247806074196061    
200 453973694165307953197296969697410619233826

Here is my Fortran code:

program cycle
    use iso_fortran_env
  implicit none
  character(200) :: str
  integer :: n
    integer(kind=int64) :: x1, result, x2, x3

    do n = 40, 500, 40
        x1 = n
        result = 1
        x2 = 0
        x3 = 1
        do 
            if (x1 > 1) then
                x2 = result
                result = result + x3 
                x3 = x2     
                x1 = x1 - 1
            else
                exit
            end if
        end do
        write(str,'(i64)')  result
        print *, n, adjustl(str)
    end do
end program cycle

Here is the sample output:

 40 165580141 
 80 37889062373143906 
120 790376311979428689  
160 9217463444206948445 
200 3721511182311577122

As you can see, it gets the first two numbers right but the rest are beyond the reach of 64 bit integers. I've looked at other questions (1) but I'm really interested in a simple way, preferably built into the language itself. In Ruby and Go I have little trouble. Am I overlooking something obvious in Fortran? Is there a better option I can use in my code?

Solution

There is no built-in "big number" support, but we can first check whether there is a larger integer kind available (as mentioned by Francescalus above and also many previous pages (e.g. this page). On my computer with gfortran-6.1, the compiler seems to support 128-bit integer kind, so I could calculate the result up to n=160 or so.

program cycle
...
integer, parameter :: verylong = selected_int_kind(32)
integer(verylong) :: x1, result, x2, x3

print *, "int32 = ", int32   !! from iso_fortran_env
print *, "int64 = ", int64
print *
print *, "kind..(16) => ", selected_int_kind(16)  !! 8
print *, "kind..(32) => ", selected_int_kind(32)  !! 16
print *, "kind..(40) => ", selected_int_kind(40)  !! -1 (not available)
print *, "kind..(64) => ", selected_int_kind(64)  !! -1 (not available)
print *
print *, "sizeof(x1)       = ", sizeof(x1), "(bytes)"       !! GNU extension
print *, "storage_size(x1) = ", storage_size(x1), "(bits)"  !! F2008
print *, "huge(x1)         = ", huge(x1)                    !! largest integer
...

Results:

 int32 =            4
 int64 =            8

 kind..(16) =>            8
 kind..(32) =>           16
 kind..(40) =>           -1
 kind..(64) =>           -1

 sizeof(x1)       =                    16 (bytes)
 storage_size(x1) =          128 (bits)
 huge(x1)         =  170141183460469231731687303715884105727

 n=          40 res= 165580141
 n=          80 res= 37889062373143906
 n=         120 res= 8670007398507948658051921
 n=         160 res= 1983924214061919432247806074196061
 n=         200 res= 37016692776042937155243383431825151522
 n=         240 res= -159769225356713774587328406036589956191
 ...

While there is no built-in "BigInt" type, it is rather straightforward to use an external library (for example, fmlib linked from this page). Since various operators and assignment are overloaded, almost no modifications are necessary to your codes.

Procedures:

1) Download the file FMfiles.zip and extract FM.95, FMZM90.f95, and FMSAVE.f95

2) Make a library file as

gfortran -c -O2 FMSAVE.f95 FMZM90.f95 FM.f95
ar rv fmlib.a FM*.o

3) Modify your code as follows (the modified parts are marked with arrows).

program cycle
    use FMZM           !<----- a module for handling big numbers
    implicit none
    character(200) :: str
    integer :: n
    type(IM) :: x1, result, x2, x3     !<----- IM = BigInt, FM = BigFloat

    do n = 40, 500, 40
        x1 = n
        result = 1
        x2 = 0
        x3 = 1
        do 
            if (x1 > 1) then
                x2 = result
                result = result + x3 
                x3 = x2     
                x1 = x1 - 1
            else
                exit
            end if
        end do
        str = IM_format( 'i200', result )   !<----- convert BigInt to string
        print *, n, trim( adjustl(str) )    !<----- print huge integers
    end do
end program cycle

4) Compile (assuming that "test.f90" is the above code):

gfortran test.f90 fmlib.a
./a.out

5) Results

   n result
  40 165580141
  80 37889062373143906
 120 8670007398507948658051921
 160 1983924214061919432247806074196061
 200 453973694165307953197296969697410619233826
 240 103881042195729914708510518382775401680142036775841
 280 23770696554372451866815101694984845480039225387896643963981
 320 5439356428629292972296177350244602806380313370817060034433662955746
 360 1244666864935793005828156005589143096022236302705537193166716344690085611761
 400 284812298108489611757988937681460995615380088782304890986477195645969271404032323901
 440 65172495098135102433647404982700073500075401759827878315356483347951218369680224170989749666
 480 14913169640232740127827512057302148063648650711209401966150219926546779697987984279570098768737999681

We can verify the result by noting that result for n is actually equal to fibonacci(n+1), so for example we have fibonacci(481) for n = 480.