Chemical Equation Balancer

# Chemical equation balancing using linear algebra
import numpy as np
from fractions import Fraction

def parse_formula(formula):
    """Parse a chemical formula like H2O into {'H': 2, 'O': 1}"""
    import re
    pattern = r"([A-Z][a-z]?)([0-9]*)"
    counts = {}
    for element, count in re.findall(pattern, formula):
        counts[element] = counts.get(element, 0) + int(count or 1)
    return counts

def balance_equation(reactants, products):
    """
    Balance a chemical equation.
    reactants: list of formula strings, e.g. ["H2", "O2"]
    products:  list of formula strings, e.g. ["H2O"]
    Returns coefficients for [reactants..., products...]
    """
    # Collect all elements
    all_compounds = reactants + products
    all_elements = set()
    for compound in all_compounds:
        all_elements.update(parse_formula(compound).keys())
    elements = sorted(all_elements)

    # Build matrix: rows=elements, cols=compounds
    # Reactants are positive, products are negative
    n = len(all_compounds)
    m = len(elements)
    matrix = [[Fraction(0)] * n for _ in range(m)]

    for j, compound in enumerate(all_compounds):
        parsed = parse_formula(compound)
        sign = 1 if j < len(reactants) else -1
        for i, element in enumerate(elements):
            matrix[i][j] = Fraction(sign * parsed.get(element, 0))

    # Gaussian elimination to find null space
    # (simplified — works for most single-reaction equations)
    # For production use, consider sympy:
    # from sympy import Matrix
    # M = Matrix(matrix)
    # coeffs = M.nullspace()[0]
    return "Use sympy for robust balancing"

# Example with sympy (recommended)
from sympy import Matrix, lcm
from functools import reduce

def balance_sympy(reactants, products):
    import re
    all_compounds = reactants + products
    elements = sorted(set(
        e for f in all_compounds
        for e, _ in re.findall(r"([A-Z][a-z]?)([0-9]*)", f)
    ))

    rows = []
    for elem in elements:
        row = []
        for i, f in enumerate(all_compounds):
            counts = {e: int(n or 1) for e, n in
                      re.findall(r"([A-Z][a-z]?)([0-9]*)", f)}
            sign = 1 if i < len(reactants) else -1
            row.append(sign * counts.get(elem, 0))
        rows.append(row)

    M = Matrix(rows)
    null = M.nullspace()[0]
    denom = reduce(lcm, [v.q for v in null])
    coeffs = [int(v * denom) for v in null]
    return coeffs

coeffs = balance_sympy(["H2", "O2"], ["H2O"])
print(coeffs)  # [2, 1, 2] -> 2H2 + O2 -> 2H2O