import math
# Attempted extrapolation of Gidney 2025

# Assumes the logical error of a surface code scales as C*(physical_error_rate/T)^((d+1)/2)
SURFACE_CODE_C = 0.14 
SURFACE_CODE_T = 0.01225

def get_distance(target_logical_erorr, physical_error):
	return math.ceil(2*math.log(target_logical_erorr/SURFACE_CODE_C)/math.log(phys_error/SURFACE_CODE_T)-1) # rearrange the usual formula


def estimate(n,l,s,f,  phys_error):
	# m, lenm as defined in Table 2
	m = math.ceil(n/2 + n/s)
	lenm = math.ceil(math.log(m,2))
	# Maximum of total qubits in Table 4
	logical_qubits = m + f + l +  lenm + max(2*f+l,3*l,2*l+2*lenm)
	# This is based on loop4 specifically
	hot_logical = logical_qubits - m
	cold_logical = m
	# Straight from section 3.2
	compute = 7*18
	target_hot_logical_error = 1e-15 # just keeping this constant
	d_hot = get_distance(target_hot_logical_error, phys_error)
	hot_phys_to_logical_ratio = 2*(d_hot+1)**2
	target_cold_logical_error = 1e-7 # based on a d=11 underlying code
	d_cold = get_distance(target_cold_logical_error,phys_error)
	cold_phys_to_logical_ratio = 2*(d_cold+1)**2 *1.5 # 1.5 is based on the 430 qubits from simulations to the 288 per surface code patch at that distance
	cold_phys_to_logical_ratio = min(cold_phys_to_logical_ratio,hot_phys_to_logical_ratio) # don't use yoked codes if they underperform
	total_physical_count = (hot_logical+compute)*hot_phys_to_logical_ratio + cold_logical*cold_phys_to_logical_ratio
	return total_physical_count, d_hot

# use the same range of s,l,f from Section 3
estimates = {}
for n in [1024,2048,4096,8192,16384]:
	for phys_error in [SURFACE_CODE_T * (m/10) for m in range(9,0,-1)] + [1e-3, 1e-4,1e-5,1e-6,1e-7,1e-8,1e-9]:
		min_est = (float('inf'), 0)
		for s in range(2,14):	
			for l in range(18,25):
				for f in range(24,59):
					if estimate(n,l,s,f, phys_error)[0] < min_est[0]:
						min_est = estimate(n,l,s,f, phys_error)
		estimates[(n,phys_error)] = min_est
		print(n,phys_error,min_est)
print(estimates)