Integer sequences

1

A Hacker News post makes me wonder which sequences should be easily available in Fortran, perhaps in stdlib – maybe primes and Fibonacci numbers to start. It’s interesting that 2,3,5,7,11,13,17,19,23 matches 228 results in the On-line Encyclopedia of Integer Sequences, not just the first N primes. The Encyclopedia provides programs in Mathematica and Maple for sequences. It would be good to have codes in non-proprietary languages such as C and Fortran. Reading the comments at HN, many people have used the Encyclopedia to discover patterns in integers they observed.

2 Likes
2

This would be a great fpm package to provide integer sequences! Feel free to start one and post here, I am sure others will help test / improve it.

1 Like
3

I often visit the site mathworld.wolfram.com for mathematical gems - not always to find something of immediate use :slight_smile:. One of the subjects that intrigues me is that of integer sequences. I would be happy to contribute to such a package!
(One sequence I encountered recently, not on this site but here: The Recaman's sequence, has a rather nice visualisation.)

2 Likes
4

Recaman sequence is a good one to highlight Fortran array capabilities and a decent comparative test case for the FINDLOC intrinsic relative to other approaches used at the Rosetta Code site.where the C++ example makes one really wonder as to where is the formula translation for:
The Recamán’s sequence {isplaystyle a_{0},a_{1},a_{2}ots } is defined as:

{isplaystyle a_{n}={egin{cases}0&&{ext{if }}n=0a_{n-1}-n&&{ext{if }}a_{n-1}-n>0{ext{ and is not already in the sequence}}a_{n-1}+n&&{ext{otherwise}}nd{cases}}}

module intseq_m
contains
   subroutine CalcRecamanSequence( a )
   ! https://en.wikipedia.org/wiki/Recam%C3%A1n%27s_sequence
      integer, intent(inout) :: a(0:)
      integer :: n, an
      logical :: InSequence
      if ( size(a) == 0 ) return
      a(0) = 0
      do n = 1, size(a)-1
         InSequence = .true.
         an = a(n-1) - n
         if ( an > 0 ) then
            InSequence = ( findloc( a(0:n-1), an, dim=1 ) /= 0 )
         end if
         if ( InSequence ) then
            a(n) = a(n-1) + n
         else
            a(n) = an
         end if
      end do
   end subroutine 
end module
   use intseq_m, only : CalcRecamanSequence
   integer :: x( 50 )
   call CalcRecamanSequence( x )
   print "('The first terms of the sequence are: ', (*(g0,1x)))", x
end

C:\Temp>ifort /standard-semantics seq.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.2.0 Build 20210228_000000
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 14.28.29337.0
Copyright (C) Microsoft Corporation. All rights reserved.

-out:seq.exe
-subsystem:console
seq.obj

C:\Temp>seq.exe
The first terms of the sequence are: 0 1 3 6 2 7 13 20 12 21 11 22 10 23 9 24 8 25 43 62 42 63 41 18 42 17 43 16 44 15 45 14 46 79 113 78 114 77 39 78 38 79 37 80 36 81 35 82 34 83

C:\Temp>

3 Likes
5

This inspired me to write a more compact version:

module integer_sequences
    implicit none

contains
function Recaman( n )
    integer, intent(in)                :: n
    integer, dimension(:), allocatable :: Recaman

    integer                            :: i

    if ( n <= 0 ) then
        allocate( Recaman(0) )
        return
    endif

    Recaman = [0]

    do i = 1,n-1
        Recaman = [ Recaman, merge( Recaman(i)-i, Recaman(i)+i, &
                                    Recaman(i)-i > 0 .and. all( Recaman /= Recaman(i)-i ) ) ]
    enddo
end function Recaman

end module integer_sequences

! test --
!     Straightforward test
!
program test_recaman
    use integer_sequences

    implicit none

    write(*,'(10g5.0)' ) Recaman(50)
end program test_recaman
1 Like
6

Now we should benchmark @FortranFan’s and @Arjen’s solutions. The complexity seems O(n^2) for both, but I would expect @FortranFan’s to be a lost faster, as there are no allocations/reallocations.

I think that must be why the Rosetta page has some solutions using a hash table, I guess then the complexity is only O(n).

7

Side remark: An even more functional-style solution would use recursion to avoid the do-loop.

If the performance data are to reliable, it would require a lot more eleemtns of the series.

But a hash table? That would used to avoid looping over the preceding elements to see if the candidate is already present. We could use a bit array for that. Perhaps a nice demonstration of @wclodius 's work?

2 Likes
8

@certik, certainly that would be the case. Most C++ implementations now have highly optimized standard containers that it’s greatly diminishing returns with trying to do better than them.

So with Fortran stdlib in mind, as I had suggested a while ago at the GitHub site, I would suggest a mixed-language approach where certain “core” backend stuff may be in C++, Assembler, etc. with suitable Fortran APIs to consume the Fortran stdlib APIs effectively.

Thus for example, one can have the Recaman Sequence calculator in C++ like so:

#include <set>
#include "ISO_Fortran_binding.h"
using namespace std;

extern "C" void CalcRecamanSequence(CFI_cdesc_t*);

extern "C" void CalcRecamanSequence(CFI_cdesc_t* desc_a) {
// https://en.wikipedia.org/wiki/Recam%C3%A1n%27s_sequence
   int size_a = (int)desc_a->dim[0].extent;
   int* a = (int*)desc_a->base_addr;
   std::set<int> unique_vals{ 0 };
   a[0] = 0;
   int n = 1;
   while (n < size_a ) {
      int an = a[n - 1] - n;
      if (an > 0 && unique_vals.find(an) == unique_vals.end()) {
         a[n] = an;
      }
      else {
         a[n] = a[n - 1] + n;
      }
      if (unique_vals.find(a[n]) == unique_vals.end()) unique_vals.insert(a[n]);
      n++;
   }
   return;
}

And a Fortran wrapper for it that can also provide the functional style API of interest to some:

module IntegerSequences_m
   use, intrinsic :: iso_c_binding, only : c_int
   interface
      subroutine CalcRecamanSequence( a ) bind(C, name="CalcRecamanSequence")
         import :: c_int
         implicit none
         integer(c_int), intent(inout) :: a(:)
      end subroutine
   end interface
contains
   function RecamanSequence( n ) result(r)
   ! Employ C++ implementation
      ! Argument list
      integer, intent(in) :: n
      ! Function result
      integer(c_int), allocatable :: r(:)
      allocate( r(max(n,0)) )
      if (n >=0 ) call CalcRecamanSequence( r )
      return 
   end function RecamanSequence
end module

With your respect to your suggestions for benchmarking, that’s always good. But if the goal is to try to do better than C++, as I indicated above, I think it will take a lot of effort with perhaps only marginal benefits. The above method I show works reasonably well, it’s nearly 100x faster than what I showed upthread using FINDLOC and about 300x faster than the one by @Arjen that suffers from repeated allocations.

C:\temp>cl /c /O2 /EHsc /W3 CalcRecaman.cpp
Microsoft (R) C/C++ Optimizing Compiler Version 19.28.29337 for x64
Copyright (C) Microsoft Corporation. All rights reserved.

CalcRecaman.cpp
C:\Program Files (x86)\Intel\oneAPI\compiler\latest\windows\compiler\include\ISO_Fortran_binding.h(156): warning C4200: nonstandard extension used: zero-sized array in struct/union
C:\Program Files (x86)\Intel\oneAPI\compiler\latest\windows\compiler\include\ISO_Fortran_binding.h(156): note: This member will be ignored by a defaulted constructor or copy/move assignment operator

C:\temp>ifort /c /standard-semantics /O3 /warn:all /stand:f18 Recaman.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.2.0 Build 20210228_000000
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.

C:\temp>link CalcRecaman.obj Recaman.obj /subsystem:console /out:Recaman.exe
Microsoft (R) Incremental Linker Version 14.28.29337.0
Copyright (C) Microsoft Corporation. All rights reserved.

C:\temp>Recaman.exe
Block 0: Calculate 50 numbers in the sequence.
to check the method.
The first terms of the sequence using the function are: 0 1 3 6 2 7 13 20 12 21 11 22 10 23 9 24 8 25 43 62 42 63 41 18 42 17 43 16 44 15 45 14 46 79 113 78 114 77 39 78 38 79 37 80 36 81 35 82 34 83
The first terms of the sequence using the C++ method are: 0 1 3 6 2 7 13 20 12 21 11 22 10 23 9 24 8 25 43 62 42 63 41 18 42 17 43 16 44 15 45 14 46 79 113 78 114 77 39 78 38 79 37 80 36 81 35 82 34 83
Both method agree? T

Block 1: Calculate 200000 numbers in the sequence using functional style.
Average CPU usage: 0.064 seconds.
a(N) = 361078
Check a? T

Block 2: Calculate 200000 numbers in the sequence using C++ “subroutine”.
Average CPU usage: 0.062 seconds.
a(N) = 361078
Check a? T

C:\temp>

1 Like
10

With respect to @certik’s point about performance benchmarking, how much better can a foreign method export from Haskell used in Fortran via C interoperability perform compared to options posted here?

12

Indeed, there are several criteria to consider:

  • correctly implements a function (close to mathematics)
  • how easy it is to write the code
  • how easy it is to read and understand the code
  • performance
  • size of the code

This particular example might not be wroth spending the effort to optimize. It’s just my habit to write code that is in general performing, even if we don’t need the performance in the particular case.

Update: I added and re-ordered some of the criteria based on comments below.

15

“correctly implements a function” is a given, the first thing the Fortran-based code examples here strive to show, one just has to look upthread.

Should this site be for Fortran language, what’s of interest far more than anything else is relevance to Fortran and how to apply the knowledge from the discourse usefully and efficiently and preferably immediately in Fortran. None of the criterion matter when there is little to no connection on offer with Fortran; there will always be other infinitely better sites on how to do something in some other programming language or paradigm.

1 Like
16

Thanks @FortranFan and @Arjen for writing Fortran versions. I think it’s very good to compare to other languages, thanks @pmk for writing the Haskell version.