multiloop_sunset_mellin.sage

# Define the Mellin transform
def Mellin(L, z_vars=None):
    """
    Define the Mellin transform for a given list of masses.
    
    Mellin = GAMMA(1+z1+...+z_{L+1})^2 * sin(Pi*(z1+...+z_{L+1})) / Pi 
             * (m1^2/t)^z1 * ... * (m_{L+1}^2/t)^z_{L+1} 
             / (GAMMA(1+z1)^2 * ... * GAMMA(1+z_{L+1})^2)
    
    Parameters:
    -----------
    L : list
        List of masses [m1, m2, ..., m_{L+1}]
    z_vars : list, optional
        List of variable names for z1, z2, ..., z_{L+1}
        If None, will create variables z1, z2, ..., z_n where n = len(L)
    
    Returns:
    --------
    Symbolic expression for the Mellin transform
    """
    from sage.symbolic.ring import SR
    from sage.functions.gamma import gamma
    from sage.functions.trig import sin
    from sage.symbolic.constants import pi
    
    n = len(L)  # Number of masses
    
    # Create symbolic variables for z1, z2, ..., z_n
    if z_vars is None:
        z_vars = [SR.var(f'z{i+1}') for i in range(n)]
    elif len(z_vars) != n:
        raise ValueError(f"Length of z_vars must match length of L (got {len(z_vars)}, expected {n})")
    
    # Create symbolic variable t
    t = SR.var('t')
    
    # Sum of all z variables
    z_sum = sum(z_vars)
    
    # Numerator: GAMMA(1 + z1 + ... + z_n)^2 * sin(Pi*(z1 + ... + z_n)) / Pi
    numerator = gamma(1 + z_sum)^2 * sin(pi * z_sum) / pi
    
    # Product terms: (m_i^2/t)^z_i for each mass
    for i in range(n):
        numerator *= (L[i]^2 / t)^z_vars[i]
    
    # Denominator: GAMMA(1 + z1)^2 * ... * GAMMA(1 + z_n)^2
    denominator = 1
    for i in range(n):
        denominator *= gamma(1 + z_vars[i])^2
    
    # Complete Mellin transform
    M = numerator / denominator
    
    return M


# Example usage:
# For masses [m1, m2, m3]
#L = [SR.var('m1'), SR.var('m3'), SR.var('m3')]
#M = Mellin(L)
#print(M)

# Or with specific numeric masses
#L_numeric = [1, 3, 5]
#M_numeric = Mellin(L_numeric)
#print(M_numeric)

# Or with custom z variable names
#z_custom = [SR.var('z'), SR.var('w'), SR.var('u')]
#M_custom = Mellin([1, 2, 3], z_vars=z_custom)
#print(M_custom)

# Compute the derivatives of the Mellin transform
def MellinFirstDerivativesEvaluate(L, n_values, z_vars=None, simplify=True):
    """
    Compute the mixed partial derivative d^n/(dz1...dz_n) Mellin(L)
    and evaluate at z1=n1, z2=n2, ..., z_n=n_n
    
    This takes the first derivative with respect to each z variable once,
    then evaluates at the specified point.
    
    Parameters:
    -----------
    L : list
        List of masses [m1, m2, ..., m_n]
    n_values : list
        List of values [n1, n2, ..., n_n] at which to evaluate
    z_vars : list, optional
        List of variable names for z1, z2, ..., z_n
    simplify : bool, optional
        Whether to simplify the result (default: True)
    
    Returns:
    --------
    d^n/(dz1...dz_n) Mellin(L) evaluated at (n1, n2, ..., n_n)
    """
    from sage.symbolic.ring import SR
    
    n = len(L)
    
    if len(n_values) != n:
        raise ValueError(f"Length of n_values must match length of L (got {len(n_values)}, expected {n})")
    
    # Create symbolic variables if not provided
    if z_vars is None:
        z_vars = [SR.var(f'z{i+1}') for i in range(n)]
    
    # Get the Mellin transform
    M = Mellin(L, z_vars=z_vars)
    
    # Take first derivative with respect to each z variable
    result = M
    for z in z_vars:
        result = diff(result, z)
    
    # Create substitution dictionary
    substitutions = {z_vars[i]: n_values[i] for i in range(n)}
    
    # Evaluate at the given point
    result = ((-1)^(sum(z_vars))*result).subs(substitutions).subs(euler_gamma==0)
    
    if simplify:
        result = result.simplify_full()
    
    return result

# Compute the truncated first derivative with respect to each variables
def MellinFirstDerivativesSumTruncated(L, Nmax, z_vars=None, simplify=True, verbose=True):
    """
    Compute the multiple sum of MellinFirstDerivativesEvaluate over all n=[n1,...,n_{L+1}]
    such that n1 + n2 + ... + n_{L+1} <= Nmax, with each n_i >= 0.
    
    Parameters:
    -----------
    L : list
        List of masses [m1, m2, ..., m_{L+1}]
    Nmax : int
        Maximum order for truncation (sum of all n_i)
    z_vars : list, optional
        List of variable names for z1, z2, ..., z_{L+1}
    simplify : bool, optional
        Whether to simplify the result (default: True)
    verbose : bool, optional
        Whether to print progress information (default: True)
    
    Returns:
    --------
    Sum of MellinFirstDerivativesEvaluate(L, n) over all valid n
    """
    from sage.symbolic.ring import SR
    
    L_len = len(L)

    if verbose:
        import logging
        logging.basicConfig()
        logging.getLogger().setLevel(logging.INFO)
        logger = logging.getLogger("Mellin_derivative")
        result_total_count = count_solutions(L_len-1, Nmax)
        logger.info(f"number of terms to compute: {result_total_count}")
        count_step= 10^floor(log(result_total_count, 10) - 1)
    
    # Create symbolic variables if not provided
    if z_vars is None:
        z_vars = [SR.var(f'z{i+1}') for i in range(L_len)]
    
    # Helper function to generate compositions
    def compositions_with_length(n, k):
        """
        Generate all compositions of n into exactly k non-negative parts.
        """
        if k == 1:
            yield [n]
            return
        
        for i in range(n + 1):
            for comp in compositions_with_length(n - i, k - 1):
                yield [i] + comp
    
    # Initialize the sum
    total_sum = 0
    count = 0
    
    # Iterate through all possible sums from 0 to Nmax
    for total in range(Nmax + 1):
        # For each total, generate all compositions of 'total' into L_len non-negative parts
        compositions = compositions_with_length(total, L_len)
        
        for n_tuple in compositions:
            try:
                contribution = -1/t*MellinFirstDerivativesEvaluate(L, list(n_tuple), z_vars=z_vars, simplify=False)
                total_sum += contribution
                count += 1
                
                if verbose and count % count_step == 0:
                    logger.info(f"Processed {count} terms...")
                    
            except Exception as e:
                if verbose:
                    logger.info(f"Warning: Error computing for n={n_tuple}: {e}")
                continue
    
    if verbose:
        logger.info(f"Total: computed {count} terms")
    
    if simplify:
        logger.info("Simplifying result...")
        total_sum = total_sum.simplify_full()
    
    return (-1)^(L_len-1)*total_sum


# Example usage:
#if __name__ == '__main__':
    #from sage.symbolic.ring import SR
    
    ## Example with symbolic masses
    #L_sym = [SR.var('m1'), SR.var('m2'), SR.var('m3')]
    #n_vals = [0, 1, 2]
    
    #result = MellinFirstDerivativesEvaluate(L_sym, n_vals)
    #print("d³/(dz1 dz2 dz3) Mellin([m1,m2,m3]) at z1=0, z2=1, z3=2:")
    #print(result)
    #print()
    
    ## Example with numeric masses at origin
    #L_num = [1, 3, 5]
    #n_vals_origin = [0, 0, 0]
    
    #result_origin = MellinFirstDerivativesEvaluate(L_num, n_vals_origin)
    #print("d³/(dz1 dz2 dz3) Mellin([1,3,5]) at z1=0, z2=0, z3=0:")
    #print(result_origin)
    #print()
    
    ## Example at a different point
    #n_vals2 = [1, 2, 3]
    #result2 = MellinFirstDerivativesEvaluate(L_num, n_vals2)
    #print("d³/(dz1 dz2 dz3) Mellin([1,3,5]) at z1=1, z2=2, z3=3:")
    #print(result2)

def c_k(k, m):
    """
    Compute c_k(n) using Sage's built-in bell_polynomial.
    Uses psi(k, x) for polygamma function (psi^{(k)}(x)).
    Assumes bell_polynomial uses variables x_1, x_2, ..., x_n (1-indexed).
    This is equation (18) of the paper
    """
    k = Integer(k)
    m = Integer(m)
    
    if k == 0:
        return SR(0)
    
    result = SR(0)
    
    for r1 in range(k+1):
        r1 = Integer(r1)
        
        # KroneckerDelta[Mod[k-r1,2]-1] means k-r1 must be ODD
        if (k - r1) % 2 != 1:
            continue
        
        # (-Pi^2)^((k-r1-1)/2)
        exponent = (k - r1 - 1) // 2
        pi_factor = (SR(-1) * SR(pi)**2)**exponent
        
        for a in range(r1+1):
            a = Integer(a)
            
            # Binomial coefficients: C(k,r1) * C(r1,a)
            binom_k_r1 = binomial(k, r1)
            binom_r1_a = binomial(r1, a)
            
            # Compute Y_a using bell_polynomial
            if a == 0:
                Y_a = SR(1)
            else:
                Y_a_poly = bell_polynomial(a)
                # Get variables and substitute
                vars_in_poly = Y_a_poly.variables()
                vars_sorted = sorted(vars_in_poly, key=str)
                
                # Substitute x_j -> psi(j-1, m+1) for j = 1, ..., a
                # vars_sorted[0] is x_1, which should become psi(0, m+1)
                subs_dict = {var: psi(idx, m+1).subs(euler_gamma==0) for idx, var in enumerate(vars_sorted)}
                Y_a = Y_a_poly.subs(subs_dict)
            
            # Compute Y_{r1-a}
            if r1 - a == 0:
                Y_r1_a = SR(1)
            else:
                Y_r1_a_poly = bell_polynomial(r1-a)
                vars_in_poly = Y_r1_a_poly.variables()
                vars_sorted = sorted(vars_in_poly, key=str)
                
                # Substitute x_j -> psi(j-1, m+1) for j = 1, ..., r1-a
                subs_dict = {var: psi(idx, m+1).subs(euler_gamma==0) for idx, var in enumerate(vars_sorted)}
                Y_r1_a = Y_r1_a_poly.subs(subs_dict)
            
            term = binom_k_r1 * binom_r1_a * pi_factor * Y_a * Y_r1_a
            result = result + term
    
    # (-1)^(m+1) factor
    result = result * ((-1)**(m+1))
    return SR(result).simplify_full()

def int_sunset_L(L, p_sq, m_list, max_r_sum=5,simplify=True, verbose=True):
    """
    Compute Intsunset^{(L)}(p^2, underline{m}^2).
    
    Formula:
    -1/p^2 * sum_{r_1,...,r_{L+1} >= 0} [(sum r_i)!^2 / prod(r_i!^2)] * prod((m_i^2/p^2)^{r_i})
    * sum_{r=1}^{L+1} c_r(sum r_i) * P_{L+1}^{L+1-r}(ell_1(r_1), ..., ell_{L+1}(r_{L+1}))
    
    where ell_i(r) = log(-m_i^2/p^2) - 2* sum_{k=1}^{r-1} 1/k
    
    Args:
        L: The loop order (number of propagators minus 1)
        p_sq: p^2 (momentum squared)
        m_list: list of m_i values (length should be L+1)
        max_r_sum: maximum sum of r_i for truncation (since infinite sum)
    
    Returns:
        The truncated sum for Intsuset^{(L)}
    """
    L = Integer(L)
    p_sq = SR(p_sq)
    n_props = L + 1  # Number of propagators
    
    if len(m_list) != n_props:
        raise ValueError(f"m_sq_list must have length {n_props} (L+1), got {len(m_list)}")
    
    m_list = [SR(m) for m in m_list]
    m_sq_list= [m**2 for m in m_list]
    
    result = SR(0)
    count=0

    if verbose:
        import logging
        logging.basicConfig()
        logging.getLogger().setLevel(logging.INFO)
        logger = logging.getLogger("Mellin_representation")
        result_total_count = count_solutions(L, max_r_sum)*(L+1)
        logger.info(f"number of terms to compute: {result_total_count}")
        count_step= 10^floor(log(result_total_count, 10) - 1)
    
    # Sum over all (r_1, ..., r_{L+1}) with 0 <= sum(r_i) <= max_r_sum
    for r_tuple in product(range(max_r_sum+1), repeat=n_props):
        r_sum = sum(r_tuple)
        if r_sum > max_r_sum:
            continue
        
        # Multinomial coefficient: (r_1 + ... + r_{L+1})!^2 / (r_1!^2 * ... * r_{L+1}!^2)
        multinomial_coef = factorial(r_sum)**2
        for r in r_tuple:
            multinomial_coef = multinomial_coef // (factorial(r)**2)
        
        # Product of (m_i^2 / p^2)^{r_i}
        mass_ratio_product = SR(1)
        for i, r in enumerate(r_tuple):
            mass_ratio_product = mass_ratio_product * (-m_sq_list[i] / p_sq)**r
        
        # Compute ell values: ell_i(r_i) for each i
        ell_values = [ell_i(i, r, m_sq_list[i], p_sq) for i, r in enumerate(r_tuple)]
        
        # Inner sum over r from 1 to L+1
        inner_sum = SR(0)
        for r in range(1, n_props+1):  # r = 1, ..., L+1
            # c_r(sum of all r_i)
            c_val = c_k(r, r_sum)
            count +=1
            
            
            # P_{L+1}^{L+1-r}: elementary symmetric polynomial of degree L+1-r
            degree = n_props - r  # L+1-r
            sym_poly = elementary_symmetric_polynomial(ell_values, degree)
                      
            inner_sum = inner_sum + c_val * sym_poly

        if verbose and count % count_step == 0:
            logger.info(f"Processed {count} terms...")
        
        term = multinomial_coef * mass_ratio_product * inner_sum
        result = result + term

    if verbose:
        logger.info(f"Total: computed {count} terms")
    
    # Overall factor: -1/p^2
    result = ((-1)^(L)*result / p_sq).subs(euler_gamma==0) 
    
    if simplify:
        if verbose:
            logger.info("Simplifying result...")
        result = result.simplify_full()
    
    return result