Python search algorithm for OEIS integer sequence A059405

OEIS_A059405.py

# MIT License
#
# Copyright (c) 2023 Neil Webber
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Search for values in the sequence: https://oeis.org/A059405

import itertools
import math
import time


def _makeparams(n, /, *, base=10, zeropass=False):
    # This is the heart of the search.
    #
    # A list of multiples, called 'mults' is made for each digit 'd'
    # in the candidate value 'n'. All of these mults lists are gathered
    # into a list of lists, called, surprisingly enough, mults_lists.
    #
    # For example, if 'n' is 384 then the mults for digit '3' will be:
    #    [3, 9, 27, 81, 243]
    #
    # HOWEVER, as an important optimization, values that do not evenly
    # divide 'n' are eliminated; that optimization reduces the mults
    # for '3' to just [3] when 'n' is 384.
    #
    # Similarly, again with 'n' 384, the mults for 8 would be
    #   [8, 64]
    # and for 4 would be
    #   [4, 16, 64]
    #
    # A further optimization occurs for any digit 'd' that is present
    # more than once in the number. Since (by definition) all numbers
    # in the sequence must be powers of ALL their digits, instead of making
    # k identical mults for a digit that occurs 'k' times, one mult is
    # made with minimum value d**k. Think of this as a regular mult for 'd'
    # but including another d**1 for each additional occurrence of 'd'
    #
    # For example, in 6144 the mults for '4' are:
    #
    #   [16, 64, 256, 1024]
    #
    # reflecting both the start value optimization (4*4 = 16 in this case)
    # and filtering out 4096 because 4064 does not divide 6144.
    #
    # No mults are made for 1 digits (because: obvious reason).
    #
    # If zeropass is False (default), any number that has any zero digits
    # is rejected (because cannot be in the sequence, by definition).
    # If zeropass is True, zero digits are ignored, which results in many
    # more numbers being in the (modified) sequence.
    #

    # first find the raw digits, with repeats
    xn = n
    rawdigs = []
    while xn > 0:
        xn, d = divmod(xn, base)
        # note that numbers containing zeros cannot be in the sequence
        # so just abandon this candidate, and 1's are multiplicative
        # identities and thus do not need to be in the search space
        if d == 0 and not zeropass:
            return None
        elif d > 1:
            # if the digit doesn't divide n evenly, n can't be in sequence.
            # Note: this "optimization" actually is slower for numbers
            # that are in sequence (especially if repeated digits because
            # this test runs multiple times in those cases) but for numbers
            # out of the sequence it short-circuits most of the code and
            # since they dominate this is much much faster this way.
            # Without this test here searching for numbers is ~2.5x slower (!)
            if (n % d) != 0:
                return None

            # this is a candidate digit. NOTE: duplicates handled below.
            rawdigs.append(d)

    # Make the mults lists as described above.
    mults_lists = []
    already_seen = set()
    for d in rawdigs:
        if d not in already_seen:
            already_seen.add(d)
            mult = d ** rawdigs.count(d)
            mults = []
            while mult <= n:
                if n % mult == 0:
                    mults.append(mult)
                    mult *= d
                else:
                    break

            # mults is now 1 or more multiples of d; put it in mults_lists
            # NOTE: mults is known not-empty at this point because it
            #       was already tested that d divides n.
            mults_lists.append(mults)

    return mults_lists or None      # change [] to None on general principles


# optimization makes simple things less clear, but this is a substantial
# improvement by avoiding the _makeparams work for numbers that have zero
# digits in them (which cannot possibly be in sequence). This avoids even
# taking any look at all at quite a few of those. This is used, when
# possible, in lieu of itertools.count()
#

def fewerzeros(start=0):

    bunch = 1000     # hardwired, empirically chosen
    g = itertools.count(start)
    if start % bunch != 0:
        while start % bunch != 0:
            yield (start := next(g))

    while start < bunch:
        yield (start := next(g))

    # precomputed. Note that everything < 100 has a zero in it and
    # can be ignored because the "bunch" is 1000 and the above made
    # sure that start is currently at least 1000

    nzbunch = [111, 112, 113, 114, 115, 116, 117, 118,
               119, 121, 122, 123, 124, 125, 126, 127, 128, 129, 131,
               132, 133, 134, 135, 136, 137, 138, 139, 141, 142, 143,
               144, 145, 146, 147, 148, 149, 151, 152, 153, 154, 155,
               156, 157, 158, 159, 161, 162, 163, 164, 165, 166, 167,
               168, 169, 171, 172, 173, 174, 175, 176, 177, 178, 179,
               181, 182, 183, 184, 185, 186, 187, 188, 189, 191, 192,
               193, 194, 195, 196, 197, 198, 199, 211, 212, 213, 214,
               215, 216, 217, 218, 219, 221, 222, 223, 224, 225, 226,
               227, 228, 229, 231, 232, 233, 234, 235, 236, 237, 238,
               239, 241, 242, 243, 244, 245, 246, 247, 248, 249, 251,
               252, 253, 254, 255, 256, 257, 258, 259, 261, 262, 263,
               264, 265, 266, 267, 268, 269, 271, 272, 273, 274, 275,
               276, 277, 278, 279, 281, 282, 283, 284, 285, 286, 287,
               288, 289, 291, 292, 293, 294, 295, 296, 297, 298, 299,
               311, 312, 313, 314, 315, 316, 317, 318, 319, 321, 322,
               323, 324, 325, 326, 327, 328, 329, 331, 332, 333, 334,
               335, 336, 337, 338, 339, 341, 342, 343, 344, 345, 346,
               347, 348, 349, 351, 352, 353, 354, 355, 356, 357, 358,
               359, 361, 362, 363, 364, 365, 366, 367, 368, 369, 371,
               372, 373, 374, 375, 376, 377, 378, 379, 381, 382, 383,
               384, 385, 386, 387, 388, 389, 391, 392, 393, 394, 395,
               396, 397, 398, 399, 411, 412, 413, 414, 415, 416, 417,
               418, 419, 421, 422, 423, 424, 425, 426, 427, 428, 429,
               431, 432, 433, 434, 435, 436, 437, 438, 439, 441, 442,
               443, 444, 445, 446, 447, 448, 449, 451, 452, 453, 454,
               455, 456, 457, 458, 459, 461, 462, 463, 464, 465, 466,
               467, 468, 469, 471, 472, 473, 474, 475, 476, 477, 478,
               479, 481, 482, 483, 484, 485, 486, 487, 488, 489, 491,
               492, 493, 494, 495, 496, 497, 498, 499, 511, 512, 513,
               514, 515, 516, 517, 518, 519, 521, 522, 523, 524, 525,
               526, 527, 528, 529, 531, 532, 533, 534, 535, 536, 537,
               538, 539, 541, 542, 543, 544, 545, 546, 547, 548, 549,
               551, 552, 553, 554, 555, 556, 557, 558, 559, 561, 562,
               563, 564, 565, 566, 567, 568, 569, 571, 572, 573, 574,
               575, 576, 577, 578, 579, 581, 582, 583, 584, 585, 586,
               587, 588, 589, 591, 592, 593, 594, 595, 596, 597, 598,
               599, 611, 612, 613, 614, 615, 616, 617, 618, 619, 621,
               622, 623, 624, 625, 626, 627, 628, 629, 631, 632, 633,
               634, 635, 636, 637, 638, 639, 641, 642, 643, 644, 645,
               646, 647, 648, 649, 651, 652, 653, 654, 655, 656, 657,
               658, 659, 661, 662, 663, 664, 665, 666, 667, 668, 669,
               671, 672, 673, 674, 675, 676, 677, 678, 679, 681, 682,
               683, 684, 685, 686, 687, 688, 689, 691, 692, 693, 694,
               695, 696, 697, 698, 699, 711, 712, 713, 714, 715, 716,
               717, 718, 719, 721, 722, 723, 724, 725, 726, 727, 728,
               729, 731, 732, 733, 734, 735, 736, 737, 738, 739, 741,
               742, 743, 744, 745, 746, 747, 748, 749, 751, 752, 753,
               754, 755, 756, 757, 758, 759, 761, 762, 763, 764, 765,
               766, 767, 768, 769, 771, 772, 773, 774, 775, 776, 777,
               778, 779, 781, 782, 783, 784, 785, 786, 787, 788, 789,
               791, 792, 793, 794, 795, 796, 797, 798, 799, 811, 812,
               813, 814, 815, 816, 817, 818, 819, 821, 822, 823, 824,
               825, 826, 827, 828, 829, 831, 832, 833, 834, 835, 836,
               837, 838, 839, 841, 842, 843, 844, 845, 846, 847, 848,
               849, 851, 852, 853, 854, 855, 856, 857, 858, 859, 861,
               862, 863, 864, 865, 866, 867, 868, 869, 871, 872, 873,
               874, 875, 876, 877, 878, 879, 881, 882, 883, 884, 885,
               886, 887, 888, 889, 891, 892, 893, 894, 895, 896, 897,
               898, 899, 911, 912, 913, 914, 915, 916, 917, 918, 919,
               921, 922, 923, 924, 925, 926, 927, 928, 929, 931, 932,
               933, 934, 935, 936, 937, 938, 939, 941, 942, 943, 944,
               945, 946, 947, 948, 949, 951, 952, 953, 954, 955, 956,
               957, 958, 959, 961, 962, 963, 964, 965, 966, 967, 968,
               969, 971, 972, 973, 974, 975, 976, 977, 978, 979, 981,
               982, 983, 984, 985, 986, 987, 988, 989, 991, 992, 993,
               994, 995, 996, 997, 998, 999]

    while True:
        for d in nzbunch:
            yield start + d
        start += bunch
        # if start is a multiple of 100*bunch or 10*bunch then even more
        # skipping is possible
        if (start % (1000*bunch)) == 0:
            start += 111*bunch
        elif (start % (100*bunch)) == 0:
            start += 11*bunch
        elif (start % (10*bunch)) == 0:
            start += bunch


def A059405(*, start=1, base=10, zeropass=False):
    """Generate numbers in the A059405 sequence.

    Keyword arguments:
    start    -- First value to search. Default 1.
    base     -- Number base, for generating digits. Default 10.
                NOTE: Minimum 3. No number in binary can be in the sequence
                      (because of how the sequence is defined).
    zeropass -- If True, zero digits are ignored, which tremendously
                increases the set of numbers in the (modified) sequence.
                Default is False, which yields correct A059405 values.
    """
    if base == 2:               # in binary there are no numbers in sequence
        return
    if base < 2:                # C'mon man
        raise ValueError("base must be 2 or more")

    countg = itertools.count(start)
    if not (zeropass or (base != 10)):
        countg = fewerzeros(start)

    if start <= 1:
        yield 1
    for n in countg:
        # the hard work is all in _makeparams; all this does is
        # iterate through the combinatoric product of all the mults
        mults_lists = _makeparams(n, base=base, zeropass=zeropass)
        if mults_lists:
            for candy in itertools.product(*mults_lists):
                if math.prod(candy) == n:
                    yield n
                    break


def is_A059405(n, base=10):
    """Test a single number for membership in the sequence."""

    # 1 is a special case
    if n == 1:
        return True

    mults_lists = _makeparams(n, base=base)
    if mults_lists:
        for candy in itertools.product(*mults_lists):
            if math.prod(candy) == n:
                return True
    return False


def find_sequence_numbers(i=0, /, *, g=None):
    """Convenient function for searching for sequence values."""
    if g is None:
        g = A059405(start=i)
    prev = 0
    while True:
        t0 = time.time()
        a = next(g)
        t1 = time.time()
        nsecs = 1000 * 1000 * (t1 - t0) / (a - prev)
        print(f"Found: {a}, elapsed={t1 - t0:.2f} ({nsecs:.2f}nsec per eval)")
        prev = a


if __name__ == "__main__":
    import unittest

    class TestMethods(unittest.TestCase):

        def test_quicksequence(self):
            # this data cut-and-pasted from https://oeis.org/A059405
            # but cut off so as not to take too long to run test
            testvec = (
                1, 2, 3, 4, 5, 6, 7, 8, 9, 128, 135, 175, 384, 432, 672,
                735, 1296, 1715, 6144, 6912, 13824, 18432, 23328, 34992,
                82944, 93312, 131712, 248832, 442368, 1492992, 2239488,
                2333772, 2612736)
            g = A059405()
            for t in testvec:
                with self.subTest(t=t):
                    self.assertEqual(t, next(g))

        def test_makeparams(self):
            # these examples were hand-computed
            # each test vector is a tuple: (n, makeparamsresults)
            testvec = ((123, None),
                       # 2232:
                       #   2,2,2 -> [8]. Note that 16 and beyond don't divide
                       #   3 -> 3, 9
                       (2232, [[8], [3, 9]]),

                       # 111132:
                       #   2 -> [2, 4]. 8 and beyond don't divide
                       #   3 -> [3, 9, 27, 81]. 243 doesn't divide

                       (111132, [[2, 4], [3, 9, 27, 81]]),

                       # 8136:
                       #   8 -> [8]. 64 does not divide.
                       #   3 -> [3, 9]. 27 does not divide.
                       #   6 -> [6, 36]. 216 does not divide.
                       (8136, [[8], [3, 9], [6, 36]]))

            for n, rslt in testvec:
                with self.subTest(n=n, rslt=rslt):
                    m = _makeparams(n)
                    # put the mults lists into canonical order otherwise
                    # testvecs would have to be constructed knowing
                    # implementation details of _makeparams
                    if m is not None:
                        m = sorted([sorted(x) for x in m])
                    if rslt is not None:
                        rslt = sorted([sorted(x) for x in rslt])
                    self.assertEqual(m, rslt)

        def test_makeparams2(self):
            # semantic tests of a makeparams result. Likely redundant
            # given the other tests, but here it is anyway. For every
            # value in a _makeparams result, make sure that it divides.
            for i in range(1000):
                p = _makeparams(i)
                if p is not None:
                    for m in p:
                        for dv in m:
                            self.assertEqual(i % dv, 0)

        def testzeropass(self):
            g = A059405(start=1250, zeropass=True)
            self.assertEqual(next(g), 1250)

            g = A059405(start=1250, zeropass=False)
            self.assertNotEqual(next(g), 1250)

            g = A059405(start=1250)
            self.assertNotEqual(next(g), 1250)

        def testmore(self):
            # this data cut-and-pasted from
            #     https://oeis.org/A059405/b059405.txt
            # and auto-transformed into this list
            #
            testvec = (
                3981312, 4128768, 4741632, 9289728, 12192768, 13226976,
                13395375, 13436928, 13716864, 21233664, 21337344, 22127616,
                22478848, 27433728, 27869184, 28449792, 49787136)
            g = A059405(start=testvec[0])
            for t in testvec:
                with self.subTest(t=t):
                    self.assertEqual(t, next(g))

        def testfz(self):
            # test the "fewer zeros" count optimization
            g1 = itertools.count(1)
            g2 = fewerzeros(1)
            for i in range(10000000):
                # every result from itertools.count that has NO zeros
                # in it must come (in order) from fewerzeros. So this first
                # bit of code gets the next "no zeros" from g1
                for a in g1:
                    if '0' not in str(a):
                        break

                # now g2 should eventually return this value without
                # going beyond it
                for b in g2:
                    if '0' not in str(b):
                        break
                    self.assertTrue(b < a)

                self.assertEqual(a, b)

    unittest.main()

This new version which tests for divisibility earlier (in determining rawdigits in _makeparams) is 2.5x faster. Also added more tests.

Even moar faster, by skipping over many numbers containing zero digits (not even calling _makeparams on them)