/* Kevin Ryde, December 2024 Usage: ./a.out [-v] with compile-time chosen N and steps LE or LRE This is some C code to find terms of A039745(n) = maximum steps LE to any permutation, "diameter" A378881 = permutations at diameter A186144(n) = how many permutations at diameter A048200(n) = number of steps LE to reverse a permutation A061545(n) = number of ways for that reversal or A186783(n) = maximum steps LRE to any permutation, "diameter" A186752(n) = number of steps LRE to reverse a permutation A378834(n) = number of ways for that reversal Steps LE or LRE are L = rotate left 1 place R = rotate right 1 place E = exchange elements at first 2 positions The method is a breadth-first explore through permutations at successive distance. Re-visits are identified by a flags array. Multiple ways to a permutation are combined by sort and add the counts of ways. Expect to peak at maybe 1.5*(n!) bytes of memory. The flags are random access so not well suited to disk based swap. Configuration ------------- The desired permutation size n and steps type are #defines in the Configuration section of the code below. #define N sequence n #define STEPS_TYPE steps LE or LRE Those are compile-time constants to give the compiler maximum opportunity to optimise when attempting n=13 or n=14. Changing the source and re-compiling is inconvenient, but only a handful of small n are feasible anyway. Memory use is kept down by bit fields in struct idx_and_ways. #define IDX_FIELD_BITS 33 #define WAYS_FIELD_BITS 31 These supplied sizes are enough for n <= 13 using 8 bytes per element. Configuration for n=14 ---------------------- To attempt n=14, try IDX_FIELD_BITS 37 WAYS_FIELD_BITS 43 The idx field must hold numbers 0 .. 14!-1 so 37 bits needed. The ways field is a guess at how big the largest number of ways might be. (Maybe 35 bits is enough?) GCC "__attribute__ ((packed))" keeps the struct to only as many bytes as needed. This reduces peak memory use. The two fields may as well add up to a multiple of 8 bits so as to use all of the bytes in the struct. Output ------ Output is a free form report with for instance, and among other information A039745(10) = 58 LE diameter A186144(10) = 1 LE num permutations at diameter: A048200(10) = 55 LE reversal A061545(10) = 512 LE reversal ways Command-line parameter "-v" prints some more verbose progress messages (which distance d is being considered). Permutations ------------ A permutation represented as a byte array with 0 based elements uint8_t perm[N] elements 0..N-1 Each permutation has an index number "idx", idx_t idx = perm_and_inv_to_idx (perm,inv) range 0 .. n!-1 inclusive The indexing of permutations is chosen so E exchange is simply "XOR 1" (flip the least significant bit). Otherwise it only matters that every permutation has a distinct index. See "Permutation Indexing" below. Data Types ---------- idx_t holds an index outside the struct and is uint64_t to allow for n=13 with 13! using 33 bits, as noted above. ways_t holds a ways count outside the struct and is uint64_t anticipating a 64 bit build and that more than 32 bits might be needed at n=14. These and their corresponding PRIidx and PRIways printf types can be changed if desired. A 32 bit build works for n <= 12. At n=13 the flags[] array is fine at 742 mbytes, but the idx_list[] malloc() peaks around 5 gbytes which exceeds 32 bit addressing. Incidentally, since flags[] is a large static array (BSS), expect an immediate segfault from the loader, with no output, if not enough memory (plus swap space) for it. malloc()s later will fail with messages. Calculation ----------- Within calculate(), idx_list[] = array of idx_and_ways at distance d num = how many of those idx_list_expand() takes each idx and "expands" by applying L,E or L,R,E steps, but only keep those not before in flats[]. This expand is in-place with re-writing early entries after they've been read, and appending after the end if necessary. Each "old" permutation makes at most 2 new. LE = both can be new and never seen before LRE = For distance d>=1, one of L,R,E is the inverse of the step which lead to the present permutation, so was seen before, so up to 2 new For distance d=0, the identity permutation had no previous step, so 3 new there. idx_list[] is realloc()ed up to 2*num+3 for this reason and for the scheme in idx_list_expand() of unconditional stores to both re-write and append locations, so both must be in the alloc. idx_list_sort() uses qsort() in-place to bring copies of the same idx adjacent. Copies occur when there are different ways to reach idx. idx_list_uniq() merges the "ways" count of adjacent same idx, and collapses to one of each idx. It also marks all those in flags[] as having been seen at what will now be a previous distance. Things Not Done --------------- This form is after attempts at more or less RAM, more or less disk, and sorting or not. Sorting to merge "ways" is a significant cost. It's only needed to sum ways from multiple occurrences of the same idx. If the ways were unwanted then expand can immediately flag newly seen permutations, leaving new idx_list all distinct (and in no particular order). The alternative to idx is to pack the elements of a permutation in nibbles of a 64 bit word (or similar). Keeping sorted lists of all which occur at a given distance allows a sort and "multi-scan" to eliminate duplicates. The passes over that data are sequential, but are also large, so impractical at n=13 for a typical mechanical hard drive. An alternative to flags[] is to record bytes distance[idx]. That eliminates any idx_list since a pass over distance[] can find idx of a given d. But this doesn't have "ways" counts and uses more memory than "flags without ways". Both diameter and reversal distances look like quadratics, through to n=13. If that could be proven then searching for them here would be unnecessary. GNU C library qsort(), around the time of writing, chooses between quick sort (in-place) and merge sort (use extra memory) depending on what physical RAM sysconf() reports, or something like that. If the large sorts here somehow use more RAM than you want then might try the undocumented _quicksort() which is always quick sort (in qsort_r() style). */ /* ------------------------------------------------------------------------ */ /* Configuration */ #define N 10 /* sequence n */ /* STEPS_TYPE_LE for A039745, A048200, etc. STEPS_TYPE_LRE for A186783, A186752, etc. */ #define STEPS_TYPE STEPS_TYPE_LE /* Set WANT_ASSERT to 1 to enable self-checks (slower). Set WANT_DEBUG to 1 for some rough development prints. */ #define WANT_ASSERT 0 #define WANT_DEBUG 0 /* ------------------------------------------------------------------------ */ #if ! WANT_ASSERT #define NDEBUG #endif #include #include #include #include #include #include #include #include #include #include #include #include /* ------------------------------------------------------------------------ */ /* Generic */ #if WANT_DEBUG #define DEBUG(expr) do { expr; } while (0) #else #define DEBUG(expr) do { } while (0) #endif #ifdef __GNUC__ #define LIKELY(cond) __builtin_expect((cond) != 0, 1) #define UNLIKELY(cond) __builtin_expect((cond) != 0, 0) #define ATTRIBUTE_PRINTF __attribute__ ((format (printf,1,2))) #define ATTRIBUTE_PACKED __attribute__ ((packed)) #else #define LIKELY(cond) (cond) #define UNLIKELY(cond) (cond) #define ATTRIBUTE_PRINTF #define ATTRIBUTE_PACKED #endif #define numberof(array) (sizeof(array)/sizeof((array)[0])) #define SIZE_T_MAX (~ (size_t) 0) #define MAX(x,y) ((x)>=(y) ? (x) : (y)) #define MIN(x,y) ((x)<=(y) ? (x) : (y)) /* evaluate to ceil(n/d), with no risk of overflow */ #define DIV_CEIL(n,d) (((n)/(d)) + (((n)%(d)) != 0)) noreturn ATTRIBUTE_PRINTF void error (const char *format, ...) { va_list ap; va_start (ap, format); vfprintf (stderr, format, ap); va_end (ap); exit(1); } static void * malloc_or_die(size_t size) { void *ptr = malloc(size); if (UNLIKELY (ptr == NULL)) error("Cannot malloc() %zu bytes\n", size); return (ptr); } /* with ptr=NULL allowed for initial malloc() too */ static void * realloc_or_die(void *ptr, size_t size) { if (size==0) { if (ptr) free(ptr); return (NULL); } else { ptr = (ptr==NULL ? malloc(size) : realloc(ptr, size)); if (UNLIKELY (ptr == NULL)) error("Cannot realloc() to %zu bytes\n", size); return (ptr); } } /* Return x+y or die if that would overflow size_t. */ static size_t size_add_or_die (size_t x, size_t y) { if (UNLIKELY (x > SIZE_T_MAX - y)) error("Add %zu + %zu overflows size_t (of %zu bytes)\n", x, y, sizeof(size_t)); return (x+y); } /* Return x*y or die if that would overflow size_t. */ static size_t size_multiply_or_die (size_t x, size_t y) { if (UNLIKELY (x > SIZE_T_MAX / y)) error("Multiply %zu * %zu overflows size_t (of %zu bytes)\n", x, y, sizeof(size_t)); return (x*y); } /* num many elements of size bytes each */ void * malloc_num_or_die(size_t num, size_t size) { return (malloc_or_die(size_multiply_or_die(num,size))); } /* realloc to num many elements of size bytes each */ void * realloc_num_or_die(void *ptr, size_t num, size_t size) { return (realloc_or_die(ptr, size_multiply_or_die(num,size))); } static void ferror_die (FILE *fp, const char *filename, const char *action) { if (ferror(fp)) error("Error %s %s: %s\n", action, filename, strerror(errno)); } /* CPU time consumed by this process so far, in seconds */ static double cputime (void) { struct timespec t; if (clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &t) != 0) error("Cannot clock_gettime() CPUTIME: %s\n", strerror(errno)); return (t.tv_sec + t.tv_nsec/1e9); } /* ------------------------------------------------------------------------ */ /* for(N=1,15, print("#elif N==",N,"\n#define N_FACTORIAL ",N!)); */ #if N==1 #define N_FACTORIAL 1 #elif N==2 #define N_FACTORIAL 2 #elif N==3 #define N_FACTORIAL 6 #elif N==4 #define N_FACTORIAL 24 #elif N==5 #define N_FACTORIAL 120 #elif N==6 #define N_FACTORIAL 720 #elif N==7 #define N_FACTORIAL 5040 #elif N==8 #define N_FACTORIAL 40320 #elif N==9 #define N_FACTORIAL 362880 #elif N==10 #define N_FACTORIAL 3628800 #elif N==11 #define N_FACTORIAL 39916800 #elif N==12 #define N_FACTORIAL 479001600 #elif N==13 #define N_FACTORIAL 6227020800 #elif N==14 #define N_FACTORIAL 87178291200 #elif N==15 #define N_FACTORIAL 1307674368000 #else #error Unrecognised N #endif typedef uint64_t idx_t; #define PRIidx PRIu64 #define IDX_T_MAX ((idx_t) (~ (idx_t) 0)) typedef uint64_t ways_t; #define PRIways PRIu64 #define WAYS_T_MAX ((ways_t) (~ (ways_t) 0)) #define WAYS_T_BITS ((int) (CHAR_BIT * sizeof(ways_t))) /* With "packed", GCC allows the bit fields to be any size and fits them in as many bytes as necessary (doesn't have to be a multiple of uint64_t). May as well have the fields sum to a multiple of 8 bits so use the whole of the bytes. */ struct idx_and_ways { #define IDX_FIELD_BITS 33 #define WAYS_FIELD_BITS 31 uint64_t idx : IDX_FIELD_BITS ATTRIBUTE_PACKED; uint64_t ways : WAYS_FIELD_BITS ATTRIBUTE_PACKED; }; #define WAYS_FIELD_MAX \ ((ways_t) (UINT64_MAX >> (64-WAYS_FIELD_BITS))) static uint8_t flags[DIV_CEIL(N_FACTORIAL,8)]; #define MAX_DISTANCE 300 static idx_t num_at_d[MAX_DISTANCE + 1]; int reversal_distance = -1; ways_t reversal_ways = 0; struct idx_and_ways *antipodes_idx_list; size_t antipodes_num; static uint8_t perm_identity[N]; /* 0..N-1 */ static uint8_t perm_reversal[N]; /* N-1..0 */ static uint8_t perm_decrement[N]; /* N-1,0,1,...,N-2 is -1 mod N */ static uint8_t perm_increment[N]; /* 1,2,...,N-1,0 is +1 mod N */ idx_t reversal_idx; /* dst->ways += add, with overflow check. */ static inline ways_t ways_addto (struct idx_and_ways *dst, ways_t add) { assert (add <= WAYS_FIELD_MAX); if (sizeof(ways_t) > sizeof(uint64_t)) error("oops, ways_t bigger than uint64_t\n"); ways_t sum = dst->ways + add; /* If ways field is full ways_t then sum < add means wrap-around, and gcc recognises as check carry flag. Otherwise smaller than ways_t so sum is always good and just check it fits the bit field stored. */ if (WAYS_FIELD_MAX == WAYS_T_MAX ? UNLIKELY(sum < add) : UNLIKELY(sum > WAYS_FIELD_MAX)) error("ways overflow\n"); dst->ways = sum; return sum; } /* ------------------------------------------------------------------------ */ /* Permutations perm[N] is a byte array permutation of 0..N-1 idx is numbering 0..N!-1 for a permutation. */ void perm_print (const uint8_t *perm) { int i; for (i=0; i < N; i++) printf("%d%s", perm[i], i==N-1 ? "" : ","); } /* Printout 1-based (rather than 0-based). */ void perm_print1 (const uint8_t *perm) { int i; for (i=0; i < N; i++) printf("%d%s", perm[i]+1, i==N-1 ? "" : ","); } int check_perm (const uint8_t *perm) { /* DEBUG(printf("check_perm() "); perm_print(perm); printf("\n")); */ int i; uint8_t seen[N]; memset(seen,'\0',sizeof(seen)); for (i=0; i < N; i++) { int t = perm[i]; if (! (0 <= t && t < N)) error("oops, perm[] i=%d bad term t=%d\n", i,t); if (seen[t]) error("oops, perm[] i=%d already seen term t=%d\n", i,t); seen[t] = 1; } for (i=0; i < N; i++) assert (seen[i]); /* all of 0..N-1 seen */ return (1); } static inline void perm_copy (uint8_t *dst, const uint8_t * restrict src) { memcpy (dst, src, N*sizeof(uint8_t)); } /* Set inv[] to inverse permutation of perm[]. inv[] and perm[] arrays must not overlap. */ static inline void perm_inverse (uint8_t *inv, const uint8_t * restrict perm) { assert (inv >= perm+N || perm >= inv+N); /* no overlap */ assert (check_perm(perm)); int i; for (i=0; i < N; i++) { assert (0 <= perm[i] && perm[i] < N); inv[perm[i]] = i; } } /* inv[] is the inverse of a permutation. inv_rotate_left() sets dst to the corresponding inverse for that permutation with its elements rotated 1 place to the left. inv_rotate_right() similarly, rotated 1 place to the right. Left means each inv[i] changed to inv[i]-1. Right means each inv[i] changed to inv[i]+1. In both cases with wrap-around, so mod N, elements 0..N-1. */ static inline void inv_rotate_left (uint8_t *dst, const uint8_t * restrict inv) { int i; for (i=0; i < N; i++) dst[i] = perm_decrement[inv[i]]; } static inline void inv_rotate_right (uint8_t *dst, const uint8_t * restrict inv) { int i; for (i=0; i < N; i++) dst[i] = perm_increment[inv[i]]; } /* ------------------------------------------------------------------------ */ /* Permutation Indexing Permutations are indexed by "idx" using swaps selected by the factorial base representation of idx, idx = digits r[1..N-1] in factorial base so idx = Sum r[i] * (i+1)! with 0 <= r[i] <= i i=1..N-1 most least significant significant radix N N-1 3 2 digit r[N-1], r[N-2], ..., r[2], r[1] This is r[i] = A301652(idx,i+1), taking terms as 0 if i+1 is beyond the end of row idx (so high 0 digits). inv_from_idx() calculates r[i] digits by successive divisions. N is a compile-time constant so the compiler can unroll the loop and use its optimisations for division by constant. GNU C does this (under -O3 optimisation level), and it might even inter-leave the 2 independent steps L,R for LRE (though hard to tell under much inlining!). Each r[i] says where to swap permutation entry perm[i] to, start perm[i] = i identity permutation for i = N-1 to 1 swap perm[i] <-> perm[r[i]] Notice these swaps are by decreasing i. This ensures that the final swap is the first two elements, or not. swap perm[1] <-> perm[0] if r[1] = 0 swap swap perm[1] <-> perm[1] if r[1] = 1 no-op The E operation exchanging perm[0] <-> perm[1] is as easy as E(idx) = idx ^ 1 XOR An inv[] inverse of the idx permutation is calculated with the same r[i] digits but i run in the opposite order, inv_from_idx(): start inv[i] = i identity permutation for i = 1 to N-1 swap inv[i] <-> inv[r[i]] The loop always has inv[i] = i since it's the first time inv[i] has been touched. inv_from_idx() uses that to implement the "swap" by two stores. idx_list_expand() works with inv[] inverses since inv_from_idx() is a tiny bit less work than a direct "perm_from_idx()" would be. The converse operation, permutation to idx, uses perm[] and inv[] together to find each r[i]. The final swap in inv_to_idx() is at i = N-1 and moves element i to position r[i]. So r[i] is recovered by finding where is i. perm[i] is exactly that, inv[r[i]] = i perm[i] = r[i] After undoing this final swap of i=N-1, in both inv[] and perm[], the same way to find r[i] holds at i = N-2, and so on. perm_and_inv_to_idx(): for i = N-1 to 1 r[i] = perm[i] swap inv[i] <-> inv[r[i]] (undo) corresponding update perm[] The undo swap sets inv[i] = i and corresponding perm[i] = i. perm_and_inv_to_idx() doesn't bother to store those since nothing later uses them. An attraction of this indexing is that swaps mean the to and from operations are a little less work than the more common lexicographic ordering (or "co-lexicographic"). */ static inline void inv_from_idx (uint8_t *inv, idx_t idx) { assert (idx < N_FACTORIAL); assert (memset(inv,'\xAA',N)); int i; inv[0] = 0; for (i=1; i < N; i++) { int m = i+1; size_t r = idx % m; idx /= m; int t = inv[r]; inv[i] = t; inv[r] = i; } assert (check_perm(inv)); } static void perm_from_idx (uint8_t *perm, idx_t idx) { assert (idx < N_FACTORIAL); uint8_t inv[N]; inv_from_idx (inv, idx); perm_inverse (perm, inv); } /* Return the index idx ranging 0 <= idx < N! of permutation perm[] and whose inverse is inv[]. Destroys the contents of perm[] and inv[]. */ static inline size_t perm_and_inv_to_idx (uint8_t *perm, uint8_t *inv) { idx_t idx = 0; int i; for (i=N-1; i>=1; i--) { int r = perm[i]; int t = inv[i]; inv[r] = t; perm[t] = r; idx = (i+1)*idx + r; } return (idx); } /* Return the index idx ranging 0 <= idx < N! of the permutation whose inverse is the given inv[]. Destroys the contents of inv[]. */ static inline size_t inv_to_idx (uint8_t *inv) { uint8_t perm[N]; perm_inverse (perm, inv); return (perm_and_inv_to_idx (perm,inv)); } /* without modifying perm[] */ static size_t perm_const_to_idx (const uint8_t *perm) { uint8_t copy[N],inv[N]; perm_copy (copy,perm); perm_inverse(inv,perm); return (perm_and_inv_to_idx (copy,inv)); } static void init_perms (void) { int i; for (i=0; i < N; i++) { perm_identity[i] = i; perm_reversal[i] = N-1 - i; perm_decrement[i] = (i+N-1) % N; perm_increment[i] = (i+1) % N; } reversal_idx = perm_const_to_idx(perm_reversal); } /* Showing permutation elements as 0..N-1. */ void perm_print_idx (idx_t idx) { uint8_t perm[N]; perm_from_idx (perm, idx); perm_print (perm); } /* Showing permutation elements as 1..N. */ void perm_print1_idx (idx_t idx) { uint8_t perm[N]; perm_from_idx(perm,idx); perm_print1 (perm); } /* Compare xp and yp by their ->idx fields numerically. No particular significance to that order, it only has to get equal idx values consecutive. */ static int idx_and_ways_cmp (const void *xp, const void *yp) { idx_t x_idx = ((const struct idx_and_ways *) xp)->idx; idx_t y_idx = ((const struct idx_and_ways *) yp)->idx; return (x_idx > y_idx ? 1 : x_idx == y_idx ? 0 : -1); } /* Compare xp an yp by their permutations in lexicographic order. This is only used on the few antipodal permutations so doesn't need to be efficient. */ static int idx_and_ways_cmp_by_perm (const void *xp, const void *yp) { idx_t x_idx = ((const struct idx_and_ways *) xp)->idx; idx_t y_idx = ((const struct idx_and_ways *) yp)->idx; uint8_t x_perm[N], y_perm[N]; perm_from_idx (x_perm, x_idx); perm_from_idx (y_perm, y_idx); return (memcmp (x_perm, y_perm, sizeof(x_perm))); } /* Compare xp an yp by their permutations in lexicographic order. This is only used on the few maximum ways permutations so doesn't need to be efficient. */ static int idx_cmp_by_perm (const void *xp, const void *yp) { idx_t x_idx = * (const idx_t *) xp; idx_t y_idx = * (const idx_t *) yp; uint8_t x_perm[N], y_perm[N]; perm_from_idx (x_perm, x_idx); perm_from_idx (y_perm, y_idx); return (memcmp (x_perm, y_perm, sizeof(x_perm))); } /* ------------------------------------------------------------------------ */ static int option_verbose = 0; static size_t peak_malloc = 0; #define STEPS_TYPE_LE 0 #define STEPS_TYPE_LRE 1 #define NUM_STEPS_TYPES 2 #define STEPS_WANT_RIGHT (STEPS_TYPE == STEPS_TYPE_LRE) #define NUM_STEP_ALTERNATIVES (2 + STEPS_WANT_RIGHT) const char *steps_type_name[NUM_STEPS_TYPES] = { "LE", "LRE" }; #define VALUE_TYPE_DIAMETER 0 #define VALUE_TYPE_NUM_AT_DIAMETER 1 #define VALUE_TYPE_REVERSAL_DISTANCE 2 #define VALUE_TYPE_REVERSAL_WAYS 3 #define NUM_VALUE_TYPES 4 /* info[STEPS_TYPE][VALUE_TYPE_...] */ static const struct info { const char *anum; const char *name; uint64_t want[16]; } info_array[NUM_STEPS_TYPES][NUM_VALUE_TYPES] = { { /* STEPS_TYPE_LE */ { "A039745", "LE diameter", { -1, /* no n=0 */ 0, 1, 2, 6, 11, 18, 25, 35, 45, 58, 71, /* n=11 */ 87, /* n=12 */ 103 /* n=13, by Ben Whitmore */ } }, { "A186144", "LE num permutations at diameter", { -1, /* no n=0 */ 1, 1, 3, 3, 2, 1, 2, 1, 2, 1, 2, 1, 2 /* n=13 */ } }, { "A048200", "LE reversal distance", { -1, /* no n=0 */ 0, 1, 2, 4, 10, 15, 23, 32, 42, 55, 67, 84, /* n=12 by Sean A. Irvine */ 98 /* n=13 */ } }, { "A061545", "LE reversal ways", { -1, /* no n=0 */ 0, /* n=1 */ 2, /* n=2 in 2 ways L or E */ 1, 1, 1, 2, 4, 16, 40, 512, 1472, 65536, /* n=11,12 */ 208896 /* n=13 */ } } }, { /* STEPS_TYPE_LRE */ { "A186783", "LRE diameter", { -1, /* no n=0 */ 0, 1, 2, 6, 10, 15, 21, 28, 36, 45, 55, /* n=11 */ 66, /* n=12 */ 78 /* n=13 by Ben Whitmore */ } }, { "", "LRE num permutations at diameter", { -1, /* no n=0 */ 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 /* n=13 */ } }, { "A186752", "LRE reversal distance", { -1, /* no n=0 */ 0, 1, 2, 4, 8, 13, 19, 26, 34, 43, 53, 64, /* n=12 */ 76 /* n=13 */ } }, { "A378834", "LRE reversal ways", { -1, /* no n=0 */ 1, 3, 2, 2, 2, 6, 6, 32, 32, /* n=8,9 */ 240, 240, /* n=10,11 */ 2304, 2304 /* n=12,13 */ } }, } }; /* want[] padded with 0s, find last non-0 for length */ int samples_length (const struct info *info) { int i; for (i = numberof(info->want)-1; i>=0; i--) if (info->want[i] != 0) return (i+1); return (0); } static void check_value (int n, int value, int value_type) { assert (0 <= STEPS_TYPE && STEPS_TYPE < NUM_STEPS_TYPES); assert (0 <= value_type && value_type < NUM_VALUE_TYPES); if (! (n >= 1)) error("oops, check_value() is for n>=1\n"); const struct info *info = &info_array[STEPS_TYPE][value_type]; int length = samples_length(info); if (n >= length) { if (option_verbose) printf("(%s %s beyond sample data)\n", info->anum, info->name); return; } if (value == info->want[n]) { printf("(%s %s good vs sample data)\n", info->anum, info->name); return; } error("oops, %s %s n=%d, got %d want %"PRIu64"\n", info->anum, info->name, n, value, info->want[n]); } /* ------------------------------------------------------------------------ */ /* maximum number of ways to reach any permutation */ static ways_t ways_max; static idx_t ways_max_num; #define WAYS_MAX_LIST_SIZE 32 static idx_t ways_max_list[WAYS_MAX_LIST_SIZE]; static inline void ways_max_accumulate (idx_t idx, ways_t ways) { if (UNLIKELY (ways >= ways_max)) { if (ways == ways_max) { ways_max_num++; } else { ways_max = ways; ways_max_num = 1; } if (ways_max_num <= numberof(ways_max_list)) ways_max_list[ways_max_num-1] = idx; } } static inline uint8_t flag_get (idx_t idx) { if (idx >= N_FACTORIAL) printf("idx %"PRIidx"\n",idx); assert (idx < N_FACTORIAL); size_t wpos = idx / 8; size_t bpos = idx % 8; assert (wpos < numberof(flags)); return ((flags[wpos] >> bpos) & 1); } static inline void flag_set (idx_t idx) { assert (idx < N_FACTORIAL); size_t wpos = idx / 8; size_t bpos = idx % 8; assert (wpos < numberof(flags)); flags[wpos] |= ((uint8_t)1 << bpos); DEBUG(uint8_t inv[N]; uint8_t perm[N]; inv_from_idx(inv,idx); perm_inverse(perm,inv); printf(" flag_set() %"PRIidx" = ",idx); perm_print(perm); printf("\n")); } void show_configuration (void) { static int done = 0; if (done) return; done = 1; if (WANT_ASSERT || option_verbose ) printf("assert()s %s\n", WANT_ASSERT ? "enabled" : "disabled"); if (option_verbose) { done = 1; #ifdef __GNUC__ printf("Have __GNUC__\n"); #else printf("Not __GNUC__\n"); #endif printf("idx_t %zu bytes = %zu bits\n", sizeof(idx_t), sizeof(idx_t)*CHAR_BIT); printf("ways_t %zu bytes = %zu bits\n", sizeof(ways_t), sizeof(ways_t)*CHAR_BIT); printf("struct idx_and_ways %zu bytes = %zu bits\n", sizeof(struct idx_and_ways), sizeof(struct idx_and_ways)*CHAR_BIT); printf(" idx field %d bits\n", IDX_FIELD_BITS); printf(" ways field %d bits\n", WAYS_FIELD_BITS); int field_bits_total = IDX_FIELD_BITS + WAYS_FIELD_BITS; printf(" total %d bits is %s full struct\n", field_bits_total, field_bits_total == sizeof(struct idx_and_ways)*CHAR_BIT ? "=" : "<"); printf("flags[] %zu bytes = %.3lf gbytes\n", sizeof(flags), sizeof(flags) / 1024.0 / 1024.0 / 1024.0); printf("\n"); } idx_t idx_t_max = (~ (idx_t) 0); if (! (idx_t_max >= N_FACTORIAL - 1)) error ("oops, idx_t too small for N=%d\n", N); } /* Total number of permutations seen at distances 0..d inclusive. */ static idx_t total_num (int d) { idx_t total = 0; int i; for (i=0; i <= d; i++) total += num_at_d[i]; return (total); } /* Total number of permutations seen at distances 0..d inclusive as a percentage of N_FACTORIAL = number of all permutations. */ static double total_num_percent (int d) { idx_t total = total_num (d); double percent = total * 100.0 / N_FACTORIAL; /* don't let something < 100% show as 100% */ if (total != N_FACTORIAL) percent = MIN(percent, 99.999); return (percent); } void show_idx_list (struct idx_and_ways *idx_list, size_t num) { size_t i; for (i=0; i < num; i++) { idx_t idx = idx_list[i].idx; ways_t ways = idx_list[i].ways; printf(" idx=%"PRIidx" perm ",idx); perm_print_idx(idx); printf(" ways=%"PRIways"\n", ways); } } /* Expand elements idx_list[0..num-1] in-place. Return new num for new idx_list[0..num-1]. */ static inline size_t idx_list_expand (int d, struct idx_and_ways *idx_list, size_t num, size_t alloc) { if (option_verbose) printf(" expand\n"); /* Elements are taken from successive *from_ptr. write_ptr is an unused location at beginning of idx_list. It always trails 1 or more places behind from_ptr. append_ptr is an unused location after the original idx_list. . +----------+----------+----------+----------+ . | | | | | . +----------+----------+----------+----------+ . write_ptr from_ptr ... append_ptr A new element, with to_idx = one of L,R,E, is stored unconditionally to *write_ptr and *append_ptr since both are unused locations. If to_idx was seen before, per flags_get(), then add=0 and neither write_ptr nor append_ptr advances. If to_idx not seen before, then add=1 and one of write_ptr or append_ptr advances. If new_write_ptr reaches from_ptr under the advance, then that means it won't point at an unused location, so advance append_ptr instead. The worst case is . . d=0 3 new writes, all to append_ptr . so alloc >= 3 = num + 3 = 2*num+2 suffices . (original num=1 from d=0) . . d>=1 2 new per old, . first 2 always to append_ptr . then 1 to write_ptr, 1 to append_ptr . so alloc >= num + num+1 = 2*num+1 suffices In both cases further +1 so append_ptr destination is always an available unused place. */ struct idx_and_ways *from_ptr = idx_list; struct idx_and_ways *end_ptr = idx_list + num; struct idx_and_ways *write_ptr = idx_list; struct idx_and_ways *append_ptr = end_ptr; while (LIKELY (from_ptr != end_ptr)) { /* struct will get new idx, but each ways unchanged */ struct idx_and_ways elem = *from_ptr++; idx_t idx = elem.idx; DEBUG(printf(" from pos=%td idx=%"PRIidx" perm ", from_ptr-idx_list, idx); perm_print_idx(idx); printf(" ways=%"PRIways"\n", (uint64_t) elem.ways)); if (UNLIKELY(idx == reversal_idx)) { assert (reversal_distance == -1); /* not seen before */ reversal_distance = d; reversal_ways = elem.ways; } uint8_t inv[N]; inv_from_idx(inv,idx); uint8_t left_inv[N]; uint8_t right_inv[N]; inv_rotate_left (left_inv, inv); inv_rotate_right (right_inv, inv); size_t exchange_idx = idx ^ 1; size_t left_idx = inv_to_idx(left_inv); size_t right_idx = inv_to_idx(right_inv); #define ACCUMULATE(to_idx) \ do { \ DEBUG(printf(" to idx=%"PRIidx" perm ", to_idx); \ perm_print_idx(to_idx); \ if (flag_get(to_idx)) printf(" seen before"); \ printf(" w=%td f=%td a=%td\n", \ write_ptr - idx_list, \ from_ptr - idx_list, \ append_ptr - end_ptr)); \ assert (write_ptr <= from_ptr \ && from_ptr <= end_ptr \ && end_ptr <= append_ptr \ && append_ptr < idx_list + alloc); \ \ elem.idx = to_idx; \ *write_ptr = elem; \ *append_ptr = elem; \ \ size_t add = 1 - flag_get(to_idx); \ struct idx_and_ways *new_write_ptr = write_ptr + add; \ struct idx_and_ways *new_append_ptr = append_ptr + add; \ if (new_write_ptr == from_ptr) \ append_ptr = new_append_ptr; \ else \ write_ptr = new_write_ptr; \ } while (0) ACCUMULATE (exchange_idx); ACCUMULATE (left_idx); if (STEPS_WANT_RIGHT) ACCUMULATE (right_idx); } if (append_ptr >= idx_list + alloc) error ("oops, append_ptr went past alloc\n"); size_t append_num = append_ptr - end_ptr; size_t write_num = write_ptr - idx_list; size_t new_num = write_num + append_num; DEBUG(printf(" final w=%td, f=%td, e=%td, a=%td, new_num %zu\n", write_ptr - idx_list, from_ptr - idx_list, end_ptr - idx_list, append_ptr - end_ptr, new_num)); memmove(write_ptr, end_ptr, append_num*sizeof(write_ptr[0])); return (new_num); } int check_sorted (struct idx_and_ways *idx_list, size_t num) { size_t i; for (i=1; i< num; i++) { idx_t idx = idx_list[i].idx; idx_t prev_idx = idx_list[i-1].idx; if (prev_idx > idx) { printf("bad %p:\n", idx_list); show_idx_list (idx_list, num); error("oops, not sorted i=%zu %"PRIidx" vs %"PRIu64"\n", i, prev_idx, idx); } } return(1); } /* Sort the elements idx_list[0..num-1]. */ static inline void idx_list_sort (struct idx_and_ways *idx_list, size_t num) { if (option_verbose) printf(" sort (%zu entries, %.3lf gbytes)\n", num, num * sizeof(idx_list[0]) / 1024.0 / 1024.0 / 1024.0); qsort (idx_list, num, sizeof(idx_list[0]), idx_and_ways_cmp); assert(check_sorted (idx_list, num)); } /* idx_list[] contains "num" entries with the first at from_ptr. Entries are sorted by idx, so same idx are adjacent. Merge same idx by summing their "ways". Store resulting entries to the start idx_list[0..new_num-1] and return "new_num". */ static inline size_t idx_list_uniq (struct idx_and_ways *idx_list, size_t num) { if (option_verbose) printf(" uniq\n"); size_t i; struct idx_and_ways *from_ptr = idx_list; struct idx_and_ways *to_ptr = idx_list; idx_t prev_idx = IDX_T_MAX; for (i=0; i < num; i++) { idx_t idx = from_ptr->idx; idx_t ways = from_ptr->ways; from_ptr++; if (UNLIKELY (idx == prev_idx)) { assert (to_ptr > idx_list); ways = ways_addto (to_ptr-1, ways); } else { to_ptr->idx = idx; to_ptr->ways = ways; to_ptr++; } ways_max_accumulate (idx, ways); prev_idx = idx; flag_set (idx); } size_t new_num = to_ptr - idx_list; assert (new_num <= num); return (new_num); } /* Calculate various global variables for this N. Return the diameter. */ int calculate (void) { struct idx_and_ways *idx_list; size_t num = 1; idx_list = malloc_num_or_die (num, sizeof(idx_list[0])); { idx_t idx = perm_const_to_idx(perm_identity); idx_list[0].idx = idx; idx_list[0].ways = 1; flag_set(idx); DEBUG(show_idx_list(idx_list,num)); } int d; for (d = 0; num != 0; d++) { if (d >= MAX_DISTANCE) error("oops, distance exceeds MAX_DISTANCE = %d\n", MAX_DISTANCE); num_at_d[d] = num; if (option_verbose) printf("at d=%d, %zu this distance (total %.3lf%%)\n", d, num, total_num_percent(d)); /* 2*num+3 suffices, per notes in idx_list_expand() */ size_t alloc = size_multiply_or_die (num, 2); alloc = size_add_or_die (alloc, 3); idx_list = realloc_num_or_die (idx_list, alloc, sizeof(idx_list[0])); peak_malloc = MAX(peak_malloc, alloc * sizeof(idx_list[0])); size_t new_num = idx_list_expand(d,idx_list,num,alloc); assert (num <= alloc); if (new_num == 0) break; num = new_num; idx_list_sort(idx_list, num); num = idx_list_uniq (idx_list, num); idx_list = realloc_num_or_die (idx_list, num, sizeof(idx_list[0])); } antipodes_idx_list = idx_list; antipodes_num = num; return (d); } void show (void) { if (! (N >= 2)) error("must have N >= 2\n"); printf("N=%d, steps %s, on N! = %"PRIidx" permutations\n", N, steps_type_name[STEPS_TYPE], (idx_t) N_FACTORIAL); double start_time = cputime(); int diameter = calculate(); double end_time = cputime(); if (option_verbose) { printf("\n"); printf("N=%d, steps %s, on N! = %"PRIidx" permutations\n", N, steps_type_name[STEPS_TYPE], (idx_t) N_FACTORIAL); } printf("CPU time %.1lf seconds\n", end_time - start_time); printf(" peak malloc() %zu bytes = %.1lf mbytes = %.2lf gbytes\n", peak_malloc, peak_malloc / 1024.0 / 1024.0, peak_malloc / 1024.0 / 1024.0 / 1024.0); { const struct info *info = &info_array[STEPS_TYPE][VALUE_TYPE_DIAMETER]; printf("%s(%d) = %d %s\n", info->anum, N, diameter, info->name); } { const struct info *info = &info_array[STEPS_TYPE][VALUE_TYPE_REVERSAL_DISTANCE]; printf("%s(%d) = %d %s\n", info->anum, N, reversal_distance, info->name); } { const struct info *info = &info_array[STEPS_TYPE][VALUE_TYPE_REVERSAL_WAYS]; printf("%s(%d) = %"PRIways" %s\n", info->anum, N, reversal_ways, info->name); } { printf("number at distance: "); size_t max_at_distance = 0; int i; for (i=0; i <= diameter; i++) { printf("%"PRIidx"%s", num_at_d[i], i==diameter ? "" : ","); max_at_distance = MAX(max_at_distance,num_at_d[i]); } printf("\n"); printf("maximum %zu at d=",max_at_distance); const char *separator = ""; for (i=0; i <= diameter; i++) { if (num_at_d[i] == max_at_distance) { printf("%s%d", separator,i); separator = ","; } } printf("\n"); } { const struct info *info = &info_array[STEPS_TYPE][VALUE_TYPE_NUM_AT_DIAMETER]; if (info->anum[0]) printf("%s(%d) = ", info->anum, N); printf("%zu %s:\n", antipodes_num, info->name); assert (antipodes_num == num_at_d[diameter]); if (STEPS_TYPE == STEPS_TYPE_LE) printf(" A378881 =\n"); /* list of permutations */ qsort(antipodes_idx_list, antipodes_num, sizeof(antipodes_idx_list[0]), idx_and_ways_cmp_by_perm); size_t i; for (i=0; i < antipodes_num; i++) { idx_t idx = antipodes_idx_list[i].idx; ways_t ways = antipodes_idx_list[i].ways; printf(" "); perm_print1_idx(idx); printf(" in %"PRIways" ways\n", ways); } } { printf("maximum ways %"PRIways", at %"PRIidx" permutations:\n", ways_max, ways_max_num); qsort(ways_max_list, ways_max_num, sizeof(ways_max_list[0]), idx_cmp_by_perm); size_t i; for (i=0; i < MIN(ways_max_num,WAYS_MAX_LIST_SIZE); i++) { printf(" "); perm_print1_idx(ways_max_list[i]); printf("\n"); } if (ways_max_num > WAYS_MAX_LIST_SIZE) printf(" ... and more\n"); } check_value (N, diameter, VALUE_TYPE_DIAMETER); check_value (N, antipodes_num, VALUE_TYPE_NUM_AT_DIAMETER); check_value (N, reversal_distance, VALUE_TYPE_REVERSAL_DISTANCE); check_value (N, reversal_ways, VALUE_TYPE_REVERSAL_WAYS); { idx_t total = total_num(diameter); if (total != N_FACTORIAL) error("oops, total at distance should be %"PRIidx" got %"PRIidx"\n", (idx_t) N_FACTORIAL, total); } ferror_die (stdout, "stdout", "writing"); } int main (int argc, char *argv[]) { DEBUG (option_verbose = 1); setbuf(stdout,NULL); init_perms(); int i; for (i = 1; i < argc; i++) { const char *arg = argv[i]; if (strcmp(arg,"-v")==0) { option_verbose++; } else { error("Unrecognised command line option: %s\n", arg); } } show_configuration(); show (); return(0); }