Riordan numbers
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 show up in several places in set theory. They are closely related to Motzkin numbers, and may be used to derive them.
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.
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
C++
#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 n) {
std::string str = num.get_str();
size_t len = str.size();
if (len > n)
str = str.substr(0, n / 2) + "..." + str.substr(len - n / 2);
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).
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
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
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
... 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
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)
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.
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)
Wren
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)