Jacobsthal numbers are an integer sequence related to Fibonacci numbers. Similar to Fibonacci, where each term is the sum of the previous two terms, each term is the sum of the previous, plus twice the one before that. Traditionally the sequence starts with the given terms 0, 1.


   J0 = 0
   J1 = 1
   Jn = Jn-1 + 2 × Jn-2


Terms may be calculated directly using one of several possible formulas:


   Jn = ( 2n - (-1)n ) / 3



Jacobsthal-Lucas numbers are very similar. They have the same recurrence relationship, the only difference is an initial starting value J0 = 2 rather than J0 = 0.

Terms may be calculated directly using one of several possible formulas:

   
   JLn = 2n + (-1)n
   

Jacobsthal oblong numbers is the sequence obtained from multiplying each Jacobsthal number Jn by its direct successor Jn+1.


Jacobsthal primes are Jacobsthal numbers that are prime.



Task
 Find and display the first 30 Jacobsthal numbers
 Find and display the first 30 Jacobsthal-Lucas numbers
 Find and display the first 20 Jacobsthal oblong numbers
 Find and display at least the first 10 Jacobsthal primes


See also
 Wikipedia: Jacobsthal number
 Numbers Aplenty - Jacobsthal number
 OEIS:A001045 - Jacobsthal sequence (or Jacobsthal numbers)
 OEIS:A014551 - Jacobsthal-Lucas numbers.
 OEIS:A084175 - Jacobsthal oblong numbers
 OEIS:A049883 - Primes in the Jacobsthal sequence
 Related task: Fibonacci sequence
 Related task: Leonardo numbers




BEGIN # find some Jacobsthal and related Numbers                             #
    INT max jacobsthal = 29;        # highest Jacobsthal number we will find #
    INT max oblong     = 20; # highest Jacobsthal oblong number we will find #
    INT max j prime    = 20;     # number of Jacobsthal prinmes we will find #
    PR precision 200 PR                 # set the precision of LONG LONG INT #
    PR read "primes.incl.a68" PR                   # include prime utilities #
    [ 0 : max jacobsthal ]LONG INT j;         # will hold Jacobsthal numbers #
    [ 0 : max jacobsthal ]LONG INT jl;  # will hold Jacobsthal-Lucas numbers #
    [ 1 : max oblong     ]LONG INT jo; # will hold Jacobsthal oblong numbers #
    # calculate the Jacobsthal Numbers and related numbers                   #
    # Jacobsthal      : J0  = 0, J1  = 1, Jn  = Jn-1  + 2 × Jn-2             #
    # Jacobsthal-Lucas: JL0 = 2, JL1 = 1, JLn = JLn-1 + 2 × JLn-2            #
    # Jacobsthal oblong: JOn = Jn x Jn-1                                     #
    j[ 0 ] := 0; j[ 1 ] := 1; jl[ 0 ] := 2; jl[ 1 ] := 1; jo[ 1 ] := 0;
    FOR n FROM 2 TO UPB j DO
        j[  n ] := j[  n - 1 ] + ( 2 * j[  n - 2 ] );
        jl[ n ] := jl[ n - 1 ] + ( 2 * jl[ n - 2 ] )
    OD;
    FOR n TO UPB jo DO
        jo[ n ] := j[ n ] * j[ n - 1 ]
    OD;
    # prints an array of numbers with the specified legend                   #
    PROC show numbers = ( STRING legend, []LONG INT numbers )VOID:
         BEGIN
            INT n count := 0;
            print( ( "First ", whole( ( UPB numbers - LWB numbers ) + 1, 0 ), " ", legend, newline ) );
            FOR n FROM LWB numbers TO UPB numbers DO
                print( ( " ", whole( numbers[ n ], -11 ) ) );
                IF ( n count +:= 1 ) MOD 5 = 0 THEN print( ( newline ) ) FI
            OD
         END # show numbers # ;
    # show the various numbers numbers                                       #
    show numbers( "Jacobsthal Numbers:",        j  );
    show numbers( "Jacobsthal-Lucas Numbers:",  jl );
    show numbers( "Jacobsthal oblong Numbers:", jo );
    # find some prime Jacobsthal numbers                                     #
    LONG LONG INT  jn1 := j[ 1 ], jn2 := j[ 0 ];
    INT  p count := 0;
    print( ( "First ", whole( max j prime, 0 ), " Jacobstal primes:", newline ) );
    print( ( "   n  Jn", newline ) );
    FOR n FROM 2 WHILE p count < max j prime DO
        LONG LONG INT jn = jn1 + ( 2 * jn2 );
        jn2        := jn1;
        jn1        := jn;
        IF is probably prime( jn ) THEN
            # have a probably prime Jacobsthal number                        #
            p count +:= 1;
            print( ( whole( n, -4 ), ": ", whole( jn, 0 ), newline ) )
        FI
    OD
END