Created
June 8, 2021 18:52
Revisions
-
hdf created this gist
Jun 8, 2021 .There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,67 @@ import sys import numpy as np import matplotlib.pyplot as plt ## Recursive divisibility histogram # How many times can we divide a number by an integer, and by how many integers? # Divisors start from 2 and increase by 1 until we reach half of the number. # Show 3D histogram. n = int(sys.argv[1]) if len(sys.argv) > 1 and sys.argv[1].isdigit() else 180 def div_hist(n): ds = {} for d in range(2, n//2): ds[d] = 0 n2 = n while(n2%d==0): n2 /= d ds[d] += 1 if ds[d] < 2: del ds[d] return ds #print(div_hist(n)) def div_hist3d(m): ds = {} for n in range(2, m+1): ds[n] = {} #print('working on', n) for d in range(2, n//2): ds[n][d] = 0 n2 = n #print(' trying to divide', n2, 'by', d) while(n2%d==0): #print(' ', n2, 'is divisible by', d) n2 //= d ds[n][d] += 1 ds[n] = [e[1] for e in ds[n].items()] ds[n] += [0]*(m-1-len(ds[n])) # zero fill return [e[1] for e in ds.items()] ds = np.array(div_hist3d(n))[::-1] #print(ds, ds.shape) fig = plt.figure() ax = fig.add_subplot(projection='3d') _x = _y = np.arange(len(ds)) _xx, _yy = np.meshgrid(_x, _y, indexing="ij") x, y = _xx[::-1].ravel(), _yy.ravel() z = ds.ravel() if np.count_nonzero(z) < 1: print('No divisors found.') exit() x, y, z = map(np.array, zip(*[(x[i], y[i], z[i]) for i in range(0, len(z)) if z[i] > 0])) ax.set_xlim(2, len(ds)+2) ax.set_ylim(2, (len(ds)+1)//2) #ax.set_zlim(0, n//2) #print(x, y, z) ax.bar3d(x+2, y+2, 0, 1-.1, 1/2-.2, z) ax.set_xlabel('Number') ax.set_ylabel('Divisor') ax.set_zlabel('Recursive divisibility') ax.set_title('Recursive divisibility of numbers up to '+str(n)) plt.show()