OEIS A109732 and (related) A261414

gistfile1.txt

# 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/A109732

import itertools


def A109732():
    """Generate numbers in the A109732 sequence."""

    # From OEIS A109732:
    #   a(1) = 1
    #   For n > 1, a(n) is the smallest number not already present
    #   which is entailed by the rules:
    #          (i) k present => 2k+1 present;
    #          (ii) 3k present => k present.

    a = {1}
    yield 1
    biggest = 1
    while True:
        # The biggest value found so far would generate (2*biggest) + 1
        # as the next value, which by definition isn't in the sequence
        # (or it would be the biggest). Therefore this is a sufficiently-large
        # starting point for the next value that would be "smallest" in
        # the sequence (and not already in the sequence), plus 1 (so that
        # if it ends up being found, it is smaller)
        smallest = (2*biggest) + 1
        for ka in a:
            k, r = divmod(ka, 3)
            if r == 0 and k not in a:
                if k < smallest:
                    smallest = k
            else:
                # make the next 2k+1 entry
                k = (ka*2) + 1
                if k not in a and k < smallest:
                    smallest = k
        a.add(smallest)
        yield smallest
        if biggest < smallest:
            biggest = smallest


def A109732_oddconjecture(n):
    # Unproven conjecture: every odd number eventually appears.
    # This tests that conjecture for all odd numbers below n.
    # Success = this function returns. Failure = "it takes forever"
    g = A109732()
    x = set(2*k + 1 for k in range(n//2))
    while len(x) > 0:
        k = next(g)
        if k in x:
            x.remove(k)


def A109732_oddconjecture_gen(*, report=100):
    # Like .._oddconjecture, but a generator that yields a result
    # each time the next "report" odd numbers are found. For example,
    # if report==100 (the default), then this yields 100 as soon as
    # all the odd numbers 1, 3, 5 ... 99 have been found, then yields
    # 200 once 101, 103, .. 199 have been found, etc.
    #
    # NOTE (OF COURSE) THAT EVENTUALLY THIS TAKES VERY LONG TO YIELD RESULTS
    #
    g = A109732()
    found = set()

    for nth in itertools.count():
        group = set(k
                    for k in range((nth*report)+1, (nth+1)*report, 2)
                    if k not in found)
        while len(group) > 0:
            k = next(g)
            if k in group:
                group.remove(k)
            if k in found:         # repeats should never happen
                raise TypeError(f"{k} duplicated!")
            found.add(k)
        yield (nth+1) * report, len(found)


# not a generator - returns JUST the n'th value of the sequence
def A261414_n(n):
    """A261414: 2^n+1 appears in A109732 at position a(n)"""

    # CAUTION: time taken to compute increases rapidly with n
    needed = (2**n) + 1
    for a732pos, c in enumerate(A109732(), start=1):
        if c == needed:
            return a732pos


if __name__ == "__main__":
    import unittest

    class TestMethods(unittest.TestCase):

        testvec = (1, 3, 7, 15, 5, 11, 23, 31, 47, 63, 21, 43, 87, 29, 59,
                   95, 119, 127, 175, 191, 239, 255, 85, 171, 57, 19, 39,
                   13, 27, 9, 55, 79, 111, 37, 75, 25, 51, 17, 35, 71, 103,
                   115, 143, 151, 159, 53, 107, 207, 69, 139, 215, 223, 231,
                   77, 155, 279, 93, 187, 287, 303, 101, 203)

        def test_seq(self):
            g = A109732()
            for t in self.testvec:
                self.assertEqual(t, next(g))

        oddtest_reasonable = 1000

        def test_odd(self):
            A109732_oddconjecture(self.oddtest_reasonable)

        def test_A261414(self):
            tv = (2, 5, 30, 38, 201, 242, 689, 1806, 7175)
            for i, t in enumerate(tv, start=1):
                self.assertEqual(t, A261414_n(i))

    unittest.main()