generate.sage

back to table · edit · history · where entries came from · files · download

2030 bytes, as of the version from 2026-08-09 09:11 (current). Recorded here, not run.

#References:
#[1] B. Connon "A recurrence relation for the Li/Keiper constants in terms of the Stieltjes constants"

import yaml
import os
from utils.utils import numbers_to_yaml
from utils.utils import real_interval_to_sage_string
from utils.utils import blur_real_interval

path = 'data/Series_expansions/Keiper-Li_coefficients/'

prec10 = 100 #relative precision in base 10
n_max = 1000
n_range = [0..n_max]

print("n_max:",n_max)

RIFprec = RealIntervalField(prec10 * 3.4 * 5)

#xi(s) = 1/2 * s*(s-1)*pi^(-s/2)*gamma(s/2)*zeta(s)

psi_n_3o2 = {}
for n in n_range:
	psi_n_3o2[n] = RIFprec(psi(n,3/2))

print("finished computing psi's")

binom = {}
for n in n_range:
	binom[n,0] = RIFprec(1)
	binom[n,n] = RIFprec(1)
for n in [1..n_max]:
	for k in [1..n-1]:
		binom[n,k] = binom[n-1,k-1] + binom[n-1,k]

print("finished computing pascal's triangle")

fact = {0: 1}
for n in n_range:
	if n>=1:
		fact[n] = n*fact[n-1]

print("finished computing factorials")

stiel = {}
for n in n_range:
	print("n:",n)
	stiel[n] = RIFprec(blur_real_interval(stieltjes(n).n(RIFprec.prec()+10)))
	#Blurring is done as stieltjes(n).n() might return an integer,
	#which is then interpreted by RIFprec as exact, 
	#which is wrong (though probably not an issue, but let's just be safe).

print("finished computing Stieltjes constants")

lamda = {}
eta = {}

for n in n_range:
	#(3.35) of [1]:
	lamda[n] = RIFprec(1/2*n*(euler_gamma-2*log(2)+2-log(pi)) + \
		sum(
			binom[n,j] / fact[j-1] *
			(1/2^j * psi_n_3o2[j-1] - fact[j-1]*eta[j-1])
			for j in range(2,n+1)
		)
	)

	#(4.5) of [1], originally due to Coffey:
	eta[n] = (-1)^(n+1)*(n+1)/fact[n] * stiel[n] + \
		(-1)^(n+1) * sum( 
			(-1)^(k+1)/fact[n-k-1]*stiel[n-k-1] * eta[k]
			for k in range(n)
		)
		
numbers = {}

for n in n_range:
	n_str = str(n)
	number = lamda[n]	
	
	numbers[n_str] = real_interval_to_sage_string(
		number,
		max_digits = prec10,
	).replace('?','')
		 

filename = os.path.join(path, 'numbers.yaml')
yaml.dump(numbers, stream = open(filename, 'w'), sort_keys = False)