question_marks.py

import numpy as np


def stern_sequence(n):
    import numpy as np
    a = np.array([0, 1], dtype=np.int64)
    for _ in range(1, n):
        a = np.vstack((a, np.roll(a, -1) + a)).reshape((-1,), order='F')
        a[-1] += 1
    return a


def bit_reversal_permutation(n):
    import numpy as np
    a = np.arange(2 ** n, dtype=np.uint64)
    r = 6
    total_dist = 2 ** r
    for k in range(0, r):
        dist = 2 ** k
        mask = (2 ** total_dist - 1) // (2 ** dist + 1)
        a = ((a >> dist) & mask) | ((a & mask) << dist)

    a >>= total_dist - n
    return a


def main():
    # P for Legendre polynomials
    # T for Chebyshev polynomials of the first kind
    # U for Chebyshev polynomials of the second kind
    polynomial_type = 'T'
    # Compose with either the box, question mark, or identity functions: box, ?, id
    compose = '?'
    out_filename = 'out.png'

    n1 = 16
    w = 2**n1

    curve_count = 2**7
    lw = 0.1

    ss = stern_sequence(n1 + 2)
    br = bit_reversal_permutation(n1 + 1)
    if compose == '?':
        x0 = ss[br[:2**n1]].astype(np.float64) / ss[br[:2**n1] + 1].astype(np.float64)
        x = np.arange(2**n1) / 2**n1
    elif compose == 'box':
        x0 = np.arange(2**n1) / 2**n1
        x = ss[br[:2**n1]].astype(np.float64) / ss[br[:2**n1] + 1].astype(np.float64)
    elif compose == 'id':
        x0 = np.arange(2**n1) / 2**n1
        x = np.arange(2**n1) / 2**n1
    else:
        raise

    pts = np.empty((3, len(x)), dtype=np.float64)
    pts[0] = 1
    if polynomial_type in {'T', 'P'}:
        pts[1] = x
    elif polynomial_type in {'U'}:
        pts[1] = 2 * x
    else:
        raise

    import matplotlib.pyplot as plt
    plt.clf()
    plt.close()
    plt.figure()
    for k in range(min(2, curve_count)):
        plt.plot(x0, pts[k], linewidth=lw)
    for k in range(2, curve_count):
        if polynomial_type in {'P'}:
            pts[-1] = ((2 * k - 1) * x * pts[1] - (k - 1) * pts[0]) / k  # Legendre
        elif polynomial_type in {'U', 'T'}:
            pts[-1] = 2 * x * pts[1] - pts[0]  # Chebyshev

        pts[:-1] = pts[1:]

        plt.plot(x0, pts[-1], linewidth=lw)
        print(f'{k/(curve_count-1):2.2%} complete')

    plt.savefig(out_filename, dpi=600)


if __name__ == '__main__':
    main()