# Integer sequences

**URL:** <https://fortran-lang.discourse.group/t/integer-sequences/1081>\
**Category:** Uncategorized\
**Created:** [April 18, 2021, 3:38pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081 "2021-04-18T15:38:53Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [April 18, 2021, 3:38pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/1 "2021-04-18T15:38:53Z")

</div>

A Hacker News [post](https://news.ycombinator.com/item?id=26848517) 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](https://oeis.org/search?q=2%2C3%2C5%2C7%2C11%2C13%2C17%2C19%2C23&sort=&language=english&go=Search) matches 228 results in the [On-line Encyclopedia of Integer Sequences](https://oeis.org/), 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.

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [April 18, 2021, 8:18pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/2 "2021-04-18T20:18:52Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\
**Post date:** [April 19, 2021, 1:19pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/3 "2021-04-19T13:19:19Z")

</div>

I often visit the site [mathworld.wolfram.com](http://mathworld.wolfram.com) for mathematical gems - not always to find something of immediate use 🙂. 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](https://wiki.tcl-lang.org/page/The+Recaman%27s+sequence), has a rather nice visualisation.)

---

<div class="post-metadata">

**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\
**Post date:** [April 19, 2021, 3:51pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/4 "2021-04-19T15:51:21Z")

</div>

> [@Arjen](#):
>
> …One sequence I encountered recently, not on this site but here: [The Recaman’s sequence](https://wiki.tcl-lang.org/page/The+Recaman%27s+sequence) …

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**](http://rosettacode.org/wiki/Recaman%27s_sequence).where the C++ example makes one really wonder as to where is the **for** mula **trans** lation for:  
The Recamán’s sequence ![{isplaystyle a_{0},a_{1},a_{2}ots }](https://wikimedia.org/api/rest_v1/media/math/render/svg/f98223b2560de036588b919037549a149aeb7eda) 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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/89ab504e66d289aeb91cdd3c9b3d88212a13f506)

```Fortran
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\>

---

<div class="post-metadata">

**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\
**Post date:** [April 20, 2021, 7:47am UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/5 "2021-04-20T07:47:16Z")

</div>

> [@FortranFan](#):
>
> `https://en.wikipedia.org/wiki/Recam%C3%A1n%27s_sequence`

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

```

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [April 21, 2021, 5:11pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/6 "2021-04-21T17:11:53Z")

</div>

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)`.

---

<div class="post-metadata">

**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\
**Post date:** [April 22, 2021, 6:55am UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/7 "2021-04-22T06:55:48Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\
**Post date:** [April 22, 2021, 4:37pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/8 "2021-04-22T16:37:49Z")

</div>

> [@certik](#):
>
> 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)` .

@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:

```C++
#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:

```Fortran
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\>

---

<div class="post-metadata">

**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\
**Post date:** [April 22, 2021, 7:05pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/10 "2021-04-22T19:05:40Z")

</div>

> [@pmk](#):
>
> I think Haskell is a better tool for these jobs

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?

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [April 22, 2021, 8:31pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/12 "2021-04-22T20:31:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\
**Post date:** [April 22, 2021, 8:57pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/15 "2021-04-22T20:57:56Z")

</div>

> [@certik](#):
>
> … there are several criteria to consider …

“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.

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [April 22, 2021, 9:14pm UTC](https://fortran-lang.discourse.group/t/integer-sequences/1081/16 "2021-04-22T21:14:18Z")

</div>

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.
