Simple logarithmic integral

logarithmic_integral.py

from decimal import Decimal, getcontext
import math
from itertools import cycle

ctx = getcontext()
ctx.prec = 256
LI_PRECISION = 300
# Constant from: http://www.plouffe.fr/simon/constants/gamma.txt
euler_mascheroni = Decimal(
	'.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495146314472498070824809605040144865428362241739976449235362535003337429373377376739427925952582470949160087352039481656708532331517766115286211995015079847937450857057400299213547861466940296043254215190587755352')
# Cache for faster
fact_mul_k_cache = tuple(math.factorial(k) * k for k in range(1, LI_PRECISION + 1))


# A dead simple algorithm for approximating the logarithmic integral
# Accuracy is near equal to sympy's li()
def li(input, li_precision=LI_PRECISION):
	assert input > 1 and li_precision <= LI_PRECISION
	global fact_mul_k_cache, euler_mascheroni
	ln_x = ctx.ln(input)
	approx_sum: Decimal = sum(
		ctx.divide(ctx.power(ln_x, k), fact_mul_k_cache[k - 1]) for k in range(1, li_precision + 1))
	return ctx.ln(ln_x) + euler_mascheroni + approx_sum


# A more complex version of li, using a different approximation method
# Claims to have a 'more rapid' convergence
# Testing at 10000000 shows 20% fewer iterations are needed to obtain similar accuracy
def li_2(input, li_precision=int(LI_PRECISION//1.25)):
	assert input > 1 and li_precision <= LI_PRECISION
	global fact_mul_k_cache, euler_mascheroni
	pos_neg_cycle = cycle((1, -1))
	ln_x = ctx.ln(input)
	approx_sum = sum(
		ctx.multiply(
			ctx.divide(
				ctx.multiply(next(pos_neg_cycle), ctx.power(ln_x, n)),
				ctx.multiply(math.factorial(n), 1 << (n - 1))
			),
			# summation k=0 -> floor((n-1)/2) requires at least 1 iteration at [0]
			# therefore end of range() must be at least 1
			sum(ctx.divide(1, 2 * k + 1) for k in range(0, (n - 1) // 2 + 1))
		)
		for n in range(1, li_precision + 1)
	)
	return ctx.ln(ln_x) + euler_mascheroni + ctx.multiply(ctx.sqrt(input), approx_sum)

if __name__ == "__main__":
	import sympy
	print(w := sympy.li(10000000).evalf(258))
	print(x := li(10000000))
	print(y := li_2(10000000))
	print(z1 := min(i for i, (x, y) in enumerate(zip(str(w), str(x))) if x != y))
	print(z2 := min(i for i, (x, y) in enumerate(zip(str(x), str(y))) if x != y))
	print(str(x) == str(y))