diff --git a/examples/exercises/algebra/polynomial_equations_examples.py b/examples/exercises/algebra/polynomial_equations_examples.py
new file mode 100644
index 00000000..ebb54d0a
--- /dev/null
+++ b/examples/exercises/algebra/polynomial_equations_examples.py
@@ -0,0 +1,105 @@
+"""Examples of generated problems from the PolynomialEquations exercise.
+
+This file demonstrates different types of polynomial equation problems that can be generated
+at various difficulty levels.
+"""
+
+from reasoning_gym.curricula.algebra.polynomial_equations_curriculum import PolynomialEquationsCurriculum
+from reasoning_gym.exercises.algebra.polynomial_equations import PolynomialEquationsExercise
+import random
+
+def main():
+    # Initialize with fixed seed for reproducibility
+    curriculum = PolynomialEquationsCurriculum()
+    exercise = PolynomialEquationsExercise()
+    curriculum.rng = random.Random(42)
+
+    print("\n========================================\n")
+
+    # Level 0: Linear equations (ax + b = 0)
+    curriculum.set_attr_level("num_terms", 0)  # 2 terms
+    curriculum.set_attr_level("coefficient_value", 0)  # Small coefficients (1-10)
+    curriculum.set_attr_level("max_degree", 0)  # Linear equations
+    curriculum.set_attr_level("operators", 0)  # Just + operator
+    curriculum.set_attr_level("sign", 0)  # No signs
+    curriculum.set_attr_level("var_name", 0)  # Basic variables (x, y, z)
+    problem = exercise.generate(curriculum)
+    print("Level 0 (Linear Equations):")
+    print(problem)
+
+    print("\n========================================\n")
+
+    # Level 1: Quadratic equations (ax² + bx + c = 0)
+    curriculum.set_attr_level("num_terms", 1)  # 3 terms
+    curriculum.set_attr_level("coefficient_value", 1)  # Medium coefficients (1-50)
+    curriculum.set_attr_level("max_degree", 1)  # Quadratic equations
+    curriculum.set_attr_level("operators", 1)  # +, - operators
+    curriculum.set_attr_level("sign", 1)  # Allow +/-
+    curriculum.set_attr_level("var_name", 0)  # Basic variables
+    problem = exercise.generate(curriculum)
+    print("Level 1 (Quadratic Equations):")
+    print(problem)
+
+    print("\n========================================\n")
+
+    # Level 2: Cubic equations with larger coefficients
+    curriculum.set_attr_level("num_terms", 2)  # 4 terms
+    curriculum.set_attr_level("coefficient_value", 2)  # Large coefficients (1-100)
+    curriculum.set_attr_level("max_degree", 2)  # Cubic equations
+    curriculum.set_attr_level("operators", 1)  # +, - operators
+    curriculum.set_attr_level("sign", 1)  # Allow +/-
+    curriculum.set_attr_level("var_name", 1)  # All ASCII letters
+    problem = exercise.generate(curriculum)
+    print("Level 2 (Cubic Equations):")
+    print(problem)
+
+    print("\n========================================\n")
+
+    # Random Examples with Different Seeds
+    print("Random Examples (Different Seeds):")
+    for seed in range(10, 15):
+        curriculum.rng = random.Random(seed)
+        # Randomly set curriculum levels
+        curriculum.set_attr_level("num_terms", random.randint(0, 2))
+        curriculum.set_attr_level("coefficient_value", random.randint(0, 2))
+        curriculum.set_attr_level("max_degree", random.randint(0, 2))
+        curriculum.set_attr_level("operators", random.randint(0, 1))
+        curriculum.set_attr_level("sign", random.randint(0, 1))
+        curriculum.set_attr_level("var_name", random.randint(0, 2))
+        problem = exercise.generate(curriculum)
+        print(f"\nRandom Example (Seed {seed}):")
+        print(problem)
+
+    print("\n========================================\n")
+
+    # Special Cases
+    print("Special Cases:")
+
+    # Case 1: Greek variable names with high degree
+    curriculum.set_attr_level("num_terms", 2)  # 4 terms
+    curriculum.set_attr_level("coefficient_value", 1)  # Medium coefficients
+    curriculum.set_attr_level("max_degree", 2)  # Cubic equations
+    curriculum.set_attr_level("var_name", 2)  # Greek letters
+    problem = exercise.generate(curriculum)
+    print("\nGreek Variables with High Degree:")
+    print(problem)
+
+    # Case 2: Maximum terms with small coefficients
+    curriculum.set_attr_level("num_terms", 2)  # Maximum terms
+    curriculum.set_attr_level("coefficient_value", 0)  # Small coefficients
+    curriculum.set_attr_level("max_degree", 1)  # Quadratic equations
+    problem = exercise.generate(curriculum)
+    print("\nMaximum Terms with Small Coefficients:")
+    print(problem)
+
+    # Case 3: Linear equation with large coefficients
+    curriculum.set_attr_level("num_terms", 0)  # 2 terms
+    curriculum.set_attr_level("coefficient_value", 2)  # Large coefficients
+    curriculum.set_attr_level("max_degree", 0)  # Linear equations
+    curriculum.set_attr_level("var_name", 0)  # Basic variables
+    problem = exercise.generate(curriculum)
+    print("\nLinear Equation with Large Coefficients:")
+    print(problem)
+
+if __name__ == "__main__":
+    main() 
\ No newline at end of file
diff --git a/reasoning_gym/algebra/__init__.py b/reasoning_gym/algebra/__init__.py
index 69d4b91e..40e22837 100644
--- a/reasoning_gym/algebra/__init__.py
+++ b/reasoning_gym/algebra/__init__.py
@@ -1,9 +1,7 @@
-from .polynomial_equations import PolynomialEquationsConfig, PolynomialEquationsDataset
-from .simple_equations import SimpleEquationsConfig, SimpleEquationsDataset
+from .polynomial_equations import PolynomialEquationsExercise
+# from .simple_equations import SimpleEquationsConfig, SimpleEquationsDataset
 
 __all__ = [
-    "SimpleEquationsDataset",
-    "SimpleEquationsConfig",
-    "PolynomialEquationsConfig",
-    "PolynomialEquationsDataset",
+    # "SimpleEquationsDataset",
+    "PolynomialEquationsExercise",
 ]
diff --git a/reasoning_gym/algebra/polynomial_equations.py b/reasoning_gym/algebra/polynomial_equations.py
index ed7e857f..fea264f9 100644
--- a/reasoning_gym/algebra/polynomial_equations.py
+++ b/reasoning_gym/algebra/polynomial_equations.py
@@ -1,150 +1,115 @@
-import random
-import string
-from dataclasses import dataclass
-from typing import Optional, Tuple
+"""
+Polynomial equation exercise that generates equations and finds their real solutions.
+"""
 
-from sympy import Eq, Symbol, expand, solve
+from typing import Dict, Any
+from sympy import Symbol, expand, solve, Eq, parse_expr
 
-from ..factory import ProceduralDataset, register_dataset
-
-
-@dataclass
-class PolynomialEquationsConfig:
+class PolynomialEquationsExercise:
     """
-    Configuration for polynomial equation task generation.
+    Generates random polynomial equations and finds their real solutions.
+    The polynomial is formed by summing random terms of the form: coeff * x^exponent.
+    Then we solve "polynomial_expr = 0" using Sympy.
     """
 
-    min_terms: int = 2  # Minimum number of polynomial terms
-    max_terms: int = 4  # Maximum number of polynomial terms
-    min_value: int = 1  # Minimum value for coefficients
-    max_value: int = 100  # Maximum value for coefficients
-    min_degree: int = 1  # Minimum polynomial degree
-    max_degree: int = 3  # Maximum polynomial degree
-    operators: Tuple[str, ...] = (
-        "+",
-        "-",
-    )  # Allowed operators between terms, Avoid adding '*' or '/' because they will affect the degree
-    seed: Optional[int] = None
-    size: int = 500
+    def __init__(self):
+        self.curriculum = None
 
-    def validate(self) -> None:
-        """Validate configuration parameters."""
-        assert self.min_terms > 0, "min_terms must be positive."
-        assert self.max_terms >= self.min_terms, "max_terms must be >= min_terms."
-
-        assert self.min_value > 0, "min_value must be positive."
-        assert self.max_value >= self.min_value, "max_value must be >= min_value."
-
-        assert self.min_degree >= 1, "min_degree must be >= 1."
-        assert self.max_degree >= self.min_degree, "max_degree must be >= min_degree."
-
-        allowed_ops = {"+", "-"}
-        assert len(self.operators) > 0, "operators tuple cannot be empty."
-        assert all(op in allowed_ops for op in self.operators), "Invalid operator found. Must be a subset of {+, -}."
-
-
-class PolynomialEquationsDataset(ProceduralDataset):
-    """
-    Generates random polynomial equations of degree in [min_degree, max_degree].
-    - The polynomial is formed by summing random terms of the form: coeff * x^exponent.
-    - Then we solve "polynomial_expr = 0" using Sympy.
-    - The solution may be real or complex; we filter real solutions by default for simplicity.
-    """
-
-    def __init__(self, config: PolynomialEquationsConfig):
-        self._prompt_templates = [
-            "Find the real value(s) of {variable} in the equation: {polynomial_expanded} = 0",
-            "Solve for real {variable}: {polynomial_expanded} = 0",
-            "Determine the real value(s) of {variable} tha satisfies: {polynomial_expanded} = 0",
-            "Solve the polynomial equation for real {variable}:\n{polynomial_expanded} = 0",
-        ]
-        super().__init__(config=config, seed=config.seed, size=config.size)
-
-    def __getitem__(self, idx: int) -> dict:
+    def generate(self, curriculum: Any) -> Dict[str, Any]:
         """
-        Generate a single polynomial equation item.
+        Generate a polynomial equation problem using the curriculum.
 
         Returns:
-            A dict with:
+            Dict containing:
                 - question: str (e.g. "Solve the polynomial equation: 2*x^2 - 3*x + 1 = 0")
                 - answer: str (the sorted list of real solutions, e.g. "[0.5, 1.0]")
-                - metadata: dict with details (polynomial_expr, degree, etc.)
+                - metadata: dict with details (polynomial_expr, symbolic_info, etc.)
         """
-        rng = random.Random(self.seed + idx)
+        self.curriculum = curriculum
+        template = curriculum.get_template(curriculum.rng)
+        return template.eval(self, curriculum.rng)
 
-        # Get variable and generate polynomial equation in standard form
-        variable = self._get_variable(rng)
-        degree = rng.randint(self.config.min_degree, self.config.max_degree)
-        polynomial_expr = self._generate_polynomial_expr(rng, variable, degree)
-        polynomial_expanded = expand(polynomial_expr)
+    def _parse_expression(self, metadata: Dict[str, Any]) -> Dict[str, Any]:
+        """
+        Parse the template metadata into structured data.
 
-        # Solve the polynomial = 0
-        # We filter real solutions only
-        solutions = solve(Eq(polynomial_expanded, 0), variable, dict=False)
-        real_solutions = []
-        for sol in solutions:
-            if sol.is_real:
-                # Evaluate symbolic solution to a floating approximation
-                real_solutions.append(float(sol.evalf()))
-        real_solutions.sort()
-        answer_str = str(real_solutions)
-
-        return {
-            "question": rng.choice(self._prompt_templates).format(
-                variable=variable,
-                polynomial_expanded=polynomial_expanded,
-            ),
-            "answer": answer_str,
-            "metadata": {
-                "polynomial_expr": str(polynomial_expanded),
-                "variable": variable,
-                "degree": degree,
-                "real_solutions": real_solutions,
+        The metadata structure is expected to be:
+        {
+            "expression": {
+                "term_0": {
+                    "sign": str,      # "" or "-"
+                    "coeff": str,     # coefficient value with "*" if needed
+                    "variable": str,  # variable name or empty
+                    "exponent": str   # "**N" for degree N > 1, or empty
+                },
+                "term_1": {...},      # Same structure as term_0
+                ...,
+                "op_0": str,         # "+" or "-" between terms
+                "op_1": str,         # More operators if needed
+                ...
             },
+            "variable": {
+                "var": str           # The variable name used in the equation
+            }
         }
 
-    def _get_variable(self, rng: random.Random) -> str:
-        """Get a random lowercase variable name"""
-        return rng.choice(string.ascii_lowercase)
-
-    def _generate_polynomial_expr(self, rng: random.Random, variable: Symbol, degree: int):
-        """
-        Randomly generate a polynomial expression of 'degree'.
-        We'll use the config parameters:
-            - min_terms, max_terms: how many total terms to combine
-            - min_value, max_value: range for coefficients
-            - operators: to decide sign flips or direct addition
-
-         Args:
-            rng: Random number generator
-            variable: Variable symbol to use in equation
-            degree: Highest degree. We ensure that there is at least one term with exponent=degree
-
+        Args:
+            metadata: Raw metadata from template evaluation
         Returns:
-            Polynomial string
+            Dictionary containing:
+                - terms: List[str] of formatted term strings
+                - operators: List[str] of operators between terms
+                - variable: str, the variable name used
         """
-        x = Symbol(variable)
+        expr_parts = metadata["expression"]
 
-        # Choose the number of terms and their respective degrees
-        num_terms = rng.randint(self.config.min_terms, self.config.max_terms)
-        # Keep track of exponents, exponents can repeat or skip but we force the highest exponent
-        chosen_exponents = [degree]
-        # Fill the rest randomly in [0, degree]
-        for _ in range(num_terms - 1):
-            exp = rng.randint(0, degree)
-            chosen_exponents.append(exp)
+        # Extract terms and operators
+        terms = []
+        operators = []
+        i = 0
+        while f"term_{i}" in expr_parts:
+            term_dict = expr_parts[f"term_{i}"]
+            terms.append("".join(term_dict[k] for k in ("sign", "coeff", "variable", "exponent")))
+            # Get operator if it exists
+            if f"op_{i}" in expr_parts:
+                operators.append(expr_parts[f"op_{i}"])
+            i += 1
 
-        # Now build the polynomial expression: sum_{term}( coeff * x^exponent ), with optional sign
-        polynomial_expr = 0
-        for exp in chosen_exponents:
-            coeff = rng.randint(self.config.min_value, self.config.max_value)
-            # If '-' in operators, we can randomly flip the sign
-            if "-" in self.config.operators and rng.random() < 0.5:
-                coeff = -coeff
-            term_expr = coeff * (x**exp)
-            polynomial_expr += term_expr
+        return {
+            "terms": terms,
+            "operators": operators,
+            "variable": metadata["variable"]["var"]
+        }
 
-        return polynomial_expr
+    def _evaluate_expression(self, parsed: Dict[str, Any]) -> str:
+        """
+        Evaluate the polynomial equation and find its real solutions.
 
+        Args:
+            parsed: Dictionary containing parsed expression data
+        Returns:
+            String representation of the sorted list of real solutions
+        """
+        # Create sympy symbol from parsed variable
+        var = Symbol(parsed["variable"])
 
-register_dataset("polynomial_equations", PolynomialEquationsDataset, PolynomialEquationsConfig)
+        # Build expression from parsed terms and operators
+        expr = parsed["terms"][0]
+        for i, op in enumerate(parsed["operators"]):
+            expr = f"{expr} {op} {parsed['terms'][i + 1]}"
+
+        try:
+            sympy_expr = parse_expr(expr, local_dict={parsed["variable"]: var})
+            expanded = expand(sympy_expr)
+            solutions = solve(Eq(expanded, 0), var, dict=False)
+
+            # Filter and sort real solutions
+            real_solutions = []
+            for sol in solutions:
+                if sol.is_real:
+                    real_solutions.append(float(sol.evalf()))
+            real_solutions.sort()
+
+            return str(real_solutions)
+        except Exception as e:
+            return f"Error evaluating expression: {expr}\nError: {str(e)}"
diff --git a/reasoning_gym/curricula/__init__.py b/reasoning_gym/curricula/__init__.py
index 6a212f27..f4c071b1 100644
--- a/reasoning_gym/curricula/__init__.py
+++ b/reasoning_gym/curricula/__init__.py
@@ -1,4 +1,4 @@
-# from .algebra import *
+from .algebra import *
 # from .algorithmic import *
 from .arithmetic import *
 # from .code import *
@@ -11,8 +11,9 @@ from .arithmetic import *
 # Re-export all Curriculum classes
 __all__ = []
 for module in [
-    arithmetic
-    # algebra, algorithmic, arithmetic, code,
+    arithmetic,
+    algebra,
+    # algorithmic, arithmetic, code,
     # cognition, games, geometry, graphs, logic
 ]:
     __all__.extend([name for name in module.__all__ if name.endswith('Curriculum')]) 
\ No newline at end of file
diff --git a/reasoning_gym/curricula/algebra/__init__.py b/reasoning_gym/curricula/algebra/__init__.py
new file mode 100644
index 00000000..4718d3cb
--- /dev/null
+++ b/reasoning_gym/curricula/algebra/__init__.py
@@ -0,0 +1,7 @@
+from .polynomial_equations_curriculum import PolynomialEquationsCurriculum
+# from .simple_equations import SimpleEquationsConfig, SimpleEquationsDataset
+
+__all__ = [
+    # "SimpleEquationsCurriculum",
+    "PolynomialEquationsCurriculum",
+]
diff --git a/reasoning_gym/curricula/algebra/polynomial_equations_curriculum.py b/reasoning_gym/curricula/algebra/polynomial_equations_curriculum.py
new file mode 100644
index 00000000..b367ec57
--- /dev/null
+++ b/reasoning_gym/curricula/algebra/polynomial_equations_curriculum.py
@@ -0,0 +1,150 @@
+"""
+Curriculum definition for polynomial equation exercises.
+"""
+
+import string
+from typing import Dict, Any
+from reasoning_gym.core.base_curriculum import BaseCurriculum
+from reasoning_gym.core.attributes import AttributeDefinition, AttributeType
+from reasoning_gym.core.template import Template
+
+
+class PolynomialEquationsCurriculum(BaseCurriculum):
+    def __init__(self):
+        super().__init__("PolynomialEquationsCurriculum")
+
+    def _init_curriculum(self) -> None:
+        """Initialize the polynomial equations curriculum configuration"""
+        # Define valid attribute types
+        self._valid_types = {
+            AttributeType.STATIC,   # For operators
+            AttributeType.UBOUND,   # For ranges like num_terms, degree, value
+            AttributeType.APPEND,    # For operators that accumulate
+            AttributeType.APPEND_LIST    # For variables that accumulate
+        }
+
+        # Define attributes
+        self._attributes = {
+            "num_terms": AttributeDefinition(
+                levels=[2, 3, 4],  # From min_terms/max_terms
+                default_level=0,
+                description="Number of polynomial terms",
+                attr_type=AttributeType.UBOUND,
+                min_value=2  # Ensure at least 2 terms
+            ),
+            "coefficient_value": AttributeDefinition(
+                levels=[10, 50, 100],  # From min_value/max_value
+                default_level=0,
+                description="Maximum value for coefficients",
+                attr_type=AttributeType.UBOUND,
+                min_value=1  # Ensure non-zero coefficients
+            ),
+            "max_degree": AttributeDefinition(
+                levels=[1, 2, 3],  # From min_degree/max_degree
+                default_level=0,
+                description="Maximum polynomial degree",
+                attr_type=AttributeType.UBOUND
+            ),
+            "operators": AttributeDefinition(
+                levels=["+", "-"],  # Keep original operators
+                default_level=0,
+                description="Allowed operators between terms",
+                attr_type=AttributeType.APPEND
+            ),
+            "sign": AttributeDefinition(
+                levels=["", "-"],  # Remove explicit + sign
+                default_level=0,
+                description="Sign of the coefficient",
+                attr_type=AttributeType.APPEND
+            ),
+            "var_name": AttributeDefinition(
+                levels=[list("xyz"), list(string.ascii_lowercase + string.ascii_uppercase), list("αβγρθφψω")],
+                default_level=0,
+                description="Variables to use in polynomials",
+                attr_type=AttributeType.APPEND_LIST
+            )
+        }
+
+        # Define templates with symbolic placeholders
+        self._templates = [
+            Template(
+                template="Find the real value(s) of {variable} in the equation: {expression} = 0",
+                parts={"expression": "polynomial_expression", "variable": "variable_name"}
+            ),
+            Template(
+                template="Solve for real {variable}: {expression} = 0",
+                parts={"expression": "polynomial_expression", "variable": "variable_name"}
+            ),
+            Template(
+                template="Determine the real value(s) of {variable} that satisfies: {expression} = 0",
+                parts={"expression": "polynomial_expression", "variable": "variable_name"}
+            ),
+            Template(
+                template="Solve the polynomial equation for real {variable}:\n{expression} = 0",
+                parts={"expression": "polynomial_expression", "variable": "variable_name"}
+            )
+        ]
+
+		# TODO: must always be at least one var
+        # Define symbolic structure
+        self._symbolic = {
+            # Define composition templates
+            "templates": {
+                # Variable name template
+                "variable_name": lambda refs: {
+                    "template": "{var}",
+                    "parts": {
+                        "var": lambda refs=refs: refs["var"](refs)
+                    }
+                },
+                # Expression structure
+                "polynomial_expression": lambda refs: (
+                    n_terms := refs["num_terms"](),
+                    {
+                        "template": "{term_0}" + "".join(f" {{op_{i}}} {{term_{i+1}}}" 
+                                                        for i in range(n_terms - 1)),
+                        "parts": {
+                            **{f"term_{i}": "polynomial_term" for i in range(n_terms)},
+                            **{f"op_{i}": lambda refs=refs: refs["operator"](refs)()
+                               for i in range(n_terms - 1)}
+                        }
+                    }
+                )[-1],
+                # Term structure
+                "polynomial_term": lambda refs: (
+                    coeff := refs["coefficient"](refs)(),
+                    deg := refs["degree"](refs)(),
+                    var := refs["var"](refs),
+                    {
+                        "template": "{sign}{coeff}{variable}{exponent}",
+                        "parts": {
+                            "sign": lambda refs=refs: refs["sign_term"](refs)(),
+                            "coeff": lambda: (
+                                f"{coeff}*" if (deg > 0 and coeff != 1) else
+                                f"{coeff}" if deg == 0 else
+                                ""  # No coefficient if 1 and has variable
+                            ),
+                            "variable": lambda: (
+                                "" if deg == 0 else
+                                f"{var}"
+                            ),
+                            "exponent": lambda: (
+                                "" if deg <= 1 else
+                                f"**{deg}"
+                            )
+                        }
+                    }
+                )[-1]
+            },
+            # Define shared variables that need to be consistent across templates
+            "shared_vars": {
+                "var": lambda refs: refs["var_name"]()
+            },
+            # Define value generators
+            "generators": {
+                "coefficient": lambda refs: lambda: refs["coefficient_value"](),
+                "degree": lambda refs: lambda: refs["max_degree"](),
+                "operator": lambda refs: lambda: refs["operators"](),
+                "sign_term": lambda refs: lambda: refs["sign"]()
+            }
+        }
\ No newline at end of file
diff --git a/reasoning_gym/exercises/__init__.py b/reasoning_gym/exercises/__init__.py
index 46d171c9..22f3c604 100644
--- a/reasoning_gym/exercises/__init__.py
+++ b/reasoning_gym/exercises/__init__.py
@@ -1,4 +1,4 @@
-# from .algebra import *
+from .algebra import *
 # from .algorithmic import *
 from .arithmetic import *
 # from .code import *
@@ -11,8 +11,9 @@ from .arithmetic import *
 # Re-export all Dataset classes
 __all__ = []
 for module in [
-    arithmetic
-    # algebra, algorithmic, arithmetic, code,
+    arithmetic,
+    algebra,
+    # algorithmic, arithmetic, code,
     # cognition, games, geometry, graphs, logic
 ]:
     __all__.extend([name for name in module.__all__ if name.endswith('Exercise')]) 
\ No newline at end of file
diff --git a/tests/test_chain_sum.py b/tests/test_chain_sum.py
index 9db63af1..e40b1435 100644
--- a/tests/test_chain_sum.py
+++ b/tests/test_chain_sum.py
@@ -397,7 +397,7 @@ class TestChainSumGenerate(unittest.TestCase):
     def test_comprehensive_random_evaluation(self):
         """Test 1000 random problems across all levels to verify correct evaluation"""
         num_samples = 1000
-        
+
         # Statistics tracking
         stats = {
             'operator_counts': {},      # Count of each operator used
diff --git a/tests/test_polynomial_equations.py b/tests/test_polynomial_equations.py
index 6e1bb0c0..8547ab3d 100644
--- a/tests/test_polynomial_equations.py
+++ b/tests/test_polynomial_equations.py
@@ -1,117 +1,577 @@
-import pytest
-from sympy import Symbol, sympify
+from reasoning_gym.curricula.algebra.polynomial_equations_curriculum import PolynomialEquationsCurriculum
+from reasoning_gym.exercises.algebra.polynomial_equations import PolynomialEquationsExercise
+import unittest
+import random
+from sympy import solve, Symbol, Eq, parse_expr
 
-from reasoning_gym import create_dataset
-from reasoning_gym.algebra.polynomial_equations import PolynomialEquationsConfig, PolynomialEquationsDataset
+class TestPolynomialEquationsParsing(unittest.TestCase):
+    """Test parsing of polynomial expressions and terms"""
 
+    def setUp(self):
+        self.exercise = PolynomialEquationsExercise()
 
-def test_polynomial_config_validation():
-    """Test that invalid configs raise appropriate errors"""
-    with pytest.raises(AssertionError):
-        PolynomialEquationsConfig(min_terms=0).validate()
+    def test_parse_expression(self):
+        """Test parsing of polynomial expressions"""
+        test_metadata = {
+            'type': 'direct',
+            'executed_parts': {
+                'terms': ['2*x**2', '3*x', '1'],
+                'operators': ['+', '+'],
+                'variable': 'x'
+            }
+        }
 
-    with pytest.raises(AssertionError):
-        PolynomialEquationsConfig(min_value=0).validate()
+        parsed = test_metadata['executed_parts']
+        self.assertEqual(parsed["terms"], ["2*x**2", "3*x", "1"])
+        self.assertEqual(parsed["operators"], ["+", "+"])
+        self.assertEqual(parsed["variable"], "x")
 
-    with pytest.raises(AssertionError):
-        PolynomialEquationsConfig(min_degree=0, max_degree=3).validate()
+    def test_parse_negative_terms(self):
+        """Test parsing of expressions with negative terms"""
+        test_metadata = {
+            'type': 'direct',
+            'executed_parts': {
+                'terms': ['-2*x**2', '4*x'],
+                'operators': ['+'],
+                'variable': 'x'
+            }
+        }
 
-    with pytest.raises(AssertionError):
-        PolynomialEquationsConfig(min_degree=4, max_degree=3).validate()
+        parsed = test_metadata['executed_parts']
+        self.assertEqual(parsed["terms"], ["-2*x**2", "4*x"])
+        self.assertEqual(parsed["operators"], ["+"])
+        self.assertEqual(parsed["variable"], "x")
 
-    with pytest.raises(AssertionError):
-        PolynomialEquationsConfig(operators=("^",)).validate()
+class TestPolynomialEquationsEvaluation(unittest.TestCase):
+    """Test evaluation of polynomial equations"""
 
+    def setUp(self):
+        self.exercise = PolynomialEquationsExercise()
 
-def test_polynomial_equations_dataset_basic():
-    """Test dataset creation and length"""
-    dataset_size = 50
-    config = PolynomialEquationsConfig(
-        min_terms=2,
-        max_terms=3,
-        min_value=1,
-        max_value=5,
-        min_degree=1,
-        max_degree=2,
-        seed=42,
-        size=dataset_size,
-    )
+    def test_quadratic_equation(self):
+        """Test evaluation of quadratic equations"""
+        parsed = {
+            "terms": ["x**2", "-5*x", "6"],
+            "operators": ["+", "+"],
+            "variable": "x"
+        }
+        result = self.exercise._evaluate_expression(parsed)
+        expected = "[2.0, 3.0]"  # x^2 - 5x + 6 = 0 has roots at x = 2 and x = 3
+        self.assertEqual(result, expected)
 
-    dataset = PolynomialEquationsDataset(config)
+    def test_linear_equation(self):
+        """Test evaluation of linear equations"""
+        parsed = {
+            "terms": ["2*x", "-4"],
+            "operators": ["+"],
+            "variable": "x"
+        }
+        result = self.exercise._evaluate_expression(parsed)
+        expected = "[2.0]"  # 2x - 4 = 0 has root at x = 2
+        self.assertEqual(result, expected)
 
-    assert len(dataset) == dataset_size
+    def test_no_real_solutions(self):
+        """Test equations with no real solutions"""
+        parsed = {
+            "terms": ["x**2", "1"],
+            "operators": ["+"],
+            "variable": "x"
+        }
+        result = self.exercise._evaluate_expression(parsed)
+        expected = "[]"  # x^2 + 1 = 0 has no real solutions
+        self.assertEqual(result, expected)
 
+class TestPolynomialEquationsGeneration(unittest.TestCase):
+    """Test problem generation"""
 
-def test_polynomial_equations_dataset_items():
-    """Test that generated items have correct structure"""
-    ds = create_dataset(
-        "polynomial_equations",
-        min_terms=2,
-        max_terms=3,
-        min_value=1,
-        max_value=5,
-        min_degree=1,
-        max_degree=2,
-        size=3,
-        seed=100,
-    )
+    def setUp(self):
+        self.curriculum = PolynomialEquationsCurriculum()
+        self.exercise = PolynomialEquationsExercise()
+        self.rng = random.Random(42)
+        self.curriculum.rng = self.rng
 
-    for item in ds:
-        assert "question" in item
-        assert "answer" in item
-        assert "metadata" in item
+    def test_problem_structure(self):
+        """Test that generated problems have the correct structure"""
+        problem = self.exercise.generate(self.curriculum)
 
-        # Check metadata
-        assert isinstance(item["metadata"]["polynomial_expr"], str)
-        assert isinstance(item["metadata"]["variable"], str)
-        assert isinstance(item["metadata"]["degree"], int)
-        assert isinstance(item["metadata"]["real_solutions"], list)
+        # Check basic structure
+        self.assertIn("question", problem)
+        self.assertIn("answer", problem)
+        self.assertIn("metadata", problem)
 
-        # Check polynomial_expr existence
-        poly_str = item["metadata"]["polynomial_expr"]
-        # Ensure it can parse with sympy
-        sympify(poly_str)
+        # Check metadata structure
+        metadata = problem["metadata"]
+        self.assertEqual(metadata["type"], "direct")
+        self.assertIn("executed_parts", metadata)
+        executed_parts = metadata["executed_parts"]
+        self.assertIn("terms", executed_parts)
+        self.assertIn("operators", executed_parts)
+        self.assertIn("variable", executed_parts)
 
+    def test_term_generation(self):
+        """Test generation of polynomial terms"""
+        # Set curriculum to basic settings
+        self.curriculum.set_attr_level("coefficient_value", 0)  # 1-10
+        self.curriculum.set_attr_level("max_degree", 0)  # degree 1
+        self.curriculum.set_attr_level("sign", 0)  # No signs
 
-def test_polynomial_equations_dataset_deterministic():
-    """Test dataset reproducibility with fixed seed."""
-    cfg = PolynomialEquationsConfig(seed=999, size=3)
-    ds1 = PolynomialEquationsDataset(cfg)
-    ds2 = PolynomialEquationsDataset(cfg)
+        problem = self.exercise.generate(self.curriculum)
+        executed_parts = problem["metadata"]["executed_parts"]
 
-    for i in range(len(ds1)):
-        assert ds1[i] == ds2[i], "Polynomial datasets with same seed should match exactly."
+        # Check we have at least one term
+        self.assertTrue(len(executed_parts["terms"]) > 0)
 
+        # Check first term format
+        first_term = executed_parts["terms"][0]
+        self.assertTrue(isinstance(first_term, str))
+        self.assertTrue(first_term.replace('*', '').replace('x', '').replace('-', '').replace('.', '').isdigit() or
+                       first_term == 'x')
 
-def test_polynomial_solutions_evaluation():
-    """Test that real_solutions satisfy the polynomial equation."""
-    ds = create_dataset(
-        "polynomial_equations",
-        min_terms=2,
-        max_terms=4,
-        min_value=1,
-        max_value=10,
-        min_degree=1,
-        max_degree=3,
-        size=5,
-        seed=42,
-    )
+    def test_operator_generation(self):
+        """Test generation of operators"""
+        self.curriculum.set_attr_level("operators", 1)  # +, -
+        self.curriculum.set_attr_level("num_terms", 0)  # 2 terms
 
-    for item in ds:
-        # Extract the polynomial expression and solutions
-        poly_str = item["metadata"]["polynomial_expr"]
-        real_solutions = item["metadata"]["real_solutions"]
-        x = Symbol(item["metadata"]["variable"])
-        # Parse the polynomial expression
-        poly_expr = sympify(poly_str)
+        problem = self.exercise.generate(self.curriculum)
+        executed_parts = problem["metadata"]["executed_parts"]
 
-        # Verify that each solution satisfies the polynomial
-        for solution in real_solutions:
-            # Evaluate the expression with the solution substituted
-            evaluated_value = poly_expr.subs(x, solution)
+        # Check we have operators for n-1 terms
+        self.assertEqual(len(executed_parts["operators"]), len(executed_parts["terms"]) - 1)
 
-            # Ensure the evaluated value is close to zero (numerical stability threshold)
-            assert abs(evaluated_value) < 1e-6, (
-                f"Solution {solution} does not satisfy the polynomial {poly_str}. "
-                f"Evaluated value: {evaluated_value}"
-            )
+        # Check operator is valid
+        if executed_parts["operators"]:
+            self.assertIn(executed_parts["operators"][0], ["+", "-"])
+
+class TestPolynomialEquationsComprehensive(unittest.TestCase):
+    """Comprehensive tests for polynomial equations"""
+
+    def setUp(self):
+        self.curriculum = PolynomialEquationsCurriculum()
+        self.exercise = PolynomialEquationsExercise()
+        self.rng = random.Random(42)
+        self.curriculum.rng = self.rng
+
+    def test_variable_consistency(self):
+        """Test that the same variable is used consistently throughout the equation"""
+        num_samples = 50
+
+        for _ in range(num_samples):
+            problem = self.exercise.generate(self.curriculum)
+            executed_parts = problem["metadata"]["executed_parts"]
+            var_name = executed_parts["variable"]
+
+            # Check variable appears in question
+            self.assertIn(var_name, problem["question"])
+
+            # Check variable is used consistently in terms
+            for term in executed_parts["terms"]:
+                if var_name in term:  # If term has a variable
+                    self.assertIn(var_name, term)
+
+    def test_coefficient_ranges(self):
+        """Test that coefficients are within expected ranges"""
+        self.curriculum.set_attr_level("coefficient_value", 0)  # 1-10
+        num_samples = 50
+
+        for _ in range(num_samples):
+            problem = self.exercise.generate(self.curriculum)
+            executed_parts = problem["metadata"]["executed_parts"]
+
+            for term in executed_parts["terms"]:
+                # Extract coefficient if term has one
+                if '*' in term:
+                    coeff = term.split('*')[0]
+                    if coeff and coeff != '-':  # Skip if empty or just a minus sign
+                        coeff = float(coeff)
+                        self.assertLessEqual(abs(coeff), 10)
+                        self.assertGreater(abs(coeff), 0)
+
+    def test_degree_constraints(self):
+        """Test that polynomial degrees respect the curriculum settings"""
+        self.curriculum.set_attr_level("max_degree", 0)  # Level 0 means max degree 1
+        num_samples = 50
+
+        for _ in range(num_samples):
+            problem = self.exercise.generate(self.curriculum)
+            executed_parts = problem["metadata"]["executed_parts"]
+
+            max_degree = 0
+            for term in executed_parts["terms"]:
+                if "**" in term:
+                    degree = int(term.split("**")[1])
+                    max_degree = max(max_degree, degree)
+                elif executed_parts["variable"] in term:  # Variable without exponent means degree 1
+                    max_degree = max(max_degree, 1)
+
+            self.assertLessEqual(max_degree, 1)
+
+    def test_solution_validity(self):
+        """Test that generated solutions are valid"""
+        num_samples = 50
+
+        for _ in range(num_samples):
+            problem = self.exercise.generate(self.curriculum)
+            executed_parts = problem["metadata"]["executed_parts"]
+
+            # Parse the answer string to get solutions
+            solutions = eval(problem["answer"])  # Safe since we control the input
+
+            if solutions:  # If there are real solutions
+                # Verify each solution satisfies the equation
+                var = Symbol(executed_parts["variable"])
+                expr = executed_parts["terms"][0]
+
+                # Reconstruct the expression
+                for i, term in enumerate(executed_parts["terms"][1:], 1):
+                    expr += f" {executed_parts['operators'][i-1]} {term}"
+
+                # Verify each solution
+                sympy_expr = parse_expr(expr)
+                for sol in solutions:
+                    result = abs(float(sympy_expr.subs(var, sol)))
+                    self.assertAlmostEqual(result, 0, places=10)
+
+    def test_comprehensive_random_evaluation(self):
+        """Test 1000 random problems across all levels to verify correct generation and evaluation"""
+        num_samples = 1000
+
+        # Statistics tracking
+        stats = {
+            'operator_counts': {},      # Count of each operator used
+            'degree_counts': {},        # Count of polynomial degrees
+            'term_counts': {},          # Distribution of number of terms
+            'variable_counts': {},      # Count of each variable used
+            'coefficient_stats': {      # Track coefficient statistics
+                'min': float('inf'),
+                'max': float('-inf'),
+                'total': 0,
+                'count': 0,
+                'unique': set()
+            },
+            'solution_stats': {         # Track solution statistics
+                'no_solutions': 0,      # Count of equations with no real solutions
+                'one_solution': 0,      # Count of equations with exactly one solution
+                'two_solutions': 0,     # Count of equations with exactly two solutions
+                'min': float('inf'),    # Minimum solution value
+                'max': float('-inf'),   # Maximum solution value
+            },
+            'level_distribution': {     # Track curriculum level usage
+                'max_degree': {},
+                'num_terms': {},
+                'coefficient_value': {},
+                'operators': {},
+                'sign': {},
+                'var_name': {}
+            }
+        }
+
+        for _ in range(num_samples):
+            # Randomly set curriculum levels
+            levels = {
+                'max_degree': self.rng.randint(0, 2),
+                'num_terms': self.rng.randint(0, 2),
+                'coefficient_value': self.rng.randint(0, 2),
+                'operators': self.rng.randint(0, 1),
+                'sign': self.rng.randint(0, 1),
+                'var_name': self.rng.randint(0, 2)
+            }
+
+            # Update level distribution stats
+            for attr, level in levels.items():
+                stats['level_distribution'][attr][level] = stats['level_distribution'][attr].get(level, 0) + 1
+
+            # Set curriculum levels
+            for attr, level in levels.items():
+                self.curriculum.set_attr_level(attr, level)
+
+            problem = self.exercise.generate(self.curriculum)
+            executed_parts = problem["metadata"]["executed_parts"]
+            terms = executed_parts["terms"]
+            operators = executed_parts["operators"]
+            variable = executed_parts["variable"]
+
+            # Update operator statistics
+            for op in operators:
+                stats['operator_counts'][op] = stats['operator_counts'].get(op, 0) + 1
+
+            # Update term count statistics
+            num_terms = len(terms)
+            stats['term_counts'][num_terms] = stats['term_counts'].get(num_terms, 0) + 1
+
+            # Update variable statistics
+            stats['variable_counts'][variable] = stats['variable_counts'].get(variable, 0) + 1
+
+            # Calculate and update degree statistics
+            max_degree = 0
+            for term in terms:
+                if "**" in term:
+                    degree = int(term.split("**")[1])
+                    max_degree = max(max_degree, degree)
+                elif variable in term:  # Variable without exponent means degree 1
+                    max_degree = max(max_degree, 1)
+            stats['degree_counts'][max_degree] = stats['degree_counts'].get(max_degree, 0) + 1
+
+            # Update coefficient statistics
+            for term in terms:
+                if '*' in term:
+                    coeff = term.split('*')[0]
+                    if coeff and coeff not in ['-', '+']:
+                        try:
+                            value = abs(float(coeff))
+                            stats['coefficient_stats']['min'] = min(stats['coefficient_stats']['min'], value)
+                            stats['coefficient_stats']['max'] = max(stats['coefficient_stats']['max'], value)
+                            stats['coefficient_stats']['total'] += value
+                            stats['coefficient_stats']['count'] += 1
+                            stats['coefficient_stats']['unique'].add(value)
+                        except ValueError:
+                            # Skip if coefficient is not a number (e.g., just a variable)
+                            continue
+
+            # Update solution statistics
+            solutions = eval(problem["answer"])  # Safe since we control the input
+            num_solutions = len(solutions)
+            if num_solutions == 0:
+                stats['solution_stats']['no_solutions'] += 1
+            elif num_solutions == 1:
+                stats['solution_stats']['one_solution'] += 1
+                stats['solution_stats']['min'] = min(stats['solution_stats']['min'], solutions[0])
+                stats['solution_stats']['max'] = max(stats['solution_stats']['max'], solutions[0])
+            elif num_solutions == 2:
+                stats['solution_stats']['two_solutions'] += 1
+                stats['solution_stats']['min'] = min(stats['solution_stats']['min'], min(solutions))
+                stats['solution_stats']['max'] = max(stats['solution_stats']['max'], max(solutions))
+
+            # Verify solution correctness
+            if solutions:  # If there are real solutions
+                var = Symbol(variable)
+                expr = terms[0]
+                for i, term in enumerate(terms[1:], 1):
+                    expr += f" {operators[i-1]} {term}"
+
+                # Create local dict with the variable symbol
+                local_dict = {variable: var}
+                sympy_expr = parse_expr(expr, local_dict=local_dict)
+                for sol in solutions:
+                    result = abs(float(sympy_expr.subs(var, sol)))
+                    self.assertAlmostEqual(result, 0, places=10)
+
+        # Print comprehensive statistics
+        print("\nComprehensive Random Evaluation Statistics:")
+        print("-" * 50)
+
+        print("\nOperator Distribution:")
+        total_ops = sum(stats['operator_counts'].values())
+        for op, count in sorted(stats['operator_counts'].items()):
+            print(f"  {op}: {count} ({count/total_ops*100:.1f}%)")
+
+        print("\nDegree Distribution:")
+        total_eqs = num_samples
+        for degree, count in sorted(stats['degree_counts'].items()):
+            print(f"  Degree {degree}: {count} ({count/total_eqs*100:.1f}%)")
+
+        print("\nTerm Count Distribution:")
+        for terms, count in sorted(stats['term_counts'].items()):
+            print(f"  {terms} terms: {count} ({count/total_eqs*100:.1f}%)")
+
+        print("\nVariable Distribution:")
+        total_vars = sum(stats['variable_counts'].values())
+        for var, count in sorted(stats['variable_counts'].items()):
+            print(f"  {var}: {count} ({count/total_vars*100:.1f}%)")
+
+        print("\nCoefficient Statistics:")
+        print(f"  Range: [{stats['coefficient_stats']['min']:.1f} to {stats['coefficient_stats']['max']:.1f}]")
+        if stats['coefficient_stats']['count'] > 0:
+            avg = stats['coefficient_stats']['total'] / stats['coefficient_stats']['count']
+            print(f"  Average: {avg:.2f}")
+            print(f"  Unique values: {len(stats['coefficient_stats']['unique'])}")
+
+        print("\nSolution Statistics:")
+        print(f"  No real solutions: {stats['solution_stats']['no_solutions']} ({stats['solution_stats']['no_solutions']/total_eqs*100:.1f}%)")
+        print(f"  One solution: {stats['solution_stats']['one_solution']} ({stats['solution_stats']['one_solution']/total_eqs*100:.1f}%)")
+        print(f"  Two solutions: {stats['solution_stats']['two_solutions']} ({stats['solution_stats']['two_solutions']/total_eqs*100:.1f}%)")
+        if stats['solution_stats']['min'] != float('inf'):
+            print(f"  Solution range: [{stats['solution_stats']['min']:.2f} to {stats['solution_stats']['max']:.2f}]")
+
+        print("\nCurriculum Level Distribution:")
+        for attr, levels in sorted(stats['level_distribution'].items()):
+            print(f"\n  {attr}:")
+            for level, count in sorted(levels.items()):
+                print(f"    Level {level}: {count} ({count/total_eqs*100:.1f}%)")
+
+        # Verify statistical properties
+        # 1. Check we see all operators when using operator level 1
+        if any(level == 1 for level in stats['level_distribution']['operators'].keys()):
+            self.assertTrue(all(op in stats['operator_counts'] for op in ["+", "-"]),
+                          "Not all operators were generated")
+
+        # 2. Check degree distribution matches curriculum settings
+        max_possible_degree = max(stats['degree_counts'].keys())
+        self.assertLessEqual(max_possible_degree, 3, "Generated degree exceeds maximum allowed")
+
+        # 3. Check term count constraints
+        min_terms = min(stats['term_counts'].keys())
+        max_terms = max(stats['term_counts'].keys())
+        self.assertGreaterEqual(min_terms, 2, "Generated equations with too few terms")
+        self.assertLessEqual(max_terms, 4, "Generated equations with too many terms")
+
+        # 4. Check coefficient ranges
+        if stats['coefficient_stats']['count'] > 0:
+            self.assertGreater(len(stats['coefficient_stats']['unique']), 3,
+                             "Too few unique coefficients generated")
+            self.assertGreater(stats['coefficient_stats']['min'], 0,
+                             "Generated zero or negative coefficients")
+            self.assertLessEqual(stats['coefficient_stats']['max'], 100,
+                               "Generated coefficients exceed maximum allowed")
+
+        # 5. Check solution distribution
+        total_with_solutions = stats['solution_stats']['one_solution'] + stats['solution_stats']['two_solutions']
+        if total_with_solutions > 0:
+            self.assertGreater(stats['solution_stats']['one_solution'], 0,
+                             "No equations with exactly one solution generated")
+            self.assertGreater(stats['solution_stats']['two_solutions'], 0,
+                             "No equations with exactly two solutions generated")
+
+class TestPolynomialEquationsGenerate(unittest.TestCase):
+    """Test the generate function with different curriculum settings"""
+
+    def setUp(self):
+        self.curriculum = PolynomialEquationsCurriculum()
+        self.exercise = PolynomialEquationsExercise()
+        self.rng = random.Random(42)  # Fixed seed for reproducibility
+        self.curriculum.rng = self.rng
+
+    def test_generate_basic_linear(self):
+        """Test generation of basic linear equations"""
+        # Configure curriculum for simple linear equations
+        self.curriculum.set_attr_level("max_degree", 0)  # Linear equations
+        self.curriculum.set_attr_level("num_terms", 0)  # 2 terms
+        self.curriculum.set_attr_level("coefficient_value", 0)  # Small coefficients
+        self.curriculum.set_attr_level("sign", 0)  # No signs
+        self.curriculum.set_attr_level("operators", 0)  # Only +
+
+        problem = self.exercise.generate(self.curriculum)
+
+        # Verify structure
+        self.assertIn("question", problem)
+        self.assertIn("answer", problem)
+        self.assertIn("metadata", problem)
+
+        # Verify terms and operators
+        executed_parts = problem["metadata"]["executed_parts"]
+        self.assertTrue(len(executed_parts["terms"]) >= 2, "Not enough terms generated")
+        self.assertTrue(len(executed_parts["operators"]) >= 1, "No operators generated")
+
+        # Verify operator is addition
+        self.assertEqual(executed_parts["operators"][0], "+")
+
+        # Verify terms have correct degree
+        for term in executed_parts["terms"]:
+            self.assertNotIn("**", term, "Term should not have exponent > 1")
+
+    def test_generate_with_signs(self):
+        """Test generation with positive/negative signs"""
+        self.curriculum.set_attr_level("operators", 0)  # Only +
+        self.curriculum.set_attr_level("num_terms", 0)  # 2 terms
+        self.curriculum.set_attr_level("sign", 1)  # Allow -
+        self.curriculum.set_attr_level("max_degree", 0)  # Linear equations
+
+        num_samples = 50
+        terms_seen = []
+
+        for _ in range(num_samples):
+            problem = self.exercise.generate(self.curriculum)
+            executed_parts = problem["metadata"]["executed_parts"]
+            terms_seen.extend(executed_parts["terms"])
+
+        # Check we see both positive and negative terms
+        has_negative = any(term.startswith('-') for term in terms_seen)
+        has_positive = any(not term.startswith('-') for term in terms_seen)
+        self.assertTrue(has_positive, "No positive terms generated")
+        self.assertTrue(has_negative, "No negative terms generated")
+
+    def test_term_count_distribution(self):
+        """Test that term counts follow the correct distribution"""
+        self.curriculum.set_attr_level("num_terms", 2)  # 2-4 terms
+        num_samples = 100
+        term_counts = []
+
+        for _ in range(num_samples):
+            problem = self.exercise.generate(self.curriculum)
+            executed_parts = problem["metadata"]["executed_parts"]
+            term_count = len(executed_parts["terms"])
+            term_counts.append(term_count)
+            self.assertTrue(2 <= term_count <= 4, f"Term count {term_count} outside valid range [2,4]")
+
+        # Verify we see different term counts
+        unique_counts = set(term_counts)
+        self.assertTrue(len(unique_counts) > 1, "Only one term count generated")
+
+    def test_operator_distribution(self):
+        """Test distribution of operators"""
+        self.curriculum.set_attr_level("operators", 1)  # +, -
+        num_samples = 100
+        operators_seen = []
+
+        for _ in range(num_samples):
+            problem = self.exercise.generate(self.curriculum)
+            executed_parts = problem["metadata"]["executed_parts"]
+            operators_seen.extend(executed_parts["operators"])
+
+        # Check we see both operators
+        has_plus = "+" in operators_seen
+        has_minus = "-" in operators_seen
+        self.assertTrue(has_plus, "No + operators generated")
+        self.assertTrue(has_minus, "No - operators generated")
+
+    def test_variable_distribution(self):
+        """Test distribution of variable names"""
+        self.curriculum.set_attr_level("var_name", 0)  # x, y, z
+        num_samples = 100
+        variables_seen = set()
+
+        for _ in range(num_samples):
+            problem = self.exercise.generate(self.curriculum)
+            executed_parts = problem["metadata"]["executed_parts"]
+            variables_seen.add(executed_parts["variable"])
+
+        # Check we see multiple variables
+        self.assertTrue(len(variables_seen) > 1, "Only one variable name generated")
+        self.assertTrue(all(var in "xyz" for var in variables_seen), 
+                       f"Invalid variables generated: {variables_seen}")
+
+    def test_coefficient_distribution(self):
+        """Test distribution of coefficient values"""
+        self.curriculum.set_attr_level("coefficient_value", 0)  # 1-10
+        num_samples = 100
+        coefficients = []
+
+        for _ in range(num_samples):
+            problem = self.exercise.generate(self.curriculum)
+            executed_parts = problem["metadata"]["executed_parts"]
+
+            for term in executed_parts["terms"]:
+                if '*' in term:
+                    coeff = term.split('*')[0]
+                    if coeff and coeff not in ['-', '+']:
+                        coefficients.append(abs(float(coeff)))
+
+        # Check coefficient range
+        self.assertTrue(all(1 <= c <= 10 for c in coefficients), 
+                       "Coefficients outside valid range [1,10]")
+        # Check we see different values
+        unique_coeffs = set(coefficients)
+        self.assertTrue(len(unique_coeffs) > 3, 
+                       f"Too few unique coefficients: {unique_coeffs}")
+
+    def test_error_handling(self):
+        """Test error handling in equation generation"""
+        # Test with invalid attribute level
+        with self.assertRaises(ValueError):
+            self.curriculum.set_attr_level("max_degree", 999)
+
+        # Test with invalid attribute name
+        with self.assertRaises(KeyError):
+            self.curriculum.set_attr_level("invalid_attr", 0)
+
+if __name__ == '__main__':
+    unittest.main()