/*****************************************************************************/ /* http://www.research.att.com/~njas/sequences/gatonore.c.txt */ /* */ /* gatonore.c -- A lean and fast C-program to count empirically */ /* the amount of distinct non-recursive gatomorphisms. */ /* */ /* See the entries A089831-A089843 in Neil Sloane's Online-Encyclopedia */ /* of Integer Sequences at: http://www.research.att.com/~njas/sequences/ */ /* */ /* Coded 2003 by Antti Karttunen, http://www.iki.fi/~kartturi/ */ /* and placed in Public Domain. */ /* First edited 7. Nov 2003. -- Last edited 5. Dec 2003. */ /* */ /* Mirrored at: */ /* http://www.iki.fi/~kartturi/matikka/Nekomorphisms/gatonore.c */ /* */ /* The idea of program in short: */ /* In principle (*), iterate over all partial_sums(A089836) compositions */ /* of binary tree "transformation clauses", and examine which ones */ /* give bijections which are not encountered before. Count them, */ /* and give the results as table A089831, etc. */ /* */ /* (* in practice we can do some obvious filtering, see the functions */ /* clause_seq_next_composition and clause_seq_next) */ /* */ /* For the essence of the fast implementation, see the functions */ /* apply_clauses_to_x and select_vectors_needed */ /* */ /* If you have any improvements/suggestions, you can mail them to me */ /* to my address Antti.Karttunen (at) iki.fi */ /* */ /* To get correct results, the command line parameters must be large enough. */ /* Both gatonore 127 8 19394489 and gatonore 127 9 15485863 */ /* produce the results that are submitted to OEIS as of 5 Dec 2003. */ /* */ /* See also http://www.research.att.com/~njas/sequences/A089840p.txt */ /* See also http://www.research.att.com/~njas/sequences/gatomorf.c.txt */ /* */ /*****************************************************************************/ #include "stdio.h" #include "stdlib.h" #define MAX_CLAUSES 20 #define MAX_CLAUSE_SIZE 20 /* With our 64 bits MAXSIZE is 31, with 32 bit machines it's only 15. (because we need one extra bit for the last leaf of binary trees) */ #ifdef ONLY32BITS #define MAXSIZE 21 #define MAXSIZE_FOR_TBBS 15 #else #define MAXSIZE 31 #define MAXSIZE_FOR_TBBS 31 #endif typedef unsigned char BYTE; typedef unsigned long int ULI; typedef unsigned short int USI; #define FILLED_BYTE ((BYTE)((~((BYTE) 0)))) #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 running on Alpha, or maybe some small 32-bit machine. */ typedef unsigned long int ULLI; /* Only one long on DEC Alpha's */ 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; typedef USI SGRANK; /* 2 bytes for "Short Global Ranks", i.e. valid up to size n=10. */ typedef USI SLRANK; /* 2 bytes for "Short Local Ranks", i.e. valid up to size n=11. */ typedef ULI PERMRANK; /* 4 bytes for permutation ranks, valid upto size n=11, i.e. 12 leaves. */ typedef ULI NORERANK; /* Ditto for the rank of distinct bijections in A089840. */ typedef BYTE SSIZE; /* One byte for short size. */ typedef ULLI LRANK; /* Local RANK, semantically different from a global rank. */ #define MAX_COMPUTATION_SIZE 10 #define MAX_BINEXP 255 /* I.e. only up to clause sequences of the total size=8. */ /**********************************************************************/ int globopt_output_cycle_lists = 0; char *globopt_list_begin = "("; char *globopt_list_delim = ""; char *globopt_list_end = ")\n"; char *globopt_elem_delim = ","; char *globopt_pair_begin = "("; char *globopt_pair_delim = " "; char *globopt_pair_end = ")"; char *globopt_comment_begin = ";;"; int globopt_HTML = 0; char *glob_author_info = "Antti Karttunen (His_Firstname.His_Surname(AT)iki.fi)"; char *glob_datestr = "Mon DD 20YY"; /* E.g. , Oct 28 2003 */ /**********************************************************************/ /* */ /* Added 14. Oct 2003: SEXP-structure (a simple planar binary tree) */ /* because most of the gatomorphisms are easier to define in such */ /* Lisp-like terms, than to work directly on the binary sequences. */ /* */ /**********************************************************************/ 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) long int glob_cons_cells_in_use = 0; /* For detecting possible leaks. */ 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); */ glob_cons_cells_in_use++; 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((NULL == reused) || 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); glob_cons_cells_in_use--; } /* Recursively free the list allocated with cons, up to the upto_n:th node. */ void free_cons_list(SEXP node,int upto_n) { int i; SEXP next_cdr; for(i=0; i < upto_n; i++) { if(!PAIR(node)) { fprintf(stderr, "free_cons_list: Fatal error: the list to be freed is only %u elems long, shorter than expected %u.\n", i,upto_n); exit(1); } next_cdr = CDR(node); free(node); glob_cons_cells_in_use--; node = next_cdr; } } /**********************************************************************/ #define LONG_ONE ((ULLI) 1) #define two_to(n) (LONG_ONE << (n)) #define A000217(n) (((n)*(n+1))>>1) #define tr2seqind(n,m) (A000217(n)+m) #define tr2seqind11(n,m) tr2seqind(n-1,m-1) /* tr2seqind(0,0) -> 0 tr2seqind(1,0) -> 1 tr2seqind(1,1) -> 2 tr2seqind(2,0) -> 3 tr2seqind(2,1) -> 4 tr2seqind(2,2) -> 5 tr2seqind(3,0) -> 6 etc. */ #define ar2seqind(row,col) (((row+col)*(row+col) + (3*row) + col)>>1) /* ar2seqind(0,0) -> 0 ar2seqind(0,1) -> 1 ar2seqind(1,0) -> 2 ar2seqind(0,2) -> 3 ar2seqind(1,1) -> 4 ar2seqind(2,0) -> 5 ar2seqind(0,3) -> 6 etc. */ int binwidth(int n) { int i=0; while(n != 0) { n >>= 1; i++; } return(i); } PERMRANK sA000142[13] = { 1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800,479001600 }; #define A000142(n) (sA000142[n]) #ifdef ONLY32BITS #define SS 20 #else #define SS 33 #endif RANK sA00108[SS] = {1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786, 208012, 742900, 2674440, 9694845, 35357670, 129644790, 477638700, 1767263190 /* The last value < 2^32 at n=19 */ #ifndef ONLY32BITS ,6564120420, 24466267020, 91482563640, 343059613650, 1289904147324, 4861946401452, 18367353072152, 69533550916004, 263747951750360, 1002242216651368, 3814986502092304, 14544636039226909, 55534064877048198 #endif }; #define A000108(n) (sA00108[n]) #define Cat(n) (sA00108[n]) #define A001246(n) (Cat(n)*Cat(n)) /* The last value < 2^32 at n=19 */ RANK sA014137[SS] = {1, 2, 4, 9, 23, 65, 197, 626, 2056, 6918, 23714, 82500, 290512, 1033412, 3707852, 13402697, 48760367, 178405157, 656043857, ((RANK)2423307047) #ifndef ONLY32BITS ,8987427467, 33453694487, 124936258127, 467995871777, 1757900019101, 6619846420553, 24987199492705, 94520750408709, 358268702159069, 1360510918810437, 5175497420902741, 19720133460129650, 75254198337177848 #endif }; #define A014137(n) (sA014137[n]) /* 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 }; /* See Frank Ruskey's thesis at: http://www.cs.uvic.ca/~fruskey/Publications/Thesis/Thesis.html */ RANK CatalanRank(SIZE n,TBBS a) { int m = -1; int y = 0; LRANK r = 0; while(a > 0) { if(0 == (a & 1)) { m++; } else { y++; r += CatTriangle(m,y); } a >>= 1; } return(r); } TBBS CatalanUnrank(SIZE n,LRANK r) { int m = n-1; int y = n; TBBS a = 0; while(m >= 0) { LRANK c = CatTriangle(m,y); a <<= 1; if(r >= c) { y--; a++; r -= c; } else { m--; } } return(a); } static int CRS_m, CRS_y; static LRANK 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,LRANK 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) { LRANK 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)); } } void print_sexp(SEXP s) { putchar('('); while(!NULLP(s)) { print_sexp(CAR(s)); s = CDR(s); } putchar(')'); } /* See the section "Number Conversion" at the end of the excerpt: http://www.iki.fi/~kartturi/matikka/kl10exmp.txt */ int fprint_ulli(FILE *fp,ULLI x) { int s = 0; if(x >= 10) { s = fprint_ulli(fp,(x/((ULLI)10))); } fputc(('0' + (x%((ULLI)10))),fp); return(s+1); } /* 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"); */ } void CheckRankings(int upto_n) /* Well, superficially... */ { int n; RANK r,r2,uplim; for(n=0; n <= upto_n; n++) { uplim = Cat(n); for(r=0; r < uplim; r++) { TBBS tbs = CatalanUnrank(n,r); if(globopt_output_cycle_lists) { fprint_ulli(stdout,tbs); printf(" "); } r2 = CatalanRank(n,tbs); if(r2 != r) { fflush(stdout); fprintf(stderr,"CatalanRank(%d,",n); fprint_ulli(stderr,tbs); fprintf(stderr,")="); fprint_ulli(stderr,r2); fprintf(stderr," != "); fprint_ulli(stderr,r); fprintf(stderr,"\n"); exit(1); } } if(globopt_output_cycle_lists) { printf("\n"); } } fprintf(stdout,"Ranking & Unranking OK upto n=%d.\n", upto_n); } void CheckSexpRankings(int upto_n) /* Well, superficially... */ { int n; RANK r,r2,uplim; SEXP old_s = NULL; for(n=0; n <= upto_n; n++) { uplim = Cat(n); for(r=0; r < uplim; r++) { SEXP s = CatalanUnrankSexp(n,r,&old_s); if(globopt_output_cycle_lists) { print_sexp(s); printf(" "); } r2 = CatalanRankSexp(n,s); if(r2 != r) { fflush(stdout); fprintf(stderr,"CatalanRankSexp(%d,s)=",n); fprint_ulli(stderr,r2); fprintf(stderr," != "); fprint_ulli(stderr,r); fprintf(stderr,"\n"); exit(1); } old_s = s; } if(globopt_output_cycle_lists) { printf("\n"); } } fprintf(stdout,"Ranking & Unranking OK upto n=%d.\n", upto_n); } /**********************************************************************/ void fprint_USI_vector(FILE *fp,USI *vec,int veclen) { int i; fprintf(fp,"%s",globopt_list_begin); for(i=0; i < veclen; i++) { if(i>0) { fprintf(fp,"%s",globopt_elem_delim); } fprintf(fp,"%u",*(vec+i)); } fprintf(fp,"%s",globopt_list_end); } void fprint_vector_of_USI_vectors(FILE *fp,USI **vec,int veclen) { while(NULL != *vec) { fprint_USI_vector(fp,*vec++,veclen); } } void fprint_ULI_vector(FILE *fp,ULI *vec,int veclen) { int i; fprintf(fp,"%s",globopt_list_begin); for(i=0; i < veclen; i++) { if(i>0) { fprintf(fp,"%s",globopt_elem_delim); } fprintf(fp,"%lu",*(vec+i)); } fprintf(fp,"%s",globopt_list_end); } #define NA_SENTINEL ((ULI) -1) /* I.e. 4294967295 with 32 bits. */ void fprint_ULI_vector_with_NA_SENTINELs(FILE *fp,ULI *vec,int veclen) { int i; fprintf(fp,"%s",globopt_list_begin); for(i=0; i < veclen; i++) { if(i>0) { fprintf(fp,"%s",globopt_elem_delim); } if(NA_SENTINEL == *(vec+i)) { fprintf(fp,"NA"); } else { fprintf(fp,"%lu",*(vec+i)); } } fprintf(fp,"%s",globopt_list_end); } void fprint_ULLI_vector(FILE *fp,ULLI *vec,int veclen) { int i; fprintf(fp,"%s",globopt_list_begin); for(i=0; i < veclen; i++) { if(i>0) { fprintf(fp,"%s",globopt_elem_delim); } fprint_ulli(fp,*(vec+i)); } fprintf(fp,"%s",globopt_list_end); } void fprint_ULLI_vector_up_to_first_0(FILE *fp,ULLI *vec,int veclen) { int i; fprintf(fp,"%s",globopt_list_begin); for(i=0; (i < veclen) && (0 != *(vec+i)); i++) { if(i>0) { fprintf(fp,"%s",globopt_elem_delim); } fprint_ulli(fp,*(vec+i)); } fprintf(fp,"%s",globopt_list_end); } char *w6d(char *tb,int n) { sprintf(tb,"%06u",n); return(tb); } void output_HTML_Alink(FILE *fp,int Anum) { char tb[81] = { 'A' }; char *Astr = (w6d(tb+1,Anum)-1); fprintf(fp,"%s",Astr,Astr); } void output_n_chars(int i, char c) { while(i--) { putchar(c); } } void output_n_spaces(int i) { output_n_chars(i,' '); } /* How about templates? */ int USIvec_find_pos_of_1st_larger_than_one(USI *vec,int veclen) { int pos_of_1st_term_gte_2 = 0; while((pos_of_1st_term_gte_2 < veclen) && (*(vec+pos_of_1st_term_gte_2) < 2)) { pos_of_1st_term_gte_2++; } if(pos_of_1st_term_gte_2 >= veclen) { pos_of_1st_term_gte_2 = 1; } /* Not found, use 1. */ else { pos_of_1st_term_gte_2++; } /* Because One-based. */ return(pos_of_1st_term_gte_2); } int ULIvec_find_pos_of_1st_larger_than_one(ULI *vec,int veclen) { int pos_of_1st_term_gte_2 = 0; while((pos_of_1st_term_gte_2 < veclen) && (*(vec+pos_of_1st_term_gte_2) < 2)) { pos_of_1st_term_gte_2++; } if(pos_of_1st_term_gte_2 >= veclen) { pos_of_1st_term_gte_2 = 1; } /* Not found, use 1. */ else { pos_of_1st_term_gte_2++; } /* Because One-based. */ return(pos_of_1st_term_gte_2); } int ULLIvec_find_pos_of_1st_larger_than_one(ULLI *vec,int veclen) { int pos_of_1st_term_gte_2 = 0; while((pos_of_1st_term_gte_2 < veclen) && (*(vec+pos_of_1st_term_gte_2) < 2)) { pos_of_1st_term_gte_2++; } if(pos_of_1st_term_gte_2 >= veclen) { pos_of_1st_term_gte_2 = 1; } /* Not found, use 1. */ else { pos_of_1st_term_gte_2++; } /* Because One-based. */ return(pos_of_1st_term_gte_2); } #define VECTYPE_USI 16 #define VECTYPE_ULI 32 #define VECTYPE_ULLI 64 #define OEIS_SEQ 0 #define OEIS_TABLE 1 #define VECTYPE_NORERANK VECTYPE_ULI #define fprint_NORERANK_vector_with_NA_SENTINELs fprint_ULI_vector_with_NA_SENTINELs typedef void (*PF_VEC_OUT)(FILE *,void *,int); void output_OEIS_sequence(FILE *fp,int Anum, void *vec,int veclen, int vectype, PF_VEC_OUT fprintvecfun, int offset, int seqtype, char *name ) { char tb[81] = { 'A' }; char *Astr = (w6d(tb+1,Anum)-1); int pos_of_1st_term_gte_2 = ((VECTYPE_ULLI == vectype) ? ULLIvec_find_pos_of_1st_larger_than_one(vec,veclen) : ((VECTYPE_USI == vectype) ? USIvec_find_pos_of_1st_larger_than_one(vec,veclen) : ULIvec_find_pos_of_1st_larger_than_one(vec,veclen) ) ); /* PF_VEC_OUT fprintvecfun = ((VECTYPE_ULLI == vectype) ? fprint_ULLI_vector : ((VECTYPE_USI == vectype) ? fprint_USI_vector : fprint_ULI_vector ) ); */ char *save_globopt_list_begin = globopt_list_begin; char *save_globopt_list_end = globopt_list_end; char *save_globopt_elem_delim = globopt_elem_delim; /* Set these for fprint_vector: */ globopt_list_begin = ""; globopt_list_end = ""; globopt_elem_delim = ","; fprintf(fp,"%%I %s\n",Astr); if(globopt_HTML) { fprintf(fp,"%%S "); output_HTML_Alink(fp,Anum); fprintf(fp," "); fprintvecfun(fp,vec,veclen); fprintf(fp, ""); } else { fprintf(fp,"%%S %s ",Astr); fprintvecfun(fp,vec,veclen); } fprintf(fp,"\n"); fprintf(fp,"%%N %s %s\n",Astr,name); /* Name should end with period. */ fprintf(fp,"%%H %s A. Karttunen, C-program for computing the initial terms of this sequence\n", Astr); fprintf(fp,"%%K %s nonn%s\n",Astr,((OEIS_TABLE == seqtype) ? ",tabl" : "")); fprintf(fp,"%%O %s %u,%u\n",Astr,offset,pos_of_1st_term_gte_2); fprintf(fp,"%%A %s %s, %s\n", Astr, glob_author_info, glob_datestr); /* E.g. Antti Karttunen (Firstname.Surname@iki.fi), Nov 11 2003 */ fprintf(fp,"\n"); fflush(fp); globopt_list_begin = save_globopt_list_begin; globopt_list_end = save_globopt_list_end; save_globopt_elem_delim = save_globopt_elem_delim; } /**********************************************************************/ int sexp_length(SEXP s) { int i = 0; while(PAIR(s)) { i++; s = CDR(s); } return(i); } /* Place the subtrees of t1 extending past the tips of t2 into the SEXP-vector subtrees. The positions in the constructed list correspond to the tips (leaves) of t2 in preorder. This version returns back immediately if A082858(t1,t2) != t2, that is, if t1 is not a super-tree of t2. An example how to write re-entrant code by using doubly indirect pointing, i.e. we do not increment p_to_subtrees itself, but a SEXP pointer to which it points at! */ int sexp_bintree_difference_aux(SEXP t1,SEXP t2,SEXP **p_to_subtrees) { if(!PAIR(t2)) { **p_to_subtrees = t1; ++*p_to_subtrees; return(1); } else if(!PAIR(t1)) { return(0); } /* In this order! */ else { if(0 == sexp_bintree_difference_aux(CAR(t1),CAR(t2),p_to_subtrees)) { return(0); } else { return(sexp_bintree_difference_aux(CDR(t1),CDR(t2),p_to_subtrees)); } } } int sexp_bintree_difference(SEXP t1,SEXP t2,SEXP *subtrees) { SEXP *copy_of_subtree_pointer = subtrees; /* With which we can safely play! */ return(sexp_bintree_difference_aux(t1,t2,©_of_subtree_pointer)); } /* Transform an arbitrarily shaped tree t1 to a "right comb", provided that A082858(t1,t2) = t2. 1 2 2 3 0 \/ 1 \/ \/ 3 --> 0 \/ \/ \/ 7 8 1 2 3 4 6 \/ \/ \/ 5 \/ \ / 4 \/ 0 \/6 7 3 \/ \/5 \/ --> 2 \/ \ \/ 8 1 \/ \ \/ 0 \/ \/ \/ NOTE: It's the responsibility of the caller also to free the main stem (tree_size cons cells) of the tree new_tree possibly constructed here! */ SEXP sexp_transform1(SEXP t1,SEXP t2,const int t2_size) { int i; SEXP tmp; SEXP stem[MAX_CLAUSE_SIZE+1]; SEXP new_tree; if(0 == sexp_bintree_difference(t1,t2,stem)) { return(NULL); } new_tree = stem[t2_size]; for(i=t2_size-1; (i >= 0); i--) { new_tree = cons(stem[i],new_tree); } return(new_tree); } /* Transform a right comb "t1" into an arbitrarily shaped tree new_tree guided by the "guiding tree" t2. Do this by destructively modifying t2, by overwriting the NILs at its tips by the corresponding subtrees of t1. I.e. we will have A082858(new_tree,t2) = t2. This is an inverse transform of sexp_transform1, i.e. we should have sexp_transform2(sexp_transform1(t1,t2),t2) = t1 provided that A082858(t1,t2) = t2. One SHOULD NOT call this if length(t1) < t2_size 2 3 1 2 1 \/ 0 \/ 0 \/ --> \/ 3 \/ \/ 7 8 6 \/ 1 2 3 4 5 \/ \/ \/ 4 \/ \ / 3 \/ 0 \/6 7 2 \/ --> \/5 \/ 1 \/ \ \/ 8 0 \/ \ \/ \/ \/ Think that some people advocate programming languages with no POINTERS! */ void sexp_transform2_aux(SEXP *t1,SEXP *t2,int *n_nodes_left) { if(!PAIR(*t1)) { fprintf(stderr,"sexp_transform2_aux: length(t1) < size(t2), premature end!\n"); exit(1); } if(!PAIR(*t2)) /* Reached one of the tips of t2. */ { if(0 == *n_nodes_left) { *t2 = CDR(*t1); } else { *t2 = CAR(*t1); if(0 != --*n_nodes_left) { *t1 = CDR(*t1); } } } else { sexp_transform2_aux(t1,&(CAR(*t2)),n_nodes_left); sexp_transform2_aux(t1,&(CDR(*t2)),n_nodes_left); } } /* This should protect the variables in the original caller: */ void sexp_transform2(SEXP t1,SEXP t2,int t2_size) /* <= t1_length */ { sexp_transform2_aux(&t1,&t2,&t2_size); } /* We swap the tp1:th and tp2:th subtree (zero-based) sprouting leftward from the main stem, unless tp2 is tree_size, in which case it is the (tp2-1):th cdr from the root. 2 3 0 3 1 \/ 1 \/ 0 \/ transposition (tp1=0, tp2=2) results 2 \/ \/ \/ transposition (tp1=1, tp2=3) results: 2 1 3 \/ 0 \/ \/ 0 1 \/ Note that tp1 < tp2 is required! */ int sexp_transpose(SEXP sexp,const int tree_size,const int tp1,const int tp2) { int i; SEXP tmp; SEXP stem[MAX_CLAUSE_SIZE+1]; for(i=0; (i < tree_size) && PAIR(sexp); ) { stem[i++] = sexp; sexp = CDR(sexp); } /* Can't do it on this tree because its main stem is not long enough. So return zero as a sign of failure. */ if(i < tree_size) { return(0); } /* No-op "transposition", return immediately with success: */ if(tp1 == tp2) { return(1); } if(tree_size == tp2) { tmp = CDR(stem[tp2-1]); SET_CDR(stem[tp2-1],CAR(stem[tp1])); } else /* tp2 < tree_size, swap two leftward leaning subtrees. */ { tmp = CAR(stem[tp2]); SET_CAR(stem[tp2],CAR(stem[tp1])); } SET_CAR(stem[tp1],tmp); return(1); } void precompute_transform1(SGRANK *vec,int upto_n,int t2_size,LRANK t2_lrank) { int t1_size; /* The size of the trees for which we are computing this transformation in the inner loop. */ LRANK t1_lrank; SEXP t1 = NULL; SEXP t2 = CatalanUnrankSexp(t2_size,t2_lrank,NULL); SEXP new_tree = NULL; int i = 0; vec[i++] = 0; /* NIL to NIL, always! */ /* t2 is the "guiding" tree, against which we compare every tree t1 generated in the loop below: */ for(t1_size=1; t1_size <= upto_n; t1_size++) { for(t1_lrank = 0; t1_lrank < Cat(t1_size); t1_lrank++) { if(t1_size < t2_size) /* The subtree does not fit yet. */ { vec[i++] = 0; } else /* Our left-comb of tree_size nodes could be a subtree of */ { /* this sexp, at least in principle. */ t1 = CatalanUnrankSexp(t1_size,t1_lrank,&t1); new_tree = sexp_transform1(t1,t2,t2_size); if(NULLP(new_tree)) { vec[i++] = 0; } /* A082858(t1,t2) != t2 */ else /* t1 is a super-tree of t2. */ { vec[i++] = CatalanRankSexpGlobal(t1_size,new_tree); free_cons_list(new_tree,t2_size); } } } } free_cons_tree(t1); free_cons_tree(t2); } /* Could be faster if we coded a modified variant of CatalanUnrankSexp that used the subtrees of t1 as its leaf-nodes. Here we do some extra unranking of t2, because it is physically modified. However, as the time taken by these precomputing routines is minimal to the run proper, we don't care now. */ void precompute_transform2(SGRANK *vec,int upto_n,int t2_size,LRANK t2_lrank) { int t1_size; /* The size of the trees for which we are computing this transformation in the inner loop. */ LRANK t1_lrank; SEXP t1 = NULL; SEXP t2 = NULL; SEXP old_t2 = NULL; int i = 0; vec[i++] = 0; /* NIL to NIL, always! */ for(t1_size=1; t1_size <= upto_n; t1_size++) { for(t1_lrank = 0; t1_lrank < Cat(t1_size); t1_lrank++) { if(t1_size < t2_size) /* The subtree does not fit yet. */ { vec[i++] = 0; } else { t1 = CatalanUnrankSexp(t1_size,t1_lrank,&t1); if(sexp_length(t1) < t2_size) { vec[i++] = 0; } else { t2 = CatalanUnrankSexp(t2_size,t2_lrank,&old_t2); free_cons_tree(old_t2); /* If any residues left... */ /* Modifies t2 by overwriting its NILs with t1's subtrees. */ sexp_transform2(t1,t2,t2_size); vec[i++] = CatalanRankSexpGlobal(t1_size,t2); /* Free the main stem of t1. The other nodes of t1 have been passed to t2, so their time will come later, when t2 is eventually recycled: */ free_cons_list(t1,t2_size); t1 = NULL; /* Avoid double deallocations. */ old_t2 = t2; } } } } free_cons_tree(t1); free_cons_tree(t2); } /* The tree on which we do these subtree transpositions is the right-leaning "comb", i.e. the lexicographically first tree of size tree_size. Note that tp1 < tp2 is required! */ void precompute_subtree_transposition(SGRANK *vec,int upto_n,int tree_size, int tp1,int tp2) { int t_size; /* The size of the trees for which we are computing this subtree transposition in the inner loop. */ SEXP sexp = NULL; LRANK r; int i = 0; vec[i++] = 0; /* NIL to NIL, always! */ for(t_size=1; t_size <= upto_n; t_size++) { for(r = 0; r < Cat(t_size); r++) { if(t_size < tree_size) /* The subtree does not fit yet. */ { vec[i++] = 0; } else /* Our left-comb of tree_size nodes could be a subtree of */ { /* this sexp, at least in principle. */ sexp = CatalanUnrankSexp(t_size,r,&sexp); if(sexp_transpose(sexp,tree_size,tp1,tp2)) { vec[i++] = CatalanRankSexpGlobal(t_size,sexp); } else { vec[i++] = 0; } } } } free_cons_tree(sexp); } /* Precompute transpositions <1,0>; <2,1>, <2,0>; <3,2>, <3,1>, <3,0>; <4,3>, <4,2>, <4,1>, <4,0>; .... , , ... to the vectors that are allocated on the fly, and stored into vec_ptrs. */ USI **precompute_transposition_vectors(int tree_size,int upto_n,int vec_size) { USI **org_vec_ptrs; USI **vec_ptrs = ((USI **) calloc(A000217(tree_size),sizeof(USI *))); int i,j; if(NULL == vec_ptrs) { fprintf(stderr, "precompute_transposition_vectors(tree_size=%u): Couldn't allocate a vector of %u pointers where to store the transposition vectors!\n", tree_size,A000217(tree_size)); exit(1); } else { org_vec_ptrs = vec_ptrs; } for(i=1; i <= tree_size; i++) { for(j=i; j > 0; j--) { USI *new_vec = ((USI *) calloc(vec_size,sizeof(USI))); *vec_ptrs++ = new_vec; if(NULL == new_vec) { fprintf(stderr, "precompute_transposition_vectors(tree_size=%u): Couldn't allocate a %uth vector of %u short integers! i=%u, j=%u\n", tree_size,(vec_ptrs-org_vec_ptrs),vec_size,i,j); exit(1); } precompute_subtree_transposition(new_vec,upto_n,tree_size,j-1,i); } } { fprintf(stderr,"precompute_transposition_vectors: precomputed %u vectors for the clause trees of size %u\n", (vec_ptrs-org_vec_ptrs),tree_size ); } if(0 != glob_cons_cells_in_use) { fprintf(stderr,"precompute_transposition_vectors(tree_size=%u): glob_cons_cells_in_use=%lu\n", tree_size,glob_cons_cells_in_use ); } return(org_vec_ptrs); } USI **precompute_transform_vectors(int tree_size,int upto_n,int vec_size) { USI **org_vec_ptrs; USI **vec_ptrs = ((USI **) calloc(2*Cat(tree_size),sizeof(USI *))); int tree_lrank; if(NULL == vec_ptrs) { fprintf(stderr, "precompute_transform_vectors(tree_size=%u): Couldn't allocate a vector of %lu pointers where to store the transposition vectors!\n", tree_size,2*Cat(tree_size)); exit(1); } else { org_vec_ptrs = vec_ptrs; } for(tree_lrank = 0; tree_lrank < Cat(tree_size); tree_lrank++) { USI *vec1 = ((USI *) calloc(vec_size,sizeof(USI))); USI *vec2; int i; *vec_ptrs++ = vec1; if(NULL == vec1) { fprintf(stderr, "precompute_transform_vectors(tree_size=%u): Couldn't allocate a %uth vector of %u short integers! tree_lrank=%u\n", tree_size,(vec_ptrs-org_vec_ptrs),vec_size,tree_lrank); exit(1); } vec2 = ((USI *) calloc(vec_size,sizeof(USI))); *vec_ptrs++ = vec2; if(NULL == vec2) { fprintf(stderr, "precompute_transform_vectors(tree_size=%u): Couldn't allocate a %uth vector of %u short integers! tree_lrank=%u\n", tree_size,(vec_ptrs-org_vec_ptrs),vec_size,tree_lrank); exit(1); } precompute_transform1(vec1,upto_n,tree_size,tree_lrank); precompute_transform2(vec2,upto_n,tree_size,tree_lrank); for(i=0; i< vec_size; i++) { if(((0 != vec1[i]) && (vec2[vec1[i]] != i)) || ((0 != vec2[i]) && (vec1[vec2[i]] != i)) ) { fprintf(stderr, "precompute_transform_vectors(tree_size=%u): Fatal error: precompute_transform1 and _transform2 didn't produce inverses of each other! Offending terms at i=%u. tree_lrank=%u. The vec1 and vec2 are printed below:\n", tree_size,i,tree_lrank ); fflush(stderr); fprint_USI_vector(stderr,vec1,vec_size); fprint_USI_vector(stderr,vec2,vec_size); exit(1); } } } { fprintf(stderr,"precompute_transform_vectors: precomputed %u vectors for the clause trees of size %u\n", (vec_ptrs-org_vec_ptrs),tree_size ); } if(0 != glob_cons_cells_in_use) { fprintf(stderr,"precompute_transform_vectors(tree_size=%u): glob_cons_cells_in_use=%lu\n", tree_size,glob_cons_cells_in_use ); } return(org_vec_ptrs); } /**********************************************************************/ /* Better to pass this by reference... */ struct evaluation_context { int upto_size_n; int vec_size; int bijectivity_broken_at_size_n; int clear_next_time_upto_size_n; ULLI checksum; SGRANK *vals_obtained; SGRANK *vals_obtained2; /* When we need compositions. */ USI **all_vecs[(2*MAX_CLAUSE_SIZE)+2]; }; typedef struct evaluation_context ECOTEX; #ifdef DEBUG #define ECOTEX_CLEAR_vals_obtained(ep)\ printf("Clearing up to size %u, %u bytes...\n",ep->clear_next_time_upto_size_n,(sizeof(SGRANK)*A014137(ep->clear_next_time_upto_size_n)));\ memset(((BYTE *)ep->vals_obtained),0,(sizeof(SGRANK)*A014137(ep->clear_next_time_upto_size_n))) #else #define ECOTEX_CLEAR_vals_obtained(ep)\ memset(((BYTE *)ep->vals_obtained),0,(sizeof(SGRANK)*A014137(ep->clear_next_time_upto_size_n))) #endif #define ECOTEX_CLEAR_all_vals_obtained(ep)\ memset(((BYTE *)ep->vals_obtained),0,(sizeof(SGRANK)*ep->vec_size)) /* Was: memset(((BYTE *)ep->vals_obtained),0,(sizeof(SGRANK)*A014137(ep->upto_size_n))) */ typedef int (*CLAUSESEQ_COMPFUN)(USI ***,int,ECOTEX *); ECOTEX *new_ecotex(int upto_size_n) { int tree_size,i; ECOTEX *ep = ((ECOTEX *) calloc(1,sizeof(ECOTEX))); if(NULL == ep) { fprintf(stderr, "new_ecotex: Couldn't allocate a new structure of %u bytes!\n", sizeof(ECOTEX) ); exit(1); } ep->upto_size_n = upto_size_n; ep->vec_size = A014137(upto_size_n); ep->bijectivity_broken_at_size_n = 0; ep->clear_next_time_upto_size_n = 0; ep->vals_obtained = ((SGRANK *) calloc(ep->vec_size,sizeof(SGRANK))); if(NULL == ep->vals_obtained) { fprintf(stderr, "new_ecotex: Couldn't allocate a value table vals_obtained of %u short integers (%u bytes). upto_size_n=%u\n", ep->vec_size,ep->vec_size*sizeof(SGRANK),upto_size_n ); exit(1); } ep->vals_obtained2 = ((SGRANK *) calloc(ep->vec_size,sizeof(SGRANK))); if(NULL == ep->vals_obtained2) { fprintf(stderr, "new_ecotex: Couldn't allocate a value table vals_obtained2 of %u short integers (%u bytes). upto_size_n=%u\n", ep->vec_size,ep->vec_size*sizeof(SGRANK),upto_size_n ); exit(1); } for(i=0, tree_size = 1; tree_size <= upto_size_n; tree_size++) { ep->all_vecs[i++] = precompute_transposition_vectors(tree_size,upto_size_n,ep->vec_size); ep->all_vecs[i++] = precompute_transform_vectors(tree_size,upto_size_n,ep->vec_size); } return(ep); } /**********************************************************************/ #define globrank(tree_size,lrank) ((0 == (tree_size)) ? 0 : A014137(tree_size-1)+(lrank)) /* Size = 2+2+4+1 = 9 bytes. */ struct s_raw_clause { SSIZE s_size; SLRANK s_lrank; SLRANK d_lrank; PERMRANK p_rank; }; /* Size = 4+4+1+1 = 10 bytes in 32-bit machines, 8+4+1+1 = 14 bytes if 64-bit pointers. */ struct s_raw_clauseseq_header { SSIZE n_clauses; BYTE binexp; /* Only upto binexp=255, but even that is enough! */ struct s_clause *next_clauseseq; /* In the hash table bucket, NULL if last one. */ NORERANK nth_distinct_bijection; /* Zero-based position in A089840. */ }; /* Clause sequence is essentially a vector of clauses (CLAUSE *) where the first allocated clause structure is however interpreted as struct_raw_clauseseq_header, instead of s_raw_clause: */ struct s_clause { union { struct s_raw_clause cl; /* In this order, so the explicit initializations think they initialize clauses, not headers. */ struct s_raw_clauseseq_header hd; } u; }; typedef struct s_clause CLAUSE; typedef CLAUSE *CLAUSESEQ; #define CLAUSESEQ_next_clauseseq(CS) ((&CS[0])->u.hd.next_clauseseq) /* In the hash bucket. */ #define CLAUSESEQ_n_clauses(CS) ((&CS[0])->u.hd.n_clauses) /* The number of clauses that follow. */ #define CLAUSESEQ_binexp(CS) ((&CS[0])->u.hd.binexp) /* The integer n whose binexp's runcounts are used for the composition of binwidth(n). */ #define CLAUSESEQ_norerank(CS) ((&CS[0])->u.hd.nth_distinct_bijection) #define SET_CLAUSESEQ_norerank(CS,NTH) (CLAUSESEQ_norerank(CS) = (NTH)) /* We should check for overflows > 4294967295 */ /* For explicit initializations. More than slightly dangerous! For now we have to omit binexp, because only the very first members, s_size and n_clauses of the structures s_raw_clause and s_raw_clauseseq_header are of the same size. */ #define CLAUSESEQ_begin(binexp,n_clauses) { (n_clauses),0,0,0 } #define CLAUSESEQ_nthclause(CS,n) (CS[n]) #define CLAUSE_size(c) ((unsigned int)((c).u.cl.s_size)) #define SET_CLAUSE_size(c,s) (((c).u.cl.s_size) = ((SSIZE) (s))) #define CLAUSE_s_lrank(c) ((c).u.cl.s_lrank) #define CLAUSE_d_lrank(c) ((c).u.cl.d_lrank) #define CLAUSE_p_rank(c) ((c).u.cl.p_rank) #define SET_CLAUSE_s_lrank(c,r) (((c).u.cl.s_lrank) = (r)) #define SET_CLAUSE_d_lrank(c,r) (((c).u.cl.d_lrank) = (r)) #define SET_CLAUSE_p_rank(c,r) (((c).u.cl.p_rank) = (r)) #define CLAUSE_s_grank(c) globrank(CLAUSE_size(c),CLAUSE_s_lrank(c)) #define CLAUSE_d_grank(c) globrank(CLAUSE_size(c),CLAUSE_d_lrank(c)) void fprint_clause_seq(FILE *fp,CLAUSESEQ clause_seq) { int cs_length = CLAUSESEQ_n_clauses(clause_seq); int i; fprintf(fp,"c%u",CLAUSESEQ_binexp(clause_seq)); for(i=1; i <= cs_length; i++) { fprintf(fp,"__%u_%u_%u_%u", CLAUSE_size(CLAUSESEQ_nthclause(clause_seq,i)), CLAUSE_s_lrank(CLAUSESEQ_nthclause(clause_seq,i)), CLAUSE_d_lrank(CLAUSESEQ_nthclause(clause_seq,i)), CLAUSE_p_rank(CLAUSESEQ_nthclause(clause_seq,i)) ); } } CLAUSESEQ clause_seq_dup(CLAUSESEQ clause_seq) { int cs_size = (1+CLAUSESEQ_n_clauses(clause_seq)); CLAUSESEQ dup_seq = ((CLAUSESEQ) calloc(cs_size,sizeof(CLAUSE))); if(NULL == dup_seq) { fprintf(stderr, "clause_seq_dup: couldn't allocate %u bytes for duplicating the clause:\n", (cs_size*sizeof(CLAUSE)) ); fprint_clause_seq(stderr,clause_seq); fprintf(stderr,"\n"); exit(1); } else { memcpy(dup_seq,clause_seq,(sizeof(CLAUSE)*cs_size)); return(dup_seq); } } typedef ULLI FILTERMASK; #define FILTER_IF_NONULTIMATE_1S two_to(0) #define FILTER_LONE_NONSTANDALONES two_to(1) #define contains(mask,bit) ((mask) & (bit)) /* Tries to increment clause's field one by one, and returns 0 if it wrapped around back to zero. */ int increment_clause(CLAUSE *cp) { int s = CLAUSE_size(*cp); if((++(CLAUSE_p_rank(*cp))) < A000142(s+1)) { return(1); } SET_CLAUSE_p_rank(*cp,0); if((++(CLAUSE_d_lrank(*cp))) < A000108(s)) { return(1); } SET_CLAUSE_d_lrank(*cp,0); if((++(CLAUSE_s_lrank(*cp))) < A000108(s)) { return(1); } SET_CLAUSE_s_lrank(*cp,0); return(0); } /* Note that this should be used only on clause sequences with enough extra CLAUSE-structures reserved in the end, because the clause sequence can here effectively grow in length. (e.g when binexp changes from 15 to 16, it grows in length from 1 to 2 clauses.) */ int clause_seq_next_composition(CLAUSESEQ clause_seq,FILTERMASK filmask) { int org_binexp,binexp,i,c,prevbit; do_it_again: org_binexp = binexp = ++(CLAUSESEQ_binexp(clause_seq)); i=0; c=0; prevbit = (binexp & 1); /* Compute the runcounts of new binexp, to get another composition (an ordered partition): */ while(0 != binexp) { c++; binexp >>= 1; if((binexp & 1) != prevbit) { if(contains(filmask,FILTER_IF_NONULTIMATE_1S) && (0 != binexp) && (1 != org_binexp) && (1 == c)) { goto do_it_again; } SET_CLAUSE_size(CLAUSESEQ_nthclause(clause_seq,++i),c); SET_CLAUSE_s_lrank(CLAUSESEQ_nthclause(clause_seq,i),0); SET_CLAUSE_d_lrank(CLAUSESEQ_nthclause(clause_seq,i),0); SET_CLAUSE_p_rank(CLAUSESEQ_nthclause(clause_seq,i),0); c = 0; prevbit = (1-prevbit); } } CLAUSESEQ_n_clauses(clause_seq) = i; return(i); } /* Increment clauses with an odometer-principle, from right to left: */ int clause_seq_next(CLAUSESEQ clause_seq,FILTERMASK filmask) { int cs_length = CLAUSESEQ_n_clauses(clause_seq); int i; for(i= cs_length; i > 0; i--) { CLAUSE *cp = &(CLAUSESEQ_nthclause(clause_seq,i)); if(increment_clause(cp)) { if(contains(filmask,FILTER_LONE_NONSTANDALONES) && (1 == cs_length) && (CLAUSE_s_lrank(*cp) != CLAUSE_d_lrank(*cp)) ) { /* When d_rank has incremented past s_rank (means p_rank is 0), then set s_rank to be equal to d_rank, if we are filtering off lone nonstandalones already at this stage: */ CLAUSE_s_lrank(*cp) = CLAUSE_d_lrank(*cp); } return(i); } /* Otherwise, a clause wrapped around, so let's continue incrementing clauses to the left. */ } /* If we fall from the above loop, it means that ALL clauses have wrapped around to zero, so now it's time to get the next composition, i.e. to increment the binexp of the clause sequence: */ clause_seq_next_composition(clause_seq,filmask); return(i); } /**********************************************************************/ /**********************************************************************/ /* ;; The following algorithm is a slight modification of unrank1 ;; algorithm as presented by W. Myrvold and F. Ruskey, in ;; Ranking and Unranking Permutations in Linear Time, ;; Inform. Process. Lett. 79 (2001), no. 6, 281-284. ;; Available on-line: http://www.cs.uvic.ca/~fruskey/Publications/RankPerm.html (define (permute-A060118 elems size permrank) (let ((p (vector-head elems size))) (let unrankA060118 ((r permrank) (i 1) ) (cond ((zero? r) p) (else (let* ((j (1+ i)) (m (modulo r j)) ) (cond ((not (zero? m)) ;; Swap at i and (i-(r mod (i+1))) (let ((org-i (vector-ref p i))) (vector-set! p i (vector-ref p (- i m))) (vector-set! p (- i m) org-i) ) ) ) (unrankA060118 (/ (- r m) j) j) ) ) ) ) ) ) */ void reverse_n_vectors(USI **vecs,int n_vecs) { int i,j; for(i=0,j=n_vecs-1; i < j; i++, j--) { USI *tmp = vecs[i]; vecs[i] = vecs[j]; vecs[j] = tmp; } } /* The following in effect uses the same algorithm A060118. Note that in the loop we set m to be the ith least significant digit (1-based) in original permrank's factorial expansion. */ int select_transposition_vectors(USI **res,USI **transp_vecs,int permrank) { int vecs_needed=0; int i,m; for(i=1; (permrank != 0); transp_vecs += i, i++) { m = (permrank % (i+1)); if(0 != m) { res[vecs_needed++] = transp_vecs[m-1]; } permrank = (permrank - m)/(i+1); } return(vecs_needed); } /* The maximum number of vectors that we might be selected here is tree_size+2. Thus there must be at least tree_size+3 pointers of space in res. Note that clause's tree_size cannot be more than upto_size_n (the second parameter of the command line), otherwise the references to all_vecs will overflow. */ int select_vectors_needed(USI **res,USI ***all_vecs,int tree_size,LRANK s_rank,LRANK d_rank,int permrank) { int vecs_needed=0; USI **transp_vecs = all_vecs[2*(tree_size-1)]; USI **transform_vecs = all_vecs[(2*(tree_size-1))+1]; /* The first optimization: use transform1 only when either s_rank is not zero, OR when ALL ranks are zeroes. The latter because we need in that case at least one vector reference to know whether the computed tree is a super-tree of the tree corresponding to s_rank and d_rank, i.e. the right comb: */ if((0 != s_rank) || ((0 == permrank) && (0 == d_rank)) ) { res[vecs_needed++] = transform_vecs[(2*s_rank)]; /* transform1 vectors on even indices. */ } vecs_needed += select_transposition_vectors(res+vecs_needed,transp_vecs,permrank); /* The second "optimization": We don't need transform2 when destination tree is the right comb. */ if(0 != d_rank) { res[vecs_needed++] = transform_vecs[(2*d_rank)+1]; /* transform2 vectors on odd indices. */ } res[vecs_needed++] = NULL; /* As an end marker. */ return(vecs_needed); /* Return one more than the number of vecs needed, because of ending NULL. */ } /* In essence, the inverse of the clause sequence ((s1,d1,p1), (s2,d2,p2), ..., (sn,dn,pn)) is the clause sequence: ((d1,s1,A060125(p1)), (d2,s2,A060125(p2)), ..., (dn,sn,A060125(pn))) (where the triplet (si,di,pi) indicates the source tree, destination tree and permutation ranks of the ith clause) but because I don't care to code A060125 in C (unless I find a more elegant algorithm for it) I instead here simply swap s and d, and reverse the order of transposition vectors selected: */ int select_vectors_for_inverse(USI **res,USI ***all_vecs,int tree_size,LRANK d_rank,LRANK s_rank,int permrank) { int vecs_needed=0, t_vecs_needed; USI **transp_vecs = all_vecs[2*(tree_size-1)]; USI **transform_vecs = all_vecs[(2*(tree_size-1))+1]; if((0 != s_rank) || ((0 == permrank) && (0 == d_rank)) ) { res[vecs_needed++] = transform_vecs[(2*s_rank)]; /* transform1 vectors on even indices. */ } t_vecs_needed = select_transposition_vectors(res+vecs_needed,transp_vecs,permrank); reverse_n_vectors(res+vecs_needed,t_vecs_needed); vecs_needed += t_vecs_needed; if(0 != d_rank) { res[vecs_needed++] = transform_vecs[(2*d_rank)+1]; /* transform2 vectors on odd indices. */ } res[vecs_needed++] = NULL; /* As an end marker. */ return(vecs_needed); /* Return one more than the number of vecs needed, because of ending NULL. */ } /* For maximum performance: Unroll the loops here, and inline this to a calling function. Note: when n_clauses is 0, then it doesn't matter what is the contents of clauseseqs_selected_vectors, but instead this will return x intact back, as an empty set of clauses correspond to a single default clause, that is, an identity function A001477. */ SGRANK apply_clauses_to_x(SGRANK x,USI ***clauseseqs_selected_vectors,int n_clauses) { int i=0; SGRANK y; USI **selected_vectors; do { if(i == n_clauses) /* No clause matched, return the x back as a default. */ { return(x); } selected_vectors = clauseseqs_selected_vectors[i++]; y = (*selected_vectors)[x]; /* Apply the first vector (transform 1) to x. */ } while(0 == y); /* until we find a clause that matches. */ while(NULL != *++selected_vectors) { y = (*selected_vectors)[y]; } return(y); } /* This compares the the clause_seqs results to another one's results (computed with the next one, to_obtain_signature_perm, and placed by it to ep->vals_obtained) and returns the first point (x) where their values differ (Note: x=0 is not checked, as it is assumed to stay always zero). and zero if they are equal. */ int to_compare_to_another_signature_perm(USI ***clauseseqs_selected_vectors,int n_clauses,ECOTEX *ep) { int t_size; LRANK t_lrank; SGRANK x; SGRANK y; for(t_size=1, x=1; t_size <= ep->upto_size_n; t_size++) { for(t_lrank = 0; t_lrank < Cat(t_size); t_lrank++, x++) { y = apply_clauses_to_x(x,clauseseqs_selected_vectors,n_clauses); if(y != ep->vals_obtained[x]) { return(x); } } } return(0); } int to_obtain_signature_perm(USI ***clauseseqs_selected_vectors,int n_clauses,ECOTEX *ep) { int t_size; LRANK t_lrank; SGRANK x; SGRANK y; ULLI checksum = 0; ep->vals_obtained[0] = 0; for(t_size=1, x=1; t_size <= ep->upto_size_n; t_size++) { for(t_lrank = 0; t_lrank < Cat(t_size); t_lrank++, x++) { y = apply_clauses_to_x(x,clauseseqs_selected_vectors,n_clauses); ep->vals_obtained[x] = y; checksum = (31*checksum)+y; /* Magic checksum algorithm... MUST be the same as in to_check_the_bijectivity */ } } ep->checksum = checksum; return(1); } /* We use this one for the latter, when we want to compute the composition of two clauseseqs, as this applies its gatomorphism onto *the old* values of ep->vals_obtained vector, instead of x directly. */ int to_obtain_composite_signature_perm(USI ***clauseseqs_selected_vectors,int n_clauses,ECOTEX *ep) { int t_size; LRANK t_lrank; SGRANK x; SGRANK y; ULLI checksum = 0; /* Make a safe copy of the previously computed ep->vals_obtained vector: */ memcpy(ep->vals_obtained2,ep->vals_obtained,(sizeof(SGRANK)*ep->vec_size)); ep->vals_obtained[0] = 0; for(t_size=1, x=1; t_size <= ep->upto_size_n; t_size++) { for(t_lrank = 0; t_lrank < Cat(t_size); t_lrank++, x++) { y = apply_clauses_to_x(ep->vals_obtained2[x],clauseseqs_selected_vectors,n_clauses); ep->vals_obtained[x] = y; checksum = (31*checksum)+y; /* Magic checksum algorithm... */ } } ep->checksum = checksum; return(1); } int to_check_the_bijectivity(USI ***clauseseqs_selected_vectors,int n_clauses,ECOTEX *ep) { int t_size; LRANK t_lrank; SGRANK x; SGRANK y; ULLI checksum = 0; SGRANK *vals_obtained = ep->vals_obtained; ECOTEX_CLEAR_vals_obtained(ep); for(t_size=1, x=1; t_size <= ep->upto_size_n; t_size++) { for(t_lrank = 0; t_lrank < Cat(t_size); t_lrank++, x++) { y = apply_clauses_to_x(x,clauseseqs_selected_vectors,n_clauses); if(0 == y) { fprintf(stderr,"to_check_the_bijectivity: Bad internal error: y=0 for x=%u\n",x); exit(1); } if(0 != (vals_obtained[y])++) /* Not the first time for this value y? */ { ep->bijectivity_broken_at_size_n = t_size; /* Because not injective. */ ep->clear_next_time_upto_size_n = t_size; return(0); } checksum = (31*checksum)+y; /* Magic checksum algorithm... */ } } ep->clear_next_time_upto_size_n = ep->upto_size_n; ep->checksum = checksum; return(1); } int compute_clause_seq(CLAUSESEQ_COMPFUN for_what,CLAUSESEQ clause_seq,ECOTEX *ep) { USI *space_for_selected_vectors[MAX_CLAUSES*(MAX_CLAUSE_SIZE+3)]; /* Vector of vectors. */ USI **ssp = space_for_selected_vectors; /* Pointer to above. */ USI **clauseseqs_selected_vectors[MAX_CLAUSES+1]; /* Vector of vector of vectors */ int n_clauses = CLAUSESEQ_n_clauses(clause_seq); int i; for(i=1; i <= n_clauses; i++) { clauseseqs_selected_vectors[i-1] = ssp; ssp += select_vectors_needed(ssp,ep->all_vecs, CLAUSE_size(CLAUSESEQ_nthclause(clause_seq,i)), CLAUSE_s_lrank(CLAUSESEQ_nthclause(clause_seq,i)), CLAUSE_d_lrank(CLAUSESEQ_nthclause(clause_seq,i)), CLAUSE_p_rank(CLAUSESEQ_nthclause(clause_seq,i)) ); #ifdef DEBUG fprint_vector_of_USI_vectors(stdout,clauseseqs_selected_vectors[i-1],A014137(ep->upto_size_n)); #endif } return(for_what(clauseseqs_selected_vectors,n_clauses,ep)); } int compute_inverse_of_clause_seq(CLAUSESEQ_COMPFUN for_what,CLAUSESEQ clause_seq,ECOTEX *ep) { USI *space_for_selected_vectors[MAX_CLAUSES*(MAX_CLAUSE_SIZE+3)]; /* Vector of vectors. */ USI **ssp = space_for_selected_vectors; /* Pointer to above. */ USI **clauseseqs_selected_vectors[MAX_CLAUSES+1]; /* Vector of vector of vectors */ int n_clauses = CLAUSESEQ_n_clauses(clause_seq); int i; for(i=1; i <= n_clauses; i++) { clauseseqs_selected_vectors[i-1] = ssp; ssp += select_vectors_for_inverse(ssp,ep->all_vecs, CLAUSE_size(CLAUSESEQ_nthclause(clause_seq,i)), CLAUSE_s_lrank(CLAUSESEQ_nthclause(clause_seq,i)), CLAUSE_d_lrank(CLAUSESEQ_nthclause(clause_seq,i)), CLAUSE_p_rank(CLAUSESEQ_nthclause(clause_seq,i)) ); #ifdef DEBUG fprint_vector_of_USI_vectors(stdout,clauseseqs_selected_vectors[i-1],A014137(ep->upto_size_n)); #endif } return(for_what(clauseseqs_selected_vectors,n_clauses,ep)); } /**********************************************************************/ CLAUSESEQ clauseseq_found_in_bucket(CLAUSESEQ clauseseqs_with_same_checksum,CLAUSESEQ *p_to_last,ECOTEX *ep) { static int max_bucklen=0; int bucklen=1; while(NULL != clauseseqs_with_same_checksum) { if(0 == compute_clause_seq(to_compare_to_another_signature_perm,clauseseqs_with_same_checksum,ep)) { return(clauseseqs_with_same_checksum); } *p_to_last = clauseseqs_with_same_checksum; clauseseqs_with_same_checksum = CLAUSESEQ_next_clauseseq(clauseseqs_with_same_checksum); bucklen++; } /* This message is not printed until after we get the first gatomorphism with the same hash table index as a previously encountered, but different gatomorphism: */ if(bucklen > max_bucklen) { max_bucklen = bucklen; fprintf(stderr,"clauseseq_found_in_bucket: maximum bucket length of the hash table will/might soon be: %u\n", max_bucklen); } return(NULL); /* Didn't find it! */ } /* Returns the first previously encountered CLAUSESEQ of which clause_seq is a duplicate of, or if it is so far unique, then add it to a hash table, and return NULL. */ CLAUSESEQ duplicate_found_in_hashtable(CLAUSESEQ clause_seq,CLAUSESEQ *hashtable,int hashsize,ECOTEX *ep,CLAUSESEQ *p_to_stored_one,int compute_signature_perm_in_any_case) { int hindex; CLAUSESEQ dup,last_in_bucket; if(compute_signature_perm_in_any_case) { compute_clause_seq(to_obtain_signature_perm,clause_seq,ep); } hindex = (ep->checksum % hashsize); if(NULL == hashtable[hindex]) /* The very first one for this index. */ { CLAUSESEQ_next_clauseseq(clause_seq) = NULL; /* Should be already! */ hashtable[hindex] = clause_seq_dup(clause_seq); if(NULL != p_to_stored_one) { *p_to_stored_one = hashtable[hindex]; } return(NULL); } /* Else, we have to check whether it produces *really* the same sequence as one of the clause sequences in this bucket. So, first, obtain the signature permutation for this clause sequence (in case it has not already been computed above, and stored in ep->vals_obtained): */ if(!compute_signature_perm_in_any_case) { compute_clause_seq(to_obtain_signature_perm,clause_seq,ep); } dup = clauseseq_found_in_bucket(hashtable[hindex],&last_in_bucket,ep); /* Didn't find a duplicate, so add this new one to the end of the list: */ if(NULL == dup) { CLAUSESEQ tmp; CLAUSESEQ_next_clauseseq(clause_seq) = NULL; /* Should be already! */ CLAUSESEQ_next_clauseseq(last_in_bucket) = tmp = clause_seq_dup(clause_seq); if(NULL != p_to_stored_one) { *p_to_stored_one = tmp; } } return(dup); } CLAUSESEQ composition_found_in_hashtable(CLAUSESEQ clause_seq1,CLAUSESEQ clause_seq2,CLAUSESEQ *hashtable,int hashsize,ECOTEX *ep) { int hindex; CLAUSESEQ last_in_bucket; /* Not used here, but clauseseq_found_in_bucket expects it anyway. */ compute_clause_seq(to_obtain_signature_perm,clause_seq1,ep); compute_clause_seq(to_obtain_composite_signature_perm,clause_seq2,ep); hindex = (ep->checksum % hashsize); if(NULL == hashtable[hindex]) /* The very first one for this index? So composition is not among the so-far computed. */ { return(NULL); } /* Else, we check whether the compositions produces really the same sequence as one of the clause sequences in this bucket: */ return(clauseseq_found_in_bucket(hashtable[hindex],&last_in_bucket,ep)); } CLAUSESEQ inverse_found_in_hashtable(CLAUSESEQ clause_seq1,CLAUSESEQ *hashtable,int hashsize,ECOTEX *ep) { int hindex; CLAUSESEQ last_in_bucket; /* Not used here, but clauseseq_found_in_bucket expects it anyway. */ compute_inverse_of_clause_seq(to_obtain_signature_perm,clause_seq1,ep); hindex = (ep->checksum % hashsize); if(NULL == hashtable[hindex]) /* The very first one for this index? So the inverse is not among the so-far computed. */ { return(NULL); } /* Else, we check that the inverse produces really the same sequence as one of the clause sequences in this bucket: */ return(clauseseq_found_in_bucket(hashtable[hindex],&last_in_bucket,ep)); } /* The main loop is here! We start from an identity gatomorphism of no clauses but the implicit default. Here are the 50 A-numbers you requested: 89831 --- 89880. */ int check_upto_binexp_n(int upto_binexp_n,ECOTEX *ep, int hashsize,CLAUSESEQ *hashtable,CLAUSESEQ *first_clauseseqs,FILTERMASK filmask, int misc_seq_size,int misc_table_size,int max_order_tested_to, ULLI *A089832,ULLI *A089833, NORERANK *A089839,SGRANK *A089840,NORERANK *misc_seq_vec) { int max_bijectivity_broken_at_size_n=0; int prev_binexp = 0, prev_n_clauses = 0, prev_tot_size = 0; int row,i; ULLI bijections=0, distinct_bijections_among_this_size=0, distinct_bijections=0, duplicates=0; CLAUSE clause_seq[MAX_CLAUSES+1]; CLAUSESEQ dupe = NULL; prev_binexp = CLAUSESEQ_n_clauses(clause_seq) = CLAUSESEQ_binexp(clause_seq) = 0; ECOTEX_CLEAR_all_vals_obtained(ep); do { /* fprint_clause_seq(stdout,clause_seq); printf(" in the top of the loop!\n"); fflush(stdout); */ if(compute_clause_seq(to_check_the_bijectivity,clause_seq,ep)) { CLAUSESEQ copied_one = NULL; /* fprint_clause_seq(stdout,clause_seq); printf(" is a bijection\n"); fflush(stdout); */ bijections++; dupe = duplicate_found_in_hashtable(clause_seq,hashtable,hashsize,ep, &copied_one, /* If was not duplicate, i.e. if dupe is NULL. */ 0); /* No need for signature perm, as long as the bucket is empty. */ if(NULL == dupe) /* equal to: (NULL != copied_one) */ { if(distinct_bijections < misc_seq_size) { first_clauseseqs[distinct_bijections] = copied_one; } SET_CLAUSESEQ_norerank(copied_one,distinct_bijections); /* Set its rank in A089840. */ /* If the fourth argument misc_table_size was specified on the command line, then we are collecting from the first misc_table_size distinct bijections their signature permutations, to the square array A089840: */ if((NULL != A089840) && (distinct_bijections < misc_table_size)) { compute_clause_seq(to_obtain_signature_perm,copied_one,ep); for(i=0; i < (misc_table_size-distinct_bijections); i++) { A089840[ar2seqind(distinct_bijections,i)] = ep->vals_obtained[i]; } ECOTEX_CLEAR_all_vals_obtained(ep); } distinct_bijections_among_this_size++; distinct_bijections++; /* In total. */ } else { duplicates++; #ifdef DEBUG fprint_clause_seq(stdout,clause_seq); fprintf(stdout," (%lx -> %u) is a duplicate of ",ep->checksum,(ep->checksum % hashsize)); fprint_clause_seq(stdout,dupe); fprintf(stdout,"\n"); #endif } } if(ep->bijectivity_broken_at_size_n > max_bijectivity_broken_at_size_n) { fprintf(stdout,"New record for bijectivity_broken_at (%u) occurred with the clause seq: ", ep->bijectivity_broken_at_size_n); fprint_clause_seq(stdout,clause_seq); fprintf(stdout,"\n"); max_bijectivity_broken_at_size_n = ep->bijectivity_broken_at_size_n; } clause_seq_next(clause_seq,filmask); /* fprint_clause_seq(stdout,clause_seq); printf(" is the next clause seq.\n"); fflush(stdout); */ if(CLAUSESEQ_binexp(clause_seq) != prev_binexp) { int i; int tot_size = binwidth(prev_binexp); printf("c%u*: ",prev_binexp); fprint_ulli(stdout,bijections); printf(" "); fprint_ulli(stdout,distinct_bijections_among_this_size); printf("\n"); fflush(stdout); A089832[tot_size] += distinct_bijections_among_this_size; if(tot_size > 0) { A089833[tr2seqind11(tot_size,prev_n_clauses)] += distinct_bijections_among_this_size; } if(tot_size != prev_tot_size) { printf("A089832[%u]=",prev_tot_size); fflush(stdout); fprint_ulli(stdout,A089832[prev_tot_size]); printf(": "); fflush(stdout); for(i=1; i <= prev_tot_size; i++) { printf(" "); fprint_ulli(stdout,A089833[tr2seqind11(prev_tot_size,i)]); } printf("\nduplicates so far: "); fprint_ulli(stdout,duplicates); printf("\n"); fflush(stdout); prev_tot_size = tot_size; } distinct_bijections_among_this_size = bijections = 0; prev_binexp = CLAUSESEQ_binexp(clause_seq); prev_n_clauses = CLAUSESEQ_n_clauses(clause_seq); } } while(CLAUSESEQ_binexp(clause_seq) <= upto_binexp_n); for(row=1; row <= binwidth(upto_binexp_n); row++) { printf("A089832[%u]=",row); fprint_ulli(stdout,A089832[row]); printf(": "); for(i=1; i <= row; i++) { if(i > 1) { printf(","); } fprint_ulli(stdout,A089833[tr2seqind11(row,i)]); } printf("\n"); } globopt_list_begin = ""; globopt_list_end = "\n"; printf("max_bijectivity_broken_at_size_n=%u\n", max_bijectivity_broken_at_size_n); printf("duplicates in total: "); fprint_ulli(stdout,duplicates); printf("\n"); fflush(stdout); output_OEIS_sequence(stdout,89831,A089833,A000217(binwidth(upto_binexp_n)), VECTYPE_ULLI,(PF_VEC_OUT)fprint_ULLI_vector,1,OEIS_TABLE, "The triangle T(n,m) read as T(1,1); T(2,1), T(2,2); T(3,1), T(3,2), T(3,3); ... Number of distinct non-recursive gatomorphisms of total n opened nodes, divided into m non-default clauses."); output_OEIS_sequence(stdout,89832,A089832,binwidth(upto_binexp_n)+1, VECTYPE_ULLI,(PF_VEC_OUT)fprint_ULLI_vector,0,OEIS_SEQ, "Number of distinct non-recursive gatomorphisms of total n opened nodes. Row sums of A089831."); /* The first column of A089831: A089833 = (Cat(n)*((n+1)!-1)) = A001813[n]-A000108[n] 0,1,10,115,1666,30198,665148,17296851,518916970,... (not in OEIS) A089834 = A089832-A089833, the row sums of A089831 excluding the first column. A089835= Cat(n)*Cat(n)*(n+1)! (= A001246[n]*A000142[n+1] = A001813[n]*A000108[n]) 1,2,24,600,23520,1270080,87816960,7420533120,742053312000 (not in EIS) = number of clauses of tree size n. (Offset 0). A089836 = INVERT(A089835) = INVERT([seq(Cat(n)*Cat(n)*(n+1)!,n=1..9)]); [2,28,704,26800,1404416,94890112,7887853568,779773444864,89407927009280,...] */ if(misc_table_size > 0) /* E.g. (NULL != A089840) */ { int antidiagonal,other; output_OEIS_sequence(stdout,89840,A089840,A000217(misc_table_size), VECTYPE_USI,(PF_VEC_OUT)fprint_USI_vector,0,OEIS_TABLE, "The signature permutations of non-recursive gatomorphisms, sorted according to the minimum number of opening nodes needed in their defining clauses." ); for(antidiagonal=0, i=0; antidiagonal < misc_table_size; antidiagonal++) { for(other=0; other <= antidiagonal; other++) { /* printf("Composition A089840[%u] o A089840[%u] = ",other,(antidiagonal-other)); */ CLAUSESEQ comp = composition_found_in_hashtable(first_clauseseqs[(antidiagonal-other)],first_clauseseqs[other],hashtable,hashsize,ep); A089839[i++] = ((NULL == comp) ? NA_SENTINEL : CLAUSESEQ_norerank(comp)); } } output_OEIS_sequence(stdout,89839,A089839,A000217(misc_table_size), VECTYPE_NORERANK,(PF_VEC_OUT)fprint_NORERANK_vector_with_NA_SENTINELs,0,OEIS_TABLE, "The composition table T(n,m): ... of the nth and mth gatomorphisms in the table A089840." ); } /* Compute the involutions A089837 & A089838. */ if(misc_seq_size > 1) { int i; for(i=0; i < misc_seq_size; i++) { CLAUSESEQ comp = composition_found_in_hashtable(first_clauseseqs[i],first_clauseseqs[1],hashtable,hashsize,ep); misc_seq_vec[i] = ((NULL == comp) ? NA_SENTINEL : CLAUSESEQ_norerank(comp)); } output_OEIS_sequence(stdout,89838,misc_seq_vec,misc_seq_size, VECTYPE_NORERANK,(PF_VEC_OUT)fprint_NORERANK_vector_with_NA_SENTINELs,0,OEIS_SEQ, "The row 1 of A089839." ); for(i=0; i < misc_seq_size; i++) { CLAUSESEQ comp = composition_found_in_hashtable(first_clauseseqs[1],first_clauseseqs[i],hashtable,hashsize,ep); misc_seq_vec[i] = ((NULL == comp) ? NA_SENTINEL : CLAUSESEQ_norerank(comp)); } output_OEIS_sequence(stdout,89837,misc_seq_vec,misc_seq_size, VECTYPE_NORERANK,(PF_VEC_OUT)fprint_NORERANK_vector_with_NA_SENTINELs,0,OEIS_SEQ, "The column 1 of A089839." ); for(i=0; i < misc_seq_size; i++) { CLAUSESEQ comp = composition_found_in_hashtable(first_clauseseqs[i],first_clauseseqs[i],hashtable,hashsize,ep); misc_seq_vec[i] = ((NULL == comp) ? NA_SENTINEL : CLAUSESEQ_norerank(comp)); } output_OEIS_sequence(stdout,89841,misc_seq_vec,misc_seq_size, VECTYPE_NORERANK,(PF_VEC_OUT)fprint_NORERANK_vector_with_NA_SENTINELs,0,OEIS_SEQ, "The main diagonal of A089839." /* Should be A089839(A089837(n),A089838(n)) or vice versa? */ ); /* We could compute the position of each one's A057163-conjugate, in similar way: */ for(i=0; i < misc_seq_size; i++) { CLAUSESEQ comp = inverse_found_in_hashtable(first_clauseseqs[i],hashtable,hashsize,ep); misc_seq_vec[i] = ((NULL == comp) ? NA_SENTINEL : CLAUSESEQ_norerank(comp)); } output_OEIS_sequence(stdout,89843,misc_seq_vec,misc_seq_size, VECTYPE_NORERANK,(PF_VEC_OUT)fprint_NORERANK_vector_with_NA_SENTINELs,0,OEIS_SEQ, "Involution of natural numbers which shows the position of the inverse of each non-recursive gatomorphism in table A089840" ); for(i=0; i < misc_seq_size; i++) { int order=1; CLAUSESEQ comp = first_clauseseqs[i]; do { if(comp == first_clauseseqs[0]) { break; } comp = composition_found_in_hashtable(comp,first_clauseseqs[i],hashtable,hashsize,ep); } while((NULL != comp) && (++order <= max_order_tested_to)); /* misc_seq_vec[i] = (NORERANK) ((NULL == comp) ? NA_SENTINEL : ((first_clauseseqs[0] == comp) ? order : 0)); */ misc_seq_vec[i] = (NORERANK) ((first_clauseseqs[0] == comp) ? order : 0); } /* Note: The order of the 15th and 21st elements (A089859 & A089863), should be 4, not infinite! I.e. we must compute this with high enough value of upto_binexp_n (127 is recommended). */ output_OEIS_sequence(stdout,89842,misc_seq_vec,misc_seq_size, VECTYPE_NORERANK,(PF_VEC_OUT)fprint_NORERANK_vector_with_NA_SENTINELs,0,OEIS_SEQ, "The order of each element in A089840, 0 if not finite." ); } return(distinct_bijections); } /**********************************************************************/ /* See notes at the definition of CLAUSESEQ_begin macro ! */ CLAUSE gmA001477[] = { CLAUSESEQ_begin(0,0) }; /* A089840[0] */ CLAUSE gmA069770[] = { CLAUSESEQ_begin(1,1), { 1, 0, 0, 1 } };/* A089840[1] */ CLAUSE gmA072796[] = { CLAUSESEQ_begin(3,1), { 2, 0, 0, 1 } };/* A089840[2] */ CLAUSE gmA089850[] = { CLAUSESEQ_begin(3,1), { 2, 0, 0, 2 } };/* A089840[3] */ CLAUSE gmA089851[] = { CLAUSESEQ_begin(3,1), { 2, 0, 0, 3 } };/* A089840[4] */ CLAUSE gmA089852[] = { CLAUSESEQ_begin(3,1), { 2, 0, 0, 4 } };/* A089840[5] */ CLAUSE gmA089853[] = { CLAUSESEQ_begin(3,1), { 2, 0, 0, 5 } };/* A089840[6] */ CLAUSE gmA089854[] = { CLAUSESEQ_begin(3,1), { 2, 1, 1, 1 } };/* A089840[7] */ CLAUSE gmA072797[] = { CLAUSESEQ_begin(3,1), { 2, 1, 1, 2 } };/* A089840[8] */ CLAUSE gmA089855[] = { CLAUSESEQ_begin(3,1), { 2, 1, 1, 3 } };/* A089840[9] */ CLAUSE gmA089856[] = { CLAUSESEQ_begin(3,1), { 2, 1, 1, 4 } };/* A089840[10] */ CLAUSE gmA089857[] = { CLAUSESEQ_begin(3,1), { 2, 1, 1, 5 } };/* A089840[11] */ CLAUSE gmA074679[] = { CLAUSESEQ_begin(4,2), { 2, 0, 1, 0,}, { 1, 0, 0, 1 } }; /* A089840[12] */ CLAUSE gmA089858[] = { CLAUSESEQ_begin(4,2), { 2, 0, 1, 1,}, { 1, 0, 0, 1 } }; /* A089840[13] */ CLAUSE gmA073269[] = { CLAUSESEQ_begin(4,2), { 2, 0, 1, 2,}, { 1, 0, 0, 1 } }; /* A089840[14] */ CLAUSE gmA089859[] = { CLAUSESEQ_begin(4,2), { 2, 0, 1, 4,}, { 1, 0, 0, 1 } }; /* A089840[15] */ CLAUSE gmA089860[] = { CLAUSESEQ_begin(4,2), { 2, 0, 1, 5,}, { 1, 0, 0, 1 } }; /* A089840[16] */ CLAUSE gmA074680[] = { CLAUSESEQ_begin(4,2), { 2, 1, 0, 0 }, { 1, 0, 0, 1 } }; /* A089840[17] */ CLAUSE gmA089861[] = { CLAUSESEQ_begin(4,2), { 2, 1, 0, 1,}, { 1, 0, 0, 1 } }; /* A089840[18] */ CLAUSE gmA073270[] = { CLAUSESEQ_begin(4,2), { 2, 1, 0, 2,}, { 1, 0, 0, 1 } }; /* A089840[19] */ CLAUSE gmA089862[] = { CLAUSESEQ_begin(4,2), { 2, 1, 0, 3,}, { 1, 0, 0, 1 } }; /* A089840[20] */ CLAUSE gmA089863[] = { CLAUSESEQ_begin(4,2), { 2, 1, 0, 4,}, { 1, 0, 0, 1 } }; /* A089840[21] */ CLAUSE gmA0joku1[] = { CLAUSESEQ_begin(4,2), { 2, 0, 1, 0,}, { 2, 1, 0, 4 } }; /* A089840[258] */ CLAUSE gmA0joku2[] = { CLAUSESEQ_begin(4,2), { 2, 1, 0, 0,}, { 2, 0, 1, 4 } }; /* A089840[264] */ CLAUSE gmA073281[] = { CLAUSESEQ_begin(115,3), { 3, 1, 1, 7,}, { 2, 0, 1, 2 }, { 2, 1, 0, 2 } }; /* A089840[1654023] */ CLAUSE gmA089864[] = { CLAUSESEQ_begin(115,3), { 3, 2, 2, 7,}, { 2, 0, 0, 2 }, { 2, 1, 1, 1 } }; /* A089840[1654694] */ CLAUSE gmA089864b[] = { CLAUSESEQ_begin(115,3), { 3, 2, 2, 7,}, { 2, 1, 1, 1 }, { 2, 0, 0, 2 } }; /* The last two clauses in either order. */ CLAUSE gmA082351[] = { CLAUSESEQ_begin(24,2), { 3, 2, 4, 0 }, { 2, 1, 1, 5 } }; /* A089840[4069] */ CLAUSE gmA082352[] = { CLAUSESEQ_begin(24,2), { 3, 4, 2, 0 }, { 2, 1, 1, 3 } }; /* A089840[4253] */ CLAUSE gmA082353[] = { CLAUSESEQ_begin(24,2), { 3, 2, 0, 0 }, { 2, 0, 0, 3 } }; /* A089840[3886] A057163-conjugates of the above ones. */ CLAUSE gmA082354[] = { CLAUSESEQ_begin(24,2), { 3, 0, 2, 0 }, { 2, 0, 0, 5 } }; /* A089840[3702] */ struct clauses_and_Anum { CLAUSESEQ clauses; int Anum; }; struct clauses_and_Anum Submitted_nr_gatomorphisms[] = { { gmA001477, 1477 }, { gmA069770, 69770 }, { gmA072796, 72796 }, { gmA089850, 89850 }, { gmA089851, 89851 }, { gmA089852, 89852 }, { gmA089853, 89853 }, { gmA089854, 89854 }, { gmA072797, 72797 }, { gmA089855, 89855 }, { gmA089856, 89856 }, { gmA089857, 89857 }, { gmA074679, 74679 }, { gmA089858, 89858 }, { gmA073269, 73269 }, { gmA089859, 89859 }, { gmA089860, 89860 }, { gmA074680, 74680 }, { gmA089861, 89861 }, { gmA073270, 73270 }, { gmA089862, 89862 }, { gmA089863, 89863 }, { gmA0joku1, 0 }, { gmA0joku2, 0 }, { gmA082351, 82351 }, { gmA082352, 82352 }, { gmA082353, 82353 }, { gmA082354, 82354 }, { gmA073281, 73281 }, { gmA089864, 89864 }, { gmA089864b, 89864 }, { NULL, 0 } }; void fprint_submitted_nr_gatomorphism(FILE *fp,struct clauses_and_Anum *gmn,CLAUSESEQ *hashtable,int hashsize,ECOTEX *ep) { int Anum = gmn->Anum; CLAUSESEQ clause_seq = gmn->clauses; CLAUSESEQ dupe = NULL; char indbuf[81]; char name[81]; dupe = duplicate_found_in_hashtable(clause_seq,hashtable,hashsize,ep,NULL,1); if(NULL == dupe) { strcpy(indbuf,"NOT-FOUND-WITH-GIVEN-LIMITS"); } else { sprintf(indbuf,"%u",CLAUSESEQ_norerank(dupe)); } sprintf(name,"Non-recursive gatomorphism A089840[%s]",indbuf); output_OEIS_sequence(fp,Anum,ep->vals_obtained,ep->vec_size, VECTYPE_USI,(PF_VEC_OUT)fprint_USI_vector,0,OEIS_SEQ, name ); } void fprint_all_submitted_nr_gatomorphisms(FILE *fp,CLAUSESEQ *hashtable,int hashsize,ECOTEX *ep) { int i=0; struct t_gatom_descr *gato_descr; struct clauses_and_Anum *gmn; while((gmn = &(Submitted_nr_gatomorphisms[i++])) && (NULL != gmn->clauses)) { fprint_submitted_nr_gatomorphism(fp,gmn,hashtable,hashsize,ep); } } int main(int argc, char **argv) { ULLI x = (ULLI) (((SLLI) -1)); char *progname = *argv; char *tharg; int upto_binexp_n,upto_size_n; char *count_or_check = NULL; ECOTEX *ep; CLAUSESEQ *hashtable; int hashsize,misc_table_size=15,misc_seq_size=120,max_order_tested_to=21; CLAUSESEQ *first_clauseseqs; ULLI *A089832,*A089833; SGRANK *A089840 = NULL; /* An array of signature-permutations. */ NORERANK *A089839 = NULL; /* A table of their compositions. */ NORERANK *misc_seq_vec = NULL; FILTERMASK filmask = FILTER_IF_NONULTIMATE_1S+FILTER_LONE_NONSTANDALONES; int opt_print_all_submitted_nr_gatomorphisms = 0; int i=0; CheckTriangle(19); /* First handle the options: */ while((tharg = *++argv) && ('-' == *tharg)) { char c; argc--; while(c=*++tharg) { switch(c) { case 'G': { opt_print_all_submitted_nr_gatomorphisms++; break; } case 'S': { char *narg = ++tharg; /* Pointing to character after S */ if('\0' == *narg) { narg = *++argv; /* The next argument, hopefully not NULL. */ argc--; } if((NULL == narg) || ('\0' == *narg) || !isdigit(*narg) || ('0' == *narg)) { fprintf(stderr,"Option -S requires a numeric, non-zero argument. Example: -S 240\n"); exit(1); } misc_seq_size = atoi(narg); goto big_wheels_keep_on_rolling; break; } case 'T': { char *narg = ++tharg; /* Pointing to character after S */ if('\0' == *narg) { narg = *++argv; /* The next argument, hopefully not NULL. */ argc--; } if((NULL == narg) || ('\0' == *narg) || !isdigit(*narg) || ('0' == *narg)) { fprintf(stderr,"Option -T requires a numeric, non-zero argument. Example: -T 45\n"); exit(1); } misc_table_size = atoi(narg); goto big_wheels_keep_on_rolling; break; } case 'm': { char *narg = ++tharg; /* Pointing to character after S */ if('\0' == *narg) { narg = *++argv; /* The next argument, hopefully not NULL. */ argc--; } if((NULL == narg) || ('\0' == *narg) || !isdigit(*narg) || ('0' == *narg)) { fprintf(stderr,"Option -m requires a numeric, non-zero argument. Example: -T 61\n"); exit(1); } max_order_tested_to = atoi(narg); goto big_wheels_keep_on_rolling; break; } case 'l': { globopt_output_cycle_lists++; break; } case 'H': { globopt_HTML = 1; break; } case 'M': /* Maple-output and Haskell/Prolog-output. */ { globopt_list_begin = "["; globopt_list_delim = ","; globopt_list_end = "]"; globopt_elem_delim = ","; globopt_pair_begin = "["; globopt_pair_delim = ","; globopt_pair_end = "]"; globopt_comment_begin = "#"; break; } case 'R': { CheckRankings(upto_size_n); CheckSexpRankings(upto_size_n); exit(0); break; } default: { fprintf(stderr,"Unknown option %s !\n",tharg); goto usage; } } } big_wheels_keep_on_rolling: ; } if(argc < 4) { usage: fprintf(stderr,"usage: %s upto_binexp_n upto_size_n hashtable_size\n", progname); exit(1); } upto_binexp_n = atoi(*(argv+0)); if(upto_binexp_n > MAX_BINEXP) { fprintf(stderr, "%s: upto_binexp_n (%d) was specified greater than MAX_BINEXP=%u\n", progname,upto_binexp_n,MAX_BINEXP ); exit(1); } upto_size_n = atoi(*(argv+1)); if((upto_size_n < binwidth(upto_binexp_n)) || (upto_size_n > MAX_COMPUTATION_SIZE)) { fprintf(stderr, "%s: upto_size_n (%d) was specified either less than binwidth(%u) (= %u, the maximum clause tree size), or more than MAX_COMPUTATION_SIZE=%u\n", progname,upto_size_n,upto_binexp_n,binwidth(upto_binexp_n),MAX_COMPUTATION_SIZE ); exit(1); } hashsize = atoi(*(argv+2)); if(NULL == (hashtable = ((CLAUSESEQ *) calloc(hashsize,sizeof(CLAUSESEQ))))) { fprintf(stderr, "%s: Couldn't allocate a hash table of %u buckets (%u bytes).\n", progname,hashsize,(hashsize*sizeof(CLAUSESEQ)) ); exit(1); } ep = new_ecotex(upto_size_n); printf("glob_cons_cells_in_use=%lu\n",glob_cons_cells_in_use); { int i = (binwidth(upto_binexp_n)+1); if(NULL == (A089832 = ((ULLI *) calloc(i,sizeof(ULLI))))) { fprintf(stderr, "%s: Couldn't allocate a A089832 table of %u longlong integers (%u bytes).\n", progname,(i),(i*sizeof(ULLI)) ); exit(1); } i = tr2seqind(i,i)+1; if(NULL == (A089833 = ((ULLI *) calloc(i,sizeof(ULLI))))) { fprintf(stderr, "%s: Couldn't allocate a A089833 table of %u longlong integers (%u bytes).\n", progname,(i),(i*sizeof(ULLI)) ); exit(1); } } { int i = tr2seqind(misc_table_size+1,misc_table_size+1)+1; if((misc_table_size < 1) || (misc_table_size > ep->vec_size)) { fprintf(stderr, "%s: misc_table_size (%d) was specified either less than 1, or more than A014137(%u)=%lu\n", progname,misc_table_size,upto_size_n,A014137(upto_size_n) ); exit(1); } if(NULL == (A089839 = ((NORERANK *) calloc(i,sizeof(NORERANK))))) { fprintf(stderr, "%s: Couldn't allocate a triangular table of %u short integers (%u bytes) for A089839.\n", progname,(i),(i*sizeof(NORERANK)) ); exit(1); } if(NULL == (A089840 = ((SGRANK *) calloc(i,sizeof(SGRANK))))) { fprintf(stderr, "%s: Couldn't allocate a triangular table of %u short integers (%u bytes) for A089840.\n", progname,(i),(i*sizeof(SGRANK)) ); exit(1); } if(misc_seq_size < misc_table_size) { fprintf(stderr,"%s: You can't specify misc_seq_size (%u) < misc_table_size (%u)!\n", progname,misc_seq_size,misc_table_size); exit(1); } if(NULL == (first_clauseseqs = ((CLAUSESEQ *) calloc(misc_seq_size,sizeof(CLAUSESEQ))))) { fprintf(stderr, "%s: Couldn't allocate a first_clauseseqs table of %u clause sequences (%u bytes).\n", progname,misc_seq_size,(misc_seq_size*sizeof(CLAUSESEQ)) ); exit(1); } if(NULL == (misc_seq_vec = ((NORERANK *) calloc(misc_seq_size,sizeof(NORERANK))))) { fprintf(stderr, "%s: Couldn't allocate a misc_seq_vec of %u NORERANKs (%u bytes).\n", progname,misc_seq_size,(misc_seq_size*sizeof(NORERANK)) ); exit(1); } } check_upto_binexp_n(upto_binexp_n,ep,hashsize,hashtable,first_clauseseqs,filmask, misc_seq_size,misc_table_size,max_order_tested_to, A089832,A089833,A089839,A089840,misc_seq_vec); if(opt_print_all_submitted_nr_gatomorphisms) { fprint_all_submitted_nr_gatomorphisms(stdout,hashtable,hashsize,ep); } exit(0); }