import numpy as np
import matplotlib.pyplot as plt
import math

def sqrt3(center, x):
	if x - center > 0:
		return (x - center) ** (1/3)
	else:
		return -((center - x) ** (1/3))

def sqrt3sqr(center, x):
	if x - center > 0:
		return (x - center) ** (-2/3)
	elif x - center < 0:
		return (center - x) ** (-2/3)
	else:
		return math.inf

def pompeiuInverse(n, x):
	q = 1.8
	f = 1
	s = 0
	for i in range(2, n+1):
		subs = 0
		for j in range(1, i):
			if math.gcd(i, j) == 1:
				subs += sqrt3(j/i, x)
		s += subs / f
		f *= q
	return s

def pompeiuInverseDerivative(n, x):
	q = 1.8
	f = 1
	s = 0
	for i in range(2, n+1):
		subs = 0
		for j in range(1, i):
			if math.gcd(i, j) == 1:
				subs += sqrt3sqr(j/i, x)
		s += subs / f
		f *= q
	return s/3

def pompeiu(n, x, eps):
	a = -5
	b = 5
	y1 = pompeiuInverse(n, a)
	y3 = pompeiuInverse(n, b)
	c = 0
	
	while (b - a)*2 >= eps:
		c = (a + b) / 2
		y2 = pompeiuInverse(n, c)
		if x < y2:
			y3 = y2
			b = c
		elif x > y2:
			y1 = y2
			a = c
		else:
			break
	return c
	
def pompeiuDerivative(n, x, eps):
	return 1/pompeiuInverseDerivative(n, pompeiu(n, x, eps))

n = 50
zerosx = []
for i in range(2, n+1):
	for j in range(1, i):
		if math.gcd(i, j) == 1:
			zerosx.append(pompeiuInverse(n, j/i))
zerosx.sort()
x = []
y = []
for i in range(0, len(zerosx) - 1):
	x.append(zerosx[i])
	y.append(0)
	d = (zerosx[i+1]-zerosx[i]) / 20
	for j in range(1, 20):
		c = zerosx[i]+j * d
		x.append(c)
		y.append(pompeiuDerivative(n, c, 0.0001))
x.append(zerosx[len(zerosx)-1])
y.append(0)
zeroy = [0 for i in x]


f, ax = plt.subplots()

ax.fill_between(x, y, zeroy)

ax.set_xlim(-3.1, 3.1)
ax.set_ylim(-0.02, 0.4)

plt.savefig("Pompeiu derivative.svg")