/* Kevin Ryde, May 2023 Usage: ./a.out [limit] Calculate terms of A362909, which is indices of record highs in A362881, which in turn is a(n+1) as length of the longest arithmetic progression of its terms ending at a(n), with both positions and values in arithmetic progression. Output is b-file style so for example A362909(10)=121 is 10 121 Optional command line parameter "limit" is the largest value to seek (how many values of A362881 to look through), ./a.out 1000000 The default is no limit, or rather the size of an "int" which is effectively no limit since the algorithm isn't fast enough to be feasible through to say 2^31. Saved Data File --------------- Terms of the underlying A362881 are written as bytes to file (first byte is A362881(1)), A362881-bytes.data This is re-read on the next run so as to continue where the last run left off. Command line parameter "restart" can remove this file to start again ./a.out restart Growing file size during a run shows how the search is progressing. File output is flushed every 1 minute (since a whole output buffer takes a while to fill at large n). Implementation Notes -------------------- The strategy is a brute force calculation of the underlying A362881 in a byte array, and watch for the next "target" value (which is the n of A362909). Longest progressions in A362881 are sought by considering each possible "step" between indices ending at n. calculate_bounds() is a slightly subtle calculation of bounds on index and value steps if they're to improve on the longest progression seen so far. They reduce the amount of search and identify the first value (v[n-step]) of progressions which decrease too fast. Speed is not particularly good as n increases due to the number of candidate steps still explicitly checked. */ /* 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 /* ------------------------------------------------------------------------ */ /* Generic Helpers */ #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))) #else #define LIKELY(cond) (cond) #define UNLIKELY(cond) (cond) #define ATTRIBUTE_PRINTF #endif #define numberof(array) (sizeof(array)/sizeof((array)[0])) static noreturn ATTRIBUTE_PRINTF void error (const char *format, ...) { va_list ap; va_start (ap, format); vfprintf (stderr, format, ap); va_end (ap); exit(1); } /* can have ptr==NULL for new malloc() */ void * realloc_or_die(void *ptr, size_t size) { ptr = (ptr==NULL ? malloc(size) : realloc(ptr, size)); if (UNLIKELY (ptr == NULL)) error("Cannot realloc() to %ld bytes\n", (long) size); return (ptr); } FILE * fopen_or_die (const char *filename, const char *mode) { FILE *fp = fopen(filename,mode); if (fp==NULL) error ("Cannot open %s: %s\n", filename, strerror(errno)); return (fp); } void fread_or_die (void *ptr, size_t size, size_t num, FILE *fp, const char *filename) { size_t got = fread(ptr,size,num,fp); if (got < 0) error("Error reading %s: %s\n", filename, strerror(errno)); if (got != num) error("Incomplete read %s, want %ld got %ld\n", filename, (long) num, (long) got); } static inline void fputc_or_die (int c, FILE *fp, const char *filename) { if (fputc(c,fp) == EOF) error("Error writing %s: %s", filename, strerror(errno)); } void fflush_or_die (FILE *fp, const char *filename) { if (fflush(fp) != 0) error("Error writing %s: %s", filename, strerror(errno)); } long ftell_or_die (FILE *fp, const char *filename) { long pos = ftell(fp); if (pos == -1) error("ftell error on %s: %s\n", filename, strerror(errno)); return pos; } void ferror_die (FILE *fp, const char *filename, const char *action) { if (ferror(fp)) error("Error %s %s\n", action, filename); } void sigaction_or_die (int signum, const struct sigaction *act, struct sigaction *oldact) { if (sigaction (signum, act, oldact) != 0) error("Cannot sigaction signum=%d: %s\n", signum, strerror(errno)); } void setitimer_or_die (int which, const struct itimerval *new_value, struct itimerval *old_value) { if (setitimer(which, new_value, old_value) != 0) error("Cannot setitimer %d: %s\n", which, strerror(errno)); } /* ------------------------------------------------------------------------ */ /* Signal Handling */ static volatile int flag_flush = 0; void sig_handler_flush (int signum) { flag_flush = 1; } void init_flush () { struct sigaction s; s.sa_handler = sig_handler_flush; sigfillset(&s.sa_mask); s.sa_flags = SA_RESTART; sigaction_or_die(SIGALRM, &s, NULL); struct itimerval t; t.it_interval.tv_sec = 60; /* 1 minute */ t.it_interval.tv_usec = 0; t.it_value = t.it_interval; setitimer_or_die(ITIMER_REAL, &t, NULL); } /* ------------------------------------------------------------------------ */ /* Expected A362909 Terms (known at the time of writing) */ static const int A362909_data[] = { -1, /* no term A362909(0) */ /* n=1 */ 1, 3, 6, 10, 22, 31, 34, 37, 40, 121, /* n=11 */ 1035, 1078, 1121, /* n=14 */ 7607, 9453, 13667, 13729, /* n=18 */ 50040, 50123, 50206, /* n=21 */ 615964, 616448, 616932, 617416, /* n=25 */ 2828280, 2851232, 2874184, 2897136, /* n=29 */ 5350614, 5583250, 5594326, /* n=32 */ 17310256, 17312210, 17994974, 18008820, /* n=36 */ 19051432, 19069162, 19086892, 19987144, /* n=40 */ 21144970, 21193360, 21740716, 21777460 }; /* check that A362909(target) = n */ void check_A362909 (int target, int n) { assert (target >= 1); if (target < numberof(A362909_data) && A362909_data[target] != n) { error ("oops, target=%d A362909_data[]=%d got %d\n", target, A362909_data[target], n); } } /* ------------------------------------------------------------------------ */ /* Terms of A362909 Found */ /* found[found_end] is the last found value */ int found[256]; /* A362909 */ int found_end = 0; static inline int found_get (int target) { assert (1 <= target && target <= found_end); return (found[target]); } void found_set (int target, int n) { found_end++; if (found_end >= numberof(found)) error("oops, found[] array length exceeded"); if (target != found_end) error("oops found_set() expected %d next, got %d\n", found_end, target); found[found_end] = n; } /* ------------------------------------------------------------------------ */ /* Calculate Bounds */ /* calculate_bounds() takes v_n = v[n] and calculates bound constraint ----- ---------- step_max step <= step_max first_value_min v[n-step] >= first_value_min To improve on "longest" progression, need longest+1 many values descending v[n], v[n]-vstep, ..., v[n]-longest*vstep. Values are all >= 1 so the end value must be v[n]-longest*vstep >= 1 vstep <= (v[n]-1)/longest vstep_max = floor( (v[n]-1)/longest ) This is enforced as a smallest permitted first value v[n-step] in the progression, first_value_min = v[n] - vstep_max The endmin value under this vstep_max downward progression of longest+1 many terms is endmin = v[n]-longest*vstep_max A value >= endmin only occurs in v[] at or after its first found A362909(endmin). This limits the index step. Index steps are indices n, n-step, ..., n-longest*step and must have n-longest*step >= A362909(endmin) step <= (n-A362909(endmin))/longest step_max = floor( (n-A362909(endmin))/longest ) By construction 1 <= endmin <= v[n] so that A362909(endmin) is available in the found[] terms of A362909 so far. (A simpler calculation would replace A362909(endmin) with just the first index 1.) */ void calculate_bounds (int v_n, int n, int longest, int *step_maxp, int *first_value_minp) { int v_n_sub_1 = v_n - 1; int vstep_max = v_n_sub_1 / longest; int remainder = v_n_sub_1 % longest; *first_value_minp = v_n - vstep_max; int endmin = remainder + 1; assert (endmin == v_n - longest*vstep_max); int step_max = (n - found_get(endmin)) / longest; *step_maxp = step_max; assert (found_get(endmin) >= 1); assert (n - longest*step_max >= found_get(endmin)); assert (n - longest*(step_max+1) < found_get(endmin)); DEBUG(printf(" bounds longest %d vstep_max %d first_value_min %d\n", longest, vstep_max, *first_value_minp)); DEBUG(printf(" endmin %d found_get %d step_max %d\n", endmin, found_get(endmin), step_max)); } /* ------------------------------------------------------------------------ */ const char data_filename[] = "A362881-bytes.data"; FILE *data_fp; /* increase size of v[] if necessary ready to store v[n] */ static inline void v_grow (unsigned char **vp, size_t *vsizep, int n) { assert (n>=0); if (UNLIKELY(n >= *vsizep)) { *vsizep = (n+1)*4/3 + 256; (*vp) = realloc_or_die(*vp, *vsizep); } } /* Set v[n] = value and write to data_fp. Increase the size of v[] if necessary. */ static inline void v_store (unsigned char **vp, size_t *vsizep, int n, int value) { DEBUG(printf("v_store() n=%d value %d ftell(data_fp)=%ld\n", n, value, ftell(data_fp))); v_grow (vp,vsizep, n); assert (n < *vsizep); (*vp)[n] = value; assert (n == ftell_or_die(data_fp,data_filename) + 1); fputc_or_die(value, data_fp, data_filename); if (flag_flush) { flag_flush = 0; fflush_or_die(data_fp, data_filename); } } /* find_longest() returns the length of the longest progression within v[1..n] ending at v[n]. Must have n>=2 so the code can assume a progression of length 2 always exists. */ static int find_longest (const unsigned char *v, int n) { DEBUG(printf("find_longest() n=%d v[n]=%d", n, v[n]); if (n <= 30) { int i; printf(" v ="); for(i=1;i<=n;i++) printf(" %d",v[i]); } printf("\n")); assert (n>=2); int v_n = v[n]; int longest = 2; int step_max, first_value_min; calculate_bounds (v_n, n, longest, &step_max, &first_value_min); int step; for (step = 1; LIKELY(step <= step_max); step++) { DEBUG(printf(" try step %d first v[n-step] = %d\n", step, v[n-step])); int len = 2; /* length of this progression */ int pos = n - step; assert (pos >= 1); int want_value = v[pos]; int value_step = v_n - want_value; if (UNLIKELY(want_value < first_value_min)) continue; assert (v_n - longest*value_step >= 1); while (LIKELY((pos -= step) >= 1)) { want_value -= value_step; DEBUG(printf(" compare pos=%d want_value %d but got %d\n", pos, want_value, v[pos])); if (LIKELY(v[pos] != want_value)) break; len++; } if (UNLIKELY(len > longest)) { DEBUG(printf(" new longest %d with step=%d\n", len, step)); longest = len; calculate_bounds (v_n, n, longest, &step_max, &first_value_min); } } DEBUG(printf(" result longest %d\n", longest)); return (longest); } void output_record_high (int target, int n) { printf ("%d %d\n", target, n); check_A362909 (target, n); } /* Read data_filename into v[], or set initial values if no existing file. */ void read_data_file (unsigned char **vp, size_t *vsizep, int *np, int *targetp) { init_flush(); data_fp = fopen_or_die(data_filename, "a+b"); fseek (data_fp, 0, SEEK_END); /* 1..ndata inclusive in file */ int ndata = ftell_or_die (data_fp, data_filename); rewind(data_fp); ferror_die (data_fp, data_filename, "seeking"); v_grow(vp,vsizep, ndata); (*vp)[0] = 255; /* no value at array v[0] */ fread_or_die ((*vp)+1, 1,ndata, data_fp, data_filename); if (ndata == 0) { printf("# A362909\n"); output_record_high(1,1); /* A362909(1) = 1 */ } else { printf("# resume from n=%d in %s\n", ndata, data_filename); } while (ndata < 2) v_store (vp,vsizep, ++ndata,1); /* A362881(1,2) = 1 */ *np = ndata; int target = 1; int n; for (n = 1; n <= ndata; n++) { if ((*vp)[n] == target) { found_set (target, n); check_A362909 (target, n); target++; } } /* target is 1 bigger than maximum in v[1..ndata] */ *targetp = target; DEBUG(printf("initial target %d\n", *targetp)); if (UNLIKELY(target > UCHAR_MAX)) error("Reached maximum in unsigned char value type"); } void find_records (int limit) { /* v[1..limit] = are terms of the underlying A362881. n = index into v[]. */ unsigned char *v = NULL; size_t vsize = 0; int n, target; read_data_file (&v,&vsize, &n,&target); while (LIKELY(n < limit)) { /* progression ending at n, to make v[n+1] */ int t = find_longest(v,n); n++; assert (t <= target); if (UNLIKELY(t == target)) { output_record_high (target,n); found_set (target,n); target++; if (UNLIKELY(target > UCHAR_MAX)) error("Reached maximum in unsigned char value type"); } v_store(&v,&vsize, n, t); } printf("# stop at limit %d\n", limit); } int main (int argc, char *argv[]) { int limit = INT_MAX; int i, end; setbuf(stdout,NULL); for (i = 1; i < argc; i++) { const char *arg = argv[i]; if (strcmp(arg,"restart")==0) { if (unlink(data_filename) != 0 && errno != ENOENT) printf("Cannot remove %s: %s\n", data_filename, strerror(errno)); } else if (sscanf(arg, "%u%n", &limit, &end) == 1 && end == strlen(arg)) { /* limit parameter */ } else { error("Unrecognised option %s\n", arg); } } find_records(limit); return(0); }