/* This work is licensed under a Creative Commons Attribution 2.5 License. Details here: http://mrob.com/cc-license.html mcsfind - Search for a sequence A description of the problem this program solves, and details of the sequence numbering system, etc. are at mrob.com/pub/math/MCS.html For usage instructions, invoke with -h as the sole argument. */ /* Limits for lengths of generated sequences; the first is for queries and the other is for the display commands -seq##, MCS## and horadam */ #define SCAN_LIMIT 64 #define ALLOC_LIMIT 256 /* Defaults for the -lc## option: 2 and 4 lines on an 80-column TTY */ #define DEFAULT_MAX_CHARS 159 #define DEFAULT_MC_CMDSEQ 319 #include #include #include #include void usage() { const char * help_text[] = { "USAGE ", "", " mcsfind [-n#] [-v] [-f#|-r] [-w#] [-lx##]... [-D##] N1 N2 N3 ...", " mcsfind -seq###", " mcsfind horadam A B R S", "", "mcsfind orders its output by progressively increasing formula complexity.", "It will search more and more complex sequences until the specified terms", "are found or other optional limits (-n#, -ls#) are reached.", "", "The following options are available:", "", "Searching options", " -w1 Original simple weighting", " -w2 Logarithmically weighted coefficients", " -n7 Match up to 7 sequences (default is 3)", " -v Show variations (different formulas that produce the same", " sequence)", " -lc300 Limit characters: display up to 300 characters of each sequence", " (default is 159, or 319 for -seqNN command)", " -lm999 Limit magnitude: Do not display (or attempt to calculate) any terms", " whose absolute value exceeds the given amount (999 in this example)", " -lo2 Limit order: Consider only defintions up to 2nd-order recurrence", " (default allows all, currently up to 3rd order)", " -ls5 Limit score: Consider only defintions of complexity less than", " or equal to the specified value (default is 12.0)", " -lp10 Limit to positive values: In order to match, a candidate sequence", " must produce a value of 2 or more within the first N terms. Legacy", " behavior of #mcs# is achieved by specifying -lp10", "", " -f0 Plain ASCII format", " -f1 ASCII with RHTF markup", " -f2 ANSI terminal (recommended on MacOS X 10.2)", " -f3 ANSI plus UTF-8 characters (subscripts, etc.)", " -r same as -f1", "", "Other mcsfind commands", " -seq440 Display sequence MCS440", " MCS440 Display sequence MCS440", " hor 0 1 2 3 Display Horadam sequence (0,1,2,3)", "", " -D27 Set debug bitmask", "", "Additional explanation and background is available at:", "", " www.mrob.com/pub/math/MCS.html", "", 0}; int i; for(i=0; help_text[i]; i++) { printf("%s\n", help_text[i]); } } /* PURPOSE This program is designed to augment MCS.per by performing a deep search for a single desired sequence. It avoids MCS's huge memory overhead by not storing all the sequences, at the cost of having to re-generate everything for every search. TTD Support matching terms with unspecified middle, e.g. "mcsfind 27 .. 143" would find sequences that contain 27 and 143 in that order, with any number of terms in between. The number of dots could specify a maximum gap length. After a few seconds, split process into multiple threads; each thread traversing the tree in a different way such that total coverage is guaranteed without overlap (example: when a parent node's complexity PC is in the range [max_cplx-MAXC .. max_cplx], explore only those child nodes whose branch number N plus PC == TN mod TT, where TN is our thread number and TT the total number of threads.) Split again each time the workload has increased fourfold. Compare wall-clock time to utime to regulate number of threads. Changes in Weights and Generation: (None of these changes should "break" existing sequence numbers, unless there turns out to be a bug in the current implementation of &seqnum) Find a better balance between fine granularity and efficiency (as it stands, "mcs -n1 388 4957" if over 5 times slower than "mcs -w1 -n1 388 4957") Make surcharge for negative a little less than it is now. Add score improvements to MCS.per and adjust its $max_cplx value accordingly Surcharge for position in sequence (for seq-howto) should be logarithmic. In general, if the score component weights are proportional to the bitlength of the corresponding Huffman symbols, the distribution of MCS numbers in the webpage listing will make more sense. &seqnum seems to have a bug -- "$b2 =~ s/[1.]$//;" allowing a ". to be deleted instead of a "1". Actually this is saved by the greedy quality of regexp phrases, but I'd like to make it more explicit huffman table needs to be expanded for larger coefficients If both Am1 and Nam1 terms exist, print formula in factored form, and less surcharge for complexity score (Similar for An and Nan) Try to do similar tests for other factored forms, general case is (K1 N + K2) (C1 + C2 An + C3 Am1) REVISION HISTORY (Brief summary of the Perl version "MCS.per"; see its revision history for more): 20041111 Began development of MCS.per 20041204 MCS.rhtf is pretty much finished (long hiatus) 20070711 Begin modifications to make it useful as a way to search for sequences. This leads to the decision to write a faster version in C. 20070921 First version. Pretty much get it generating the same sequences that MCS.per generates. {To test this: invoke these two commands and compare their output: ORIGINAL COMMAND SYNTAX circa 20070922: ./mcsfind -n100 ./MCS.per -n 100 CURRENT SYNTAX as of 20091130: ./b-mcsfind -r -w1 -lp10 -n100 -20091201} 20070922 Pass numerical arguments to search for a subsequence. It will start with very simple formulas, then progressively more complex ones. Add coeff and print2. Improve print2 formatting to agree with MCS.per. 20070923 Change default max_match to 3. Add -r option to select RHTF output (default now uses ANSI terminal color styles). Add some UTF-8 encoded superscripts and subscripts. 20070924 Add -f option to select plain ASCII, RHTF or UTF-8 in addition to the (default) ANSI. "-r" still selects RHTF (same as "-f 1"). 20070926 Add higher-granularity weight option (specified with "-w 2") with logarithmic weights. Tweak it so that the following sequences are given in the order shown: seqnum score huffman formula MCS 220 1.0 1.101110.0'100 N MCS 884 1.7 1.101110.100'100 N+1 MCS 882 2.0 1.101110.0.10 10 2 N MCS 440 2.3 1.101110.0.0'100 N^2 MCS 15520 2.7 1.1110.0.1010.0.0.0.0'100 2 A[n] (In the original weighting, all of these had the same score, 2 points.) 20071005 More changes to high-granularity weights to make things more consistent. 20071006 WEIGHT_LOG is now the default. Detect and suppress alternate definitions ("variations") of the same sequence; add -v option to get original behavior. 20071007 Change argument format ("-n 10" becomes "-n10") to make it a little less confusing when arguments and sequence members are combined in the same command line. Add enhuff() and supporting functions. Output still agrees with old perl script, using the following method to compare: ./MCS.per -n 1000 > outMCS mcs -w1 -n1000 -v -f1 > outmcsfind diff outMCS outmcsfind 20071008 Improve efficiency of search with new weighting by adding the "range" parameter to scan1. It's still pretty bad though; I'd like to find a better way without losing the benefits of higher granularity. 20071010 Add -D (debug) option. Use arrays to cache map1, map2, etc. functions. 20080206 Add --help option. 20090617 Fix minor errors in help text; add some comments. 20091122 Add -seqNNN option. Add -lcNNN option. Increase SCAN_LIMIT to 32. This option affects search speed; to test it I am using commands like "time ./mcsfind -n27 3 5 7 8 9 10". Increasing SCAN_LIMIT from 24 to 32 increased the time of this command from 25.746s to 28.987s. Add 128-bit versions of generate() and print1() to support more output from the new command -seqNNN. "./mcsfind -lc1000 MCS30992" now prints the factorials up to 31!=8.22e33. Spend quite a while diagnosing the output of MCS206182547, which stopped mysteriously after the first negative term; I discover a bug in the 128-bit multiply. (MCS206182547 is fairly interesting, an exponentially increasing oscillating sequence with a period of about 9 terms per oscillation) 20091123 Add -lo option, to find lower order recurrence relations and other things that can't be found as easily otherwise. Now "MCS 3 7 27 -lo0" finds pure polynomial formulas (in this case 2 N^3 + 2 N^2 + 3) and "MCS -lo2 2 4 6 7 8" finds complex oscillations like MCS3616484370 that more typically show up as 3rd-order recurrences. {Now with the negative sequence changes, you should add "-lp10" and "-n10" to see some oscillating sequences. Another pattern, "mcs -n7 -lo2 2 4 6 7 6 -lp10" gives 7 sequences with different oscillation periods: MCS445983762, MCS1784118290, MCS884446988, MCS7204960786, MCS1769335200, MCS57871556818, MCS57640999442 20091126 Make it possible to search for and match negative values, mainly to support my work in recreating the Floretion sequences and things like "mcs 1 0 -1 0 1 0 -1 0" (the real component of the powers of i). Add -lp option, with "-lp10" making it act as it did before today's changes. When a sequence matches, subsequent negative terms are printed out so long as it satisfied the -lpNN criterion (an example of this is "mcs -lo2 2 4 6 7 6 -lp10", which finds MCS445983762; before today this sequence was found but no terms after the "-49" were not shown, and you have to do "mcs MCS445983762" to see them.) This change has some drawbacks, in particular regarding the uniqueness matching. If you do "mcs 2 5 9", now you get MCS1972, MCS62017 and MCS250433 -- three sequences that differ only in initial 0 and -1 terms. 20091128 Add matchmulti and addseq. These implement a new, more sophisticated CRC system. Each sequence is CRC'd multiple times based on different subsequences all starting with an initial start term. (If a sequence is reported and is CRC'd out to the 40th term, and another identical sequence is found later, the new identical sequence might only go out to the 38th term because of an extra 2 terms at the beginning. That means that the already-reported sequence now needs to be CRC'd out to 38 terms to effectively match against this one. I implement this by computing all CRCs when the sequence is found, even though many of the shorter ones might not actually get used.) The result of this change is to restore the old behavior regarding 0 and 1 terms: It now no longer reports sequences as unique that only differ in initial -1, 0 or 1 terms (as before, you can use the -v flag if you want these). 20091130 Increase ALLOC_LIMIT to 256 20091201 Add -lm### and -old-negprint options; increase SCAN_LIMIT to 64. These changes are to support the #compare# script and to help ensure that my changes to the MCS.html webpage remain consistent with recent changes to mcsfind.c (most importantly, increasing the number of displayed terms, and including sequences with negative terms). 20091202 Allow commas with arbitrary spacing in argument list, e.g. "mcs 2,3,5,8,13,21,34" or "mcs 5,25 ,125, 625, 3125 , 15625". This makes it more compatible with the #DBIS# script; I frequently use both to search fr the same sequence. Add -ls option, with default 12: #mcs# will no longer search forever if you give it an impossible problem (specify -ls99 if you still want that behavior) Major change: After discovering that I cannot find the MCS number corresponding to A077966, I decide to remove the tests that were designed to eliminate certain types of constant and trivial interleaved sequences. 20091203 fflush after writing each result to stdout 20091208 Remove requirement that mode-0 sequences must have a non-negative initial term. 20091216 Add horadam command option and cmd.horadam(). Add print3 and print_d; display sequences in "A[N]=..." form whenever possible and otherwise use K instead of N. 20101029 Fix a few compiler warnings for GCC 4.2.1 20101030 Add seqnum_load(), seqnum_def() and fmt_seqnum(); RHTF format labels now include a specific label in the link, when available (using ~/data/MCS-links.txt generated by MCS.per). 20101106 Use string array in usage() 20120107 Make the "Extra arguments were ignored" message more useful. TO DO and FUTURE IMPROVEMENTS Certain search terms seem to confuse the algorithm. "mcs 1 1 2 3 5" and "mcs 2 3 5 8" do not find MCS840, but "mcs 2 3 5 8 13" and "mcs 0 1 1 2 3" work. Three different options for handling terms prior to the match pattern for purposes of determining uniqueness. Given a search for "2 5 9", the three options would cause the following three different behaviors: * MCS1972, MCS62017 and MCS250433 are distinct (current "mcs -v 2 5 9") * MCS1972, MCS62017 and MCS250433 are "identical" but MCS1972, MCS2012290 and MCS13781018 are distinct (current "mcs 2 5 9") * MCS1972 and MCS2012290 are "identical"; MCS13781018 is still distinct This can be implemented through the current system of CRCs for each reported sequence, but starting the CRC at the first matched term, rather than at the first "significant" term. 20100108 Add crc32_init. */ #include "../include/rpmtypes.h" #include #define FORMAT_PLAIN 0 // Purely plain ASCII and no RHTF #define FORMAT_RHTF 1 // Plain ASCII; styles expressed via RHTF #define FORMAT_ANSI 2 // Use ANSI terminal colors, etc. to display some // additional info #define FORMAT_UTF8 3 // Use ANSI plus UTF-8 (works in Tiger, not in Jaguar) #define WEIGHT_ORIG 1 // original weights #define WEIGHT_LOG 2 // finer weights used to give coefficients a logarithmic // distribution int o_format; int o_rhtf_simpler; int o_max_match; int o_weight; int o_variations; int o_debug; int o_max_charlen; int o_cmdseq_lc; int o_max_order; int o_negokay; int o_poswindow; int o_legacy_negprint; int o_cplx_limit; int o_horadam; s64 g_too_big; s128 g_toobig128; /* %%% for completeness we want to define huffman table values for all coefficients and then we can set maxc = max_cplx; */ #define MINC -8 #define MAXC 8 /* ----------------------------- prototypes ------------------------------ */ void crc32_init(void); void crc32_start(u32 *crc); void crc32_add8(u32 *crc, u8 data); u32 crcseq(s64 *s, int n); int crc_match(u32 c, int siglen); char * coeff(int x, int one); void init_format(void); char * coeff_hsym(int x); void enhuff(int *a, int n, char *s, int maxlen); int dehuff(char *s, int *a, int maxlen); s64 str10_s64(char * bits); int s64_str10(s64 x, char * bits, int maxlen); s64 coeff_seqnum(int mode, int cK, int cAn, int cN, int cAm1, int cNan, int cN2, int cAm2, int cNam1, int cN3, int d2, int d1, int i0); void coeff_expand(int * ca, int * mode, int * cK, int * cAn, int * cN, int * cAm1, int * cNan, int * cN2, int * cAm2, int * cNam1, int * cN3, int * d2, int * d1, int * i0); void print2(int mode, int cK, int cN, int cN2, int cN3, int i0, int cAn, int cNan, int i1, int cAm1, int cNam1, int i2, int cAm2); void print3(int mode, int i0, int cAn, int i1, int cAm1); void print_d(int mode, int cK, int cN, int cN2, int cN3, int i0, int cAn, int cNan, int i1, int cAm1, int cNam1, int i2, int cAm2); int generate(s64 * A, int mode, int cK, int cN, int cN2, int cN3, int i0, int cAn, int cNan, int i1, int cAm1, int cNam1, int i2, int cAm2, int limit, int negokay, int * p_final_i); int gen_128(s128 * A, int mode, int pcK, int pcN, int pcN2, int pcN3, int pi0, int pcAn, int pcNan, int pi1, int pcAm1, int pcNam1, int pi2, int pcAm2, int limit, int * p_final_i); void print128(s128 * A, int max_cl, int final_i, int match_where); void print1(s64 * A, int max_cl, int final_i, int match_where); int iabs(int x); void find_sig(s64 * A, int final_i, int *firstsig, int *siglen); int ilw(int x); int map1(int x); int map2(int x, int surcharge); int map3(int x); int map4(int x); int match1(s64 * A, int final_i, int * p_match_where); int matchmulti(s64 * s, int siglen); void addseq(s64 * s, int siglen); void scan1(int min_cplx, int range); s64 str_s64(char *s, char * *r); void cmd_seq(s64 p1); void cache_weights(void); /* ---------------------------- seqnum cache ----------------------------- */ /* These routines keep track of the seqnums that were generated by MCS.per which we use to figure out what type of link to emit in RHTF output mode */ int g_snc_init = 0; long g_snc_nseq; s64 * g_snc_def; /* Contains S if seqnum S is defined */ s64 * g_snc_ref; /* Seqnum of seq this is an alias for, else equal to S */ void seqnum_load(void) { char * hd; char pn[200]; FILE * IN; g_snc_nseq = 0; hd = getenv("HOME"); if (hd) { sprintf(pn, "%s/data/MCS-links.txt", hd); // printf("pn is %s\n", pn); IN = fopen(pn, "r"); if (IN) { char * res; int na; long sn1, sn2, nseq; nseq = 0; res = hd; while (res) { res = fgets(pn, 200, IN); if (res) { switch(pn[0]) { case 'N' : /* Number of sequences to allocate */ na = sscanf(pn, "N %ld", &g_snc_nseq); if (na) { // printf("N %ld\n", g_snc_nseq); g_snc_def = (s64 *) malloc((size_t) g_snc_nseq * sizeof(s64)); g_snc_ref = (s64 *) malloc((size_t) g_snc_nseq * sizeof(s64)); if (g_snc_def && g_snc_ref) { // printf("malloc successful\n"); nseq = 0; /* This is redundant but meant to cover the case of multiple N lines in the file */ } else { g_snc_nseq = 0; } } break; case 'a' : /* Normal sequence definition */ na = sscanf(pn, "a %ld", &sn1); if (na && (nseq < g_snc_nseq)) { // printf("a %ld\n", sn1); g_snc_def[nseq] = sn1; g_snc_ref[nseq++] = sn1; } break; case 'b' : /* Alias sequence definition */ na = sscanf(pn, "b %ld %ld", &sn1, &sn2); if (na && (nseq < g_snc_nseq)) { // printf("b %ld %ld\n", sn1, sn2); g_snc_def[nseq] = sn1; g_snc_ref[nseq++] = sn2; } break; default: break; } } } fclose(IN); if (nseq) { /* Set g_snc_nseq to the actual number of sequences we got */ g_snc_nseq = nseq; } } } } /* seqnum_def tells if a seqnum is defined, and if so, which seqnum it is an alias to. Pass in the seqnum S you are interested in. If the return value is: 0 The sequence is not defined S The sequence is defined and not an alias > 0 The sequence is an alias to this other seqnum */ s64 seqnum_def(s64 seqnum) { long i; if (g_snc_init == 0) { seqnum_load(); g_snc_init = 1; } for(i=0; i> 24L) ^ data; *crc = (*crc << 8L) ^ crc32_conv_8[index]; } /* Given a pointer to a sequence, and the number of terms, returns a CRC of that sequence. */ u32 crcseq(s64 *s, int n) { u32 rv; s64 t; int i, j; crc32_start(&rv); for(i=0; i>= 8; } } return rv; } int crc_match(u32 c, int siglen) { int rv; crcnode * l; rv = 0; l = crclists[siglen]; while(l) { if (l->crc == c) { return 1; } l = l->next; } return 0; } /* ---------------- coefficient handling and conversion ------------------ */ char coeff_scratch[40]; int do_plus; char * coeff(int x, int one) { coeff_scratch[0] = 0; if (x < -1) { sprintf(coeff_scratch, " - %d", -x); } else if (x == -1) { if (one) { sprintf(coeff_scratch, " - 1"); } else { sprintf(coeff_scratch, " -"); } } else if (x == 1) { if (do_plus) { if (one) { sprintf(coeff_scratch, " + 1"); } else { sprintf(coeff_scratch, " +"); } } else if (one) { sprintf(coeff_scratch, " 1"); } } else { if (do_plus) { sprintf(coeff_scratch, " + %d", x); } else { sprintf(coeff_scratch, " %d", x); } } if (x != 0) { do_plus = 1; } return coeff_scratch; } /* ------------------------- output formatting --------------------------- */ char ks_MCS[20]; int g_rhtf_labtest; char ital_on[20]; char ital_off[20]; char bold_on[20]; char bold_off[20]; char small_on[20]; char small_off[20]; char sub_on[20]; char sub_off[20]; char sup_on[20]; char sup_off[20]; char text_N[10]; char text_subN[10]; char text_K[10]; char text_subK[10]; char text_subplus[10]; char text_subminus[10]; char text_sub0[10]; char text_sub1[10]; char text_sub2[10]; char text_sub3[10]; char text_sup2[10]; char text_sup3[10]; void init_format(void) { /* Defaults for everything */ strcpy(ks_MCS, "MCS"); g_rhtf_labtest = 0; strcpy(ital_on, ""); strcpy(ital_off, ""); strcpy(bold_on, ""); strcpy(bold_off, ""); strcpy(small_on, ""); strcpy(small_off, ""); strcpy(sub_on, "["); strcpy(sub_off, "]"); strcpy(sup_on, "^"); strcpy(sup_off, ""); strcpy(text_N, "N"); strcpy(text_subN, "N"); strcpy(text_K, "K"); strcpy(text_subK, "K"); strcpy(text_subplus, "+"); strcpy(text_subminus, "-"); strcpy(text_sub0, "0"); strcpy(text_sub1, "1"); strcpy(text_sub2, "2"); strcpy(text_sub3, "3"); strcpy(text_sup2, "2"); strcpy(text_sup3, "3"); switch(o_format) { case FORMAT_PLAIN: default: break; case FORMAT_RHTF: if (o_rhtf_simpler) { strcpy(ks_MCS, "MCS"); } else { strcpy(ks_MCS, "[MCS|+mcs]"); g_rhtf_labtest = 1; } strcpy(ital_on, "*"); strcpy(ital_off, "*"); strcpy(bold_on, "#"); strcpy(bold_off, "#"); strcpy(small_on, "<$:"); strcpy(small_off, ":>"); strcpy(sub_on, "v{"); strcpy(sub_off, "}"); strcpy(sup_on, "^{"); strcpy(sup_off, "}"); break; case FORMAT_UTF8: strcpy(sub_on, ""); strcpy(sub_off, ""); strcpy(sup_on, ""); strcpy(text_N, "X"); /* We use the Greek letter chi subscript, which looks like an "X", in place of N and K */ strcpy(text_K, "X"); sprintf(text_subN, "%c%c%c", 0xe1, 0xb5, 0xaa); sprintf(text_subK, "%c%c%c", 0xe1, 0xb5, 0xaa); sprintf(text_subplus, "%c%c%c", 0xe2, 0x82, 0x8a); sprintf(text_subminus, "%c%c%c", 0xe2, 0x82, 0x8b); sprintf(text_sub0, "%c%c%c", 0xe2, 0x82, 0x80); sprintf(text_sub1, "%c%c%c", 0xe2, 0x82, 0x81); sprintf(text_sub2, "%c%c%c", 0xe2, 0x82, 0x82); sprintf(text_sub3, "%c%c%c", 0xe2, 0x82, 0x83); sprintf(text_sup2, "%c%c", 0xc2, 0xb2); sprintf(text_sup3, "%c%c", 0xc2, 0xb3); // fall through to ANSI case FORMAT_ANSI: strcpy(ital_on, "\033[1m"); strcpy(ital_off, "\033[22m"); strcpy(bold_on, "\033[1;4;94m"); strcpy(bold_off, "\033[22;24;30m"); strcpy(small_on, "\033[32m"); strcpy(small_off, "\033[0m"); break; } } #define N_HSYM 18 const char *huff_sym[N_HSYM] = {"", "0", /* 0 */ "100", /* 1 */ "110", /* -1 */ "1010", /* 2 */ "1110", /* 3 */ "10110", /* -2 */ "101110", /* 4 */ "111100", /* 5 */ "111110", /* -3 */ "1011110", /* -4 */ "1011111", /* -5 Last symbol from original code tree */ "1111010", /* 6 First expansion symbol */ "1111011", /* -6 */ "111111000", /* 7 */ "111111100", /* -7 */ "111111001", /* 8 */ "111111101"};/* -8 */ /* More could be added here if desired in the future */ int huff_name[20] = { 999, 0, 1, -1, 2, 3, -2, 4, 5, -3, -4, -5, 6, -6, 7, -7, 8, -8}; /* coeff_hsym takes a coefficient (like 3 or -7) and returns a pointer to a string of 1's and 0's representing its huffman symbol */ char * coeff_hsym(int x) { int i; for(i=0; i 0) && (bits[i-1] == '0')) { i--; } /* Remove a single trailing '1' */ if (i>=2) { i--; } bits[i] = 0; //printf(" after: %s\n", bits); rv = str_bindec(bits); if (i>=63) { fprintf(stderr, "Seqnum huffman-encoding %s is too big for s64\n", bits); exit(-1); } return(rv); } /* Converts a 64-bit integer into a character string of 1's and 0's. The reverse of this operation is performed in str_bindec The equivalent function in MCS.per is &decbin */ int s64_str10(s64 x, char * bits, int maxlen) { char *s; s64 n; int bitlen; int i; i = 0; /* Count number of significant bits in x */ bitlen = 0; n = x; while(n) { n >>= 1; bitlen++; } /* Restore value of n for use in the conversion */ n = x; s = bits; if (maxlen <= 0) { return(0); } if (maxlen == 1) { s[i] = 0; return(0); } if (n == 0) { s[i++] = '0'; s[i] = 0; return(1); } /* Get the first '1' bit into the sign position */ while(n >= 0) { n <<= 1; } /* Get all the 1 bits into the string */ while ((n != 0) && (i+1 < maxlen)) { s[i++] = ((n<0) ? '1' : '0'); n <<= 1; } /* Check for string too short. */ if (i+1 >= maxlen) { s[i] = 0; fprintf(stderr, "s64_str10: err1 x=%lld\n", x); exit(-1); } /* Add any more needed 0's */ while((i+1 < maxlen) && (i= maxlen) || (i= 1) { ca[i++] = cAn; } ca[i++] = cN; if (mode >= 2) { ca[i++] = cAm1; } if (mode >= 1) { ca[i++] = cNan; } ca[i++] = cN2; if (mode >= 3) { ca[i++] = cAm2; } if (mode == 3) { ca[i++] = cNam1; } ca[i++] = cN3; if (mode == 2) { ca[i++] = cNam1; } if (mode >= 3) { ca[i++] = ci2; } if (mode >= 2) { ca[i++] = ci1; } if (mode >= 1) { ca[i++] = i0; } enhuff(ca, i, bits, 100); if (o_debug & 2) { printf(" bits: %s\n", bits); } return(str10_s64(bits)); } /* Given a pointer to a list of tokens (coefficients), fills in all the values of the coefficient variables. The inverse of this function is performed by coeff_seqnum. */ void coeff_expand(int * ca, int * mode, int * cK, int * cAn, int * cN, int * cAm1, int * cNan, int * cN2, int * cAm2, int * cNam1, int * cN3, int * d2, int * d1, int * i0) { int i, m, ci1, ci2; i = 0; m = 4-ca[i++]; *mode = m; *cK = ca[i++]; if (m >= 1) { *cAn = ca[i++]; } else { *cAn = 0; } *cN = ca[i++]; if (m >= 2) { *cAm1 = ca[i++]; } else { *cAm1 = 0; } if (m >= 1) { *cNan = ca[i++]; } else { *cNan = 0; } *cN2 = ca[i++]; if (m >= 3) { *cAm2 = ca[i++]; } else { *cAm2 = 0; } if (m == 3) { *cNam1 = ca[i++]; } else { *cNam1 = 0; } *cN3 = ca[i++]; if (m == 2) { *cNam1 = ca[i++]; } else { *cNam1 = 0; } if (m >= 3) { ci2 = ca[i++]; } else { ci2 = 0; } *d2 = ci2 + 1; if (m >= 2) { ci1 = ca[i++]; } else { ci1 = 0; } *d1 = ci1 + 1; if (m >= 1) { *i0 = ca[i++]; } else { *i0 = 0; } } void fmt_seqnum(s64 sn); void fmt_seqnum(s64 sn) { s64 snd; snd = 0; if (g_rhtf_labtest) { snd = seqnum_def(sn); } if (snd) { /* We are printing fancy RHTF labels and the sequence is defined in MCS.rhtf */ printf("[MCS%lld|+mcs:l%lld] : ", sn, sn); } else { printf("%s%lld : ", ks_MCS, sn); } } /* Print out the sequence's definition */ void print2(int mode, int cK, int cN, int cN2, int cN3, int i0, int cAn, int cNan, int i1, int cAm1, int cNam1, int i2, int cAm2) { s64 sn; sn = coeff_seqnum(mode, cK, cAn, cN, cAm1, cNan, cN2, cAm2, cNam1, cN3, i2-i1, i1-i0, i0); printf(" "); fmt_seqnum(sn); if (mode >= 1) { printf("%sA%s%s%s%s = %d; ", ital_on, ital_off, sub_on, text_sub0, sub_off, i0); } if (mode >= 2) { printf("%sA%s%s%s%s = %d; ", ital_on, ital_off, sub_on, text_sub1, sub_off, i1); } if (mode >= 3) { printf("%sA%s%s%s%s = %d; ", ital_on, ital_off, sub_on, text_sub2, sub_off, i2); } if (mode == 0) { printf("%sA%s%s%s%s =", ital_on, sub_on, text_subK, sub_off, ital_off); } else { printf("%sA%s%s%s%s%s%s%s%s =", ital_on, ital_off, sub_on, ital_on, text_subK, ital_off, text_subplus, text_sub1, sub_off); } do_plus = 0; if (cNan) { printf("%s %s%s%s %sA%s%s%s%s", coeff(cNan, 0), ital_on, text_K, ital_off, ital_on, sub_on, text_subK, sub_off, ital_off); } if (cAn) { printf("%s %sA%s%s%s%s", coeff(cAn, 0), ital_on, sub_on, text_subK, sub_off, ital_off); } if (cNam1) { printf("%s %s%s%s %sA%s%s%s%s%s%s%s%s", coeff(cNam1, 0), ital_on, text_K, ital_off, ital_on, ital_off, sub_on, ital_on, text_subK, ital_off, text_subminus, text_sub1, sub_off); } if (cAm1) { printf("%s %sA%s%s%s%s%s%s%s%s", coeff(cAm1, 0), ital_on, ital_off, sub_on, ital_on, text_subK, ital_off, text_subminus, text_sub1, sub_off); } if (cAm2) { printf("%s %sA%s%s%s%s%s%s%s%s", coeff(cAm2, 0), ital_on, ital_off, sub_on, ital_on, text_subK, ital_off, text_subminus, text_sub2, sub_off); } if (cN3) { printf("%s %s%s%s%s%s%s", coeff(cN3, 0), ital_on, text_K, ital_off, sup_on, text_sup3, sup_off); } if (cN2) { printf("%s %s%s%s%s%s%s", coeff(cN2, 0), ital_on, text_K, ital_off, sup_on, text_sup2, sup_off); } if (cN) { printf("%s %s%s%s", coeff(cN, 0), ital_on, text_K, ital_off); } if (cK) { printf("%s", coeff(cK, 1)); } } /* This is like print2, but only displays formulas that fit the standard "A[N]=..." form. NOTE: You should call this routine as print3(mode, i0, i1, i2, cAn, cAm1, cAm2, cK) even though the arguments are internally called "cAm1" through "aAm3" */ void print3(int mode, int i0, int i1, int i2, int cAm1, int cAm2, int cAm3, int cK) { s64 sn; sn = coeff_seqnum(mode, cK, cAm1, 0, cAm2, 0, 0, cAm3, 0, 0, i2-i1, i1-i0, i0); printf(" "); fmt_seqnum(sn); if (mode >= 1) { printf("%sA%s%s%s%s = %d; ", ital_on, ital_off, sub_on, text_sub0, sub_off, i0); } if (mode >= 2) { printf("%sA%s%s%s%s = %d; ", ital_on, ital_off, sub_on, text_sub1, sub_off, i1); } if (mode >= 3) { printf("%sA%s%s%s%s = %d; ", ital_on, ital_off, sub_on, text_sub2, sub_off, i2); } printf("%sA%s%s%s%s =", ital_on, sub_on, text_subN, sub_off, ital_off); do_plus = 0; if (cAm1) { printf("%s %sA%s%s%s%s%s%s%s%s", coeff(cAm1, 0), ital_on, ital_off, sub_on, ital_on, text_subN, ital_off, text_subminus, text_sub1, sub_off); } if (cAm2) { printf("%s %sA%s%s%s%s%s%s%s%s", coeff(cAm2, 0), ital_on, ital_off, sub_on, ital_on, text_subN, ital_off, text_subminus, text_sub2, sub_off); } if (cAm3) { printf("%s %sA%s%s%s%s%s%s%s%s", coeff(cAm3, 0), ital_on, ital_off, sub_on, ital_on, text_subN, ital_off, text_subminus, text_sub3, sub_off); } if (cK) { printf("%s", coeff(cK, 1)); } } /* Dispatches to print2 or print3 depending on which is most appropriate */ void print_d(int mode, int cK, int cN, int cN2, int cN3, int i0, int cAn, int cNan, int i1, int cAm1, int cNam1, int i2, int cAm2) { if ((cN==0) && (cN2==0) && (cN3==0) && (cNan==0) && (cNam1==0)) { print3(mode, i0, i1, i2, cAn, cAm1, cAm2, cK); } else { print2(mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2); } } /* ----------------------- scanning and generation ----------------------- */ int generate(s64 * A, int mode, int cK, int cN, int cN2, int cN3, int i0, int cAn, int cNan, int i1, int cAm1, int cNam1, int i2, int cAm2, int limit, int o_negokay, int * p_final_i) { int i, lv, lvs, N, got2, stop; s64 An, Am1, Am2, Ap1; int final_i; An = Am1 = Am2 = 0; /* To avoid warnings */ i = 0; if (mode == 0) { N = 0; A[i++] = cK; } if (mode >= 1) { A[i++] = i0; } if (mode >= 2) { A[i++] = i1; } if (mode >= 3) { A[i++] = i2; } final_i = i; if (o_negokay) { got2 = mode ? mode : 1; } else if ((mode == 0) && (A[0] < 0)) { got2 = 0; // got2 = -1; } else if ((mode < 2) && (A[0] > 1)) { got2 = 1; } else if ((mode == 2) && (A[0] > 1)) { got2 = 2; } else if ((mode == 3) && (A[0] > 1)) { got2 = 3; } else if ((mode == 2) && (A[1] > 1)) { got2 = 1; } else if ((mode == 3) && (A[1] > 1)) { got2 = 2; } else if ((mode == 3) && (A[2] > 1)) { got2 = 1; } else { got2 = 0; } if (mode == 1) { An = A[0]; } if (mode == 2) { Am1 = A[0]; An = A[1]; } if (mode == 3) { Am2 = A[0]; Am1 = A[1]; An = A[2]; } stop = 0; if (mode == 0) { lvs = 1; } else { lvs = mode; } for (lv = lvs; (got2 >=0) && (lv <= limit); lv++) { if (mode == 0) { N = lv; } else { N = lv - 1; } Ap1 = cK + (cN * N) + (cN2 * N * N) + (cN3 * N * N * N) + (cAn * An) + (cNan * N * An) + (cAm1 * Am1) + (cNam1 * N * Am1) + (cAm2 * Am2); if (i g_too_big) { stop++; } else if (o_negokay) { got2++; } else if ((o_negokay == 0) && (Ap1 < 0)) { /* Want to add another condition here, like possibly: if ((Ap1 < 0) && (final_i < 11)) to make the proscription against negative terms optional */ got2 = -1; } else if (got2 < 0) { } else if (got2 > 0) { got2++; } else if (Ap1 > 1) { got2 = 1; } Am2 = Am1; Am1 = An; An = Ap1; } *p_final_i = final_i; return got2; } int gen_128(s128 * A, int mode, int pcK, int pcN, int pcN2, int pcN3, int pi0, int pcAn, int pcNan, int pi1, int pcAm1, int pcNam1, int pi2, int pcAm2, int limit, int * p_final_i) { int i, lv, lvs, N, got2, stop; s128 An, Am1, Am2, Ap1; s128 cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2; int final_i; cK = pcK; cN = pcN; cN2 = pcN2; cN3 = pcN3; i0 = pi0; cAn = pcAn; cNan = pcNan; i1 = pi1; cAm1 = pcAm1; cNam1 = pcNam1; i2 = pi2; cAm2 = pcAm2; An = Am1 = Am2 = 0; /* To avoid warnings */ i = 0; if (mode == 0) { N = 0; A[i++] = cK; } if (mode >= 1) { A[i++] = i0; } if (mode >= 2) { A[i++] = i1; } if (mode >= 3) { A[i++] = i2; } final_i = i; if ((mode == 0) && (A[0] < 0)) { // got2 = -1; } else if ((mode < 2) && (A[0] > 1)) { got2 = 1; } else if ((mode == 2) && (A[0] > 1)) { got2 = 2; } else if ((mode == 3) && (A[0] > 1)) { got2 = 3; } else if ((mode == 2) && (A[1] > 1)) { got2 = 1; } else if ((mode == 3) && (A[1] > 1)) { got2 = 2; } else if ((mode == 3) && (A[2] > 1)) { got2 = 1; } else { got2 = 0; } if (mode == 1) { An = A[0]; } if (mode == 2) { Am1 = A[0]; An = A[1]; } if (mode == 3) { Am2 = A[0]; Am1 = A[1]; An = A[2]; } stop = 0; if (mode == 0) { lvs = 1; } else { lvs = mode; } for (lv = lvs; (got2 >=0) && (lv <= limit); lv++) { if (mode == 0) { N = lv; } else { N = lv - 1; } Ap1 = cK + (cN * N) + (cN2 * N * N) + (cN3 * N * N * N) + (cAn * An) + (cNan * N * An) + (cAm1 * Am1) + (cNam1 * N * Am1) + (cAm2 * Am2); if (i g_toobig128) { stop++; // } else if (Ap1 < 0) { We no longer want to ignore negatives // got2 = -1; } else if (got2 < 0) { } else if (got2 > 0) { got2++; } else if (Ap1 > 1) { got2 = 1; } Am2 = Am1; Am1 = An; An = Ap1; } *p_final_i = final_i; return got2; } s64 match[ALLOC_LIMIT]; int g_match_len; int g_matches_found; int g_num_crc; int x_match_where; /* Print out the sequence */ void print128(s128 * A, int max_cl, int final_i, int match_where) { int i, stop, charlen; char ai[100]; stop = 0; charlen = 3; // for the final "..." for(i=0; (stop == 0) && (i<=final_i) && (charlen < max_cl); i++) { s128_str(A[i], ai); charlen = charlen + strlen(ai) + 2; if (charlen > max_cl) { stop = 1; } else if (A[i] > g_toobig128) { stop = 1; } else { if ((i >= match_where) && (i < match_where + g_match_len)) { printf("%s%s%s, ", bold_on, ai, bold_off); } else { printf("%s, ", ai); } } } printf("...\n"); } /* Print out the sequence */ void print1(s64 * A, int max_cl, int final_i, int match_where) { int i, stop, charlen; char ai[100]; stop = 0; charlen = 3; // for the final "..." for(i=0; (stop == 0) && (i<=final_i) && (charlen < max_cl); i++) { sprintf(ai, "%lld", A[i]); charlen = charlen + strlen(ai) + 2; if (charlen > max_cl) { stop = 1; } else if (A[i] > g_too_big) { stop = 1; } else { if ((i >= match_where) && (i < match_where + g_match_len)) { printf("%s%s%s, ", bold_on, ai, bold_off); } else { printf("%s, ", ai); } } } printf("...\n"); } int iabs(int x) { if (x>=0) { return x; } return -x; } /* Given a sequence, find the first "significant" term (defined as the first term whose absolute value is at least 2) and the number of terms we have calculated starting from that term. */ void find_sig(s64 * A, int final_i, int *firstsig, int *siglen) { int i, stop, Ai; int i_firstsig, i_toobig; i_firstsig = -1; i_toobig = -1; stop = 0; for(i=0; i<=final_i; i++) { Ai = iabs(A[i]); if (stop) { } else if (Ai > g_too_big) { i_toobig = i; stop++; } else { if ((i_firstsig < 0) && (Ai > 1)) { i_firstsig = i; } } } if (i_toobig < 0) { i_toobig = i; } *firstsig = i_firstsig; *siglen = i_toobig - i_firstsig; } int a_ilw[33] = {0, 10, 20, 26, 30, 33, 36, 38, 40, 42, 43, 45, 46, 47, 48, 49, 50, 51, 52, 52, 53, 54, 55, 55, 56, 56, 57, 58, 58, 59, 59, 60, 60}; // Integer logarithmic weight. Functions as iabs if using original weights. int ilw(int x) { int ora; // over-range adjust int nsc; // negative surcharge nsc = 0; if (x<0) { x = -x; nsc = 1; } if (o_weight == WEIGHT_ORIG) { return x; } ora = 0; while (x>32) { ora += 10; x >>= 1; } return a_ilw[x] + ora + nsc; } // Maps -2 -1 0 1 2 ... onto 2 1 0 0 1 ... int map1(int x) { if (x<0) { return ilw(-x); } else if (x <= 0) { return 0; } else if (x <= 1) { if (o_weight == WEIGHT_ORIG) { return 0; } else { return 1; // 0.1 surcharge } } return ilw(x-1); } // Absolute value plus a surcharge if nonzero. For example, if the // surcharge is 1, this maps -2 -1 0 1 2 ... onto 3 2 0 2 3 ... int map2(int x, int surcharge) { if (o_weight == WEIGHT_ORIG) { surcharge /= 10; } if (x != 0) { return (surcharge + ilw(x)); } return 0; } // Simple absolute-value mapping. int map3(int x) { return ilw(x); } // Like map1, but not logarithmically skewed. int map4(int x) { if (o_weight == WEIGHT_LOG) { if (x <= 0) { return (-10 * x); } else { return (10 * (x-1)); } } return (map1(x)); } /* Cached versions of the map functions */ int a_map1[1 + 2*MAXC]; int a_map2_n3[1 + 2*MAXC]; int a_map2_6[1 + 2*MAXC]; int a_map2_13[1 + 2*MAXC]; int a_map2_16[1 + 2*MAXC]; int a_map2_20[1 + 2*MAXC]; int a_map3[1 + 2*MAXC]; int a_map4[1 + 2*MAXC]; /* See if the given sequence matches the search terms. */ int match1(s64 * A, int final_i, int * p_match_where) { int i, j, rv, gg; rv = 0; for(i=0; i<=final_i; i++) { gg = 1; for(j=0; (gg > 0) && (j final_i) { gg = -1; /* ran off end of sequence */ } else if (A[i+j] != match[j]) { gg = 0; /* no match */ } } if (gg > 0) { /* we made it through all the match terms! */ *p_match_where = i; return 1; } } *p_match_where = 0; return rv; } /* CRC a sequence at its current length and all shorter lengths, looking for a match. It stops CRCing as soon as it finds a match. */ int matchmulti(s64 * s, int siglen) { int sl; u32 crc; int unique; unique = 1; sl = siglen; while((sl > 1) && unique) { crc = crcseq(s, sl); if (crc_match(crc, sl)) { unique = 0; } sl--; } return(unique); } /* Remember this CRC, and re-checksum the sequence at each shorter length */ void addseq(s64 * s, int siglen) { int sl; u32 crc; crcnode * n; sl = siglen; while (sl > 1) { crc = crcseq(s, sl); n = (crcnode *) malloc(sizeof(crcnode)); if (n) { n->crc = crc; n->next = crclists[sl]; crclists[sl] = n; } else { fprintf(stderr, "could not allocate %dth crcnode\n", g_num_crc); exit(-1); } g_num_crc++; sl--; } } s64 g_sc1, g_sc2, thispass, thishit, maxpass; /* X + 1 cK surcharge -0.3 2 X cN surcharge 0 X^2 cN2 surcharge 1.3 2 A[n] cAn surcharge 0.6 Cost of a "2" coefficient for each term: K N N^2 N^3 mode 17 20 33 40 0 i0 An N An 10 26 36 10 d1 Am1 N Am1 10 26 36 20 d2 Am2 10 26 30 */ void scan1(int min_cplx, int range) { s64 A[ALLOC_LIMIT]; int max_cplx; int mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2, limit; int i0lo, i0hi, d1, d1lo, d1hi, d2, d2lo, d2hi; int cpl1, cpl2, cpl3, cpl4, cpl5, cpl6, cpl7, cpl8, cpl9, cp10, cp11, cplx; int valid, doit, final_i, match_where; max_cplx = min_cplx + range - 1; if (o_debug & 1) { if (min_cplx == max_cplx) { printf("scan1 [%d]: ", max_cplx); } else { printf("scan1 [%d..%d]: ", min_cplx, max_cplx); } } thispass = thishit = 0; limit = SCAN_LIMIT; for (cAn = MINC; cAn <= MAXC; cAn++) { cpl1 = a_map2_6[cAn + MAXC]; // map2(cAn, 6); for (cNan = MINC; (cNan <= MAXC) && (cpl1 <= max_cplx); cNan++) { cpl2 = cpl1 + a_map2_16[cNan + MAXC]; // map2(cNan, 16); for (cAm1 = MINC; (cAm1 <= MAXC) && (cpl2 <= max_cplx); cAm1++) { cpl3 = cpl2 + a_map2_6[cAm1 + MAXC]; // map2(cAm1, 6); for (cNam1 = MINC; (cNam1 <= MAXC) && (cpl3 <= max_cplx); cNam1++) { cpl4 = cpl3 + a_map2_16[cNam1 + MAXC]; // map2(cNam1, 16); for (cAm2 = MINC; (cAm2 <= MAXC) && (cpl4 <= max_cplx); cAm2++) { cpl5 = cpl4 + a_map2_6[cAm2 + MAXC]; // map2(cAm2, 6); /* %%% add complexity points for mode */ if ((cAn == 0) && (cNan == 0) && (cAm1 == 0) && (cNam1 == 0) && (cAm2 == 0)) { mode = 0; i0lo = 0; i0hi = 0; d1lo = 0; d1hi = 0; d2lo = 0; d2hi = 0; } else if ((cAm1 == 0) && (cNam1 == 0) && (cAm2 == 0)) { mode = 1; i0lo = -1; i0hi = MAXC; d1lo = 0; d1hi = 0; d2lo = 0; d2hi = 0; } else if (cAm2 == 0) { mode = 2; i0lo = -1; i0hi = MAXC; d1lo = -1; d1hi = MAXC; d2lo = 0; d2hi = 0; } else { mode = 3; i0lo = -1; i0hi = MAXC; d1lo = -1; d1hi = MAXC; d2lo = -1; d2hi = MAXC; } cpl5 += a_map4[mode + MAXC]; // map4(mode); /* Prune this branch if it doesn't fit the -lo option */ if (mode > o_max_order) { cpl5 = max_cplx + 1; } for (cK = MINC; (cK <= MAXC) && (cpl5 <= max_cplx); cK++) { cpl6 = cpl5 + a_map2_n3[cK + MAXC]; // map2(cK, -3); for (cN = MINC; (cN <= MAXC) && (cpl6 <= max_cplx); cN++) { cpl7 = cpl6 + a_map3[cN + MAXC]; // map3(cN); for (cN2 = MINC; (cN2 <= MAXC) && (cpl7 <= max_cplx); cN2++) { cpl8 = cpl7 + a_map2_13[cN2 + MAXC]; // map2(cN2, 13); for (cN3 = MINC; (cN3 <= MAXC) && (cpl8 <= max_cplx); cN3++) { cpl9 = cpl8 + a_map2_20[cN3 + MAXC]; // map2(cN3, 20); if ( (cN == 0) && (cN2 == 0) && (cN3 == 0) ) { if( (mode == 0) /* constant sequence */ || ( (mode == 1) /* constant iterative sequence? */ && (cAn == 1) && (cNan == 0) && (cK == 0)) || ( (mode == 2) /* 2 or 3 interleaved unrelated sequences? */ && (cAn == 0) && (cNan == 0) && (cNam1 == 0) ) ) { // cpl9 = max_cplx + 1; // skip following loops } } for(i0 = i0lo; (i0 <= i0hi) && (cpl9 <= max_cplx); i0++) { cp10 = cpl9 + a_map1[i0 + MAXC]; // map1(i0); for(d1 = d1lo; (d1 <= d1hi) && (cp10 <= max_cplx); d1++) { cp11 = cp10 + a_map1[d1 + MAXC]; // map1(d1); i1 = i0 + d1; for(d2 = d2lo; (d2 <= d2hi) && (cp11 <= max_cplx); d2++) { cplx = cp11 + a_map1[d2 + MAXC]; // map1(d2); i2 = i1 + d2; doit = 1; g_sc1++; thispass++; if (doit && (cplx >= min_cplx) && (cplx <= max_cplx)) { thishit++; g_sc2++; /* Create the first o_poswindow terms of the sequence, and see if any of the first 10 terms is greater than 2 */ valid = generate(A, mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2, o_poswindow, o_negokay, &final_i); if (valid == 0 && o_negokay) { fprintf(stderr, "negative okay but not valid?\n"); print1(A, DEFAULT_MAX_CHARS, final_i, match_where); exit(-1); } if (o_negokay || (valid > 0)) { /* "Limit + o_poswindow - valid" here is to ensure that the number of "valid" terms generated is always the same in each sequence we report */ generate(A, mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2, limit + (o_poswindow - valid), (o_legacy_negprint == 0), &final_i); if (match1(A, final_i, &match_where)) { int i_firstsig; int siglen; int doit; doit = 0; if (o_variations) { doit = 1; } else { find_sig(A, final_i, &i_firstsig, &siglen); doit = matchmulti(A + i_firstsig, siglen); if (doit) { addseq(A + i_firstsig, siglen); } // printf("CRC: %u (%d, %d)\n", crc, i_firstsig, siglen); } if (doit) { print1(A, o_max_charlen, final_i, match_where); if (o_debug & 2) { printf(" mode=%d", mode); printf(" cK=%d cN=%d cN2=%d cN3=%d", cK, cN, cN2, cN3); if (mode>=1) { printf(" i0=%d cAn=%d cNan=%d", i0, cAn, cNan); } if (mode>=2) { printf(" i1=%d cAm1=%d cNam1=%d", i1, cAm1, cNam1); } if (mode>=3) { printf(" i2=%d cAm2=%d", i2, cAm2); } printf("\n"); } print_d(mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2); if (o_weight == WEIGHT_ORIG) { printf(" %s(score: %d)%s\n", small_on, cplx, small_off); } else { printf(" %s(score: %d.%d)%s\n", small_on, cplx/10, cplx%10, small_off); } g_matches_found++; fflush(stdout); if (g_matches_found >= o_max_match) { goto exit1; // return; } } } } } } } } } } } } } } } } } if (thispass > maxpass) { maxpass = thispass; } exit1: ; if (o_debug & 1) { printf(" %9lld for %lld", thishit, thispass); if (thispass) { printf(", %.3g %%", 100.0*((double) thishit) / ((double) thispass)); } printf("\n"); } } /* Given a string s, returns an integer (with optional - sign) scanned from the beginning of the string. If the parameter r is provided, it is set to point to the first character after the end of the scanned integer. */ s64 str_s64(char *s, char * *r) { s64 rv; char *c; s64 sign; sign = 1; rv = 0; c = s; if (*c == '-') { sign = -1; c++; } while(*c) { if (((*c) >= '0') && ((*c) <= '9')) { rv *= 10LL; rv += ((s64) ((*c) - '0')); } else { goto str_s64_done; } c++; } str_s64_done: ; if (r) { *r = c; } return (sign * rv); } /* testing with -seq1808378250 mode=2 cK=5 cN=-1 cN2=-1 cN3=0 i0=1 cAn=1 cNan=0 i1=4 cAm1=2 cNam1=0 bits: 1101011110010011010100110001010100 A[0] = 1; A[1] = 4; A[N+1] = A[N] + 2 A[N-1] - N^2 - N + 5 1, 4, 9, 16, 27, 44, 73, 124, 219, 400, 753, 1448, 2827, 5572, 11049, ... */ /* Does the -seqNNN command, which takes an MCS number and prints out the definition and some terms of the sequence. */ #define CMD_SEQ_LEN 100 void cmd_seq(s64 p1) { s128 A[ALLOC_LIMIT]; char bits[CMD_SEQ_LEN]; int coeff[CMD_SEQ_LEN]; int i, d1, d2; int mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2; int limit, final_i; /* these globals affect print1() */ g_match_len = g_matches_found = 0; limit = ALLOC_LIMIT; s64_str10(p1, bits, CMD_SEQ_LEN); // printf("seq %lld, bits %s\ndehuff: ", p1, bits); dehuff(bits, coeff, CMD_SEQ_LEN); // for(i=0; i<35; i++) { printf("%d ", coeff[i]); } printf("\n"); /* Convert coeff[] into the parameters for generate() */ coeff_expand(coeff, &mode, &cK, &cAn, &cN, &cAm1, &cNan, &cN2, &cAm2, &cNam1, &cN3, &d2, &d1, &i0); i1 = i0 + d1; i2 = i1 + d2; if (o_debug & 2) { printf(" mode=%d", mode); printf(" cK=%d cN=%d cN2=%d cN3=%d", cK, cN, cN2, cN3); if (mode>=1) { printf(" i0=%d cAn=%d cNan=%d", i0, cAn, cNan); } if (mode>=2) { printf(" i1=%d cAm1=%d cNam1=%d", i1, cAm1, cNam1); } if (mode>=3) { printf(" i2=%d cAm2=%d", i2, cAm2); } printf("\n"); } gen_128(A, mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2, limit, &final_i); print128(A, o_cmdseq_lc, final_i, 0); print_d(mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2); printf("\n"); } /* Like the -seqNNN command, but creates a Horadam sequence given its (p,q,r,s) */ void cmd_horadam(int p, int q, int r, int s) { s64 p1; int mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2; int d1, d2; int limit, final_i; s128 A[ALLOC_LIMIT]; /* Set all the coefficients */ mode = 2; cK = cN = cN2 = cN3 = 0; i0 = p; cAn = s; cNan = 0; i1 = q; cAm1 = r; cNam1 = 0; d1 = i1 - i0; i2 = 0; cAm2 = 0; d2 = i2 - i1; /* these globals affect print1() */ g_match_len = g_matches_found = 0; limit = ALLOC_LIMIT; p1 = coeff_seqnum(mode, cK, cAn, cN, cAm1, cNan, cN2, cAm2, cNam1, cN3, d2, d1, i0); if (o_debug & 2) { printf(" mode=2"); printf(" cK=%d cN=%d cN2=%d cN3=%d", cK, cN, cN2, cN3); printf(" i0=%d cAn=%d cNan=%d", i0, cAn, cNan); printf(" i1=%d cAm1=%d cNam1=%d", i1, cAm1, cNam1); printf("\n"); } gen_128(A, mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2, limit, &final_i); print128(A, o_cmdseq_lc, final_i, 0); print_d(mode, cK, cN, cN2, cN3, i0, cAn, cNan, i1, cAm1, cNam1, i2, cAm2); printf("\n"); } void cache_weights(void) { int i; for(i=MINC; i<=MAXC; i++) { a_map1[i+MAXC] = map1(i); a_map2_n3[i+MAXC] = map2(i, -3); a_map2_6[i+MAXC] = map2(i, 6); a_map2_13[i+MAXC] = map2(i, 13); a_map2_16[i+MAXC] = map2(i, 16); a_map2_20[i+MAXC] = map2(i, 20); a_map3[i+MAXC] = map3(i); a_map4[i+MAXC] = map4(i); } } int main(int argc, char** argv) { int i, np; s64 mt; /* match term */ s64 n, p1; int i1, i2, i3, i4; char * arg_warn; char *s; crclists = (crcnode * *) malloc(sizeof(crcnode *) * (1+SCAN_LIMIT)); if (crclists == 0) { fprintf(stderr, "could not malloc for crclists\n"); exit(-1); } for(i=0; i<=SCAN_LIMIT; i++) { crclists[i] = 0; } crc32_init(); init_format(); cache_weights(); /* compute MAXINT */ g_too_big = 1; while (g_too_big<<1 > 0) { g_too_big <<= 1; } /* At any step the value can increase by about 2N*MAXC, due to the (cNan * N * An) and (cNam1 * N * Am1) terms. Therefore, we stop trying to calculate a sequence when a term A[n] gets bigger than MAXINT / (2*MAXC*MAXT) where MAXT is the max number of terms. */ g_too_big = g_too_big / ((s64) (2 * MAXC * SCAN_LIMIT)); // g_too_big = 999999999999ll; g_toobig128 = 1; while (g_toobig128+g_toobig128 > 0) { g_toobig128 = g_toobig128+g_toobig128; } /* %%% When I add the AnAm1 term, this will need to be a lot smaller */ g_toobig128 = g_toobig128 / ((s128) (2 * MAXC * ALLOC_LIMIT)); if (0) { char s[100]; s128_str(g_toobig128, s); printf("g_toobig128 is %s\n", s); } /* Options defaults */ o_max_match = 3; o_format = 2; o_rhtf_simpler = 0; o_weight = WEIGHT_LOG; o_variations = 0; o_max_charlen = DEFAULT_MAX_CHARS; o_cmdseq_lc = DEFAULT_MC_CMDSEQ; o_max_order = 10; // Effectively no limit o_negokay = 1; o_poswindow = 10; o_legacy_negprint = 0; o_cplx_limit = 12; o_horadam = 0; g_match_len = 0; arg_warn = 0; for(i=1; i= FORMAT_PLAIN) && (n <= FORMAT_UTF8)) { o_format = n; init_format(); // printf("opt format %d\n", o_format); } } else if (strncmp(argv[i], "-lc", 3) == 0) { if (np = sscanf(argv[i], "-lc%d", &i1), np == 1) { if (i1 > 10) { o_max_charlen = i1; o_cmdseq_lc = i1; } else { fprintf(stderr, "-lc value should be at least 11\n"); usage(); exit(-1); } } else { fprintf(stderr, "-lc needs a number, e.g. '-lc500'.\n"); usage(); exit(-1); } } else if (np = sscanf(argv[i], "-lo%d", &i1), np == 1) { if (i1 >= 0) { o_max_order = i1; } else { fprintf(stderr, "-lo value should be at least 0\n"); usage(); exit(-1); } } else if (np = sscanf(argv[i], "-lm%lld", &p1), np == 1) { /* Do not display (or attempt to calculate) any terms whose absolute value exceeds the given amount */ if ((p1 >= 1) && (p1 <= g_too_big)) { g_too_big = p1; g_toobig128 = p1; } else { fprintf(stderr, "-lm value must be in range [1..%lld]\n", g_too_big); usage(); exit(-1); } } else if (np = sscanf(argv[i], "-lp%d", &i1), np == 1) { /* Limit to positive values: In order to match, a candidate sequence must produce a value of 2 or more within the first N terms. Legacy behavior of #mcs# is achieved by specifying -lp10 */ if (i1 >= (SCAN_LIMIT/2)) { fprintf(stderr, "-lp value cannot be greater than %d\n", (SCAN_LIMIT/2)-1); usage(); exit(-1); } if (i1 > 0) { o_negokay = 0; o_poswindow = i1; } else { fprintf(stderr, "-lp value should be at least 1\n"); usage(); exit(-1); } } else if (np = sscanf(argv[i], "-ls%d", &i1), np == 1) { if (i1 > 1) { o_cplx_limit = i1; } else { fprintf(stderr, "-ls value should be at least 1\n"); usage(); exit(-1); } } else if (strncmp(argv[i], "-old-negprint", 2) == 0) { o_legacy_negprint = 1; } else if (np = sscanf(argv[i], "-n%d", &i1), np == 1) { if (i1 <= 0) { fprintf(stderr, "-n argument should be positive\n"); usage(); exit(-1); } else { o_max_match = i1; // printf("opt max_match %d\n", o_max_match); } } else if (np = sscanf(argv[i], "-seq%lld", &p1), np == 1) { cmd_seq(p1); arg_warn = argv[i]; } else if (np = sscanf(argv[i], "MCS%lld", &p1), np == 1) { cmd_seq(p1); arg_warn = argv[i]; } else if (strncmp(argv[i], "horadam", strlen(argv[i])) == 0) { o_horadam = 1; } else if (strncmp(argv[i], "-w", 2) == 0) { s = argv[i]; n = str_s64(s+2, 0); if ((n >= WEIGHT_ORIG) && (n <= WEIGHT_LOG)) { o_weight = n; cache_weights(); // printf("opt weight %d\n", o_weight); } } else if (strncmp(argv[i], "-D", 2) == 0) { s = argv[i]; o_debug = str_s64(s+2, 0); printf("opt debug %d\n", o_debug); } else if (strcmp(argv[i], "-r") == 0) { o_format = FORMAT_RHTF; init_format(); } else if (strcmp(argv[i], "-r1") == 0) { o_format = FORMAT_RHTF; o_rhtf_simpler = 1; init_format(); } else if (strcmp(argv[i], "-v") == 0) { o_variations = 1; } else if ( ((argv[i])[0] == '-') || ((argv[i])[0] == ',') || isdigit((argv[i])[0]) ) { char *s; char *r; s = argv[i]; r = 0; if (*s == ',') { s++; } while(s && (*s) && (g_match_len < SCAN_LIMIT)) { mt = str_s64(s, &r); if (mt < 0) { o_negokay = 1; } if (o_negokay || (mt > 0)) { match[g_match_len++] = mt; /* printf("got term %lld leaving %s\n", mt, r); */ } s = 0; if (r && (*r == ',')) { r++; if ((*r == '-') || (isdigit(*r))) { s = r; } } } } else if (strcmp(argv[i], ",") == 0) { /* We also allow bare commas */ } else { fprintf(stderr, "Unrecognized argument '%s'\n", argv[i]); exit(-1); } } if (arg_warn) { exit(0); } if (o_horadam) { if (g_match_len >= 4) { printf("Horadam sequence Wn(%lld, %lld; %lld, %lld) :\n", match[0], match[1], match[2], match[3]); cmd_horadam(match[0], match[1], match[2], match[3]); exit(0); } else { fprintf(stderr, "Horadam requires 4 integer arguments\n"); exit(-1); } } if (o_debug & 1) { printf("Searching for: "); for(i=0; i 60)) { range = 1 + (i-60)/7; } else { range = 1; } //printf("scan %d range %d\n", i, range); scan1(i, range); i = i + range; } if (g_matches_found < o_max_match) { printf("Search ended after %d matches, because complexity", g_matches_found); if (o_weight == WEIGHT_LOG) { printf(" %d.%d exceeded limit %d.%d.\n", i/10, i%10, o_cplx_limit/10, o_cplx_limit%10); } else { printf(" %d exceeded limit %d.\n", i, o_cplx_limit); } } if (o_debug & 1) { printf("Total: %lld for %lld", g_sc2, g_sc1); if (g_sc1) { printf(" (%3.2g %%)", 100.0 * ((double) g_sc2) / ((double) (g_sc1))); } printf(" maxpass == %lld\n", maxpass); } return 0; } /* end of mcsfind.c */