This site is supported by donations to The OEIS Foundation.
Source code for Catalan ranking and unranking functions
From OeisWiki
This is the latest
This page gives implementations for Catalan ranking and unranking functions in various programming languages. For each language, at least ranking and unranking functions (both "local" and "global") between and totally balanced binary sequences (i.e. A014486) will be presented. Also, provided that the language has S-expressions (or any similar data structure), then also versions converting directly there, and optionally, in other cases (C, Haskell?), versions converting directly to/from user-defined data-structures corresponding to plane binary trees.
The algorithms have been adapted from Frank Ruskey's thesis and in some cases from D. L. Kreher's and D. R. Stinson's book, Combinatorial Algorithms, Generation, Enumeration and Search (CAGES), CRC Press (1999), Section 3.4 "Catalan families", pp. 95--104, which gives essentially the same algorithm as the former source for totally balanced binary sequences, although requiring a few more "mirroring steps" to conform to the order adopted in OEIS (A014486).
Contents
Implementation in C
The following functions are used in [1]
/* Not complete yet! */
#ifdef REQUIRES_LONG_LONG
typedef unsigned long long int ULLI; /* Two longs on 64-bit Suns */
typedef long long int SLLI; /* Ditto. */
#else /* We are on some old, small 32-bit machine? */
typedef unsigned long int ULLI; /* Only one long. */
typedef long int SLLI; /* Ditto. */
#endif
typedef ULLI TBBS; /* TBBS = totally balanced binary sequence. With 32 bits only upto size n=15, with 64 upto n=31. */
typedef SLLI RANK; /* With 32 bits we can go upto size n=19, with 64 much further. */
typedef int SIZE;
RANK CatalanRank(SIZE n,TBBS a)
{
int m = -1;
int y = 0;
RANK r = 0;
while(a > 0)
{
if(0 == (a & 1))
{
m++;
}
else
{
y++;
r += CatTriangle(m,y);
}
a >>= 1;
}
return(r);
}
TBBS CatalanUnrank(SIZE n,RANK r)
{
int m = n-1;
int y = n;
TBBS a = 0;
while(m >= 0)
{
RANK c = CatTriangle(m,y);
a <<= 1;
if(r >= c)
{
y--;
a++;
r -= c;
}
else
{
m--;
}
}
return(a);
}
/* Note: both rows and columns start from -1 */
/* Entry at CatTriangle(r,m) = CatTriangle(r,m-1) + CatTriangle(r-1,m) */
#define CatTriangle(r,c) (tA009766[(r+2)][(c+1)])
#ifdef ONLY32BITS
#define TR 21
#else
#define TR 34
#endif
RANK tA009766[][TR+1] = /* 34 rows in full version. */
{
{0},
{0, 0},
{0, 1, 0},
{0, 1, 1, 0},
{0, 1, 2, 2, 0},
{0, 1, 3, 5, 5, 0},
{0, 1, 4, 9, 14, 14, 0},
{0, 1, 5, 14, 28, 42, 42, 0},
{0, 1, 6, 20, 48, 90, 132, 132, 0},
{0, 1, 7, 27, 75, 165, 297, 429, 429, 0},
{0, 1, 8, 35, 110, 275, 572, 1001, 1430, 1430, 0},
{0, 1, 9, 44, 154, 429, 1001, 2002, 3432, 4862, 4862, 0},
{0, 1, 10, 54, 208, 637, 1638, 3640, 7072, 11934, 16796, 16796, 0},
{0, 1, 11, 65, 273, 910, 2548, 6188, 13260, 25194, 41990, 58786, 58786, 0},
{0, 1, 12, 77, 350, 1260, 3808, 9996, 23256, 48450, 90440, 149226, 208012,
208012, 0},
{0, 1, 13, 90, 440, 1700, 5508, 15504, 38760, 87210, 177650, 326876, 534888,
742900, 742900, 0},
{0, 1, 14, 104, 544, 2244, 7752, 23256, 62016, 149226, 326876, 653752,
1188640, 1931540, 2674440, 2674440, 0},
{0, 1, 15, 119, 663, 2907, 10659, 33915, 95931, 245157, 572033, 1225785,
2414425, 4345965, 7020405, 9694845, 9694845, 0},
{0, 1, 16, 135, 798, 3705, 14364, 48279, 144210, 389367, 961400, 2187185,
4601610, 8947575, 15967980, 25662825, 35357670, 35357670, 0},
{0, 1, 17, 152, 950, 4655, 19019, 67298, 211508, 600875, 1562275, 3749460,
8351070, 17298645, 33266625, 58929450, 94287120, 129644790, 129644790, 0},
{0, 1, 18, 170, 1120, 5775, 24794, 92092, 303600, 904475, 2466750, 6216210,
14567280, 31865925, 65132550, 124062000, 218349120, 347993910, 477638700,
477638700, 0},
{0, 1, 19, 189, 1309, 7084, 31878, 123970, 427570, 1332045, 3798795,
10015005, 24582285, 56448210, 121580760, 245642760, 463991880, 811985790,
1289624490, 1767263190, 1767263190, 0}
#ifndef ONLY32BITS
/* Then follows the first row with a value > 4294967295. */
,{0, 1, 20, 209, 1518, 8602, 40480, 164450, 592020, 1924065, 5722860,
15737865, 40320150, 96768360, 218349120, 463991880, 927983760, 1739969550,
3029594040, 4796857230, 6564120420, 6564120420, 0},
{0, 1, 21, 230, 1748, 10350, 50830, 215280, 807300, 2731365, 8454225,
24192090, 64512240, 161280600, 379629720, 843621600, 1771605360,
3511574910, 6541168950, 11338026180, 17902146600, 24466267020,
24466267020, 0},
{0, 1, 22, 252, 2000, 12350, 63180, 278460, 1085760, 3817125, 12271350,
36463440, 100975680, 262256280, 641886000, 1485507600, 3257112960,
6768687870, 13309856820, 24647883000, 42550029600, 67016296620,
91482563640, 91482563640, 0},
{0, 1, 23, 275, 2275, 14625, 77805, 356265, 1442025, 5259150, 17530500,
53993940, 154969620, 417225900, 1059111900, 2544619500, 5801732460,
12570420330, 25880277150, 50528160150, 93078189750, 160094486370,
251577050010, 343059613650, 343059613650, 0},
{0, 1, 24, 299, 2574, 17199, 95004, 451269, 1893294, 7152444, 24682944,
78676884, 233646504, 650872404, 1709984304, 4254603804, 10056336264,
22626756594, 48507033744, 99035193894, 192113383644, 352207870014,
603784920024, 946844533674, 1289904147324, 1289904147324, 0},
{0, 1, 25, 324, 2898, 20097, 115101, 566370, 2459664, 9612108, 34295052,
112971936, 346618440, 997490844, 2707475148, 6962078952, 17018415216,
39645171810, 88152205554, 187187399448, 379300783092, 731508653106,
1335293573130, 2282138106804, 3572042254128, 4861946401452,
4861946401452, 0},
{0, 1, 26, 350, 3248, 23345, 138446, 704816, 3164480, 12776588, 47071640,
160043576, 506662016, 1504152860, 4211628008, 11173706960, 28192122176,
67837293986, 155989499540, 343176898988, 722477682080, 1453986335186,
2789279908316, 5071418015120, 8643460269248, 13505406670700,
18367353072152, 18367353072152, 0},
{0, 1, 27, 377, 3625, 26970, 165416, 870232, 4034712, 16811300, 63882940,
223926516, 730588532, 2234741392, 6446369400, 17620076360, 45812198536,
113649492522, 269638992062, 612815891050, 1335293573130, 2789279908316,
5578559816632, 10649977831752, 19293438101000, 32798844771700,
51166197843852, 69533550916004, 69533550916004, 0},
{0, 1, 28, 405, 4030, 31000, 196416, 1066648, 5101360, 21912660, 85795600,
309722116, 1040310648, 3275052040, 9721421440, 27341497800, 73153696336,
186803188858, 456442180920, 1069258071970, 2404551645100, 5193831553416,
10772391370048, 21422369201800, 40715807302800, 73514652074500,
124680849918352, 194214400834356, 263747951750360, 263747951750360, 0},
{0, 1, 29, 434, 4464, 35464, 231880, 1298528, 6399888, 28312548, 114108148,
423830264, 1464140912, 4739192952, 14460614392, 41802112192, 114955808528,
301758997386, 758201178306, 1827459250276, 4232010895376, 9425842448792,
20198233818840, 41620603020640, 82336410323440, 155851062397940,
280531912316292, 474746313150648, 738494264901008, 1002242216651368,
1002242216651368, 0},
{0, 1, 30, 464, 4928, 40392, 272272, 1570800, 7970688, 36283236, 150391384,
574221648, 2038362560, 6777555512, 21238169904, 63040282096, 177996090624,
479755088010, 1237956266316, 3065415516592, 7297426411968, 16723268860760,
36921502679600, 78542105700240, 160878516023680, 316729578421620,
597261490737912, 1072007803888560, 1810502068789568, 2812744285440936,
3814986502092304, 3814986502092304, 0},
{0, 1, 31, 495, 5423, 45815, 318087, 1888887, 9859575, 46142811, 196534195,
770755843, 2809118403, 9586673915, 30824843819, 93865125915, 271861216539,
751616304549, 1989572570865, 5054988087457, 12352414499425, 29075683360185,
65997186039785, 144539291740025, 305417807763705, 622147386185325,
1219408876923237, 2291416680811797, 4101918749601365, 6914663035042301,
10729649537134605, 14544636039226909, 14544636039226909, 0},
{0, 1, 32, 527, 5950, 51765, 369852, 2258739, 12118314, 58261125, 254795320,
1025551163, 3834669566, 13421343481, 44246187300, 138111313215,
409972529754, 1161588834303, 3151161405168, 8206149492625, 20558563992050,
49634247352235, 115631433392020, 260170725132045, 565588532895750,
1187735919081075, 2407144796004312, 4698561476816109, 8800480226417474,
15715143261459775, 26444792798594380, 40989428837821289, 55534064877048198,
55534064877048198, 0}
#endif
};
/* We could as well compute it on runtime, of course... */
void CheckTriangle(int upto_n)
{
int r,m;
for(r=0; r <= upto_n; r++)
{
for(m=0; m < r; m++)
{
if((CatTriangle(r,m-1) + CatTriangle(r-1,m)) != CatTriangle(r,m))
{
fprintf(stderr,"(CatTriangle(%d,%d) + CatTriangle(%d,%d)) = ",
r,(m-1),(r-1),m);
fprint_ulli(stderr,(CatTriangle(r,m-1) + CatTriangle(r-1,m)));
fprintf(stderr," differs from CatTriangle(%d,%d) = ",
r,m);
fprint_ulli(stderr,CatTriangle(r,m));
fprintf(stderr,"\n");
exit(1);
}
}
if((r > 0) && (CatTriangle(r,m) != CatTriangle(r,m-1)))
{
fprintf(stderr,"(CatTriangle(%d,%d) = ",r,m);
fprint_ulli(stderr,CatTriangle(r,m));
fprintf(stderr," differs from CatTriangle(%d,%d) = ",r,(m-1));
fprint_ulli(stderr,CatTriangle(r,m-1));
fprintf(stderr,"\n");
exit(1);
}
}
/* fprintf(stderr,"Triangle OK!\n"); */
}
Implementation in C for S-expressions
Here are versions that rank and unrank directly from/to "nihilistic" S-expressions, which must be explicitly defined as a new type in C, with struct keyword. Also, certain auxiliary functions, like (re)cons has to be defined.
struct s_exp
{
struct s_exp *s_car; /* I.e. the left branch of the binary tree. */
struct s_exp *s_cdr; /* I.e. the right branch of the ---- "" ---- */
};
typedef struct s_exp SEXP_cell;
typedef SEXP_cell *SEXP;
#define NULLP(s) ((s) == NULL)
#define PAIR(s) (!NULLP(s))
#define CAR(s) ((s)->s_car)
#define CDR(s) ((s)->s_cdr)
#define SET_CAR(s,x) (((s)->s_car) = x)
#define SET_CDR(s,x) (((s)->s_cdr) = x)
/* This is of course not re-entrant... Code the other way,
if you e.g. need a parallel version.
*/
static int CRS_m, CRS_y;
static RANK CRS_r;
void CatalanRankSexpAux(SEXP node)
{
if(NULLP(node)) { CRS_m--; }
else
{
CRS_r += CatTriangle(CRS_m,CRS_y);
CRS_y--;
CatalanRankSexpAux(CAR(node));
CatalanRankSexpAux(CDR(node));
}
}
RANK CatalanRankSexp(SIZE n,SEXP node)
{
CRS_m = n-1;
CRS_y = n;
CRS_r = 0;
CatalanRankSexpAux(node);
return(CRS_r);
}
SEXP CatalanUnrankSexp(SIZE n,RANK r,SEXP *reused)
{
int m = n-1;
int y = n;
int sp = 0;
int rightson = 0;
SEXP root = NULL;
SEXP sonstack[MAXSIZE+1];
while(m >= 0)
{
RANK c = CatTriangle(m,y);
if(r >= c)
{
SEXP newbranch = recons(NULL,NULL,reused);
if(NULLP(root)) { root = newbranch; }
else
{
if(rightson) { SET_CDR(sonstack[sp],newbranch); }
else { SET_CAR(sonstack[sp],newbranch); sp++; }
}
sonstack[sp] = newbranch;
y--;
r -= c;
rightson = 0; /* The next one is a left son. */
}
else
{
m--;
sp -= rightson;
rightson = 1; /* The next one is a right son. */
}
}
return(root);
}
RANK CatalanRankGlobal(SIZE n,TBBS a)
{
if(0 == n) { return(0); }
else { return(A014137(n-1)+CatalanRank(n,a)); }
}
RANK CatalanRankSexpGlobal(SIZE n,SEXP s)
{
if(0 == n) { return(0); }
else { return(A014137(n-1)+CatalanRankSexp(n,s)); }
}
SEXP cons(SEXP carside,SEXP cdrside)
{
SEXP newcell = ((SEXP) calloc(1,sizeof(SEXP_cell)));
if(NULL == newcell)
{
fprintf(stderr,
"cons: Couldn't allocate a chunk of %u bytes to store one cons cell!\n",
sizeof(SEXP_cell));
exit(1);
}
/* fflush(stdout); fprintf(stderr,"cons called!\n"); fflush(stderr); */
SET_CAR(newcell,carside);
SET_CDR(newcell,cdrside);
return(newcell);
}
/* Taking plunge on car-side, find the first cons-cell
with either its CAR or CDR-side NULL,
and return that cell as a result,
to be reused by recons
*/
SEXP degraft(SEXP *subroot)
{
SEXP z;
if(NULLP(CAR(*subroot)))
{
z = *subroot;
*subroot = CDR(z);
return(z);
}
if(NULLP(CDR(*subroot)))
{
z = *subroot;
*subroot = CAR(z);
return(z);
}
return(degraft(&(CAR(*subroot))));
}
/* Here an all-out (not incremental like here!) deconstruction of
the previous S-expression to a separate free-list
(which just means flattening it!)
might be faster on larger S-expressions:
*/
SEXP recons(SEXP carside,SEXP cdrside,SEXP *reused)
{
if(NULLP(*reused)) { return(cons(carside,cdrside)); }
else
{
SEXP z = degraft(reused);
SET_CAR(z,carside);
SET_CDR(z,cdrside);
return(z);
}
}
/* Recursively free the tree allocated with cons. */
void free_cons_tree(SEXP node)
{
if(NULLP(node)) { return; }
if(!NULLP(CAR(node))) { free_cons_tree(CAR(node)); }
if(!NULLP(CDR(node))) { free_cons_tree(CDR(node)); }
free(node);
}
Implementation in Haskell
Implementation in Maple
# This version for global ranking function is from Peter Luschny, Aug 10 2012: A215406 := proc(n) local m, a, y, t, x, u, v; m := iquo(A070939(n), 2); a := A030101(n); y := 0; t := 1; for x from 0 to 2*m-2 do if irem(a, 2) = 1 then y := y + 1 else u := 2*m - x; v := m-1 - iquo(x+y, 2); t := t + A037012(u, v); y := y - 1 fi; a := iquo(a, 2) od; A014137(m) - t end: # seq(A215406(i), i=0..199);
The rest is extracted mostly from Antti Karttunen's old Maple-code collection.
# However, this version of A014486 was written in Aug 19 2012, and HAS NOT yet been tested with Maple: # Cf. A213704 A014486 := n -> `if`(0 = n,0,CatalanUnrank(A072643(n),(n-A014137(A072643(n)-1)))); # The characteristic function of A014486. A080116 := proc(n) local c,lev; lev := 0; c := n; while(c > 0) do lev := lev + (-1)^c; c := floor(c/2); if(lev < 0) then RETURN(0); fi; od; if(lev > 0) then RETURN(0); else RETURN(1); fi; end; A080300 := n -> `if`((0 = n) or (0 = A080116(n)),0,A014137(((A000523(n)+1)/2)-1)+A080301(n)); A080301 := n -> `if`(0 = A080116(n),-1,CatalanRank((A000523(n)+1)/2,n)); CatalanRank := proc(n,aa) local y,r,lo,a; a := aa; r := 0; y := -1; lo := 0; while (a > 0) do if(0 = (a mod 2)) then r := r+1; lo := lo + A009766(r,y); else y := y+1; fi; a := floor(a/2); od; RETURN((binomial(2*n,n)/(n+1))-(lo+1)); end; CatalanUnrank := proc(n,rr) local t,y,lo,r,a,m; r := (binomial(2*n,n)/(n+1))-(rr+1); a := 0; lo := 0; t := n; y := n-1; while(t > 0) do m := A009766(t,y); if(r < (lo+m)) then y := y-1; a := 2*a + 1; else lo := lo+m; t := t-1; a := 2*a; fi; od; RETURN(a); end; A009766 := proc(r,m) option remember; if(m < 0) then RETURN(0); fi; if(r < 0) then RETURN(0); fi; if(m > r) then RETURN(0); fi; if(0 = m) then RETURN(1); fi; RETURN(A009766(r,m-1) + A009766(r-1,m)); end; # Alternative variant: CatTrianglDirect := (r,m) -> `if`((m < 0),0,((r-m+1)*(r+m)!)/(r! * m! * (r+1))); CatalanRankGlobal := proc(a) local n; n := floor(A070939(a)/2); RETURN(add((binomial(2*j,j)/(j+1)),j=0..(n-1))+CatalanRank(n,a)); end; # Gives A014486 (the old way to produce it as a list up to a certain size): CatalanSequences := proc(upto_n) local n,a,r; a := []; for n from 0 to upto_n do for r from 0 to (binomial(2*n,n)/(n+1))-1 do a := [op(a),CatalanUnrank(n,r)]; od; od; RETURN(a); end;
Some of the old Maple-implementations of certain Catalan bijections use also following functions:
# This returns the first element of the list, and in the fashion of # some very old Lisp-variants (Interlisp), returns an empty list for # an empty list: car:=proc(l) if 0 = nops(l) then ([]) else (op(1,l)): fi: end: # This returns the rest (tail) of the list (everything else except the # first element), and in the fashion of Interlisp, returns an empty list # for an empty list: cdr:=proc(l) if 0 = nops(l) then ([]) else (l[2..nops(l)]): fi: end: # pars2binexp converts the nested list (parenthesization) structure # to the binary codes used in sequence A014486. # Examples: A014486(n) A063171(n) (same in binary) # pars2binexp([]); -> 0 0 # pars2binexp([ [] ]); -> 2 10 # pars2binexp([ [],[] ]); -> 10 1010 # pars2binexp([ [[]] ]); -> 12 1100 # pars2binexp([ [],[],[] ]); -> 42 101010 # pars2binexp([ [],[[]] ]); -> 44 101100 # pars2binexp([ [[]],[] ]); -> 50 110010 # pars2binexp([ [[],[]] ]); -> 52 110100 # pars2binexp([ [[[]]] ]); -> 56 111000 pars2binexp := proc(p) local e,s,w,x; if(0 = nops(p)) then RETURN(0); fi; e := 0; for s in p do x := pars2binexp(s); w := A000523(x); e := e * 2^(w+3) + 2^(w+2) + 2*x; od; RETURN(e); end; # This is its inverse: binexp2pars := proc(n) option remember; `if`((0 = n),[],binexp2parsR(binrev(n))); end; binexp2parsR := n -> [binexp2pars(PeelNextBalSubSeq(n)), op(binexp2pars(RestBalSubSeq(n)))]; # Like NextSubBinTree, but leaves the extra zero off: # Analogous to Lisp CAR NextBalSubSeq := proc(nn) local n,z,c; n := nn; c := 0; z := 0; while(1 = 1) do z := 2*z + (n mod 2); c := c + (-1)^n; n := floor(n/2); if(c >= 0) then RETURN(z); fi; od; end; PeelNextBalSubSeq := proc(nn) local n,z,c; if(0 = nn) then RETURN(0); fi; n := nn; c := 0; z := 0; while(1 = 1) do z := 2*z + (n mod 2); c := c + (-1)^n; n := floor(n/2); if(c >= 0) then RETURN((z - 2^(A000523(z)))/2); fi; od; end; # RestBalSubSeq := proc(nn) local n,z,c; n := nn; c := 0; while(1 = 1) do c := c + (-1)^n; n := floor(n/2); if(c >= 0) then break; fi; od; z := 0; c := -1; while(1 = 1) do z := 2*z + (n mod 2); c := c + (-1)^n; n := floor(n/2); if(c >= 0) then RETURN(z/2); fi; od; end; # Starting from the bit-0 (supposed to be 1), # this gives the totally balanced subsequence # of 1's and 0's (contained by that "root"), # followed by additional zero. # Count c contains the count of surplus 0's, each # 1 will decrement, and each 0 increment it by one. # When it gets to 1 we stop. NextSubBinTree := proc(nn) local n,z,c; n := nn; c := 0; z := 0; while(c < 1) do z := 2*z + (n mod 2); c := c + (-1)^n; n := floor(n/2); od; RETURN(z); end; A000523 := proc(n) local nn,i; # Was: floor_log_2 nn := n; for i from -1 to n do if(0 = nn) then RETURN(i); fi; nn := floor(nn/2); od; end;
Implementation in Mathematica
Implementation in Python/SAGE
The following implementations for the characteristic function of totally balanced binary sequences (A080116) and for the global ranking function (A215406) are from Peter Luschny, Aug 09 2012.
def A080116(n) : lev = 0 while n > 0 : lev += (-1)^n if lev < 0: return 0 n = n//2 return 0 if lev > 0 else 1 def A215406(n) : # CatalanRankGlobal(n) m = A070939(n)//2 a = A030101(n) y = 0; t = 1 for x in (1..2*m-1) : u = 2*m - x; v = m - (x+y+1)/2 mn = binomial(u, v) - binomial(u, v-1) t += mn*(1 - a%2) y -= (-1)^a a = a//2 return A014137(m) - t
Implementation in Scheme
The functions below have been tested to work with MIT/GNU Scheme. The code is from http://www.iki.fi/kartturi/matikka/Nekomorphisms/gatorank.scm.txt, http://www.iki.fi/kartturi/matikka/Schemuli/intfun_a.scm.txt and http://www.iki.fi/kartturi/matikka/Schemuli/intfun_b.scm.txt
For the implementation of definec see memoization.
(definec (A014486 n) (if (zero? n) 0 (A213704bi (A072643 n) (- n (A014137 (- (A072643 n) 1)))) ) ) (define (A213704bi row col) (cond ((zero? row) 0) ;; The top row all zeros. ((>= col (A000108 row)) 0) ;; On other rows, give zeros after the C(n) totally balanced binary seqs. (else (CatalanUnrank row col)) ;; But before that, give the totally balanced binary seqs. ) ) (define (CatalanUnrank size rank) (let loop ((a 0) (m (-1+ size)) ;; The row on A009766 (y size) ;; The position on row m of A009766 (rank rank) (c (A009766tr (-1+ size) size)) ) (if (negative? m) a (if (>= rank c) (loop (1+ (* 2 a)) ;; Up the mountain high m (-1+ y) (- rank c) (A009766tr m (-1+ y)) ) (loop (* 2 a) ;; Down to the valley low (-1+ m) y rank (A009766tr (-1+ m) y) ) ) ) ) ) ;; Inverse function for A014486: (define (A080300 n) (cond ((zero? (A080116 n)) 0) (else (CatalanRankGlobal n)))) (define (CatalanRankGlobal a) (if (zero? a) 0 (let ((w/2 (/ (binwidth a) 2))) (+ (A014137 (-1+ w/2)) (CatalanRank w/2 a) ) ) ) ) (define (CatalanRank w/2 a) (let loop ((a a) ;; The totally balanced binary expansion (r 0) (lo 0) (y -1) ) (if (zero? a) (- (/ (C (* 2 w/2) w/2) (1+ w/2)) (1+ lo) ) (if (zero? (modulo a 2)) (loop ;; Down to the valley (floor->exact (/ a 2)) ;; Was: (fix:lsh a -1) ;; a >>= 1 (1+ r) (+ lo (A009766tr (1+ r) y)) y ) (loop ;; Upto the mountain high. (floor->exact (/ a 2)) ;; Was: (fix:lsh a -1) r lo (1+ y) ) ) ) ) ) (define (A009766tr r m) (if (or (> m r) (< m 0)) 0 ;; Maybe we should raise an error instead?! (/ (* (1+ (- r m)) (! (+ r m))) (* (! r) (! m) (1+ r)) ) ) ) (definec (! n) (if (zero? n) 1 (* n (! (-1+ n))))) ;; A000142 (define A000108 (EIGEN-CONVOLUTION 1 *)) (define A014137 (PARTIALSUMS 0 0 A000108)) (define (A072643 n) (first_pos_with_funs_val_gte A014137 (1+ n))) (define (first_pos_with_funs_val_gte fun n) (let loop ((i 0)) (if (>= (fun i) n) i (loop (1+ i)) ) ) ) (define (A080116 c) ;; Characteristic function of A014486 (let loop ((c c) (lev 0)) (cond ((zero? c) (if (zero? lev) 1 0)) ((< lev 0) 0) (else (loop (floor->exact (/ c 2)) (+ lev (- 1 (* 2 (modulo c 2))))) ) ) ) ) (define (binwidth n) ;; = A029837(n+1) (let loop ((n n) (i 0)) (if (zero? n) i (loop (floor->exact (/ n 2)) (1+ i)) ) ) )
Implementation in Scheme for S-expressions
The following variants rank & unrank directly from/to S-expressions.
;; Rank a symbolless S-expression directly. ;; See Frank Ruskey's thesis at: ;; http://webhome.cs.uvic.ca/~ruskey/Publications/Thesis/Thesis.html ;; especially the page 19: ;; http://webhome.cs.uvic.ca/~ruskey/Publications/Thesis/ThesisPage19.png ;; This one added 14. October 2003. ;; Thanks Frank, this makes it all much leaner! (define (CatalanRankSexpAux size node) (let ((m (-1+ size)) ;; The row on A009766 (y size) ;; The position on row m of A009766 (rank 0) ) (let TreeRank ((node node)) (cond ((not (pair? node)) (set! m (-1+ m))) (else (set! rank (+ rank (A009766tr m y))) (set! y (-1+ y)) (TreeRank (car node)) (TreeRank (cdr node)) ) ) ) rank ) ) (define (CatalanRankSexp s) (let ((size (count-pars s))) (if (zero? size) 0 (+ (A014137 (-1+ size)) (CatalanRankSexpAux size s)) ) ) ) ;; Unrank from a rank (an integer) to a symbolless S-expression directly. ;; See Frank Ruskey's thesis at: ;; http://webhome.cs.uvic.ca/~ruskey/Publications/Thesis/Thesis.html ;; especially the page 24: ;; http://webhome.cs.uvic.ca/~ruskey/Publications/Thesis/ThesisPage24corrected.png (define (CatalanUnrankSexpAux size rank) (let ((sonstack (make-vector size)) (root (list)) ) (let loop ((m (-1+ size)) ;; The row on A009766 (y size) ;; The position on row m of A009766 (rank rank) (c (A009766tr (-1+ size) size)) (rightson? #f) (sp 0) ) (if (negative? m) root (cond ((>= rank c) (let ((newbranch (cons (list) (list)))) (cond ((null? root) (set! root newbranch)) (rightson? (set-cdr! (vector-ref sonstack sp) newbranch)) (else (set-car! (vector-ref sonstack sp) newbranch) (set! sp (1+ sp)) ) ) ;; cond (vector-set! sonstack sp newbranch) (loop m (-1+ y) (- rank c) (A009766tr m (-1+ y)) #f ;; Next time we have a left son. sp ;; sp already incremented above if needed. ) ) ;; let ) (else (loop (-1+ m) y rank (A009766tr (-1+ m) y) #t ;; Next time we have a right son. (- sp (if rightson? 1 0)) ) ) ) ) ) ) ) (define (CatalanUnrankSexp n) (let ((size (A072643 n))) (CatalanUnrankSexpAux size (if (zero? n) 0 (- n (A014137 (-1+ size))))) ) ) (define (count-pars a) (cond ((not (pair? a)) 0) (else (+ 1 (count-pars (car a)) (count-pars (cdr a)))) ) )
