Riordan numbers
Riordan numbers show up in several places in set theory. They are closely related to Motzkin numbers, and may be used to derive them.
Riordan numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
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
ALGOL 68
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. <lang algol68>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</lang>
- 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
F#
<lang fsharp> // 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) </lang>
- 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
J
Sequence extender:<lang J>riordanext=: (, (<: % >:)@# * 3 2 +/ .* _2&{.)</lang> Task example:<lang J> 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</lang>Stretch:<lang J> #":(1e3-1){riordanext^:(1e3) x:1 0 472
#":(1e4-1){riordanext^:(1e4) x:1 0
4765</lang>
PL/I
<lang pli>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;</lang>
- 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
... 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.
<lang pli>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 B FROM A DISCARDING */
/* THE RESULT ABD RETURNING THE CARRY FLAG */
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</lang>
- 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
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)
Raku
<lang perl6>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</lang>
- 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)
Wren
<lang ecmascript>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)
}</lang>
- 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)