In [3]:
%matplotlib inline
from functools import cache
import numpy as np
import matplotlib.pyplot as plt
import itertools
import math
from decimal import Decimal, getcontext
In [99]:
def make_pqzo(p1, p2, p3):
    @cache
    def p(i, k):
        if k < 0:
            raise ValueError("p called with k < 0")
        if k == 0:
            return [p1, p2, p3][i - 1]
        return p(i, k - 1) * (p(i, k - 1) + 2 * p(i % 3 + 1, k - 1))
    
    @cache
    def q_odd(i, k):
        return q(i, k - 1) * (q(i, k - 1) + 2 * q(i % 3 + 1, k - 1))
    
    @cache
    def z(i, k):
        return p(i, k) * (p(i, k) + 2 * p((i+1) % 3 + 1, k))
    
    @cache
    def q(i, k):
        if k < 0:
            raise ValueError("q called with k < 0")
        if k == 0:
            return p(i, 0) * (p(i, 0) + 2 * p((i + 1) % 3 + 1, 0))
        return (
            q_odd(i, k) * (z(i, k) + z(i % 3 + 1, k))
            + q_odd(i % 3 + 1, k) * z(i, k)
        )
    
    return (p, q, z, q_odd)
In [2]:
def make_pq(p1, p2, p3):
    @cache
    def p(i, k):
        if k < 0:
            raise ValueError("p called with k < 0")
        if k == 0:
            return [p1, p2, p3][i - 1]
        return p(i, k - 1) * (p(i, k - 1) + 2 * p(i % 3 + 1, k - 1))
    
    @cache
    def q(i, k):
        if k < 0:
            raise ValueError("q called with k < 0")
        if k == 0:
            return p(i, 0) * (p(i, 0) + 2 * p((i + 1) % 3 + 1, 0))
        k -= 1
        return (
          (q(i, k) * (q(i, k) + 2*q(i%3 + 1, k))) *
            (
              p(i, k + 1) * (p(i, k + 1) + 2*p((i + 1)%3 + 1, k + 1))
              + p(i%3 + 1, k + 1) * (p(i%3 + 1, k + 1) + 2*p(i, k + 1))
            )
          +
          (q(i%3 + 1, k) * (q(i%3 + 1, k) + 2*q((i + 1)%3 + 1, k))) *
            (p(i, k + 1) * (p(i, k + 1) + 2*p((i + 1)%3 + 1, k + 1)))
        )
    
    return (p, q)
In [46]:
getcontext().prec = 1000
a = Decimal(1) / Decimal(10)
b = Decimal(8) / Decimal(10)
c = a
p1, q1 = make_pq(a, b, c)
p2, q2, z2, o = make_pqzo(a, b, c)
In [11]:
p1(1, 10) - p2(1, 10)
Out[11]:
Decimal('0E-28')
In [49]:
def plot_q_v_q(q1, q2, k):
    fig, ax = plt.subplots()
    colors = ["red", "blue", "green", "pink", "aqua", "lime"]
    color_index = 0
    for fun, name in [[q1, "q1"], [q2, "q2"]]:
        for i in [1, 2, 3]:
            data = []
            for k_val in range(k):
                data.append(fun(i, k_val))
            ax.plot(data, label=f"${name}_{i}$", color=colors[color_index])
            color_index += 1
    ax.legend()
    # ax.set_title(f"$p_1^{{(0)}}={p1}$, $p_2^{{(0)}}={p2}$, $p_3^{{(0)}}={p3}$, $k \\in [0, {k}]$" if label is None else label)
    # fig.savefig(f"./de2/{name}_p1_{p1}-p2_{p2}-p3_{p3}-k_{k}_{datetime.now().strftime('%y-%m-%d-%H-%M')}.png")
    fig.show(warn=False)
In [50]:
plot_q_v_q(q1, q2, 500)
In [60]:
def plot_pqzo(p, q, z, o, k):
    fig, ax = plt.subplots()
    colors = []
    for high, low in [(1, 0), (.9, .1), (.7, .2), (.6, .2)]:
        colors.extend([(high, low, low), (low, high, low), (low, low, high)])
    color_index = 0
    for fun, name in [[p, "p"], [z, "z"], [q, "q"], [o, "o"]]:
        for i in [1, 2, 3]:
            data = []
            for k_val in (range(k) if fun != o else range(1, k)):
                data.append(fun(i, k_val))
            ax.plot(data, label=f"${name}_{i}$", color=colors[color_index])
            color_index += 1
    ax.legend()
    # ax.set_title(f"$p_1^{{(0)}}={p1}$, $p_2^{{(0)}}={p2}$, $p_3^{{(0)}}={p3}$, $k \\in [0, {k}]$" if label is None else label)
    # fig.savefig(f"./de2/{name}_p1_{p1}-p2_{p2}-p3_{p3}-k_{k}_{datetime.now().strftime('%y-%m-%d-%H-%M')}.png")
    fig.show(warn=False)
In [61]:
a = Decimal(1) / Decimal(10)
plot_pqzo(*make_pqzo(a, 1 - 2*a, a), 10)
In [70]:
def plot_pqz(p, q, z, o, k):
    fig, ax = plt.subplots()
    colors = []
    for high, low in [(1, 0), (.9, .1), (.7, .2), (.6, .2)]:
        colors.extend([(high, low, low), (low, high, low), (low, low, high)])
    color_index = 0
    for fun, name in [[p, "p"], [z, "z"]]:
        for i in [1, 2, 3]:
            data = []
            for k_val in range(k):
                data.append(fun(i, k_val))
            ax.plot(data, label=f"${name}_{i}$", color=colors[color_index])
            color_index += 1
    for i in [1, 2, 3]:
        x_axis = []
        y_axis = []
        k_val = Decimal(0)
        for k_val in range(k):
            x_axis.append(k_val)
            x_axis.append(k_val + .5)
            y_axis.append(q(i, k_val))
            y_axis.append(o(i, k_val + 1))
        ax.plot(x_axis, y_axis, label=f"$q_{i}$", color=colors[color_index])
        color_index += 1

    ax.legend()
    # ax.set_title(f"$p_1^{{(0)}}={p1}$, $p_2^{{(0)}}={p2}$, $p_3^{{(0)}}={p3}$, $k \\in [0, {k}]$" if label is None else label)
    # fig.savefig(f"./de2/{name}_p1_{p1}-p2_{p2}-p3_{p3}-k_{k}_{datetime.now().strftime('%y-%m-%d-%H-%M')}.png")
    fig.show(warn=False)
In [71]:
a = Decimal(1) / Decimal(10)
plot_pqz(*make_pqzo(a, 1 - 2*a, a), 10)
In [105]:
def plot_pq(p, q, z, o, k):
    plot_pq_range(p, q, z, o, 0, k)

def plot_pq_range(p, q, z, o, k_low, k_high, log_scale=False):
    fig, ax = plt.subplots()
    if log_scale:
        ax.set_yscale("log", base=2)
    colors = []
    for high, low in [(1, 0), (.6, .2)]:
        colors.extend([(high, low, low), (low, high, low), (low, low, high)])
    color_index = 0
    for i in [1, 2, 3]:
        data = []
        x_vals = []
        for k_val in range(k_low, k_high):
            x_vals.append(k_val)
            data.append(p(i, k_val))
        ax.plot(x_vals, data, label=f"$p_{i}$", color=colors[color_index])
        color_index += 1
    for i in [1, 2, 3]:
        x_axis = []
        y_axis = []
        k_val = Decimal(0)
        for k_val in range(k_low, k_high):
            x_axis.append(k_val)
            x_axis.append(k_val + .5)
            y_axis.append(q(i, k_val))
            y_axis.append(o(i, k_val + 1))
        ax.plot(x_axis, y_axis, label=f"$q_{i}$", color=colors[color_index])
        color_index += 1

    ax.legend()
    # ax.set_title(f"$p_1^{{(0)}}={p1}$, $p_2^{{(0)}}={p2}$, $p_3^{{(0)}}={p3}$, $k \\in [0, {k}]$" if label is None else label)
    # fig.savefig(f"./de2/{name}_p1_{p1}-p2_{p2}-p3_{p3}-k_{k}_{datetime.now().strftime('%y-%m-%d-%H-%M')}.png")
    fig.show(warn=False)
In [74]:
a = Decimal(1) / Decimal(10)
pqzo = make_pqzo(a, 1 - 2*a, a)
for k in range(10, 101, 10):
    plot_pq(*pqzo, k)
In [76]:
getcontext().prec = 2000
for i in range(1, 32, 2):
    a = Decimal(i) / Decimal(100)
    plot_pq(*make_pqzo(a, 1 - 2*a, a), 100)
In [77]:
getcontext().prec = 2000
epsilon = Decimal(i) / Decimal(2**10)
c = epsilon ** 2
a = epsilon * 2
b = 1 - a - c
plot_pq(*make_pqzo(a, b, c), 30) # note p1 high in range 10–30
In [85]:
getcontext().prec = 10000
epsilon = Decimal(i) / Decimal(2**75)
c = epsilon ** 2
a = epsilon * 2
b = 1 - a - c
plot_pq(*make_pqzo(a, b, c), 3 * 75)
plot_pq(*make_pqzo(a, b, c), 3 * 75 + 32*75)
In [98]:
getcontext().prec = 10000
epsilon = Decimal(i) / Decimal(2**75)
c = epsilon
a = epsilon
b = 1 - a - c
pqzo = make_pqzo(a, b, c)
plot_pq(*pqzo, 3 * 75)
plot_pq(*pqzo, 3 * 75 + 32*75)
plot_pq_range(*pqzo, 50, 100)
plot_pq_range(*pqzo, 2500, 3 * 75 + 32*75)
plot_pq_range(*pqzo, 60, 75)
plot_pq_range(*pqzo, 2550, 2565)
In [107]:
getcontext().prec = 10000
epsilon = Decimal(i) / Decimal(2**80)
c = epsilon
a = epsilon
b = 1 - a - c
#pqzo = make_pqzo(a, b, c)
plot_pq(*pqzo, 3 * 80)
plot_pq(*pqzo, 3 * 80 + 32*80)
plot_pq_range(*pqzo, 65, 80)
plot_pq_range(*pqzo, 0, 3 * 80, log_scale=True)
plot_pq_range(*pqzo, 0, 3 * 80 + 32*80, log_scale=True)
plot_pq_range(*pqzo, 65, 80, log_scale=True)
In [113]:
getcontext().prec = 10000
epsilon = Decimal(i) / Decimal(2**100)
c = epsilon
a = epsilon
b = 1 - a - c
#pqzo = make_pqzo(a, b, c)
plot_pq(*pqzo, 3 * 100)
plot_pq(*pqzo, 3 * 100 + 32*100)
plot_pq_range(*pqzo, 85, 100)
plot_pq_range(*pqzo, 93, 96)