Riordan numbers show up in several places in set theory. They are closely related to Motzkin numbers, and may be used to derive them.

Task
Riordan numbers
You are encouraged to solve this task according to the task description, using any language you may know.

Riordan numbers comprise the sequence a where:

   a(0) = 1, a(1) = 0, for subsequent terms, a(n) = (n-1)*(2*a(n-1) + 3*a(n-2))/(n+1)

There are other generating functions, and you are free to use one most convenient for your language.


Task
  • Find and display the first 32 Riordan numbers.


Stretch
  • Find and display the digit count of the 1,000th Riordan number.
  • Find and display the digit count of the 10,000th Riordan number.


See also



11l

Translation of: Python
F riordan(nn)
   V a = [BigInt(1), 0, 1]
   L(n) 3 .< nn
      a.append((n - 1) * (2 * a[n - 1] + 3 * a[n - 2]) I/ (n + 1))
   R a

V rios = riordan(10'000)

L(i) 32
   print(f:‘{commatize(rios[i]):18}’, end' I (i + 1) % 4 == 0 {"\n"} E ‘’)

print(‘The 1,000th Riordan has ’String(rios[999]).len‘ digits.’)
print(‘The 10,000th Rirdan has ’String(rios[9999]).len‘ digits.’)
Output:
                 1                 0                 1                 1
                 3                 6                15                36
                91               232               603             1,585
             4,213            11,298            30,537            83,097
           227,475           625,992         1,730,787         4,805,595
        13,393,689        37,458,330       105,089,229       295,673,994
       834,086,421     2,358,641,376     6,684,761,125    18,985,057,351
    54,022,715,451   154,000,562,758   439,742,222,071 1,257,643,249,140
The 1,000th Riordan has 472 digits.
The 10,000th Rirdan has 4765 digits.

Action!

Translation of: ALGOL W

Finds the first 13 Riordan numbers as Action! is limited to 16 bit integers (signed and unsiged).

;;; Find some Riordan numbers - limited to the first 13 as the largest integer
;;;                             Action! supports is unsigned 16-bit

;;; sets a to the riordan numbers 0 .. n, a must have n elements
PROC riordan( CARD n CARD ARRAY a )
  CARD i

  IF n >= 0 THEN
    a( 0 ) = 1
    IF n >= 1 THEN
      a( 1 ) = 0
      FOR i = 2 TO n DO
        a( i ) = ( ( i - 1 )
                 * ( ( 2 * a( i - 1 ) )
                   + ( 3 * a( i - 2 ) )
                   )
                 )
               / ( i + 1 )
      OD
    FI
  FI
RETURN

PROC Main()
  CARD  ARRAY r( 13 )
  CARD i

  riordan( 13, r )
  FOR i = 0 TO 12 DO
    Put( '  )
    PrintC( r( i ) )
  OD
RETURN
Output:
 1 0 1 1 3 6 15 36 91 232 603 1585 4213

ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32 and 3.0.3

Uses ALGOL 68G's LONG LONG INT which has programmer-specifiable precision. ALGOL 68G version 3 issues a warning that precision 5000 will impact performance but it still executes this program somewhat faster than version 2 does.

BEGIN # find some Riordan numbers #
    # returns a table of the Riordan numbers 0 .. n #
    OP   RIORDAN = ( INT n )[]LONG INT:
         BEGIN
            [ 0 : n ]LONG INT a;
            IF n >= 0 THEN
                a[ 0 ] := 1;
                IF n >= 1 THEN
                    a[ 1 ] := 0;
                    FOR i FROM 2 TO UPB a DO
                        a[ i ] := ( ( i - 1 )
                                  * ( ( 2 * a[ i - 1 ] )
                                    + ( 3 * a[ i - 2 ] )
                                    )
                                  )
                             OVER ( i + 1 )
                    OD
                FI
            FI;
            a
         END # RIORDAN # ;
    # returns a string representation of n with commas                       #
    PROC commatise = ( STRING unformatted )STRING:
         BEGIN
            STRING result      := "";
            INT    ch count    := 0;
            FOR c FROM UPB unformatted BY -1 TO LWB unformatted DO
                IF  ch count <= 2 THEN
                    ch count +:= 1
                ELSE
                    ch count  := 1;
                    IF unformatted[ c ] = " " THEN " " ELSE "," FI +=: result
                FI;
                unformatted[ c ] +=: result
            OD;
            result
         END; # commatise #
    # returns the length of s                                                #
    OP LENGTH = ( STRING s )INT: ( UPB s - LWB s ) + 1;
    BEGIN # show some Riordann numbers                                       #
        []LONG INT r = RIORDAN 31;
        INT shown := 0;
        FOR i FROM LWB r TO UPB r DO
            print( ( commatise( whole( r[ i ], -15 ) ) ) );
            IF ( shown +:= 1 ) = 4 THEN
                print( ( newline ) );
                shown := 0
            FI
        OD
    END;
    BEGIN # calculate the length of the 1 000th and 10 000th Riordan numbers #
        PR precision 5000 PR # allow up to 5 000 digits for LONG LONG INT    #
        LONG LONG INT r2 := -1, r1 := 1, r := 0;
        print( ( newline ) );
        FOR i FROM 2 TO 9 999 DO
            r2 := r1;
            r1 := r;
            r := ( ( i - 1 )
                 * ( ( 2 * r1 )
                   + ( 3 * r2 )
                   )
                 )
              OVER ( i + 1 );
            IF i = 999 OR i = 9 999 THEN
                STRING rs = whole( r, 0 )[ @ 1 ];
                print( ( "The ", whole( i + 1, -6 ), "th number is: "
                       , rs[ 1 : 20 ], "...", rs[ LENGTH rs - 19 : ]
                       , " with ", whole( LENGTH rs, -5 ), " digits"
                       , newline
                       )
                     )
            FI
        OD
    END
END
Output:
                  1                  0                  1                  1
                  3                  6                 15                 36
                 91                232                603              1,585
              4,213             11,298             30,537             83,097
            227,475            625,992          1,730,787          4,805,595
         13,393,689         37,458,330        105,089,229        295,673,994
        834,086,421      2,358,641,376      6,684,761,125     18,985,057,351
     54,022,715,451    154,000,562,758    439,742,222,071  1,257,643,249,140

The   1000th number is: 51077756867821111314...79942013897484633052 with   472 digits
The  10000th number is: 19927418577260688844...71395322020211157137 with  4765 digits

ALGOL W

Finds the first 22 Riordan numbers as Algol W is limited to signed 32 bit integers.

begin % -- find some Riordan numbers                                          %
    % -- sets a to the Riordan numbers 0 .. n - a must have bounds 0 :: n     %
    procedure riordan ( integer value n; integer array a ( * ) ) ;
        if n >= 0 then begin
            a( 0 ) := 1;
            if n >= 1 then begin
                a( 1 ) := 0;
                for i := 2 until n do begin
                    a( i ) := ( ( i - 1 )
                              * ( ( 2 * a( i - 1 ) )
                                + ( 3 * a( i - 2 ) )
                                )
                              )
                          div ( i + 1 )
                end for_i
            end if_n_ge_1
        end riordan ;
    begin % -- show some Riordann numbers                                     %
        integer array r ( 0 :: 21 );
        integer shown;
        riordan( 21, r );
        shown := 0;
        for i := 0 until 21 do begin
            writeon( i_w := 9, s_w := 0, " ", r( i ) );
            shown := shown + 1;
            if shown = 4 then begin
                write();
                shown := 0
            end if_shown_eq_4
        end for_i
    end;
end.
Output:
         1         0         1         1
         3         6        15        36
        91       232       603      1585
      4213     11298     30537     83097
    227475    625992   1730787   4805595
  13393689  37458330

AppleScript

on riordanNumbers(n)
    set a to {1, 0}    
    repeat with n from 3 to n
        set {an2, an1} to a's items -2 thru -1
        set a's end to (n - 2) * (an1 + an1 + 3 * an2) div n
    end repeat
    
    return a
end riordanNumbers

on intToText(int, separator)
    set groups to {}
    repeat while (int > 999)
        set groups's beginning to ((1000 + (int mod 1000 as integer)) as text)'s text 2 thru 4
        set int to int div 1000
    end repeat
    set groups's beginning to int as integer
    return join(groups, separator)
end intToText

on join(lst, delim)
    set astid to AppleScript's text item delimiters
    set AppleScript's text item delimiters to delim
    set txt to lst as text
    set AppleScript's text item delimiters to astid
    return txt
end join

on task()
    set riordans to riordanNumbers(32)
    set columnWidth to (count intToText(riordans's item 32, ",")) + 2
    set output to {"1st 32 Riordan numbers:"}
    set row to {}
    repeat with i from 1 to 32 by 4
        repeat with j from i to (i + 3)
            set end of row to text -columnWidth thru -1 of ¬
                ("                    " & intToText(riordans's item j, ","))
        end repeat
        set end of output to join(row, "")
        set row to {}
    end repeat
    return join(output, linefeed)
end task

task()
Output:
"1st 32 Riordan numbers:
                  1                  0                  1                  1
                  3                  6                 15                 36
                 91                232                603              1,585
              4,213             11,298             30,537             83,097
            227,475            625,992          1,730,787          4,805,595
         13,393,689         37,458,330        105,089,229        295,673,994
        834,086,421      2,358,641,376      6,684,761,125     18,985,057,351
     54,022,715,451    154,000,562,758    439,742,222,071  1,257,643,249,140"

Arturo

riordan: function [n].memoize[
    if zero? n -> return 1
    if one? n -> return 0

    return ((n-1) * ((2*riordan n-1) + 3*riordan n-2)) / n+1
]

riordans: map 0..31 => riordan

loop split.every: 4 riordans 'x ->
    print map x 's -> pad to :string s 13
Output:
            1             0             1             1 
            3             6            15            36 
           91           232           603          1585 
         4213         11298         30537         83097 
       227475        625992       1730787       4805595 
     13393689      37458330     105089229     295673994 
    834086421    2358641376    6684761125   18985057351 
  54022715451  154000562758  439742222071 1257643249140

BASIC

QBasic

CONST limit = 31

DIM r(0 TO limit)
PRINT "First 32 Riordan numbers:"
CALL Riordan(limit, r())
FOR i = 0 TO limit
    PRINT USING "  #############"; r(i);
    cont = cont + 1
    IF cont MOD 4 = 0 THEN PRINT
NEXT i
END

SUB Riordan (n, a())
    IF n >= 0 THEN
        a(0) = 1
        IF n >= 1 THEN
            a(1) = 0
            FOR i = 2 TO n
                a(i) = ((i - 1) * ((2 * a(i - 1)) + (3 * a(i - 2)))) / (i + 1)
            NEXT i
        END IF
    END IF
END SUB

True BASIC

Translation of: FreeBASIC
SUB riordan (n,a())
    IF n >= 0 THEN
        LET a(0) = 1
        IF n >= 1 THEN
            LET a(1) = 0
            FOR i = 2 TO n
                LET a(i) = ((i-1)*((2*a(i-1))+(3*a(i-2))))/(i+1)
            NEXT i
        END IF
    END IF
END SUB

LET limit = 31
LET cont = 0
DIM r(0)
MAT REDIM r(0 TO limit)

PRINT "First 32 Riordan numbers:"
CALL riordan(limit, r())
FOR i = 0 TO limit
    PRINT  USING "  #############": r(i);
    LET count = count + 1
    IF MOD(count, 4) = 0 THEN PRINT
NEXT i
END
Output:
Same as FreeBASIC entry.

uBasic/4tH

Translation of: FreeBASIC
l = 31
c = 0

Dim @r(l)

Print "First 32 Riordan numbers:"
Proc _Riordan(l)

For i = 0 To l
  Print Using "  ____________#"; @r(i);
  c = c + 1
  If c % 4 = 0 Then Print
Next

End

_Riordan
  Param (1)
  Local (1)
  
  If (a@ < 0) = 0 Then
    @r(0) = 1
    If (a@ < 1) = 0 Then
      @r(1) = 0
      For b@ = 2 To a@
        @r(b@) = ((b@-1) * ((2 * @r(b@-1)) + (3 * @r(b@-2)))) / (b@+1)
      Next
    EndIf
  EndIf
Return

Yabasic

Translation of: FreeBASIC
limit = 31
dim r(limit)
print "First 32 Riordan numbers:"
Riordan(limit, r())
for i = 0 to 23
    print r(i) using("#############");
    cont = cont + 1
    if mod(cont, 4) = 0  print
next i
for i = 24 to limit
	print "   ", str$(r(i));
    cont = cont + 1
    if mod(cont, 4) = 0  print
next i
end
sub Riordan (n, a())
    local i
    if n >= 0 then
        a(0) = 1
        if n >= 1 then
            a(1) = 0
            for i = 2 to n
                a(i) = ((i-1) * ((2 * a(i-1)) + (3 * a(i-2)))) / (i+1)
            next i
        fi
    fi
end sub

Bracmat

( ( a
  =
    .   !arg:0&1
      | !arg:1&0
      | !(!arg$A):>0
      |     (!arg+-1)
          * (2*a$(!arg+-1)+3*a$(!arg+-2))
          * (!arg+1)^-1
        : ?(!arg$A)
  )
& "Create a table A that is big enough to memoize 10000 values. All values are initialized to 0."
& "Table elements are selected by using the syntax !index$tableName"
& tbl$(A,10000)
& -1:?n
&   whl
  ' (!n+1:~>32:?n&out$(a$!n))
& "Theoretically, one could just ask for a(10000), but without first computing a(n) for all n < 10000 there is a risk a risk of stack overflow."
& whl'(!n+1:~>10000:?n&a$!n)
& "Apply string pattern matching @(subject:pattern) the the value of a(n) with the special pattern [?variable to find the total length of the string."
& @(a$999:? [?l1000)
& @(a$9999:? [?l10000)
& out$(str$("a(999) has " !l1000 " digits."))
& out$(str$("a(9999) has " !l10000 " digits."))
)
Output:
1
0
1
1
3
6
15
36
91
232
603
1585
4213
11298
30537
83097
227475
625992
1730787
4805595
13393689
37458330
105089229
295673994
834086421
2358641376
6684761125
18985057351
54022715451
154000562758
439742222071
1257643249140
3602118427251
a(999) has 472 digits.
a(9999) has 4765 digits.

C++

Library: GMP
#include <iomanip>
#include <iostream>

#include <gmpxx.h>

using big_int = mpz_class;

class riordan_number_generator {
public:
    big_int next();

private:
    big_int a0_ = 1;
    big_int a1_ = 0;
    int n_ = 0;
};

big_int riordan_number_generator::next() {
    int n = n_++;
    if (n == 0)
        return a0_;
    if (n == 1)
        return a1_;
    big_int a = (n - 1) * (2 * a1_ + 3 * a0_) / (n + 1);
    a0_ = a1_;
    a1_ = a;
    return a;
}

std::string to_string(const big_int& num, size_t max_digits) {
    std::string str = num.get_str();
    size_t len = str.size();
    if (len > max_digits)
        str.replace(max_digits / 2, len - max_digits, "...");
    return str;
}

int main() {
    riordan_number_generator rng;
    std::cout << "First 32 Riordan numbers:\n";
    int i = 1;
    for (; i <= 32; ++i) {
        std::cout << std::setw(14) << rng.next()
                  << (i % 4 == 0 ? '\n' : ' ');
    }
    for (; i < 1000; ++i)
        rng.next();
    auto num = rng.next();
    ++i;
    std::cout << "\nThe 1000th is " << to_string(num, 40) << " ("
              << num.get_str().size() << " digits).\n";
    for (; i < 10000; ++i)
        rng.next();
    num = rng.next();
    std::cout << "The 10000th is " << to_string(num, 40) << " ("
              << num.get_str().size() << " digits).\n";
}
Output:
First 32 Riordan numbers:
             1              0              1              1
             3              6             15             36
            91            232            603           1585
          4213          11298          30537          83097
        227475         625992        1730787        4805595
      13393689       37458330      105089229      295673994
     834086421     2358641376     6684761125    18985057351
   54022715451   154000562758   439742222071  1257643249140

The 1000th is 51077756867821111314...79942013897484633052 (472 digits).
The 10000th is 19927418577260688844...71395322020211157137 (4765 digits).

EasyLang

cnt = 32
print "First " & cnt & " Riordan numbers:"
app = 1 ; ap = 0
print app ; print ap
for n = 2 to cnt - 1
   a = (n - 1) * (2 * ap + 3 * app) / (n + 1)
   print a
   app = ap
   ap = a
.


F#

// Riordan numbers. Nigel Galloway: August 19th., 2022
let r()=seq{yield 1I; yield 0I; yield! Seq.unfold(fun(n,n1,n2)->let r=(n-1I)*(2I*n1+3I*n2)/(n+1I) in Some(r,(n+1I,r,n1)))(2I,0I,1I)}
let n=r()|>Seq.take 10000|>Array.ofSeq in n|>Array.take 32|>Seq.iter(printf "%A "); printfn "\nr[999] has %d digits\nr[9999] has %d digits" ((string n.[999]).Length) ((string n.[9999]).Length)
Output:
1 0 1 1 3 6 15 36 91 232 603 1585 4213 11298 30537 83097 227475 625992 1730787 4805595 13393689 37458330 105089229 95673994 834086421 2358641376 6684761125 18985057351 54022715451 154000562758 439742222071 1257643249140
r[999] has 472 digits
r[9999] has 4765 digits

FreeBASIC

Const limit = 31

Sub Riordan (n As Integer, a() As Integer)
    If n >= 0 Then
        a(0) = 1
        If n >= 1 Then
            a(1) = 0
            For i As Integer = 2 To n
                a(i) = ((i-1) * ((2 * a(i-1)) + (3 * a(i-2)))) / (i+1)
            Next i
        End If
    End If
End Sub

Dim As Integer r(0 To limit)
Dim As Byte cont = 0
Print "First 32 Riordan numbers:"
Riordan(limit, r())
For i As Integer = 0 To limit
    Print Using "  #############"; r(i);
    cont += 1
    If cont Mod 4 = 0 Then Print
Next i
Sleep
Output:
First 32 Riordan numbers:
              1              0              1              1
              3              6             15             36
             91            232            603           1585
           4213          11298          30537          83097
         227475         625992        1730787        4805595
       13393689       37458330      105089229      295673994
      834086421     2358641376     6684761125    18985057351
    54022715451   154000562758   439742222071  1257643249140


FutureBasic

_limit = 31

local fn Riordan( n as long, a(_limit) as long )
  long i
  if ( n >= 0 )
    a(0) = 1
    if ( n >= 1 )
      a(1) = 0
      for i = 2 to n
        a(i) = ( ( i - 1 ) * ( ( 2 * a(i-1) ) + ( 3 * a(i-2) ) ) ) / ( i + 1 )
      next
    end if
  end if
end fn

long i, count = 0, r(_limit)
printf @"First 32 Riordan numbers:"
fn Riordan( _limit, r(0) )

for i = 0 to 23
  printf @"%16ld\b", r(i)
  count++
  if count mod 4 == 0 then print
next

count = 0
for i = 24 to _limit
  printf @"%16s\b", fn StringUTF8String( Str(r(i)) )
  count++
  if count mod 4 == 0 then print
next

HandleEvents
Output:
First 32 Riordan numbers:
             1               0               1               1
             3               6              15              36
            91             232             603            1585
          4213           11298           30537           83097
        227475          625992         1730787         4805595
      13393689        37458330       105089229       295673994
     834086421      2358641376      6684761125     18985057351
   54022715451    154000562758    439742222071   1257643249140

Haskell

--------------------- RIORDAN NUMBERS --------------------

riordans :: [Integer]
riordans =
  1 :
  0 :
  zipWith
    div
    ( zipWith
        (*)
        [1 ..]
        ( zipWith
            (+)
            ((2 *) <$> tail riordans)
            ((3 *) <$> riordans)
        )
    )
    [3 ..]

-------------------------- TESTS -------------------------
main :: IO ()
main =
  putStrLn "First 32 Riordan terms:"
    >> mapM_ print (take 32 riordans)
    >> mapM_
      ( \x ->
          putStrLn $
            concat
              [ "\nDigit count of ",
                show x,
                "th Riordan term:\n",
                (show . length . show)
                  (riordans !! pred x)
              ]
      )
      [1000, 10000]
Output:
First 32 Riordan terms:
1
0
1
1
3
6
15
36
91
232
603
1585
4213
11298
30537
83097
227475
625992
1730787
4805595
13393689
37458330
105089229
295673994
834086421
2358641376
6684761125
18985057351
54022715451
154000562758
439742222071
1257643249140

Digit count of 1000th Riordan term:
472

Digit count of 10000th Riordan term:
4765

J

Sequence extender:
riordanext=: (, (<: % >:)@# * 3 2 +/ .* _2&{.)
Task example:
   riordanext^:(30) 1x 0
1 0 1 1 3 6 15 36 91 232 603 1585 4213 11298 30537 83097 227475 625992 1730787 4805595 13393689 37458330 105089229 295673994 834086421 2358641376 6684761125 18985057351 54022715451 154000562758 439742222071 1257643249140
Stretch:
   #":(1e3-1){riordanext^:(1e3) x:1 0
472
   #":(1e4-1){riordanext^:(1e4) x:1 0
4765

Java

import java.math.BigInteger;
import java.util.List;

public final class RiordanNumbers {

	public static void main(String[] args) {
		final int limit = 10_000;
		final BigInteger THREE = BigInteger.valueOf(3);
		
		BigInteger[] riordans = new BigInteger[limit];
		riordans[0] = BigInteger.ONE;
		riordans[1] = BigInteger.ZERO;
		for ( int n = 2; n < limit; n++ ) {
			BigInteger term = BigInteger.TWO.multiply(riordans[n - 1]).add(THREE.multiply(riordans[n - 2]));
			riordans[n] = BigInteger.valueOf(n - 1).multiply(term).divide(BigInteger.valueOf(n + 1));
		}
		
		System.out.println("The first 32 Riordan numbers:");
	    for ( int i = 0; i < 32; i++ ) {
	    	System.out.print(String.format("%14d%s", riordans[i], ( i % 4 == 3 ? "\n" :" " )));
	    }
	    System.out.println();
	    
	    for ( int count : List.of( 1_000, 10_000 ) ) {
	    	int length = riordans[count - 1].toString().length();
	    	System.out.println("The " + count + "th Riordan number has " + length + " digits");
	    }
	}		

}
Output:
The first 32 Riordan numbers:
             1              0              1              1
             3              6             15             36
            91            232            603           1585
          4213          11298          30537          83097
        227475         625992        1730787        4805595
      13393689       37458330      105089229      295673994
     834086421     2358641376     6684761125    18985057351
   54022715451   154000562758   439742222071  1257643249140

The 1000th Riordan number has 472 digits
The 10000th Riordan number has 4765 digits

jq

The C implementation of jq has sufficient arithmetic accuracy for the first task, but because of the stretch task, the Go implementation has been used as gojq's integer arithmetic has unbounded accuracy.

Using the program below to calculate the first 100,000 Riordan numbers, gojq takes about 4.7 seconds on a 3GHz machine.

def riordan:
  {ai: 1, a1: 0}
  | ., foreach range(1; infinite) as $i (.;
         {ai: ( ($i-1) * (2*.ai + 3*.a1) / ($i+1)),
          a1: .ai } )
  | .ai ;

def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;

def snip($n):
  tostring|lpad(6)
  + ($n | tostring | "th: \(.[:10]) .. \(.[-10:]) (\(length) digits)" );

"First 32 Riordan numbers:",
foreach limit(100000; riordan) as $riordan (0; .+1;
    if . <= 32 then $riordan
    elif . == 1000  or . == 10000 or . == 100000 then snip($riordan)
    else empty end)

Invocation: jq -nr -f riordan.jq

Output:
First 32 Riordan numbers:
1
0
1
1
3
6
15
36
91
232
603
1585
4213
11298
30537
83097
227475
625992
1730787
4805595
13393689
37458330
105089229
295673994
834086421
2358641376
6684761125
18985057351
54022715451
154000562758
439742222071
1257643249140
  1000th: 5107775686 .. 7484633052 (472 digits)
 10000th: 1992741857 .. 0211157137 (4765 digits)
100000th: 5156659846 .. 4709713332 (47704 digits)

Julia

Translation of: Python
"""  julia example for rosettacode.org/wiki/Riordan_number """

using Formatting

const riordans = zeros(BigInt, 10000)
riordans[begin] = 1

for i in firstindex(riordans)+1:lastindex(riordans)-1
    riordans[i + 1] = (i - 1) * (2 * riordans[i] + 3 * riordans[i - 1]) ÷ (i + 1)
end

for i in 0:31
    print(rpad(format(riordans[begin+i], commas = true), 18), (i + 1) % 4 == 0 ? "\n" : "")
end
println("\nThe 1,000th Riordan has $(length(string(riordans[1000]))) digits.")
println("The 10,000th Riordan has $(length(string(riordans[10_000]))) digits.")
Output:
1                 0                 1                 1
3                 6                 15                36
91                232               603               1,585
4,213             11,298            30,537            83,097
227,475           625,992           1,730,787         4,805,595
13,393,689        37,458,330        105,089,229       295,673,994
834,086,421       2,358,641,376     6,684,761,125     18,985,057,351
54,022,715,451    154,000,562,758   439,742,222,071   1,257,643,249,140

The 1,000th Riordan has 472 digits.
The 10,000th Riordan has 4765 digits.

Lua

Translation of: ALGOL 68 – basic task only, as Lua integers are (usually) limited to 2^53
do -- Riordan numbers

    local function riordan( n ) -- returns a table of the Riordan numbers 0 .. n
        local  a = {}
        if n >= 0 then
            a[ 0 ] = 1
            if n >= 1 then
                a[ 1 ] = 0
                for i = 2, n do
                    a[ i ] = math.floor( ( ( i - 1 )
                                         * ( ( 2 * a[ i - 1 ] )
                                           + ( 3 * a[ i - 2 ] )
                                           )
                                         )
                                       / ( i + 1 )
                                       )
                end
            end
        end
        return a
    end
    local function commatise( unformatted ) -- returns a string representation of n with commas
        local result, chCount = "", 0
        for c = #unformatted, 1, -1 do
            if chCount <= 2 then
                chCount = chCount + 1
            else
                chCount = 1
                result = ( unformatted:sub( c, c ) == " " and " " or "," )..result
            end
            result = unformatted:sub( c, c )..result
        end
        return result
    end

    do -- show the first 32 Riordann numbers
        local r, shown = riordan( 31 ), 0
        for i = 0, #r do
            shown = ( shown + 1 ) % 4
            io.write( commatise( string.format( "%15d", r[ i ] ) )
                    , ( shown == 0 and "\n" or "" )
                    )
        end
    end
end
Output:
                  1                  0                  1                  1
                  3                  6                 15                 36
                 91                232                603              1,585
              4,213             11,298             30,537             83,097
            227,475            625,992          1,730,787          4,805,595
         13,393,689         37,458,330        105,089,229        295,673,994
        834,086,421      2,358,641,376      6,684,761,125     18,985,057,351
     54,022,715,451    154,000,562,758    439,742,222,071  1,257,643,249,140

Mathematica/Wolfram Language

Riordan[N_] := 
 Module[{a = {1, 0, 1}}, 
  Do[AppendTo[a, ((n - 1) (2 a[[n]] + 3 a[[n - 1]])/(n + 1))], {n, 3, 
    N}];
  a]

rios = Riordan[10000];

Do[Print[ToString@NumberForm[rios[[i]], DigitBlock -> 3]], {i, 32}]

Print["The 1,000th Riordan number has ", IntegerLength[rios[[1000]]], 
  " digits."];
Print["The 10,000th Riordan number has ", 
  IntegerLength[rios[[10000]]], " digits."];
Output:
1
0
1
1
3
6
15
36
91
232
603
1,585
4,213
11,298
30,537
83,097
227,475
625,992
1,730,787
4,805,595
13,393,689
37,458,330
105,089,229
295,673,994
834,086,421
2,358,641,376
6,684,761,125
18,985,057,351
54,022,715,451
154,000,562,758
439,742,222,071
1,257,643,249,140
The 1,000th Riordan number has 472 digits.
The 10,000th Riordan number has 4765 digits.

Nim

Task

import std/strformat

iterator riordan(): int =
  var prev = 1
  var curr = 0
  yield prev
  var n = 1
  while true:
    yield curr
    inc n
    let next = (n - 1) * (2 * curr + 3 * prev) div (n + 1)
    prev = curr
    curr = next

echo &"First 32 Riordan numbers:"
var count = 0
for n in riordan():
  inc count
  stdout.write &"{n:<13}"
  stdout.write if count mod 4 == 0: '\n' else: ' '
  if count == 32: break
Output:
First 32 Riordan numbers:
1             0             1             1            
3             6             15            36           
91            232           603           1585         
4213          11298         30537         83097        
227475        625992        1730787       4805595      
13393689      37458330      105089229     295673994    
834086421     2358641376    6684761125    18985057351  
54022715451   154000562758  439742222071  1257643249140

Stretch task

Library: Nim-Integers
import std/strformat
import integers

iterator riordan(): Integer =
  var prev = newInteger(1)
  var curr = newInteger(0)
  yield prev
  var n = 1
  while true:
    yield curr
    inc n
    let next = (n - 1) * (2 * curr + 3 * prev) div (n + 1)
    prev = curr
    curr = next

var count = 0
for n in riordan():
  inc count
  if count in [1000, 10000]:
    echo &"The {count}th Riordan number has {len($n)} digits."
    if count == 10000: break
Output:
The 1000th Riordan number has 472 digits.
The 10000th Riordan number has 4765 digits.

PARI/GP

\\ Increase the stack size if necessary
default(parisize, "32M"); \\ Increase to 32MB, adjust if necessary

Riordan(N) = {
    my(a = vector(N));
    a[1] = 1; a[2] = 0; a[3] = 1;
    for (n = 3, N-1,
        a[n+1] = ((n - 1) * (2 * a[n] + 3 * a[n-1]) \ (n + 1)); \\ Integer division
    );
    return(a);
}

rios = Riordan(10000);

\\ Now print the first 32 elements in the desired format
for (i = 1, 32,{
    print1(rios[i]," ")
});
print("")

\\ Print the number of digits for the 1000th and 10000th Riordan numbers
print("The 1,000th Riordan has ", #digits(rios[1000]), " digits.");
print("The 10,000th Riordan has ", #digits(rios[10000]), " digits.");
Output:
1 0 1 1 3 6 15 36 91 232 603 1585 4213 11298 30537 83097 227475 625992 1730787 4805595 13393689 37458330 105089229 295673994 834086421 2358641376 6684761125 18985057351 54022715451 154000562758 439742222071 1257643249140 
The 1,000th Riordan has 472 digits.
The 10,000th Riordan has 4765 digits.

Perl

use v5.36;
use bigint try => 'GMP';
use experimental <builtin for_list>;
use List::Util 'max';
use List::Lazy 'lazy_list';
use Lingua::EN::Numbers qw(num2en_ordinal);

sub abbr ($d) { my $l = length $d; $l < 41 ? $d : substr($d,0,20) . '..' . substr($d,-20) . " ($l digits)" }
sub comma { reverse ((reverse shift) =~ s/(.{3})/$1,/gr) =~ s/^,//r }
sub table ($c, @V) { my $t = $c * (my $w = 2 + max map { length } @V); ( sprintf( ('%'.$w.'s')x@V, @V) ) =~ s/.{1,$t}\K/\n/gr }

my @riordan;
my $riordan_lazy = lazy_list { state @r = (1,0); state $n = 1; $n++; push @r, ($n-1) * (2*$r[1] + 3*$r[0]) / ($n+1) ; shift @r };
push @riordan, $riordan_lazy->next() for 1..1e4;

say 'First thirty-two Riordan numbers:';
say table 4, map { comma $_ } @riordan[0..31];
say 'The ' . num2en_ordinal($_) . ': ' . abbr $riordan[$_ - 1] for 1e3, 1e4;
Output:
First thirty-two Riordan numbers:
                  1                  0                  1                  1
                  3                  6                 15                 36
                 91                232                603              1,585
              4,213             11,298             30,537             83,097
            227,475            625,992          1,730,787          4,805,595
         13,393,689         37,458,330        105,089,229        295,673,994
        834,086,421      2,358,641,376      6,684,761,125     18,985,057,351
     54,022,715,451    154,000,562,758    439,742,222,071  1,257,643,249,140

The one thousandth: 51077756867821111314..79942013897484633052 (472 digits)
The ten thousandth: 19927418577260688844..71395322020211157137 (4765 digits)

Phix

with javascript_semantics
requires("1.0.2") -- mpz_get_str(comma_fill) was not working [!!!]
include mpfr.e
constant limit = 10000
sequence a = {mpz_init(1),mpz_init(0)}
for n=2 to limit do
    mpz an = mpz_init()
    mpz_mul_si(an,a[n],2)
    mpz_addmul_si(an,a[n-1],3)
    mpz_mul_si(an,an,n-1)
    assert(mpz_fdiv_q_ui(an,an,n+1)=0)
    a &= an
end for
printf(1,"First 32 Riordan numbers:\n%s\n",
 {join_by(apply(true,mpz_get_str,{a[1..32],10,true}),1,4," ",fmt:="%17s")})
for i in {1e3, 1e4} do
    printf(1,"The %6s: %s\n", {ordinal(i), mpz_get_short_str(a[i])})
end for
Output:
First 32 Riordan numbers:
                1                 0                 1                 1
                3                 6                15                36
               91               232               603             1,585
            4,213            11,298            30,537            83,097
          227,475           625,992         1,730,787         4,805,595
       13,393,689        37,458,330       105,089,229       295,673,994
      834,086,421     2,358,641,376     6,684,761,125    18,985,057,351
   54,022,715,451   154,000,562,758   439,742,222,071 1,257,643,249,140

The one thousandth: 51077756867821111314...79942013897484633052 (472 digits)
The ten thousandth: 19927418577260688844...71395322020211157137 (4,765 digits)

PL/I

showRiordan: procedure options( main ); /* find some Riordan numbers         */

   %replace maxRiordan by 32;

   /* sets a to the first n riordan numbers a must have bounds 1 : n         */
   /*      so the first number has index 1, not 0                            */
   riordan: procedure( n, a );
      declare n binary( 15 )fixed;
      declare a ( maxRiordan )decimal( 14 );

      declare ( r2, r1, ri, i ) decimal( 14 );

      if n >= 1 then do;
         r2 = 1;
         a( 1 ) = r2;
         if n >= 2 then do;
            r1 = 0;
            a( 2 ) = r1;
            do i = 2 to n;
               ri = ( ( i - 1 )
                    * ( ( 2 * r1 )
                      + ( 3 * r2 )
                      )
                    )
                  / ( i + 1 );
               a( i + 1 ) = ri;
               r2 = r1;
               r1 = ri;
            end;
         end;
      end;
   end riordan ;

   declare r ( maxRiordan )decimal( 14 ), i binary( 15 )fixed;

   call riordan( maxRiordan, r );
   do i = 1 to maxRiordan;
      put list( ' ', r( i ) );
      if mod( i, 4 ) = 0 then put skip;
   end;

end showRiordan;
Output:
                  1                   0                   1                   1
                  3                   6                  15                  36
                 91                 232                 603                1585
               4213               11298               30537               83097
             227475              625992             1730787             4805595
           13393689            37458330           105089229           295673994
          834086421          2358641376          6684761125         18985057351
        54022715451        154000562758        439742222071       1257643249140

PL/M

Works with: 8080 PL/M Compiler
... under CP/M (or an emulator)

PL/M only handles 8 and 16 bit unsigned integers but also provides two-digit BCD addition/subtraction with carry.
This sample uses the BCD facility to implement 16-digit arithmetic and solve the basic task. Ethiopian multiplication and Egyptian division are used, hence the length of the sample.

100H: /* FIND SOME RIORDAN NUMBERS                                           */

   DECLARE FALSE LITERALLY '0';
   DECLARE TRUE  LITERALLY '0FFH';

   /* CP/M SYSTEM CALL AND I/O ROUTINES                                      */
   BDOS:      PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
   PR$CHAR:   PROCEDURE( C ); DECLARE C BYTE;    CALL BDOS( 2, C );  END;
   PR$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S );  END;
   PR$NL:     PROCEDURE;   CALL PR$CHAR( 0DH ); CALL PR$CHAR( 0AH ); END;
   PR$NUMBER: PROCEDURE( N ); /* PRINTS A NUMBER IN THE MINIMUN FIELD WIDTH  */
      DECLARE N ADDRESS;
      DECLARE V ADDRESS, N$STR ( 6 )BYTE, W BYTE;
      V = N;
      W = LAST( N$STR );
      N$STR( W ) = '$';
      N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
      DO WHILE( ( V := V / 10 ) > 0 );
         N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
      END;
      CALL PR$STRING( .N$STR( W ) );
   END PR$NUMBER;

   DECLARE DEC$LAST LITERALLY '7';           /* SUBSCRIPT OF LAST DIGIT PAIR */
   DECLARE DEC$LEN  LITERALLY '8';        /* LENGTH OF A 16-DIGIT BCD NUMBER */
   DECLARE DEC$16   LITERALLY '( DEC$LEN )BYTE';    /* TYPE DECLARATION OF A */
                                            /* 16-DIGIT BCD NUMBER - 8 BYTES */

   PR$DEC: PROCEDURE( A$PTR );      /* PRINT AN UNSIGNED 16-DIGIT BCD NUMBER */
      DECLARE A$PTR ADDRESS;
      DECLARE A BASED A$PTR DEC$16;
      DECLARE ( D, ZERO$CHAR, I, V ) BYTE;
      ZERO$CHAR = ' ';
      DO I = 0 TO DEC$LAST - 1;
         V = A( I );
         D = SHR( V AND 0F0H, 4 );
         IF D = 0 THEN CALL PR$CHAR( ZERO$CHAR );
                  ELSE CALL PR$CHAR( D + ( ZERO$CHAR := '0' ) );
         D = V AND 0FH;
         IF D = 0 THEN CALL PR$CHAR( ZERO$CHAR );
                  ELSE CALL PR$CHAR( D + ( ZERO$CHAR := '0' ) );
      END;
      V = A( DEC$LAST );
      D = SHR( V AND 0F0H, 4 );
      IF D = 0 THEN CALL PR$CHAR( ZERO$CHAR );
               ELSE CALL PR$CHAR( D + '0' );
      D = V AND 0FH;
      CALL PR$CHAR( D + '0' );
   END PR$DEC ;

   /* SETS THE 16-DIGIT BCD VALUE IN A TO 0                                  */
   INIT$DEC: PROCEDURE( A$PTR );
      DECLARE A$PTR ADDRESS;
      DECLARE A BASED A$PTR ( 0 )BYTE;
      DECLARE I BYTE;
      DO I = 0 TO DEC$LAST;
         A( I ) = 0;
      END;
   END INIT$DEC ;

   /* SETS THE 16-DIGIT BCD VALUE IN A TO B                                  */
   SET$DEC: PROCEDURE( A$PTR, B );
      DECLARE A$PTR ADDRESS, B BYTE;
      DECLARE A BASED A$PTR DEC$16;
      DECLARE ( I, P, V, D1, D2 ) BYTE;
      V = B;
      P = DEC$LAST;
      DO I = 0 TO DEC$LAST;
         IF V = 0
         THEN A( P ) = 0;
         ELSE DO;  
            D1 = V MOD 10;
            D2 = ( V := V / 10 ) MOD 10;
            A( P ) = SHL( D2, 4 ) OR D1;
            V = V / 10;
         END;
         P = P - 1;
      END;
   END SET$DEC ;

   /* ASSIGN THE 16-DIGIT BCD VALUD IN B TO A                                */
   MOV$DEC: PROCEDURE( A$PTR, B$PTR );
      DECLARE ( A$PTR, B$PTR ) ADDRESS;
      DECLARE A BASED A$PTR DEC$16, B BASED B$PTR DEC$16;
      DECLARE I BYTE;
      DO I = 0 TO DEC$LAST;
         A( I ) = B( I );
      END;
   END MOV$DEC ;

   /* BCD ADDITION - ADDS B TO A, STORING THE RESULT IN A                    */
   /*     A AND B MUST HAVE 16 DIGITS                                        */
   ADD$DEC: PROCEDURE( A$PTR, B$PTR );
      DECLARE ( A$PTR, B$PTR ) ADDRESS;
      DECLARE A BASED A$PTR DEC$16, B BASED B$PTR DEC$16;
      DECLARE ( A0, A1, A2, A3, A4, A5, A6, A7 ) BYTE;
      DECLARE ( B0, B1, B2, B3, B4, B5, B6, B7 ) BYTE;
      /* SEPARATE THE DIGIT PAIRS                                            */
      A0 = A( 0 ); A1 = A( 1 ); A2 = A( 2 ); A3 = A( 3 );
      A4 = A( 4 ); A5 = A( 5 ); A6 = A( 6 ); A7 = A( 7 );
      B0 = B( 0 ); B1 = B( 1 ); B2 = B( 2 ); B3 = B( 3 );
      B4 = B( 4 ); B5 = B( 5 ); B6 = B( 6 ); B7 = B( 7 );
      /* DO THE ADDITIONS                                                    */
      A7 = DEC( A7   +  B7 );
      A6 = DEC( A6 PLUS B6 );
      A5 = DEC( A5 PLUS B5 );
      A4 = DEC( A4 PLUS B4 );
      A3 = DEC( A3 PLUS B3 );
      A2 = DEC( A2 PLUS B2 );
      A1 = DEC( A1 PLUS B1 );
      A0 = DEC( A0 PLUS B0 );
      /* RETURN THE RESULT                                                   */
      A( 0 ) = A0; A( 1 ) = A1; A( 2 ) = A2; A( 3 ) = A3;
      A( 4 ) = A4; A( 5 ) = A5; A( 6 ) = A6; A( 7 ) = A7;
   END ADD$DEC;



   /* RETURNS TRUE IF THE 16-DIGIT BCD NUMBER A IS <= B                      */
   /*    USING BCD SUBTRACTION WITH CARRY - SUBTRACTS A FROM B DISCARDING    */
   /*    THE RESULT ABD RETURNING TRUE IF THE CARRY FLAG IS CLEAR            */
   DEC$LE: PROCEDURE( A$PTR, B$PTR )BYTE;
      DECLARE ( A$PTR, B$PTR,C$PTR ) ADDRESS;
      DECLARE A BASED A$PTR DEC$16, B BASED B$PTR DEC$16, C BASED C$PTR DEC$16;
      DECLARE ( A0, A1, A2, A3, A4, A5, A6, A7 ) BYTE;
      DECLARE ( B0, B1, B2, B3, B4, B5, B6, B7 ) BYTE;
      DECLARE ( CFLAG, I ) BYTE;
      /* SEPARATE THE DIGIT PAIRS                                           */
      A0 = A( 0 ); A1 = A( 1 ); A2 = A( 2 ); A3 = A( 3 );
      A4 = A( 4 ); A5 = A( 5 ); A6 = A( 6 ); A7 = A( 7 );
      B0 = B( 0 ); B1 = B( 1 ); B2 = B( 2 ); B3 = B( 3 );
      B4 = B( 4 ); B5 = B( 5 ); B6 = B( 6 ); B7 = B( 7 );
      /* SUBTRACTION A FROM B                                               */
      CFLAG = DEC( B7   -   A7 );
      CFLAG = DEC( B6 MINUS A6 );
      CFLAG = DEC( B5 MINUS A5 );
      CFLAG = DEC( B4 MINUS A4 );
      CFLAG = DEC( B3 MINUS A3 );
      CFLAG = DEC( B2 MINUS A2 );
      CFLAG = DEC( B1 MINUS A1 );
      CFLAG = DEC( B0 MINUS A0 );
      CFLAG = CARRY; /* IF THERE'S NO CARRY, B IS > A AND SO A <= B         */
      RETURN CFLAG = 0; 
   END DEC$LE;

   /* BCD MULTIPLICATION BY AN UNSIGNED INTEGER VIA ETHIOPIAN MULTIPLICATION */
   /*     MULTIPLIES A BY B, STORES THE RESULT IN A - A MUST HAVE 16 DIGITS  */
   MUL$DEC: PROCEDURE( A$PTR, B );
      DECLARE A$PTR ADDRESS, B BYTE;
      DECLARE V BYTE, R DEC$16, ACCUMULATOR DEC$16;
      CALL MOV$DEC( .R, A$PTR );
      V = B;
      CALL INIT$DEC( .ACCUMULATOR );
      DO WHILE( V > 0 );
         IF ( V AND 1 ) = 1 THEN DO;
            CALL ADD$DEC( .ACCUMULATOR, .R );
         END;
         V = SHR( V, 1 );
         CALL ADD$DEC( .R, .R );
      END;
      CALL MOV$DEC( A$PTR, .ACCUMULATOR );
   END MUL$DEC ;

   /* POWERS OF 2 TABLE FOR THE DIVISION ROUTINE                             */
   /* 2^54 IS LARGER THAN A 10^16                                            */
   DECLARE POWERS$OF$2 (  54 /* 16 POINTERS TO 16 DIGIT BCD NUMBERS */
                       )ADDRESS;
   DECLARE POWER$DATA  ( 864 /* 54 16-DIGIT BCD NUMBERS */ )BYTE;
   DO;
      DECLARE ( P, P$POS ) ADDRESS;
      DO P = 0 TO LAST( POWER$DATA ); POWER$DATA( P ) = 0; END;
      POWER$DATA( DEC$LAST ) = 01H;  /* SET LAST DIGIT OF THE 1ST POWER TO 1 */
      P$POS = 0;
      DO P = 0 TO LAST( POWERS$OF$2 );
         POWERS$OF$2( P ) = .POWER$DATA( P$POS );
         P$POS = P$POS + DEC$LEN;
      END;
      DO P = 1 TO LAST( POWERS$OF$2 );
         CALL MOV$DEC( POWERS$OF$2( P ), POWERS$OF$2( P - 1 ) ); /* NEXT... */
         CALL ADD$DEC( POWERS$OF$2( P ), POWERS$OF$2( P     ) ); /* POWER   */
      END;
   END;

   /* BCD DIVISION BY AN UNSIGNED INTEGER VIA EGYPTIAN DIVISION              */
   /*     DIVIDES A BY B, STORES THE RESULT IN A - A MUST HAVE 16 DIGITS     */
   DIV$DEC: PROCEDURE( A$PTR, B );
      DECLARE A$PTR ADDRESS, B BYTE;
      DECLARE DOUBLINGS   (  54 /* 16 POINTERS TO 16 DIGIT BCD NUMBERS */
                          )ADDRESS;
      DECLARE DOUBLE$DATA ( 864 /* 54 16-DIGIT BCD NUMBERS */ )BYTE;
      DECLARE ( D, D$POS ) ADDRESS;
      DECLARE ACCUMULATOR DEC$16, QUOTIENT DEC$16, ACC$PLUS$$DOUBLING DEC$16;
      DECLARE MORE$DOUBLINGS BYTE;
      /* CONSTRUCT THE DOUBLINGS TABLE - A*1, A*2, A*3, ETC.                 */
      CALL SET$DEC( .DOUBLE$DATA, B );
      DOUBLINGS( 0 ) = .DOUBLE$DATA( 0 );
      D$POS          = 0;            /* START OF THE FIRST DOUBLINGS ELEMENT */
      D              = 0;
      MORE$DOUBLINGS = TRUE;
      DO WHILE( MORE$DOUBLINGS );
         D           = D + 1;
         D$POS = D$POS + DEC$LEN;            /* POSITION TO THE NEXT ELEMENT */
         DOUBLINGS( D ) = .DOUBLE$DATA( D$POS );
         CALL MOV$DEC( DOUBLINGS( D ), DOUBLINGS( D - 1 ) );
         CALL ADD$DEC( DOUBLINGS( D ), DOUBLINGS( D     ) );
         MORE$DOUBLINGS = DEC$LE( DOUBLINGS( D ), A$PTR )
                      AND D < LAST( DOUBLINGS );
      END;
      /* CONSTRUCT THE ACCUMULATOR AND QUOTIEMT                              */
      CALL INIT$DEC( .ACCUMULATOR );
      CALL INIT$DEC( .QUOTIENT    );
      D = D + 1;
      DO WHILE( D >= 1 );
         D = D - 1;
         CALL MOV$DEC( .ACC$PLUS$DOUBLING, .ACCUMULATOR );
         CALL ADD$DEC( .ACC$PLUS$DOUBLING, DOUBLINGS( D ) );
         IF DEC$LE( .ACC$PLUS$DOUBLING, A$PTR ) THEN DO;
            CALL MOV$DEC( .ACCUMULATOR, .ACC$PLUS$DOUBLING );
            CALL ADD$DEC( .QUOTIENT,    POWERS$OF$2( D ) );
         END;
      END;
      CALL MOV$DEC( A$PTR, .QUOTIENT );
   END DIV$DEC ;

   /* TASK                                                                   */

   /* SETS A TO THE RIORDAN NUMBERS 0 .. N - LAST(A) MUST BE N               */
   RIORDAN: PROCEDURE( N, A$PTR );
      DECLARE ( N, A$PTR ) ADDRESS;
      DECLARE A BASED A$PTR ( 0 )ADDRESS;
      DECLARE I ADDRESS;
      DECLARE R2 DEC$16, R1 DEC$16;
      DECLARE TWO$R1 DEC$16, THREE$R2 DEC$16;
      CALL INIT$DEC( .R2 );
      CALL INIT$DEC( .R1 );
      IF N >= 0 THEN DO;
         R2( LAST( R2 ) ) = 01H;  /* SET LAST DIGIT OF R2 TO 1, I.E., R2 = 1 */
         CALL MOV$DEC( A( 0 ), .R2 );
         IF N >= 1 THEN DO;
            CALL MOV$DEC( A( 1 ), .R1 );
            DO I = 2 TO N;
               CALL MOV$DEC( .TWO$R1, .R1 );        /* TWO$R1   = R1 ...     */
               CALL ADD$DEC( .TWO$R1, .R1 );        /*          * 2          */
               CALL MOV$DEC( .THREE$R2, .R2 );      /* THREE$R2  = R2 ...    */
               CALL ADD$DEC( .THREE$R2, .R2 );      /* THREE$R2 += R2 ...    */
               CALL ADD$DEC( .THREE$R2, .R2 );      /* THREE$R2 += R2 ...    */
               CALL ADD$DEC( .TWO$R1, .THREE$R2 );  /* TWO$R2 += THREE$R2    */
               CALL MUL$DEC( .TWO$R1, I - 1 );      /* TWO$R2 *= ( I - 1 )   */
               CALL DIV$DEC( .TWO$R1, I + 1 );      /* TWO$R1 /= ( I + 1 )   */
               CALL MOV$DEC( A( I ), .TWO$R1 );     /* A( I )  = TWO$R1      */
               CALL MOV$DEC( .R2, .R1 );            /* R2 = R1               */
               CALL MOV$DEC( .R1, A( I ) );         /* R1 = A( I )           */
            END;
         END;
      END;
   END RIORDAN ;

   /* CONSTRUCT AN ARRAY OF 16 DIGIT BCD NUMBERS                             */
   DECLARE R ( 32 )ADDRESS;                  /* THE ARRAY OF RIORDAN NUMBERS */
   DECLARE R$DATA ( 256 /* 32 * 8 */ )BYTE;   /* THE RIORDAN NUMBER'S DIGITS */
   DECLARE ( I, D$POS ) ADDRESS;
   D$POS = 0;
   DO I = 0 TO LAST( R );
      R( I ) = .R$DATA( D$POS );
      D$POS  = D$POS + DEC$LEN;
   END;
   DO I = 0 TO LAST( R$DATA ); R$DATA( I ) = 0; END;

   /* GET AND PRINT THE RIORDAN NUMBERS                                      */
   CALL RIORDAN( LAST( R ), .R );
   DO I = 0 TO LAST( R );
      CALL PR$CHAR( ' ' );
      CALL PR$DEC( R( I ) );
      IF ( I + 1 ) MOD 4 = 0 THEN CALL PR$NL;
   END;

EOF
Output:
                1                0                1                1
                3                6               15               36
               91              232              603             1585
             4213            11298            30537            83097
           227475           625992          1730787          4805595
         13393689         37458330        105089229        295673994
        834086421       2358641376       6684761125      18985057351
      54022715451     154000562758     439742222071    1257643249140

Python

def Riordan(N):
    a = [1, 0, 1]
    for n in range(3, N):
        a.append((n - 1) * (2 * a[n - 1] + 3 * a[n - 2]) // (n + 1))
    return a

rios = Riordan(10_000)

for i in range(32):
    print(f'{rios[i] : 18,}', end='\n' if (i + 1) % 4 == 0 else '')

print(f'The 1,000th Riordan has {len(str(rios[999]))} digits.')
print(f'The 10,000th Rirdan has {len(str(rios[9999]))} digits.')
Output:
                 1                 0                 1                 1
                 3                 6                15                36
                91               232               603             1,585
             4,213            11,298            30,537            83,097
           227,475           625,992         1,730,787         4,805,595
        13,393,689        37,458,330       105,089,229       295,673,994
       834,086,421     2,358,641,376     6,684,761,125    18,985,057,351
    54,022,715,451   154,000,562,758   439,742,222,071 1,257,643,249,140
The 1,000th Riordan has 472 digits.
The 10,000th Rirdan has 4765 digits.

Quackery

  [ dup -1 peek 2 *
    over -2 peek 3 * +
    over size tuck 1 - *
    swap 1+ /
    join ]               is nextterm ( [ --> [ )


  say "first 32 Riordan numbers: "
  ' [ 1 0 ]
  30 times nextterm
  witheach [ echo sp ]
  cr cr
  say "1000th Riordan number has "
  ' [ 1 0 ]
  1000 2 - times nextterm
  -1 peek number$ size echo
  say " digits"
  cr cr
  say "10000th Riordan number has "
  ' [ 1 0 ]
  10000 2 - times nextterm
  -1 peek number$ size echo
  say " digits"
Output:
first 32 Riordan numbers: 1 0 1 1 3 6 15 36 91 232 603 1585 4213 11298 30537 83097 227475 625992 1730787 4805595 13393689 37458330 105089229 295673994 834086421 2358641376 6684761125 18985057351 54022715451 154000562758 439742222071 1257643249140 

1000th Riordan number has 472 digits

10000th Riordan number has 4765 digits

Raku

use Lingua::EN::Numbers;

my @riordan = 1, 0, { state $n = 1; (++$n - 1) / ($n + 1) × (3 × $^a + 2 × $^b) } … *;

my $upto = 32;
say "First {$upto.&cardinal} Riordan numbers:\n" ~ @riordan[^$upto]».&comma».fmt("%17s").batch(4).join("\n") ~ "\n";

sub abr ($_) { .chars < 41 ?? $_ !! .substr(0,20) ~ '..' ~ .substr(*-20) ~ " ({.chars} digits)" }

say "The {.Int.&ordinal}: " ~ abr @riordan[$_ - 1] for 1e3, 1e4
Output:
First thirty-two Riordan numbers:
                1                 0                 1                 1
                3                 6                15                36
               91               232               603             1,585
            4,213            11,298            30,537            83,097
          227,475           625,992         1,730,787         4,805,595
       13,393,689        37,458,330       105,089,229       295,673,994
      834,086,421     2,358,641,376     6,684,761,125    18,985,057,351
   54,022,715,451   154,000,562,758   439,742,222,071 1,257,643,249,140

The one thousandth: 51077756867821111314..79942013897484633052 (472 digits)
The ten thousandth: 19927418577260688844..71395322020211157137 (4765 digits)

RPL

Works with: HP version 28

Needed to use unsigned integers to avoid the loss of the last digit of a(31).

≪ { #1 #0 } 
   WHILE DUP2 SIZE 1 - > REPEAT
      DUP SIZE 1 -
      DUP2 GETI 3 * 4 ROLLD
      GET 2 * ROT + 
      OVER * SWAP 2 + /
      +
   END SWAP DROP
≫ 'RIORDAN' STO
31 RIORDAN
Output:
1: { #1d #0d #1d #1d #3d #6d #15d #36d #91d #232d #603d #1585d #4213d 11298d #30537d #83097d #227475d #625992d #1730787d #4805595d #13393689d #37458330d #105089229d #295673994d #834086421d #2358641376d #6684761125d #18985057351d #54022715451d #154000562758d #439742222071d #1257643249140d }

Scala

Translation of: Java
import java.math.BigInteger
import scala.collection.mutable.ArrayBuffer

object RiordanNumbers extends App {a
  val limit = 10000
  val THREE = BigInteger.valueOf(3)

  val riordans: ArrayBuffer[BigInteger] = ArrayBuffer.fill(limit)(BigInteger.ZERO)
  riordans(0) = BigInteger.ONE
  riordans(1) = BigInteger.ZERO

  for (n <- 2 until limit) {
    val term = BigInteger.TWO.multiply(riordans(n - 1)).add(THREE.multiply(riordans(n - 2)))
    riordans(n) = BigInteger.valueOf(n - 1).multiply(term).divide(BigInteger.valueOf(n + 1))
  }

  println("The first 32 Riordan numbers:")
  for (i <- 0 until 32) {
    print(f"${riordans(i)}%14d")
    if (i % 4 == 3) println()
    else print(" ")
  }
  println()

  List(1000, 10000).foreach { count =>
    val length = riordans(count - 1).toString.length
    println(s"The ${count}th Riordan number has $length digits")
  }
}
Output:
The first 32 Riordan numbers:
             1              0              1              1
             3              6             15             36
            91            232            603           1585
          4213          11298          30537          83097
        227475         625992        1730787        4805595
      13393689       37458330      105089229      295673994
     834086421     2358641376     6684761125    18985057351
   54022715451   154000562758   439742222071  1257643249140

The 1000th Riordan number has 472 digits
The 10000th Riordan number has 4765 digits


SETL

program riordan;
    a := {[0, 1], [1, 0]};

    loop for n in [2..9999] do
        a(n) := (n-1)*(2*a(n-1) + 3*a(n-2)) div (n+1);
    end loop;

    loop for n in [0..31] do
        putchar(lpad(str a(n), 15));
        if n mod 4=3 then print; end if;
    end loop;

    loop for n in [999, 9999] do
        print("The", str (n+1)+"th Riordan number has", #str a(n), "digits.");
    end loop;
end program;
Output:
              1              0              1              1
              3              6             15             36
             91            232            603           1585
           4213          11298          30537          83097
         227475         625992        1730787        4805595
       13393689       37458330      105089229      295673994
      834086421     2358641376     6684761125    18985057351
    54022715451   154000562758   439742222071  1257643249140
The 1000th Riordan number has 472 digits.
The 10000th Riordan number has 4765 digits.

Sidef

func riordan(n) is cached {
    return 1 if (n == 0)
    return 0 if (n == 1)
    (n-1) * (2*__FUNC__(n-1) + 3*__FUNC__(n-2)) / (n+1)
}

say 32.of(riordan)

for n in (1e3, 1e4) {
    var s = Str(riordan(n-1))
    say "#{'%6s' % n.commify}th term: #{s.first(20)}..#{s.last(20)} (#{s.len} digits)"
}
Output:
[1, 0, 1, 1, 3, 6, 15, 36, 91, 232, 603, 1585, 4213, 11298, 30537, 83097, 227475, 625992, 1730787, 4805595, 13393689, 37458330, 105089229, 295673994, 834086421, 2358641376, 6684761125, 18985057351, 54022715451, 154000562758, 439742222071, 1257643249140]
 1,000th term: 51077756867821111314..79942013897484633052 (472 digits)
10,000th term: 19927418577260688844..71395322020211157137 (4765 digits)

Wren

Library: Wren-gmp
Library: Wren-fmt
import "./gmp" for Mpz
import "./fmt" for Fmt

var limit = 10000
var a = List.filled(limit, null)
a[0] = Mpz.one
a[1] = Mpz.zero
for (n in 2...limit) {
    a[n] = (a[n-1] * 2 + a[n-2] * 3) * (n-1) / (n+1)
}
System.print("First 32 Riordan numbers:")
Fmt.tprint("$,17i", a[0..31], 4)
System.print()
for (i in [1e3, 1e4]) {
   Fmt.print("$,8r: $20a ($,d digits)", i, a[i-1], a[i-1].toString.count)
}
Output:
First 32 Riordan numbers:
                1                 0                 1                 1 
                3                 6                15                36 
               91               232               603             1,585 
            4,213            11,298            30,537            83,097 
          227,475           625,992         1,730,787         4,805,595 
       13,393,689        37,458,330       105,089,229       295,673,994 
      834,086,421     2,358,641,376     6,684,761,125    18,985,057,351 
   54,022,715,451   154,000,562,758   439,742,222,071 1,257,643,249,140 

 1,000th: 51077756867821111314...79942013897484633052 (472 digits)
10,000th: 19927418577260688844...71395322020211157137 (4,765 digits)