// C++ program by Bert Dobbelaere for computation of the following OEIS sequences: // // A304241 - Indices where Mertens function A002321 reaches record amplitudes between zeros. // A304242 - Increasingly larger (in absolute value) extrema of the Mertens function A002321 between subsequent zeros. // // Note: 64 bit integer support is required. // // Operating principles: // In order to compute the sequence values above, we need all the Mertens's function values (A002321) up to the next // zero crossing beyond the highest value of A304241 reported. This program uses a single-pass approach that uses the // Mertens's function values once in ascending order of its arguments. This means that it is sufficient to compute the // Moebius function values (A008683) in ascending order and just keep a running sum. // The whole point then becomes: how to quickly compute the Moebius values ? // Factoring all its arguments blindly is very slow. // This program uses a method that is similar to the Sieve of Eratosthenes, only we store more information // per number. We can use an array up to some max value (MAX_K) for the Moebius function argument, and initialize // arr[k]=1 for (k=1..MAX_K). // Starting from 2, we do the following for all primes p so that p*p<=MAX_K: // for all k that are multiples of p (starting from 2p), so that k<=MAX_K: // arr[k]:= -arr[k]*p // _except_ for all multiples of p^2: // arr[k]:= 0 // All square-free elements for non-prime values of k are now containing the product of their prime factors up to sqrt(MAX_K), // and the sign of arr[k] shows how many primes were multiplied (mod 2). As there can be at most one factor that is // larger than sqrt(MAX_K), we just can detect it using the value of the resulting arr[k]: abs(arr[k])= BASE_TABLE_SIZE #define DO_SANITY_CHECK // If this is defined, report values of A084237 on the fly. // A084237(n) = M(10^n), where M(n) is Mertens's function. This should reduce the risk of incorrect // output due to software or hardware bugs, electrical glitches, alpha particles... // Output with defaults above: // // n A304241(n) A304242(n) Sanity check : A084237(0) must be 1 // 1 1 1 Sanity check : A084237(1) must be -1 // 2 31 -4 Sanity check : A084237(2) must be 1 // 3 114 -6 // 4 199 -8 // 5 443 -9 // 6 665 -12 Sanity check : A084237(3) must be 2 // 7 1109 -15 // 8 1637 -16 // 9 2803 -25 // 10 7021 -29 // 11 8511 35 Sanity check : A084237(4) must be -23 // 12 9861 -43 // 13 19291 51 // 14 24185 -72 // 15 31990 73 // 16 42961 -88 // 17 48433 96 // 18 59577 -113 Sanity check : A084237(5) must be -48 // 19 96014 -132 // 20 141869 -134 // 21 230399 -154 // 22 300551 240 // 23 355733 -258 // 24 603151 -278 // 25 926265 -368 Sanity check : A084237(6) must be 212 // 26 1066854 432 // 27 1793918 550 // 28 3239797 -683 // 29 5343761 -847 // 30 6481601 1060 // 31 7109110 -1078 Sanity check : A084237(7) must be 1037 // 32 10194458 1240 // 33 12874814 -1447 // 34 30919091 -2573 // 35 61913863 2845 // 36 76015339 -3448 Sanity check : A084237(8) must be 1928 // 37 110103729 -4610 // 38 179919749 -6226 // 39 360718458 6695 // 40 456877618 -8565 // 41 532476623 9132 // 42 903087703 10246 Sanity check : A084237(9) must be -222 // 43 1109331447 -15335 // 44 2597217086 -17334 // 45 3314385678 21777 // 46 3773166681 -25071 // 47 6631245058 -31206 // 48 7766842813 50286 Sanity check : A084237(10) must be -33722 // 49 15578669387 -51116 // 50 19890188718 60442 // 51 22867694771 -62880 // 52 38066335279 -81220 // 53 56808201767 -87995 Sanity check : A084237(11) must be -87856 // 54 101246135617 -129332 // 55 108924543546 170358 // 56 217309283735 -190936 // 57 297193839495 207478 // 58 330508686218 -294816 // 59 661066575037 331302 Sanity check : A084237(12) must be 62366 // 60 1252680108413 -345947 // 61 1440355022306 -368527 // 62 1653435193541 546666 // 63 2087416003490 -625681 // 64 3270926424607 -635558 // 65 3787401700373 661998 #include #include #include #include #include #if __cplusplus < 201103L // Types below work with GCC. We need 64 bit signed and unsigned integers. typedef unsigned long long u64; typedef unsigned u32; typedef long long s64; #else // C++11 has standard types #include typedef uint64_t u64; typedef uint32_t u32; typedef int64_t s64; #endif #define MAX_K ((u64)BASE_TABLE_SIZE*(u64)BASE_TABLE_SIZE) // Evaluate Mertens's function up to this argument struct BaseEntry{ u64 p; // prime u64 pmult; // current multiple of p (mod p) u64 idx; // current offset of prime multiple in overlay }; u64 n_base_primes; // Number of primes in base table // Table pointers: BaseEntry *mu_base_table; s64 *mu_overlay_table; // Current overlay boundaries: u64 overlay_window_start = 0u; u64 overlay_window_stop = 0u; // Shift the overlay window and fill it void update_overlay() { overlay_window_start = overlay_window_stop; overlay_window_stop = overlay_window_start + OVERLAY_SIZE; for(u64 k=0;k=p) sqrcount-=p; } b.pmult = sqrcount; // Update base entry for next call b.idx = idx-OVERLAY_SIZE; } } // A008683 Moebius function mu(n) - fast version taylored to our needs. // It is designed for sequential access (every incrementing k) int mu(u64 k) { if(k>=overlay_window_stop) // Out of values { update_overlay(); // Shift the overlay and fill it } s64 el=mu_overlay_table[k - overlay_window_start]; if(el<0) { if(-el < (s64)k) // one remaining prime return 1; else return -1; } if(el>0) { if(el<(s64)k) // one remaining prime return -1; else return 1; } return 0; } // Initialize Moebius function tables void init_mu_table() { // Allocate both tables. 1st alloc is too big (never mind), we only need it to store // an entry for each prime <= BASE_TABLE_SIZE mu_base_table=(BaseEntry*)calloc(BASE_TABLE_SIZE/2,sizeof(s64)); mu_overlay_table=(s64*)calloc(OVERLAY_SIZE,sizeof(s64)); assert(mu_base_table!=NULL); // Sufficient memory ? assert(mu_overlay_table!=NULL); // Sufficient memory ? // Initialize base table n_base_primes=1; mu_base_table[0].p=2; // Treat 2 separately mu_base_table[0].pmult=0; // 2%2=0 mu_base_table[0].idx=2*2; for(u64 k=3;k<=BASE_TABLE_SIZE;k+=2) { // Rudimentary primality check - only needed during init bool isprime=true; for(u64 d=3;d*d<=k;d+=2) { if(k%d==0) { isprime=false; break; } } if(isprime) { mu_base_table[n_base_primes].p=k; mu_base_table[n_base_primes].pmult=2; mu_base_table[n_base_primes].idx=2*k; n_base_primes++; } } // Overlay table will be filled when we first need it } int main() { int n=0; // Sequence index of output int mertens_k=0; // Tracks A002321 - Mertens's function int abs_max_rpt=0; // maximum abs already reported u64 extr_idx=0; // Index of most extreme value so far int extr_val=0; // Most extreme value int abs_max=0; // Most extreme value (absolute) assert(sizeof(u64)==8); // Just be sure: 64 bit assert(sizeof(s64)==8); // Just be sure: 64 bit assert(OVERLAY_SIZE>=BASE_TABLE_SIZE); assert(u64(BASE_TABLE_SIZE)<3037000000ull); // Sanity check: square must fit in s64 type init_mu_table(); // Create base table and first overlay std::cout << " n A304241(n) A304242(n)"; #ifdef DO_SANITY_CHECK u64 next_pow10=1; int next_log=0; #endif for(u64 k=1;k<=MAX_K;k++) { mertens_k+=mu(k); // Update Mertens's function with new Moebius value int a=mertens_k<0 ? -mertens_k : mertens_k; if(a>abs_max) { abs_max=a; extr_val=mertens_k; extr_idx=k; } if( (mertens_k==0) && (abs_max>abs_max_rpt)) // Report record reached before this zero crossing { n+=1; std::cout << std::endl << std::setw(3) << n << std::setw(22) << extr_idx <