""" A001414.py.txt Copyright (C) 2020 Steve Witham ess_doubleyou_at_tiac_notthis_dot_net BSD license https://opensource.org/licenses/BSD-2-Clause """ from math import log, exp from collections import Counter import matplotlib.pyplot as plt from matplotlib import colors from factor import gen_factors_powers def A001414(n): return sum(base * exp for (base, exp) in gen_factors_powers(n)) # print([A001414(n) for n in range(1, 11)]) N = 800 # 10**5 NA = [(n, A001414(n)) for n in range(2, N + 1)] NLNALA = [ (n, log(n), a, log(a)) for (n, a) in NA] HISTO = list(Counter(a for (n, a) in NA).items()) LHISTO = sorted((log(ha), hc) for (ha, hc) in HISTO) # xs, ys = zip(*NA) # plt.plot(xs, ys, '.b', label="linear", alpha=.5) ns, lns, As, lAs = zip(*NLNALA) has, hcs = zip(*LHISTO) # plt.plot(xs, ys, '.b', label="log-log") # Thus we can cut the plotting window in several hexbins # nbins = 20 # axes[1].set_title('Hexbin') # axes[1].hexbin(x, y, gridsize=nbins, cmap=plt.cm.BuGn_r) def plots1(): if False: plt.figure(figsize=(12, 12), dpi=75) plt.gca().set_title('A001414 linear-log') plt.plot(ns, lAs, '.k', markersize=.25, label="log-log") # colors.LogNorm() plt.show() if False: # 2D Histogram plt.figure(figsize=(12, 12), dpi=75) plt.gca().set_title('A001414 linear-log') plt.hist2d(ns, lAs, bins=700, cmap="Greys", norm=colors.PowerNorm(.05), label="log-log") # colors.LogNorm() plt.show() if False: plt.figure(figsize=(12, 12), dpi=75) plt.gca().set_title('A001414 log-log' % n_plot) plt.plot(lns, lAs, '.k', markersize=.25, label="log-log") # colors.LogNorm() plt.show() if False: plt.figure(figsize=(12, 12), dpi=75) plt.gca().set_title('A001414 log-log') plt.hist2d(lns, lAs, bins=700, cmap="Greys", norm=colors.PowerNorm(.05)) # colors.LogNorm() plt.show() for n_plot, lim, msz in (100, 100, 4), (200, 150, 3), (400, 225, 2), (800, 250, 1.5): plt.figure(figsize=(12, 12), dpi=75) plt.gca().set_title('A001414 linear-linear N = %d' % n_plot) pns = [n for (n, A) in zip(ns[:n_plot], As[:n_plot]) if A <= lim] pAs = [A for (n, A) in zip(ns[:n_plot], As[:n_plot]) if A <= lim] plt.plot(pns, pAs, '.k', markersize=msz, label="lin-lin") plt.show() def table(): for n in range(2, 28): print(n, A001414(n), 3*log(n, 3), 2 * log(n, 2)) print(exp(4.66), exp(4.75)) def plot_hist(): plt.figure(figsize=(6, 6), dpi=75) plt.gca().set_title("A001414 histogram N=%5.0e" % N) plt.plot(hcs, has, color=(.75, .75, .75)) plt.plot(hcs, has, '.k', markersize=2) plt.show() plots1()