number_theory_routines.sage

# Various routines

def extract_instanton(order, a_values, max_n):
    """
    Extract c(n) from a(n) = sum_{d|n} c(n/d)/d^order
    
    a_values: list or dict of a(n) values (1-indexed)
    max_n: maximum n to compute
    """
    # Initialize c array
    c = [0] * (max_n + 1)
    
    for n in range(1, max_n + 1):
        # Start with a(n)
        c[n] = a_values[n]
        
        # Subtract contributions from all proper divisors
        # a(n) = c(n)/1 + sum_{d|n, d>1} c(n/d)/d^order
        # So: c(n) = a(n) - sum_{d|n, d>1} c(n/d)/d^order
        
        divs = divisors(n)
        for d in divs:
            if d > 1:  # Skip d=1
                c[n] -= c[n//d] / (d^order)
    
    return c

def count_solutions(L, n, return_points=False):
    """
    Counts integer solutions to: r_0 + r_1 + ... + r_L <= n
    where r_i >= 0.
    """
    # Number of variables is L + 1
    num_vars = L + 1
    
    # Define inequalities in Sage format: [constant, coeff_r0, coeff_r1, ...]
    # 1. The sum constraint: n - sum(r_i) >= 0
    sum_constraint = [n] + [-1] * num_vars
    
    # 2. Non-negativity: r_i >= 0
    non_negative_constraints = []
    for i in range(num_vars):
        row = [0] * (num_vars + 1)
        row[i + 1] = 1
        non_negative_constraints.append(row)
    
    # Combine constraints
    ieqs = [sum_constraint] + non_negative_constraints
    
    # Create the Polyhedron object
    P = Polyhedron(ieqs=ieqs, base_ring=ZZ)
    
    if return_points:
        return P.integral_points()
    else:
        return P.integral_points_count()

# Example usage:
# result = count_solution(L=2, n=5)
# print(f"number of terms to compute: {result}")


def harmonic_number(n):
    """Compute H_n = sum_{k=1}^n 1/k for n >= 1, return 0 for n=0."""
    n = Integer(n)
    if n <= 0:
        return SR(0)
    return sum(QQ(1)/k for k in range(1, n+1))