# dis a resource, basically for derivative when f'(L)=1. this is Version 56, lets hope i dont need to make more versions.
import math
from fractions import Fraction
def bellpepper(n, k):
# Generates partial Bell polynomial(i forgot how this works again, if it works dont touch it)
if n == 0 and k == 0: return [(1, ())]
if n == 0 or k == 0: return []
res = []
for m in range(1, n - k + 2):
for coeff, indices in bellpepper(n - m, k - 1):
res.append((coeff * math.comb(n - 1, m - 1), (m,) + indices))
return res
def evaluate_sequence(n_max):
# Simulates the true numerical steps of the derivative histat a fixed point where f(L) L and f'(L) = 1, if not then it simulates something atleast
hist = []
# KEEP IT CLEAN. KEEP IT CLEAN. KEEP IT CLEAN. REMEMBER YOU CANT READ THE CODE WHEN DEBUGGING
h0 = {1: {(): Fraction(1)}}
for i in range(2, n_max + 1): h0[i] = {}
hist.append(h0)
for x in range(1, n_max + 2):
prev = hist[-1]
curr = {1: {(): Fraction(1)}}
for n in range(2, n_max + 1):
curr_n = {}
for a in range(1, n + 1):
bell_terms = bellpepper(n, a)
for b_coeff, indices in bell_terms:
term_profile = {(): Fraction(b_coeff)}
for idx in indices:
new_profile = {}
for p1, c1 in term_profile.items():
for p2, c2 in prev[idx].items():
m = max(len(p1), len(p2))
e1 = p1 + (0,) * (m - len(p1))
e2 = p2 + (0,) * (m - len(p2))
f_res = tuple(a_val + b_val for a_val, b_val in zip(e1, e2))
new_profile[f_res] = new_profile.get(f_res, 0) + c1 * c2
term_profile = new_profile
for p, c in term_profile.items():
if a >= 2:
m = max(len(p), a - 1)
e = p + (0,) * (m - len(p))
e_list = list(e)
e_list[a - 2] += 1
final_p = tuple(e_list)
else:
final_p = p
curr_n[final_p] = curr_n.get(final_p, 0) + c
curr[n] = {k: v for k, v in curr_n.items() if v != 0}
hist.append(curr)
return hist
def get_nth_derivative_expansion(n_max):
#Performs Langrange interpolation, or atleast tries to
hist = evaluate_sequence(n_max)
formulas = {}
x_points = list(range(len(hist)))
for n in range(1, n_max + 1):
all_profiles = set()
for h in hist:
for p in h[n].keys(): all_profiles.add(p)
formulas[n] = {}
for p in all_profiles:
y_points = [h[n].get(p, Fraction(0)) for h in hist]
poly = {0: Fraction(0)}
for i, y in enumerate(y_points):
if y == 0: continue
term = {0: Fraction(1)}
denom = 1
for j in range(len(x_points)):
if j == i: continue
denom *= (i - j)
new_term = {}
for pow, coeff in term.items():
new_term[pow + 1] = new_term.get(pow + 1, 0) + coeff
new_term[pow] = new_term.get(pow, 0) - coeff * j
term = new_term
scale = Fraction(y, denom)
for pow, coeff in term.items():
poly[pow] = poly.get(pow, 0) + coeff * scale
formulas[n][p] = {k: v for k, v in poly.items() if v != 0}
return formulas
def print_nth_expansion(target_n):
# does something i forgot, im writing these comments 3 months later
if target_n < 1: return
all_formulas = get_nth_derivative_expansion(target_n)
terms = []
for p, x_poly in all_formulas[target_n].items():
f_parts = []
for idx, power in enumerate(p):
if power == 0: continue
deriv_order = idx + 2
base_str = "f''" if deriv_order == 2 else "f'''" if deriv_order == 3 else f"f^({deriv_order})"
f_parts.append(f"({base_str}(L))^{power}" if power > 1 else f"{base_str}(L)")
f_label = " * ".join(f_parts) if f_parts else "1"
x_parts = []
for pow, coeff in sorted(x_poly.items()):
coeff_str = str(coeff) if coeff.denominator == 1 else f"({coeff.numerator}/{coeff.denominator})"
if pow == 0: x_parts.append(f"{coeff_str}")
elif pow == 1: x_parts.append(f"{coeff_str}*x")
else: x_parts.append(f"{coeff_str}*x^{pow}")
x_str = " + ".join(x_parts)
terms.append(f"[{f_label}] * ({x_str})")
print(f"c(x, {target_n}) = " + " + ".join(terms))
if __name__ == "__main__":
# change n to get derivative. yay
target_derivative_order = 5
print_nth_expansion(target_derivative_order)