#include <iostream> #include <cstdlib> #include <unistd.h> #include <vector> #include <gmpxx.h> #undef TRACE_SEARCH #undef TRACE_MODS #undef TRACE_TESTS struct msearch_data { size_t ndigits; unsigned long base; bool include_pseudo; bool reverse_meertens; std::vector<mpz_class> base_powers; std::vector<mpz_class> work_limits; std::vector<unsigned long> primes; std::vector<std::vector<mpz_class> > prime_powers; size_t node_count; size_t leaf_count; size_t test_count; std::vector<size_t> depth_counts; size_t status_interval; size_t next_status; msearch_data (unsigned long _base) : ndigits (0), base(_base), include_pseudo(false), reverse_meertens(false), status_interval(0) { base_powers.push_back (mpz_class(1)); work_limits.push_back (mpz_class(1)); primes.push_back (2); prime_powers.push_back (gen_powers (primes.back())); reset_counts(); } void set_ndigits (size_t _ndigits) { ndigits = _ndigits; while (base_powers.size() <= ndigits) { base_powers.push_back (base_powers.back() * base); work_limits.push_back ((4*base_powers.back())/base); } while (primes.size() < ndigits) { primes.push_back (nextprime ()); prime_powers.push_back (gen_powers (primes.back())); } reset_counts(); } void reset_counts() { depth_counts.resize (ndigits+1); node_count = 0; leaf_count = 0; test_count = 0; next_status = status_interval; for (std::vector<size_t>::iterator i = depth_counts.begin(); i != depth_counts.end(); ++i) { *i = 0; } } void report_counts (size_t nprefix, const mpz_class& n) { std::cout << "nodes:" << node_count << " leaves:" << leaf_count << " tests:" << test_count << " n<" << nprefix << ">:" << n << "\n"; std::cout << " depths:"; for (std::vector<size_t>::const_iterator i = depth_counts.begin(); i != depth_counts.end(); ++i) { std::cout << " " << *i; } std::cout << std::endl; } void report_counts () { std::cout << "TOTAL " << ndigits << " digits: nodes:" << node_count << " leaves:" << leaf_count << " tests:" << test_count << "\n"; std::cout << " depths:"; for (std::vector<size_t>::const_iterator i = depth_counts.begin(); i != depth_counts.end(); ++i) { std::cout << " " << *i; } std::cout << std::endl; } unsigned long nextprime () const { if (primes.size() == 0) { return 2; } else if (primes.size() == 1) { return 3; } else { unsigned long n = primes.back() + 2; while (!is_prime(n)) n += 2; return n; } } bool is_prime (unsigned long n) const { for (std::vector<size_t>::const_iterator i = primes.begin(); i != primes.end() && (*i)*(*i) <= n; ++i) { if (n % (*i) == 0) return false; } return true; } std::vector<mpz_class> gen_powers (unsigned long p) { std::vector<mpz_class> powers; mpz_class m = 1; powers.push_back (m); for (unsigned long i=1; i<base; i++) { m *= p; powers.push_back (m); } return powers; } // godel(n) = godel encoding of n, so if // n = d_1 d_2 d_3 ... d_k in base b, // godel(n) = 2^(d_1) 3^(d_2) 5^(d_3) ... p_k^(d_k) mpz_class godel(const mpz_class& _n) { mpz_class n = _n; std::vector<unsigned long> digits; mpz_class q; mpz_class r; while (n > 0) { unsigned long d = mpz_fdiv_qr_ui (q.get_mpz_t(), r.get_mpz_t(), n.get_mpz_t(), base); digits.push_back (d); n.swap(q); } mpz_class g = 1; size_t k=0; for (std::vector<unsigned long>::const_reverse_iterator i = digits.rbegin(); i != digits.rend(); ++i) { g *= prime_powers[k][*i]; ++k; } return g; } // revgodel(n) = godel encoding of n, so if // n = sum d_i b^i (0 <= d_i < b) // revgodel(n) = prod p_{i+1}^(d_i) // revgodel(n,offset) = revgodel(n*b^offset) mpz_class revgodel(const mpz_class& _n, size_t offset = 0) { mpz_class n = _n; mpz_class q; mpz_class r; mpz_class g = 1; size_t k=offset; while (n > 0) { unsigned long d = mpz_fdiv_qr_ui (q.get_mpz_t(), r.get_mpz_t(), n.get_mpz_t(), base); g *= prime_powers[k][d]; k++; n.swap(q); } return g; } // godel(n,k) = godel encoding of n as an k-digit number // (that is, perhaps with zeros prepended) mpz_class godel(const mpz_class& _n, size_t ndigits) { mpz_class n = _n; mpz_class q; mpz_class r; mpz_class g = 1; size_t k = ndigits-1; for (size_t d=0; d<ndigits && n > 0; d++) { unsigned long d = mpz_fdiv_qr_ui (q.get_mpz_t(), r.get_mpz_t(), n.get_mpz_t(), base); n.swap(q); g *= prime_powers[k][d]; --k; } return g; } // the key to the meertens searches are that if we divide a // ndigit meertens number n into a nprefix-digit prefix x and // a nsuffix-digit suffix y (so nprefix+nsuffix==ndigit), then // we have several properties: // 0 <= x < b^nprefix (because x has nprefix digits) // 0 <= y < b^nsuffix (because y has nsuffix digits) // n = x * b^nsuffix + y (because x and y are the two halves of n) // n = godel(n) (because n is a meertens number) // godel(n) = godel(x) * godel(y,ndigits) // (from the definition of godel()) // x * b^nsuffix + y = godel(x) * godel(y,ndigits) // (combining the last three equations) // and thus // x * b^nsuffix + y === 0 (mod godel(x)) // x * b^nsuffix + y === 0 (mod godel(y,ndigits)) // meertens_pref searches through prefixes x, and uses the congruence // modulo godel(x) to list possible candidate y's. // meertens_suff searches through suffixes y, and uses the congruence // modulo godel(y,ndigits) to list possible candidate x's // searching for reverse meertens numbers uses the same prefix/suffix // search structure. the only difference is that they're using the // revgodel() encodings // gx = godel(x) (or if reverse, gx = revgodel(x,nsuffix)) void meertens_pref (size_t nprefix, const mpz_class& x, const mpz_class& gx) { size_t nsuffix = ndigits - nprefix; #ifdef TRACE_SEARCH std::cout << "search nprefix " << nprefix << " x " << x << " gx " << gx << std::endl; #endif /* TRACE_SEARCH */ if (status_interval > 0 && node_count >= next_status) { report_counts(nprefix, x); next_status += status_interval; } node_count++; if (nsuffix == 0) { depth_counts[nprefix]++; if (x == gx) { std::cout << "FOUND " << x << std::endl; } } else if (gx > work_limits[nsuffix]) { leaf_count++; depth_counts[nprefix]++; // x * b^(nsuffix) + y = gx * godel(y,ndigits) // so x * b^(nsuffix) + y === 0 (mod gx) // or y === -x * b^(nsuffix) (mod gx) // and 0 <= y < b^(nsuffix) mpz_class bigx = x * base_powers[nsuffix]; mpz_class t = -bigx; mpz_class y; mpz_mod(y.get_mpz_t(), t.get_mpz_t(), gx.get_mpz_t()); #ifdef TRACE_MODS std::cout << "x " << x << " nprefix " << nprefix << " gx " << gx << " y " << y << " limit " << base_powers[nsuffix] << "\n"; #endif /* TRACE_MODS */ while (y < base_powers[nsuffix]) { test_count++; mpz_class gy = reverse_meertens ? revgodel(y) : godel(y, ndigits); #ifdef TRACE_TESTS std::cout << "test bigx " << bigx << " y " << y << " gx " << gx << " gy " << gy << "\n"; #endif /* TRACE_TESTS */ if (bigx + y == gx * gy) { std::cout << "FOUND " << (bigx + y) << std::endl; } y += gx; } } else { mpz_class newx = x * base; mpz_class newgx = gx; unsigned long p = reverse_meertens ? primes[nsuffix-1] : primes[nprefix]; for (unsigned long i=0; i<base && newgx < base_powers[ndigits]; ++i) { if (nprefix > 0 || i > 0 || include_pseudo) { meertens_pref (nprefix+1, newx, newgx); } newx++; newgx *= p; } } } // gy = godel(y,ndigits) (or if reverse_meertens, gy = revgodel(y)) void meertens_suff (size_t nsuffix, const mpz_class& y, const mpz_class& gy) { #ifdef TRACE_SEARCH std::cout << "search nsuffix " << nsuffix << " y " << y << " gy " << gy << std::endl; #endif /* TRACE_SEARCH */ size_t nprefix = ndigits - nsuffix; if (status_interval > 0 && node_count >= next_status) { report_counts(nsuffix, y); next_status += status_interval; } node_count++; if (nprefix == 0) { depth_counts[nsuffix]++; if (y == gy && (y >= base_powers[ndigits-1] || include_pseudo)) { std::cout << "FOUND " << y << std::endl; } } else if (gy > work_limits[nprefix]) { leaf_count++; depth_counts[nsuffix]++; // x * b^(nsuffix) + y = godel(x) * gy // so x * b^(nsuffix) + y === 0 (mod gy) // or x * b^(nsuffix) === -y (mod gy) // and 0 <= x < b^(nprefix) mpz_class gcd; mpz_class binv; mpz_gcdext (gcd.get_mpz_t(), binv.get_mpz_t(), NULL, base_powers[nsuffix].get_mpz_t(), gy.get_mpz_t()); mpz_class q; mpz_class r; mpz_fdiv_qr(q.get_mpz_t(),r.get_mpz_t(),y.get_mpz_t(),gcd.get_mpz_t()); // if r != 0, there are no solutions to the congruence if (r == 0) { mpz_class t = -q * binv; mpz_class gyg = gy/gcd; mpz_class x; mpz_mod(x.get_mpz_t(), t.get_mpz_t(), gyg.get_mpz_t()); #ifdef TRACE_MODS std::cout << "y " << y << " nsuffix " << nsuffix << " gy " << gy << " x " << x << " limit " << base_powers[nprefix] << "\n"; #endif /* TRACE_MODS */ while (x < base_powers[nprefix]) { if (include_pseudo || x >= base_powers[nprefix-1]) { test_count++; mpz_class gx = reverse_meertens ? revgodel(x, nsuffix) : godel(x, nprefix); #ifdef TRACE_TESTS std::cout << "test y " << y << " x " << x << " gy " << gy << " gyg " << gyg << " gx " << gx << "\n"; #endif /* TRACE_TESTS */ if (x * base_powers[nsuffix] + y == gx * gy) { std::cout << "FOUND " << (x * base_powers[nsuffix] + y) << std::endl; } } x += gyg; } } } else { mpz_class newy = y; mpz_class newgy = gy; unsigned long p = reverse_meertens ? primes[nsuffix] : primes[nprefix-1]; for (unsigned long i=0; i<base && newgy < base_powers[ndigits]; ++i) { meertens_suff (nsuffix+1, newy, newgy); newy += base_powers[nsuffix]; newgy *= p; } } } }; void usage (const char *fname) { std::cerr << "Usage: " << fname << " [flags below...]\n" << " -b base (default=10)\n" << " -p [use prefix search]\n" << " -s [use suffix search (default)]\n" << " -r [reverse-meertens]\n" << " -P [include pseudo-meertens numbers]\n" << " -N d [only search d-digit numbers]\n" << " -n d [only search <= d-digit numbers]\n" << " -l d [only search >= d-digit numbers]\n" << " -S n [report status every n nodes (default 10^7)]\n"; exit (1); } int main (int ac, char **av) { extern int optind; extern char *optarg; bool pref_search = false; bool suff_search = false; bool reverse_meertens = false; bool include_pseudo = false; unsigned long base = 10; size_t ndigits_hi = 0; size_t ndigits_lo = 1; size_t status_interval = 10000000; int c; while ((c = getopt (ac, av, "b:psPN:n:l:S:r")) != EOF) { switch (c) { case 'b': base = atoi(optarg); break; case 'p': pref_search = true; break; case 's': suff_search = true; break; case 'r': reverse_meertens = true; break; case 'P': include_pseudo = true; break; case 'N': ndigits_lo = ndigits_hi = atoi(optarg); break; case 'n': ndigits_hi = atoi(optarg); break; case 'l': ndigits_lo = atoi(optarg); break; case 'S': status_interval = atoi(optarg); break; default: usage (av[0]); } } if (optind < ac) { usage (av[0]); } msearch_data md(base); md.status_interval = status_interval; md.include_pseudo = include_pseudo; md.reverse_meertens = reverse_meertens; if (pref_search == suff_search) { // either both false or both true! pref_search = reverse_meertens; suff_search = !reverse_meertens; } for (size_t ndigits=ndigits_lo; ndigits_hi == 0 || ndigits <= ndigits_hi; ++ndigits) { std::cout << "Searching for " << ndigits << " digit base " << base << (reverse_meertens ? " reverse meertens" : " meertens") << " numbers..." << std::endl; md.set_ndigits(ndigits); mpz_class n = 0; mpz_class gn = 1; if (suff_search) { md.meertens_suff (0, n, gn); } else { md.meertens_pref (0, n, gn); } md.report_counts(); } return 0; }