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.
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
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!
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
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
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
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
ReturnYabasic
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 subBracmat
( ( 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++
#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).
Delphi
See #Pascal.
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
.
EDSAC order code
Addition only
Uses the same algorithm as the Pascal solution that calculates Riordan numbers by addition only. With EDSAC's signed 35-bit integers, output is limited to the first 27 Riordan numbers.
[Riordan numbers - Rosetta Code
EDSAC program (Initial Orders 2)
Calculates Riordan numbers by addition only.]
[Arrange the storage]
T45K P56F [H, subroutine to calculate Riordan numbers]
T47K P128F [M, main routine]
T51K P192F [G, subroutine to print 35-bit integer]
[The next 2 are examples of 'preset parameters'.(Wilkes, Wheeler and Gill)]
T46K P27F [N, how many Riordan numbers are required (in address field)]
T55K P256F [V, address of array of Riordan numbers (in address field)]
[Library subroutine M3, prints header at load time and is then overwritten.
This header leaves teleprinter on figures; does not end with CRLF.]
PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF
*RIORDAN!NUMBERS!WITHIN!EDSAC!RANGE#..PZ
[============ H parameter: Subroutine for Riordan numbers ==================
Returns the Riordan numbers R{0}, ..., R{N-1}, where It's assumed that N > 0.
Uses the preset N and V parameters (see above).
The subroutine is equivalent to the following code:
R{0} := 1;
for j := 1 to N - 1 do R{j} := 0;
for k := 2 to N - 1 do
for j := Min( 2*(k - 1), N - 1) downto k do
R{j} := R{j} + R{j-1} + R{j-2});
Note that each array element occupies 2 EDSAC addresses.]
E25K TH GK
A3F T57@ [plant return link as usual]
YF YF [add rounding twice, acc := 35-bit 1]
[4] T#V [R{0} := 1]
A58@ LD [acc := 2*N]
A4@ T4F [4F := T order for R{N}]
A4@ [acc := T order for R{0}]
[Loop to clear R{j} for j = 1, ..., N-1]
[10] A61@ [inc address in T order by 2, i.e. inc(j)]
S4F E18@ [test for j = N, jump out if so]
[13] A4F [restore acc after test]
T15@
[15] TD [(planted) clear R{j}]
A15@ E10@ [always loop bsck for nrxt j]
[18]
[End of loop to initialize array R. Here with acc = 0.
Location 15 holds T order for R{N-1}. This is used in the code below.]
A13@ T25@ [initialize 2-state order (below) to A4F]
A4@ U4F [4F := A order for R{0}]
A61@ [acc := A order for R{1}]
[Outer loop, equivalent to 'for k := 2 to N - 1']
[23] A61@ T5F [5F := T order for R{k}, 2 <= k <= N-1]
[When 2*k - 2 becomes greater than or equal to N - 1, we shall have
Min(2*k - 2, N - 1) = N - 1 for all further k. When that happens,
the following order is changed to jump over further comparisons.]
[25] XF [planted: initially A4F, later jump to location 33.]
A62@ U4F [4F := T order for R{2*k-2}]
S15@ [compare wuth T order for R{N-1}]
G33@ [jump if 2*k - 2 < N - 1]
[2*k - 2 >= N - 1, time to change order at location 25]
TF A59@ T25@ [change 2-state order]
[33] A15@ [acc := T order for R{N - 1}]
[Inner loop, equivalent to 'for j := Min( 2*(k - 1), N - 1) downto k']
[34] U44@ [plant T order for R{j}]
A60@ U41@ [plant A order for R{j}]
S61@ U42@ [plant A order for R{j-1}]
S61@ T43@ [plant A order for R{j-2}]
[41] AD [(planted) acc := R{j}]
[42] AD [(planted) add R{j-1}]
[43] AD [(planted) add R{j-2}]
[44] TD [(planted) update R{j}]
A44@ S61@ [dec(j) in the preceding T order]
S5F G51@ [if j < k then done inner loop, jump to outer loop]
A5F E34@ [restore acc after test and loop back]
[Here at end of inner loop]
[51] TF [clear acc]
A5F S15@ E57@ [jump to return link if k = N - 1]
A15@ E23@ [else restore acc after test and loop back]
[57] ZF [(planted) jump back to caller]
[58] PN [how many Riordan numbers are required (i9n address field)]
[59] E33@ [jump order for 2-state order above]
[60] MF [add to T order to make A order, same address]
[61] P2F [to change address by 2]
[62] P4F [to change address by 4]
[====================== M parameter: Main routine ======================]
E25K TM GK
[0] PF [negative count of Riordan numbers]
[1] P4F [how many per line in printout]
[2] PF [positive count of items per line]
[3] A#V [A order for first Riordan number]
[4] PN [how many Riordan numbers to print]
[5] P2F [to inc array index by 1 (inc EDSAC address by 2)]
[6] @F [carriage return]
[7] &F [line feed]
[8] !F [space]
[9] K4096F [teleprinter null]
[Enter here with acc = 0]
[10] A10@ GH [call subroutine to calculate Riordan numbers]
A1@ T2@ [number of items in line := max]
A3@ T30@ [initialize A order to start of array]
S4@ [acc := negative count of Riordan numbers]
E41@ [done if count = 0]
[18] T@ [update negative vount of Riordan numbers]
A2@ S1@ [is number of items in line = max?]
E24@ [if so, jump to print CR LF]
A1@ E26@ [else restore acc after test and join common code]
[24] O6@ O7@ [here acc = 0; print CR LF]
[26] A2F T2@ [inc number of items in line]
O8@ O8@ [print 2 spaces]
[30] AD TD [(planted) 0D := Riordan number from array]
A32@ GG !10F [print it]
A30@ A5@ T30@ [inc array index]
A@ A2F [inc negative count of Riordan numbeers]
G18@ [if not printed all, loop back]
[41] O9@ [done: print null to flush teleprinter buffer]
ZF [halt the machine]
[============================= G parameter ==================================]
[Modified from library subroutine P7; prints integer N, for 0 <= N < 10^10.]
[Prints 0 correctly, and allows caller to specify a minimum width (up to 10).]
[Input : 0D = integer (not preserved)]
[17-bit print control word follows subroutine call,]
[e.g. !5F means minimum width 5, pad on left with space.]
[49 locations; even address; workspace: 0F, 1F, 4D, 6F, 7F]
E25K TG
GKA48@U4@A27@T35@AFU4@L256FS45@T7FH46#@NDYFLDT4DS45@T1FH21@
S21@A45@G21@U4@TFV4DA1FG36@S1FLDU1FO1FF1FS1FL4FT4DAFG18@ZF
S1FL8FT4DAFA7FG43@O4@T6FE33@T46#ZPFT45ZP1024FP610D@524DP2F
[========================= M parameter again ================================]
E25K TM GK
E10Z [define entry point]
PF [acc = 0 on entry]- Output:
RIORDAN NUMBERS WITHIN EDSAC RANGE
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
Large argument
Solution to the stretch task, using the fact that the digit count of a positive integer m is 1 more than the integer part of log_10(m). To estimate log_10(R(n)) for large n, we use the asymptotic formula in OEIS:A005043 (linked from the task description). By using exact values of R(n) derived from the recurrence in the task description, the asymptotic formula is extended as far as terms in 1/n^4; see the program for details.
[Riordan numbers - Rosetta Code
EDSAC program (Initial Orders 2)
Estimates log_10 of Riordan numbers for large argument (up to 2^34 - 3).]
[From OEIS A005043, with changed notation. Asymptotic formula for
the Riordan number R(n) as n tends to infinity:
R(n) = (3^(n+2)/(sqrt(3*n*pi)*(8*n)))*(1 - 21/(16*n) + O(1/n^2)).
Using log_e(1 - x) = -x + O(x^2) for small x, we get
log_10(R(n)) = (n + 2)*log_10(3) - log_10(8*sqrt(3*pi)) - (3/2)*log_10(n)
- log_10(e)*(21/(16*n) + O(1/n^2))
The logarithm subroutines in the EDSAC library all return logs to base 2.
We use L1, which returns L1(x) = (1/32)*log_2(x), where 2^-32 <= x < 1.
Recalling that EDSAC represents a 35-bit integer n as n/(2^34), we get
log_10(R(n)) = (n + 2)*log_10(3) - A - 48*log_10(2)*L1(n/(2^34))
- log_10(e)*f(n)
where A = log_10(8*sqrt(3*pi)) + 51*log_10(2) = 16.742755...
and f(n) = 21/(16*n) + O(1/n^2). By applying this formula with values of n
for which R(n) is known exactly, we can empirically improve f(n) to
21/(16*n) - 61/(128*n^2) + (273/1024*n^3) - 4891/(16384*n^4) + O(1/n^5).]
[Arrange the storage]
T45K P56F [H: subroutine to calculate Riordan number]
T54K P104F [C: 35-bit constants, even address]
T47K P130F [M: main routine]
T48K P192F [& (Delta): library subroutine D6 for division]
T49K P228F [L: library subroutine L1 for logarithm]
T46K P266F [N: library subroutine P7 to print 35-bit integer]
T50K P302F [X: library subroutine P1 to print fraction]
[Library subroutine R2: reads constants at load time and is then ovewritten.
Since library subroutine R5, for decimal fractions, seems to be unavailable
(lost?), R2 is used to read in values as 35-bit integers (i.e. times 2^34).
WWG says maximum is 10 digits, but 11 is OK provided input is less than 2^34.]
T56K
GK T20F VD L8F A40D UD TF I40F A40F S39F G@ S2F G23F A5@ T5@ E4@ E13Z
T#C [tell R2 where to store constants]
2F8196880740F17003659824F16F12760439398F15514967838F
9792723132F3555691137F1989146886F2227316259F999#
TZ
[0] [35-bit integer 2]
[2] [log_10(3), high 35 bits]
[4] [log_10(3), low 35 bits]
[6] [constant A (above), integer part]
[8] [constant A (above), fractional part]
[10] [3*log_10(2)]
[12] [log_10(e)*(21/16)]
[14] [log_10(e)*(61/128)]
[16] [log_10(e)*(273/1024)]
[18] [log_10(e)*(4891/16384)]
[Library subroutine M3, prints header at load time and is then overwritten.
This header leaves teleprinter on figures.]
PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF
!!!!!!!!!*N!!!!!!!EST#M*!LOG#10K*R#K*N#LL@&#..PZ
[============ H parameter: Subroutine for log_10 of Riordan number ============
Input: 0D = index n as 35-bit integer, 4 <= n <= 2^34 - 3
Output: 4D = integer part of log_10(Riordan(n))
6D = fractional part]
E25K TH GK
A3F T46@ [plant return link as usual]
AD U10D [save n in 10D]
U6D [pass n to logarithm subroutine in 6D]
A#C T12D [save n + 2 in 12D]
A7@ GL [call logarithm subroutine, result in 0D]
AD T14D [save value of L1 for use be;ow]
[acc := f(n), 4th degree polynomial in 1/n, Horner's method]
A#C RD TD [division: 0D := 1]
A10D T4D [4D := n]
A16@ G& [call division routine, 0D := 1/n]
HD [mult reg := 1/n]
V18#C S16#C T10D [re-use 10D as workspace]
V10D A14#C T10D
V10D S12#C T10D
V10D
H12D [mult reg := n + 2]
S8#C [subtract fractional part of constant A]
V4#C [add (n + 2) times low 34 bits of log_10(3)]
[Now we want to (1) subtract 48*log_10(2)*L1, (2) shift 34 right.
But 48*log_10(2) is outside EDSAC range, so to get same result we
(1) shift 4 right, (2) subtract 3*log_10(2)*L1, (3) shift 30 right.]
R4F [shift acc 4 bits right]
H10#C [mult reg := 3*log_10(2)]
N14D [subtract 3*log_10(2)*L1 (where L1 was saved above)]
RF RF [shift 15 + 15 bits right]
S6#C [subtract high part of constant]
H12D [mult reg := n + 2]
V2#C [add (n + 2) times high 34 bits of log_10(3)]
U4D S4D [return integer part, and remove it]
LF LF L64F [shift 13 + 13 + 8 = 34 bits left]
T6D [return fractional part]
[46] ZF [(planted) jump back to caller]
[====================== M parameter: Main routine ======================]
E25K TM GK
T#Z PF [load time: clear sandwich bit in 35-bit constant]
TZ [load time: resume normal loading]
[0] P4D PF [35-bit constant 9]
[2] PF PF
[4] PF PF
[6] P10F [number of arguments to do (in address field)]
[7] PF
[8] MF [dot (in figures mode)]
[9] !F [space]
[10] @F [carriage return]
[11] &F [line feed]
[12] K4086F
[Enter with acc = 0]
[13] T2#@ [argument := 0]
S6@
[15] T7@
A2#@ L1F A2#@ LD [arg*10]
A#@ U2#@
TD [pass arg to print subroutine]
A23@ GN [print arg]
O9@ O9@
A2#@
TD [pass in 0D]
A29@ GH [call Riordan subroutine]
A6D T4#@ [save frac part of result]
A4D TD [pass int part in 0D]
A35@ GN O8@ [print it, followed by dot]
A4#@ TD [pass frac part in 0D]
A40@ GX P10F [print it to 10 decimals]
O10@ O11@
A7@ A2F [inc negative count of arguments]
G15@ [loop until count = 0]
O12@ [print null to flush printer buffer]
ZF [halt the machine]
[============================= L parameter ================================
Library subroutine L1, copied from WWG, 1951, pp 139-140.
38 locations. 0D := (1/32)*log_2(6D).]
E25K TL
GKA3FT33@E11@IFP1024FP512FA3@LDT6DADS4@TDS3@A6DG6@T8FS5@
T4DH6DV6DS3@E34@A3@LDYFT6DA4DADTDA4DRDYFG17@EFA3@YFT6DE29@
[============================= & parameter =============================]
[Library subroutine D6. Sets 0D := 0D/4D, divisor not 0 or -1]
E25K T&
GK A3F T34@ S4D E13@ T4D SD TD E2@ T4D AD LD TD A4D LD E8@ RD U4D LD A35@ T6D E25@ U8D N8D A6D T6D H6D S6D N4D A4D YF G21@ SD VD TD EF W1526D
[============================= N parameter =============================]
[Library subroutine P7; prints integer N in range 0 < N < 10^10.
Number is passed in 0D; output is 10 characters padded left with spaces.]
E25K TN
GKA3FT26@H28#@NDYFLDT4DS27@TFH8@S8@T1FV4DAFG31@SFLDUFOFFFSF
L4FT4DA1FA27@G11@T28#ZPFT27ZP1024FP610D@524D!FO30@SFL8FE22@
[============================= X parameter =============================]
[Library subroutine P1: print non-negative fraction in 0D, without '0.'
Number of decimals follows subroutine call, e.g. P9F for 9 decimals.]
E25K TX GK
GKA18@U17@S20@T5@H19@PFT5@VDUFOFFFSFL4FTDA5@A2FG6@EFU3FJFM1F
[======================= M parameter again =============================]
E25K TM GK
E13Z [define entry point]
PF [acc = 0 on entry]- Output:
N EST. LOG10(R(N))
9 2.3658259759 [2.365487984891]
99 43.7998316909 [43.799831691182]
999 471.7082318158 [471.708231816522]
9999 4764.2994510429 [4764.299451043293]
99999 47703.7123684837 [47703.712368484562]
999999 477111.3416154477 [477111.341615447888]
9999999 4771201.1340923364
99999999 47712112.5588619480
999999999 477121240.3065581407
9999999999 4771212531.2835200768
[Values in brackets have been added manually as a check. They are
based on exact values of R[n], derived from the recurrence relation.]
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
""" 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
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
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.
Pascal
Recurrence
Uses the recurrence in the task description, with 64-bit unsigned integers. By tweaking the algorithm it's possible to generate 46 Riordan numbers before integer overflow occurs.
program RiordanR;
{
Console program to calculate Riordan numbers, using the recurrence relation
in the Rosetta Code task description.
Command line is e.g. RiordanR 32 for first 32 Riordan numbers.
}
{$IFNDEF FPC} // if not Free Pascal Compiler
{$APPTYPE CONSOLE} // this is needed for Delphi
{$ENDIF}
uses Math, SysUtils;
const N_MAX = 46; // same for Free Pascal and Delphi
type TUint64Array = array of Uint64;
// Store Riordan numbers in the passed-in array.
procedure StoreRiordan( var R : TUint64Array);
var
S, T : Uint64;
h, k : integer;
begin
h := High(R); // srrsy is R[0..h], with h+1 elements
R[0] := 1;
R[1] := 0;
for k := 2 to h do begin
S := 3*R[k - 2] + 2*R[k - 1];
{
To get a few more values before UInt64 overflow, we avoid forming
the product (k-1)*S. Given that (k-1)*S/(k+1) is an integer, we note:
(1) if k is odd then (k-1)/2 and (k+1)/2 are coprime, so (k+1)/2 divides S;
(2) if k is even then k-1 and k+1 are coprime, so k+1 divides S.
}
if Odd(k) then
T := S div ((k+1) div 2)
else
T := 2*(S div (k+1));
R[k] := S - T;
end;
end;
// Main routine.
var
Riordan : TUint64Array;
N, k : integer;
begin
if (not SysUtils.TryStrToInt( ParamStr(1), {out} N))
or (N < 2) or (N > N_MAX) then begin
WriteLn( 'Enter RiordanR N (2 <N <= ', N_MAX,
' for first N Riordan numbers');
exit;
end;
// Call zubroutine to store first N Riordan numbers
SetLength( Riordan, N);
StoreRiordan( Riordan);
// Display Riordan numbers on console
Write( 'First ', N, ' Riordan numbers');
for k := 0 to N - 1 do begin
if k mod 3 = 0 then WriteLn
else Write(' ');
Write( Riordan[k]:24);
end;
WriteLn;
end.
- 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
Addition only
The demo program shows an algorithm to calculate Riordan numbers by addition only. The number of additions to calculate the first N Riordan numbers is about (N^2)/2, so for large N the addition method is less efficient than the recurrence in the task description. It may however be of interest for languages such as EDSAC that don't have built-in division.
Delphi note: The UInt64 type (unsigned 64-bit integer) was introduced in Delphi 7, but wasn't properly implemented until D2007. Before D2007, UInt64 is treated as a signed integer, so that only 46 (instead of 47) Riordan numbers can be generated before UInt64 overflows.
program RiordanA;
{
Console program to calculate Riordan numbers, using addition only.
Command line is e.g. RiordanA 32 for first 32 Riordan numbers.
}
{$IFNDEF FPC} // if not Free Pascal Compiler
{$APPTYPE CONSOLE} // this is needed for Delphi
{$ENDIF}
uses Math, SysUtils;
{$IFDEF FPC}
const N_MAX = 47;
{$ELSE}
const N_MAX = 46; // 47 should be OK for D2007 or later (not tested)
{$ENDIF}
type TUint64Array = array of Uint64;
// Store Riordan numbers in the passed-in array.
procedure StoreRiordan( var R : TUint64Array);
var
h, j, k : integer;
begin
h := High(R); // array is R[0..h], with h+1 elements
R[0] := 1;
for j := 1 to h do R[j] := 0;
for k := 2 to h do
for j := Math.Min( 2*(k - 1), h) downto k do
inc( R[j], R[j-1] + R[j-2]);
end;
// Main routine
var
Riordan : TUint64Array;
N, k : integer;
begin
if (not SysUtils.TryStrToInt( ParamStr(1), {out} N))
or (N < 2) or (N > N_MAX) then begin
WriteLn( 'Enter RiordanA N (2 <N <= ', N_MAX,
' for first N Riordan numbers');
exit;
end;
// Call zubroutine to store first N Riordan numbers
SetLength( Riordan, N);
StoreRiordan( Riordan);
// Display Riordan numbers on console
Write( 'First ', N, ' Riordan numbers');
for k := 0 to N - 1 do begin
if k mod 3 = 0 then WriteLn
else Write(' ');
Write( Riordan[k]:24);
end;
WriteLn;
end.
- Output:
[Same as for recurrence]
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
... 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
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
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
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)