%matplotlib inline
from functools import cache
import numpy as np
import matplotlib.pyplot as plt
import itertools
import math
from decimal import Decimal, getcontext
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)
def plot_pq(p, q, z, o, k, **kwargs):
plot_pq_range(p, q, z, o, 0, k, **kwargs)
def plot_pq_range(p, q, z, o, k_low, k_high, simplex=True, basic=True, log_scale=False):
if simplex:
fig, ax = plt.subplots()
X = np.array([1, -1, 0])
Y = np.array([0, 0, 1])
for fun, name in [[p, 'p'], [q, 'q']]:
x_axis = []
y_axis = []
for k in range(k_low, k_high):
point = np.array([fun(1, k), fun(2, k), fun(3, k)])
x_axis.append(point @ X)
y_axis.append(point @ Y)
ax.plot(x_axis, y_axis, '->', label=name)
#ax.plot(x_axis, y_axis, '-4', label=name)
ax.set_xlim(-1.1, 1.1)
ax.set_ylim(-.1, 1.1)
ax.legend()
fig.show(warn=False)
if basic:
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 = []
for k in range(k_low, k_high):
x_axis.append(k)
x_axis.append(k + .5)
y_axis.append(q(i, k))
y_axis.append(o(i, k + 1))
ax.plot(x_axis, y_axis, label=f"$q_{i}$", color=colors[color_index])
color_index += 1
ax.legend()
# ax.set_title
# fig.savefig
fig.show(warn=False)
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)
<ipython-input-18-c914d1ec050e>:6: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). fig, ax = plt.subplots()
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
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)
getcontext().prec = 10000
epsilon = Decimal(1) / 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)
getcontext().prec = 10000
epsilon = Decimal(1) / Decimal(2**100)
c = epsilon
a = epsilon
b = 1 - a - c
pqzo = make_pqzo(a, b, c)
plot_pq(*pqzo, 3 * 100 + 64*100)
getcontext().prec = 20000
epsilon = Decimal(1) / Decimal(2**50)
c = epsilon
a = epsilon
b = 1 - a - c
pqzo = make_pqzo(a, b, c)
plot_pq(*pqzo, 3 * 100, basic=False)
plot_pq(*pqzo, 6 * 100, basic=False)
plot_pq(*pqzo, 9 * 100, basic=False)
plot_pq(*pqzo, 20 * 100, basic=False)
plot_pq(*pqzo, 40 * 100, basic=False)
plot_pq(*pqzo, 60 * 100, basic=False)
plot_pq(*pqzo, 80 * 100, basic=False)
plot_pq(*pqzo, 100 * 100, basic=False)
plot_pq(*pqzo, 120 * 100, basic=False)
plot_pq(*pqzo, 140 * 100, basic=False)
plot_pq(*pqzo, 160 * 100, basic=False)
plot_pq(*pqzo, 180 * 100, basic=False)
plot_pq(*pqzo, 200 * 100, basic=False)
# plot_pq(*pqzo, 400 * 100, basic=False)
plot_pq(*pqzo, 800 * 100, basic=False)
# plot_pq(*pqzo, 1000 * 100, basic=False)
# plot_pq(*pqzo, 1500 * 100, basic=False)
# plot_pq(*pqzo, 2000 * 100, basic=False)
def plot_p_simplex(p, k_low, k_high):
fig, ax = plt.subplots()
X = np.array([1, -1, 0])
Y = np.array([0, 0, 1])
x_axis = []
y_axis = []
for k in range(k_low, k_high):
point = np.array([p(1, k), p(2, k), p(3, k)])
x_axis.append(point @ X)
y_axis.append(point @ Y)
ax.plot(x_axis, y_axis, '->')
ax.set_xlim(-1.1, 1.1)
ax.set_ylim(-.1, 1.1)
ax.legend()
fig.show(warn=False)
getcontext().prec = 2000
p = make_pqzo(Decimal(".3"), Decimal(".4"), Decimal(".3"))[0]
plot_p_simplex(p, 0, 100)
plot_p_simplex(p, 0, 200)
No handles with labels found to put in legend. No handles with labels found to put in legend.
getcontext().prec = 2000
p = make_pqzo(Decimal(".4"), Decimal(".3"), Decimal(".3"))[0]
plot_p_simplex(p, 0, 100)
plot_p_simplex(p, 0, 200)
No handles with labels found to put in legend. No handles with labels found to put in legend.
getcontext().prec = 2000
p = make_pqzo(Decimal(".6"), Decimal(".2"), Decimal(".2"))[0]
plot_p_simplex(p, 0, 100)
plot_p_simplex(p, 0, 200)
No handles with labels found to put in legend. No handles with labels found to put in legend.
def plot_trajectory(abc, k):
a, b, c = (Decimal(x) for x in abc.split())
pqzo = make_pqzo(a, b, c)
plot_pq(*pqzo, k, basic=False)
plot_trajectory(".3 .4 .3", 100)
plot_trajectory(".1 .8 .1", 100)
import random
random.seed(0)
for i in range(5):
a = random.randrange(1000)
b = random.randrange(1000 - a)
print(f".{a:03} .{b:03} .{(1000-a-b):03}")
plot_trajectory(f".{a:03} .{b:03} .{(1000-a-b):03}", 100)
.864 .098 .038 .776 .107 .117 .041 .265 .694 .988 .008 .004 .497 .207 .296
