AllLoopSunset/system_solver.sage at main · pierrevanhove/AllLoopSunset · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
# This files provides the routines for extracting the coefficients of the sunset graphs

# Helper function to extract and series expand log coefficients
def process_log_coefficients(expr, max_order):
    """Extract log(t) coefficients and series expand around u=0"""
    coeffs = expr.subs(log(-1/t)==log(-1)-log(t)).coefficients(log(t))
    return [[(x[0].subs(t=1/u)).series(u==0, max_order).truncate(), x[1]]
            for x in coeffs]

def generate_log_series(Persunset, Qsunset, Rsunset, Nmax, max_log_power):
    """
    Generate all series for different logarithmic powers.

    Parameters:
    -----------
    Lsunset : function
        The L operator for sunset diagrams
    Persunset : expression
        Period expression
    Qsunset : expression
        Q expression
    Rsunset : expression
        R expression (logarithmic term)
    Nmax : int
        Maximum number of terms
    max_log_power : int
        Maximum logarithmic power

    Returns:
    --------
    dict : log_series[log_power][n] contains the series for R^log_power * Q^(n+1)
    """
    log_series = {}

    for log_pow in range(max_log_power + 1):
        log_series[log_pow] = []

        for n in range(1, Nmax):
            if log_pow == 0:
                # No logs: just Q^n terms
                #expr = Lsunset(Persunset * Qsunset^n)
                expr = Persunset * Qsunset^n
                series = (expr.subs(t=1/u)).series(u==0, Nmax).truncate()
                log_series[log_pow].append(series)
            else:
                # With logs: R^log_pow * Q^n terms
                #expr = Lsunset(Persunset * Rsunset^log_pow * Qsunset^n)
                expr = (Persunset * Rsunset^log_pow * Qsunset^n)
                # Use custom processing function
                coefflog_ser = process_log_coefficients(expr,Nmax)
                log_series[log_pow].append(coefflog_ser)

    return log_series


# Example usage:
# log_series = generate_log_series( Persunset, Qsunset, Rsunset,
#                                   Nmax=10, max_log_power=3,
#                                   process_log_coefficients=my_process_func)

def build_equation_system(Persunset, Qsunset, Rsunset, Frobenius_basis, Mellin_loop, var_groups,  Nmax, max_log_power, return_expressions=False):
    """
    Build system of equations for each logarithmic power.

    Parameters:
    -----------
    Persunset: function
        Holomorphic period
    Qsunsent: function
        exponential of the first integral of the holomorphic period
    Rsunset: function
        First integral of the holomorphic period
    Mellin_loop: function
        Mellin_expansion
    var_groups : dict
        Dictionary of symbolic variable groups indexed by log power
    Nmax : int
        Maximum number of terms
    max_log_power : int
        Maximum logarithmic power
    return_expressions : bool, optional
        If True, return both systems and expressions dict

    Returns:
    --------
    list or tuple : System of equations, optionally with expressions dictionary
    """
    systems = []
    expressions = {} if return_expressions else None

    LeadingTerm=Persunset * (-(Loop+1) * (-Rsunset)^Loop) -Mellin_loop-Frobenius_basis
    leading_log_ser=process_log_coefficients(LeadingTerm,OrderExpansion)
    log_series=generate_log_series(Persunset, Qsunset, Rsunset,  OrderExpansion, Loop)
    for target_log_pow in range(max_log_power):
        # Start with the constant term contribution
        if target_log_pow < len(leading_log_ser):
            expr = leading_log_ser[target_log_pow][0]
        else:
            expr = 0

        # Add contributions from each variable group
        for var_log_pow in range(max_log_power):
            vars_list = var_groups[var_log_pow]

            # Determine if this variable group contributes to this log power
            # A variables (log_pow=0) only contribute to target_log_pow=0
            # B variables (log_pow=1) contribute to target_log_pow=0,1
            # etc.
            if var_log_pow >= target_log_pow:
                for n in range(Nmax - 1):
                    if var_log_pow == 0:
                        # No-log terms only contribute to no-log equation
                        if target_log_pow == 0:
                            expr += vars_list[n] * log_series[var_log_pow][n]
                    else:
                        # Extract the appropriate log coefficient
                        series_data = log_series[var_log_pow][n]
                        # Find the entry with the right log power
                        for coeff_data in series_data:
                            if coeff_data[1] == target_log_pow:
                                expr += vars_list[n] * coeff_data[0]
                                break

        if return_expressions:
            expressions[target_log_pow] = expr

        # Extract coefficients of u to form equations
        equations = [x[0] for x in expr.coefficients(u)]
        systems.extend(equations)

    if return_expressions:
        return systems, expressions
    else:
        return systems


# Example usage:
# systems = build_equation_system(log_series, var_groups, Lsumsunset0_coefflog_ser,
#                                  Nmax=10, max_log_power=3)
#
# # Or with expressions for debugging:
# systems, exprs = build_equation_system(log_series, var_groups, Lsumsunset0_coefflog_ser,
#                                         Nmax=10, max_log_power=3, return_expressions=True)