SigmaGeometry/augmented_rref.py at main · TeamSigmaGeometry/SigmaGeometry · 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
"""
augmented_rref.py
Pure Python Incremental RREF.
Optimized for Homogeneous "Virtual PI" coordinates.
"""

from fractions import Fraction
from typing import Dict, List, Optional, Tuple

class IncrementalAugmentedRREF:
    def __init__(self, initial_num_vars: int = 0):
        # Basis stores: (pivot_index, row_vector, provenance_vector)
        self.num_vars = initial_num_vars
        self.basis: List[Tuple[int, List[Fraction], Dict[int, Fraction]]] = []
        self.num_eqs = 0

    def _ensure_num_vars(self, n: int):
        if n > self.num_vars:
            for i in range(len(self.basis)):
                pivot, row, prov = self.basis[i]
                row.extend([Fraction(0)] * (n - self.num_vars))
            self.num_vars = n

    def _format_provenance(self, prov_dict: Dict[int, Fraction]) -> List[Fraction]:
        res = [Fraction(0)] * self.num_eqs
        for k, v in prov_dict.items():
            if k < self.num_eqs:
                res[k] = v
        return res

    def add_equation(self, coeffs: Dict[int, object]) -> Tuple[str, Optional[List[Fraction]]]:
        """Adds equation: sum(coeffs) = 0."""
        max_var = max(coeffs.keys()) if coeffs else -1
        self._ensure_num_vars(max_var + 1)

        row = [Fraction(0)] * self.num_vars
        for idx, val in coeffs.items():
            row[idx] = Fraction(val)

        prov = {self.num_eqs: Fraction(1)}

        # Gaussian Elimination
        for b_pivot, b_row, b_prov in self.basis:
            val = row[b_pivot]
            if val != 0:
                factor = val
                for i in range(b_pivot, len(b_row)):
                    row[i] -= factor * b_row[i]
                for src_id, src_coef in b_prov.items():
                    prov[src_id] = prov.get(src_id, Fraction(0)) - factor * src_coef

        # Find pivot
        pivot = -1
        for i in range(self.num_vars):
            if row[i] != 0:
                pivot = i
                break

        if pivot == -1:
            # Dependent
            if self.num_eqs in prov: del prov[self.num_eqs]
            final_prov = {k: -v for k, v in prov.items()}
            self.num_eqs += 1
            return "dependent", self._format_provenance(final_prov)
        else:
            # Independent
            scale = Fraction(1) / row[pivot]
            row = [c * scale for c in row]
            prov = {k: v * scale for k, v in prov.items()}

            insert_idx = len(self.basis)
            for i, (p, _, _) in enumerate(self.basis):
                if p > pivot:
                    insert_idx = i
                    break
            self.basis.insert(insert_idx, (pivot, row, prov))

            self.num_eqs += 1
            return "independent", None

    def reduce(self, coeffs: Dict[int, object]) -> Tuple[List[Fraction], Dict[int, Fraction]]:
        """Reduces a vector against the basis without modifying the basis."""
        max_var = max(coeffs.keys()) if coeffs else -1
        target_size = max(self.num_vars, max_var + 1)
        row = [Fraction(0)] * target_size
        for idx, val in coeffs.items():
            row[idx] = Fraction(val)

        prov = {}
        for b_pivot, b_row, b_prov in self.basis:
            if b_pivot >= len(row):
                if len(row) < len(b_row):
                     row.extend([Fraction(0)]*(len(b_row)-len(row)))

            if b_pivot < len(row):
                val = row[b_pivot]
                if val != 0:
                    factor = val
                    for i in range(b_pivot, len(b_row)):
                        row[i] -= factor * b_row[i]
                    for pid, pval in b_prov.items():
                        prov[pid] = prov.get(pid, Fraction(0)) + factor * pval

        return row, prov