Python implementation of the A133058 sequence of integers https://oeis.org/A133058

fly-straight-dammit.py

"""
Python implementation of the A133058 sequence of integers
https://oeis.org/A133058
"""

import functools

_PLOT = False
try:
    import matplotlib.pyplot as plt

    _PLOT = True
except ModuleNotFoundError:
    print(
        "matplotlib not installed, can't show graph. Run command: pip install matplotlib"
    )


INT_COUNT = 1000


def prime_factors(number):
    if number < 4:
        return [number]

    limit = int(number / 2) + 1
    factors = []
    for i in range(2, limit + 1):
        while number % i == 0:
            factors.append(i)
            number = int(number / i)

    if number > 1:
        return [number]

    return factors


def get_gcd(a, b):
    common_factors = [1]
    b_factors = prime_factors(b)
    for factor in prime_factors(a):
        if factor in b_factors:
            common_factors.append(factor)
            b_factors.remove(factor)
    return functools.reduce(lambda acc, i: acc * i, common_factors)


def fly_straight_dammit():
    """
    Returns a generator that sends out a series of tuples (n, a[n])
    where a[n] is the nth integer in the sequence A133058
    """
    yield 0, 1
    yield 1, 1

    i = 2
    next_num = 1
    while True:
        gcd = get_gcd(i, next_num)
        if gcd == 1:
            next_num += i + 1
        else:
            next_num = int(next_num / gcd)

        yield i, next_num
        i += 1


if __name__ == "__main__":
    x = []
    y = []

    for n in fly_straight_dammit():
        if n[0] > INT_COUNT:
            break

        x.append(n[0])
        y.append(n[1])

    if _PLOT:
        plt.scatter(x, y, s=1)
        plt.show()
    else:
        from pprint import pprint

        pprint(y[610:660])