Recursive divisibility histogram

div_hist.py

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()