diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ
index 1f41cdca..9fb6e543 100644
--- a/docs/paper/reductions.typ
+++ b/docs/paper/reductions.typ
@@ -178,6 +178,7 @@
   "PrecedenceConstrainedScheduling": [Precedence Constrained Scheduling],
   "PrimeAttributeName": [Prime Attribute Name],
   "QuadraticAssignment": [Quadratic Assignment],
+  "QuadraticDiophantineEquations": [Quadratic Diophantine Equations],
   "QuantifiedBooleanFormulas": [Quantified Boolean Formulas (QBF)],
   "RectilinearPictureCompression": [Rectilinear Picture Compression],
   "ResourceConstrainedScheduling": [Resource Constrained Scheduling],
@@ -3257,6 +3258,53 @@ A classical NP-complete problem from Garey and Johnson @garey1979[Ch.~3, p.~76],
   ]
 }
 
+#{
+  let x = load-model-example("QuadraticDiophantineEquations")
+  let a = x.instance.a
+  let b = x.instance.b
+  let c = x.instance.c
+  let config = x.optimal_config
+  let xval = config.at(0) + 1
+  let yval = int((c - a * xval * xval) / b)
+  // Enumerate all valid x values for the table
+  let max-x = calc.floor(calc.sqrt(c / a))
+  let rows = range(1, max-x + 1).map(xi => {
+    let rem = c - a * xi * xi
+    let feasible = rem > 0 and calc.rem(rem, b) == 0
+    let yi = if feasible { int(rem / b) } else { none }
+    (xi, rem, feasible, yi)
+  })
+  [
+    #problem-def("QuadraticDiophantineEquations")[
+      Given positive integers $a$, $b$, $c$, determine whether there exist positive integers $x$, $y$ such that $a x^2 + b y = c$.
+    ][
+      Quadratic Diophantine equations of the form $a x^2 + b y = c$ form one of the simplest families of mixed-degree Diophantine problems. The variable $y$ is entirely determined by $x$ via $y = (c - a x^2) slash b$, so the decision problem reduces to checking whether any $x in {1, dots, floor(sqrt(c slash a))}$ yields a positive integer $y$. This can be done in $O(sqrt(c))$ time by trial#footnote[No algorithm improving on brute-force trial of all candidate $x$ values is known; the registered complexity `sqrt(c)` reflects this direct enumeration bound.].
+
+      *Example.* Let $a = #a$, $b = #b$, $c = #c$. Then $x$ ranges over $1, dots, #max-x$:
+
+      #pred-commands(
+        "pred create --example QuadraticDiophantineEquations -o qde.json",
+        "pred solve qde.json --solver brute-force",
+        "pred evaluate qde.json --config " + config.map(str).join(","),
+      )
+
+      #align(center, table(
+        columns: 4,
+        align: center,
+        table.header([$x$], [$c - a x^2$], [Divisible by $b$?], [$y$]),
+        ..rows.map(((xi, rem, ok, yi)) => (
+          [$#xi$],
+          [$#rem$],
+          [#if ok [Yes] else [No]],
+          [#if yi != none [$#yi$] else [$dash$]],
+        )).flatten(),
+      ))
+
+      The instance is satisfiable: $x = #xval, y = #yval$ gives $#a dot #xval^2 + #b dot #yval = #c$.
+    ]
+  ]
+}
+
 #{
   let x = load-model-example("ClosestVectorProblem")
   let basis = x.instance.basis
diff --git a/problemreductions-cli/src/cli.rs b/problemreductions-cli/src/cli.rs
index 1a99531c..35484abd 100644
--- a/problemreductions-cli/src/cli.rs
+++ b/problemreductions-cli/src/cli.rs
@@ -249,6 +249,7 @@ Flags by problem type:
   ProductionPlanning             --num-periods, --demands, --capacities, --setup-costs, --production-costs, --inventory-costs, --cost-bound
   SubsetSum                       --sizes, --target
   ThreePartition                  --sizes, --bound
+  QuadraticDiophantineEquations    --coeff-a, --coeff-b, --coeff-c
   SumOfSquaresPartition           --sizes, --num-groups
   ExpectedRetrievalCost           --probabilities, --num-sectors
   PaintShop                       --sequence
@@ -754,6 +755,15 @@ pub struct CreateArgs {
     /// Target string for StringToStringCorrection (comma-separated symbol indices, e.g., "0,1,3,2")
     #[arg(long)]
     pub target_string: Option<String>,
+    /// Coefficient a for QuadraticDiophantineEquations (coefficient of x²)
+    #[arg(long)]
+    pub coeff_a: Option<u64>,
+    /// Coefficient b for QuadraticDiophantineEquations (coefficient of y)
+    #[arg(long)]
+    pub coeff_b: Option<u64>,
+    /// Constant c for QuadraticDiophantineEquations (right-hand side of ax² + by = c)
+    #[arg(long)]
+    pub coeff_c: Option<u64>,
 }
 
 #[derive(clap::Args)]
diff --git a/problemreductions-cli/src/commands/create.rs b/problemreductions-cli/src/commands/create.rs
index b171f550..fb22dba6 100644
--- a/problemreductions-cli/src/commands/create.rs
+++ b/problemreductions-cli/src/commands/create.rs
@@ -9,7 +9,8 @@ use anyhow::{bail, Context, Result};
 use problemreductions::export::{ModelExample, ProblemRef, ProblemSide, RuleExample};
 use problemreductions::models::algebraic::{
     ClosestVectorProblem, ConsecutiveBlockMinimization, ConsecutiveOnesMatrixAugmentation,
-    ConsecutiveOnesSubmatrix, FeasibleBasisExtension, SparseMatrixCompression, BMF,
+    ConsecutiveOnesSubmatrix, FeasibleBasisExtension, QuadraticDiophantineEquations,
+    SparseMatrixCompression, BMF,
 };
 use problemreductions::models::formula::Quantifier;
 use problemreductions::models::graph::{
@@ -193,7 +194,10 @@ fn all_data_flags_empty(args: &CreateArgs) -> bool {
         && args.conjuncts_spec.is_none()
         && args.deps.is_none()
         && args.query.is_none()
+        && args.coeff_a.is_none()
+        && args.coeff_b.is_none()
         && args.rhs.is_none()
+        && args.coeff_c.is_none()
         && args.required_columns.is_none()
 }
 
@@ -728,6 +732,7 @@ fn example_for(canonical: &str, graph_type: Option<&str>) -> &'static str {
         "IntegerKnapsack" => "--sizes 3,4,5,2,7 --values 4,5,7,3,9 --capacity 15",
         "SubsetSum" => "--sizes 3,7,1,8,2,4 --target 11",
         "ThreePartition" => "--sizes 4,5,6,4,6,5 --bound 15",
+        "QuadraticDiophantineEquations" => "--coeff-a 3 --coeff-b 5 --coeff-c 53",
         "BoyceCoddNormalFormViolation" => {
             "--n 6 --sets \"0,1:2;2:3;3,4:5\" --target 0,1,2,3,4,5"
         }
@@ -2418,6 +2423,32 @@ pub fn create(args: &CreateArgs, out: &OutputConfig) -> Result<()> {
             )
         }
 
+        // QuadraticDiophantineEquations
+        "QuadraticDiophantineEquations" => {
+            let a = args.coeff_a.ok_or_else(|| {
+                anyhow::anyhow!(
+                    "QuadraticDiophantineEquations requires --coeff-a, --coeff-b, and --coeff-c\n\n\
+                     Usage: pred create QuadraticDiophantineEquations --coeff-a 3 --coeff-b 5 --coeff-c 53"
+                )
+            })?;
+            let b = args.coeff_b.ok_or_else(|| {
+                anyhow::anyhow!(
+                    "QuadraticDiophantineEquations requires --coeff-b\n\n\
+                     Usage: pred create QuadraticDiophantineEquations --coeff-a 3 --coeff-b 5 --coeff-c 53"
+                )
+            })?;
+            let c = args.coeff_c.ok_or_else(|| {
+                anyhow::anyhow!(
+                    "QuadraticDiophantineEquations requires --coeff-c\n\n\
+                     Usage: pred create QuadraticDiophantineEquations --coeff-a 3 --coeff-b 5 --coeff-c 53"
+                )
+            })?;
+            (
+                ser(QuadraticDiophantineEquations::try_new(a, b, c).map_err(anyhow::Error::msg)?)?,
+                resolved_variant.clone(),
+            )
+        }
+
         // SumOfSquaresPartition
         "SumOfSquaresPartition" => {
             let sizes_str = args.sizes.as_deref().ok_or_else(|| {
@@ -7703,7 +7734,10 @@ mod tests {
             storage: None,
             quantifiers: None,
             homologous_pairs: None,
+            coeff_a: None,
+            coeff_b: None,
             rhs: None,
+            coeff_c: None,
             required_columns: None,
         }
     }
diff --git a/src/models/algebraic/mod.rs b/src/models/algebraic/mod.rs
index 134662fd..c3694736 100644
--- a/src/models/algebraic/mod.rs
+++ b/src/models/algebraic/mod.rs
@@ -8,6 +8,7 @@
 //! - [`ConsecutiveBlockMinimization`]: Consecutive Block Minimization
 //! - [`ConsecutiveOnesSubmatrix`]: Consecutive Ones Submatrix (column selection with C1P)
 //! - [`QuadraticAssignment`]: Quadratic Assignment Problem
+//! - [`QuadraticDiophantineEquations`]: Decide ax² + by = c in positive integers
 //! - [`SparseMatrixCompression`]: Sparse Matrix Compression by row overlay
 
 pub(crate) mod bmf;
@@ -18,6 +19,7 @@ pub(crate) mod consecutive_ones_submatrix;
 pub(crate) mod feasible_basis_extension;
 pub(crate) mod ilp;
 pub(crate) mod quadratic_assignment;
+pub(crate) mod quadratic_diophantine_equations;
 pub(crate) mod qubo;
 pub(crate) mod sparse_matrix_compression;
 
@@ -29,6 +31,7 @@ pub use consecutive_ones_submatrix::ConsecutiveOnesSubmatrix;
 pub use feasible_basis_extension::FeasibleBasisExtension;
 pub use ilp::{Comparison, LinearConstraint, ObjectiveSense, VariableDomain, ILP};
 pub use quadratic_assignment::QuadraticAssignment;
+pub use quadratic_diophantine_equations::QuadraticDiophantineEquations;
 pub use qubo::QUBO;
 pub use sparse_matrix_compression::SparseMatrixCompression;
 
@@ -44,6 +47,7 @@ pub(crate) fn canonical_model_example_specs() -> Vec<crate::example_db::specs::M
     specs.extend(consecutive_ones_submatrix::canonical_model_example_specs());
     specs.extend(feasible_basis_extension::canonical_model_example_specs());
     specs.extend(quadratic_assignment::canonical_model_example_specs());
+    specs.extend(quadratic_diophantine_equations::canonical_model_example_specs());
     specs.extend(sparse_matrix_compression::canonical_model_example_specs());
     specs
 }
diff --git a/src/models/algebraic/quadratic_diophantine_equations.rs b/src/models/algebraic/quadratic_diophantine_equations.rs
new file mode 100644
index 00000000..039ac21a
--- /dev/null
+++ b/src/models/algebraic/quadratic_diophantine_equations.rs
@@ -0,0 +1,219 @@
+//! Quadratic Diophantine Equations problem implementation.
+//!
+//! Given positive integers a, b, c, determine whether there exist
+//! positive integers x, y such that ax² + by = c.
+
+use crate::registry::{FieldInfo, ProblemSchemaEntry, ProblemSizeFieldEntry};
+use crate::traits::Problem;
+use crate::types::Or;
+use serde::de::Error as _;
+use serde::{Deserialize, Deserializer, Serialize};
+
+inventory::submit! {
+    ProblemSchemaEntry {
+        name: "QuadraticDiophantineEquations",
+        display_name: "Quadratic Diophantine Equations",
+        aliases: &["QDE"],
+        dimensions: &[],
+        module_path: module_path!(),
+        description: "Decide whether ax² + by = c has a solution in positive integers x, y",
+        fields: &[
+            FieldInfo { name: "a", type_name: "u64", description: "Coefficient of x²" },
+            FieldInfo { name: "b", type_name: "u64", description: "Coefficient of y" },
+            FieldInfo { name: "c", type_name: "u64", description: "Right-hand side constant" },
+        ],
+    }
+}
+
+inventory::submit! {
+    ProblemSizeFieldEntry {
+        name: "QuadraticDiophantineEquations",
+        fields: &["c"],
+    }
+}
+
+/// Quadratic Diophantine Equations problem.
+///
+/// Given positive integers a, b, c, determine whether there exist
+/// positive integers x, y such that ax² + by = c.
+///
+/// # Example
+///
+/// ```
+/// use problemreductions::models::algebraic::QuadraticDiophantineEquations;
+/// use problemreductions::{Problem, Solver, BruteForce};
+///
+/// // a=3, b=5, c=53: x=1 gives y=10, x=4 gives y=1
+/// let problem = QuadraticDiophantineEquations::new(3, 5, 53);
+/// let solver = BruteForce::new();
+/// let witness = solver.find_witness(&problem);
+/// assert!(witness.is_some());
+/// ```
+#[derive(Debug, Clone, Serialize)]
+pub struct QuadraticDiophantineEquations {
+    /// Coefficient of x².
+    a: u64,
+    /// Coefficient of y.
+    b: u64,
+    /// Right-hand side constant.
+    c: u64,
+}
+
+impl QuadraticDiophantineEquations {
+    fn validate_inputs(a: u64, b: u64, c: u64) -> Result<(), String> {
+        if a == 0 {
+            return Err("Coefficient a must be positive".to_string());
+        }
+        if b == 0 {
+            return Err("Coefficient b must be positive".to_string());
+        }
+        if c == 0 {
+            return Err("Right-hand side c must be positive".to_string());
+        }
+        Ok(())
+    }
+
+    /// Create a new QuadraticDiophantineEquations instance, returning an error
+    /// instead of panicking when inputs are invalid.
+    pub fn try_new(a: u64, b: u64, c: u64) -> Result<Self, String> {
+        Self::validate_inputs(a, b, c)?;
+        Ok(Self { a, b, c })
+    }
+
+    /// Create a new QuadraticDiophantineEquations instance.
+    ///
+    /// # Panics
+    ///
+    /// Panics if any of a, b, c is zero.
+    pub fn new(a: u64, b: u64, c: u64) -> Self {
+        Self::try_new(a, b, c).unwrap_or_else(|msg| panic!("{msg}"))
+    }
+
+    /// Get the coefficient a (coefficient of x²).
+    pub fn a(&self) -> u64 {
+        self.a
+    }
+
+    /// Get the coefficient b (coefficient of y).
+    pub fn b(&self) -> u64 {
+        self.b
+    }
+
+    /// Get the right-hand side constant c.
+    pub fn c(&self) -> u64 {
+        self.c
+    }
+
+    /// Compute the integer square root of n (floor(sqrt(n))).
+    fn isqrt(n: u64) -> u64 {
+        if n == 0 {
+            return 0;
+        }
+        let mut x = (n as f64).sqrt() as u64;
+        // Correct for floating-point imprecision.
+        while x.saturating_mul(x) > n {
+            x -= 1;
+        }
+        while (x + 1).saturating_mul(x + 1) <= n {
+            x += 1;
+        }
+        x
+    }
+
+    /// Compute the maximum value of x (floor(sqrt(c/a))).
+    /// Returns 0 if c < a (no positive x is possible since x >= 1 requires a*1 <= c).
+    fn max_x(&self) -> u64 {
+        if self.c < self.a {
+            return 0;
+        }
+        Self::isqrt(self.c / self.a)
+    }
+
+    /// Check whether a given x yields a valid positive integer y.
+    ///
+    /// Returns Some(y) if y is a positive integer, None otherwise.
+    pub fn check_x(&self, x: u64) -> Option<u64> {
+        if x == 0 {
+            return None;
+        }
+        let ax2 = self.a.checked_mul(x)?.checked_mul(x)?;
+        if ax2 >= self.c {
+            return None;
+        }
+        let remainder = self.c - ax2;
+        if !remainder.is_multiple_of(self.b) {
+            return None;
+        }
+        let y = remainder / self.b;
+        if y == 0 {
+            return None;
+        }
+        Some(y)
+    }
+}
+
+#[derive(Deserialize)]
+struct QuadraticDiophantineEquationsData {
+    a: u64,
+    b: u64,
+    c: u64,
+}
+
+impl<'de> Deserialize<'de> for QuadraticDiophantineEquations {
+    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
+    where
+        D: Deserializer<'de>,
+    {
+        let data = QuadraticDiophantineEquationsData::deserialize(deserializer)?;
+        Self::try_new(data.a, data.b, data.c).map_err(D::Error::custom)
+    }
+}
+
+impl Problem for QuadraticDiophantineEquations {
+    const NAME: &'static str = "QuadraticDiophantineEquations";
+    type Value = Or;
+
+    fn variant() -> Vec<(&'static str, &'static str)> {
+        crate::variant_params![]
+    }
+
+    fn dims(&self) -> Vec<usize> {
+        let max_x = self.max_x() as usize;
+        if max_x == 0 {
+            // No valid x exists; return empty config space.
+            return vec![0];
+        }
+        // One variable: x ranges over {1, ..., max_x}.
+        // config[0] in {0, ..., max_x - 1} maps to x = config[0] + 1.
+        vec![max_x]
+    }
+
+    fn evaluate(&self, config: &[usize]) -> Or {
+        Or({
+            if config.len() != 1 {
+                return Or(false);
+            }
+            let x = (config[0] as u64) + 1; // 1-indexed
+            self.check_x(x).is_some()
+        })
+    }
+}
+
+crate::declare_variants! {
+    default QuadraticDiophantineEquations => "sqrt(c)",
+}
+
+#[cfg(feature = "example-db")]
+pub(crate) fn canonical_model_example_specs() -> Vec<crate::example_db::specs::ModelExampleSpec> {
+    vec![crate::example_db::specs::ModelExampleSpec {
+        id: "quadratic_diophantine_equations",
+        instance: Box::new(QuadraticDiophantineEquations::new(3, 5, 53)),
+        // x=1 (config[0]=0) gives y=10: 3*1 + 5*10 = 53
+        optimal_config: vec![0],
+        optimal_value: serde_json::json!(true),
+    }]
+}
+
+#[cfg(test)]
+#[path = "../../unit_tests/models/algebraic/quadratic_diophantine_equations.rs"]
+mod tests;
diff --git a/src/models/mod.rs b/src/models/mod.rs
index 031f17c8..f19f71ae 100644
--- a/src/models/mod.rs
+++ b/src/models/mod.rs
@@ -11,8 +11,8 @@ pub mod set;
 // Re-export commonly used types
 pub use algebraic::{
     ClosestVectorProblem, ConsecutiveBlockMinimization, ConsecutiveOnesMatrixAugmentation,
-    ConsecutiveOnesSubmatrix, FeasibleBasisExtension, QuadraticAssignment, SparseMatrixCompression,
-    BMF, ILP, QUBO,
+    ConsecutiveOnesSubmatrix, FeasibleBasisExtension, QuadraticAssignment,
+    QuadraticDiophantineEquations, SparseMatrixCompression, BMF, ILP, QUBO,
 };
 pub use formula::{
     CNFClause, CircuitSAT, KSatisfiability, NAESatisfiability, QuantifiedBooleanFormulas,
diff --git a/src/unit_tests/models/algebraic/quadratic_diophantine_equations.rs b/src/unit_tests/models/algebraic/quadratic_diophantine_equations.rs
new file mode 100644
index 00000000..8f727c75
--- /dev/null
+++ b/src/unit_tests/models/algebraic/quadratic_diophantine_equations.rs
@@ -0,0 +1,173 @@
+use crate::models::algebraic::QuadraticDiophantineEquations;
+use crate::solvers::BruteForce;
+use crate::traits::Problem;
+use crate::types::Or;
+
+fn yes_problem() -> QuadraticDiophantineEquations {
+    // a=3, b=5, c=53: x=1 gives y=10, x=4 gives y=1
+    QuadraticDiophantineEquations::new(3, 5, 53)
+}
+
+fn no_problem() -> QuadraticDiophantineEquations {
+    // a=3, b=5, c=10: x=1 gives 5y=7, not integer
+    QuadraticDiophantineEquations::new(3, 5, 10)
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_basic() {
+    let problem = yes_problem();
+    assert_eq!(problem.a(), 3);
+    assert_eq!(problem.b(), 5);
+    assert_eq!(problem.c(), 53);
+    // max_x = isqrt(53/3) = isqrt(17) = 4
+    assert_eq!(problem.dims(), vec![4]);
+    assert_eq!(problem.num_variables(), 1);
+    assert_eq!(
+        <QuadraticDiophantineEquations as Problem>::NAME,
+        "QuadraticDiophantineEquations"
+    );
+    assert_eq!(
+        <QuadraticDiophantineEquations as Problem>::variant(),
+        vec![]
+    );
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_evaluate_yes() {
+    let problem = yes_problem();
+    // config[0]=0 -> x=1: 3*1 + 5y = 53, y=10
+    assert_eq!(problem.evaluate(&[0]), Or(true));
+    // config[0]=1 -> x=2: 3*4 + 5y = 53, 5y=41, not integer
+    assert_eq!(problem.evaluate(&[1]), Or(false));
+    // config[0]=2 -> x=3: 3*9 + 5y = 53, 5y=26, not integer
+    assert_eq!(problem.evaluate(&[2]), Or(false));
+    // config[0]=3 -> x=4: 3*16 + 5y = 53, 5y=5, y=1
+    assert_eq!(problem.evaluate(&[3]), Or(true));
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_evaluate_no() {
+    let problem = no_problem();
+    // max_x = isqrt(10/3) = isqrt(3) = 1
+    assert_eq!(problem.dims(), vec![1]);
+    // config[0]=0 -> x=1: 3*1 + 5y = 10, 5y=7, not integer
+    assert_eq!(problem.evaluate(&[0]), Or(false));
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_evaluate_invalid_config() {
+    let problem = yes_problem();
+    // Wrong length
+    assert_eq!(problem.evaluate(&[]), Or(false));
+    assert_eq!(problem.evaluate(&[0, 1]), Or(false));
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_solver_finds_witness() {
+    let problem = yes_problem();
+    let solver = BruteForce::new();
+    let witness = solver.find_witness(&problem).unwrap();
+    assert_eq!(problem.evaluate(&witness), Or(true));
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_solver_finds_all_witnesses() {
+    let problem = yes_problem();
+    let solver = BruteForce::new();
+    let all = solver.find_all_witnesses(&problem);
+    // Two solutions: x=1 (config[0]=0) and x=4 (config[0]=3)
+    assert_eq!(all.len(), 2);
+    assert!(all.iter().all(|sol| problem.evaluate(sol) == Or(true)));
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_solver_no_witness() {
+    let problem = no_problem();
+    let solver = BruteForce::new();
+    assert!(solver.find_witness(&problem).is_none());
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_serialization() {
+    let problem = yes_problem();
+    let json = serde_json::to_value(&problem).unwrap();
+    assert_eq!(
+        json,
+        serde_json::json!({
+            "a": 3,
+            "b": 5,
+            "c": 53,
+        })
+    );
+
+    let restored: QuadraticDiophantineEquations = serde_json::from_value(json).unwrap();
+    assert_eq!(restored.a(), problem.a());
+    assert_eq!(restored.b(), problem.b());
+    assert_eq!(restored.c(), problem.c());
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_deserialization_rejects_invalid() {
+    // a=0
+    let result: Result<QuadraticDiophantineEquations, _> =
+        serde_json::from_value(serde_json::json!({"a": 0, "b": 5, "c": 53}));
+    assert!(result.is_err());
+    // b=0
+    let result: Result<QuadraticDiophantineEquations, _> =
+        serde_json::from_value(serde_json::json!({"a": 3, "b": 0, "c": 53}));
+    assert!(result.is_err());
+    // c=0
+    let result: Result<QuadraticDiophantineEquations, _> =
+        serde_json::from_value(serde_json::json!({"a": 3, "b": 5, "c": 0}));
+    assert!(result.is_err());
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_check_x() {
+    let problem = yes_problem();
+    assert_eq!(problem.check_x(1), Some(10)); // 3 + 50 = 53
+    assert_eq!(problem.check_x(2), None); // 12 + 5y = 53, 41/5 not integer
+    assert_eq!(problem.check_x(3), None); // 27 + 5y = 53, 26/5 not integer
+    assert_eq!(problem.check_x(4), Some(1)); // 48 + 5 = 53
+    assert_eq!(problem.check_x(5), None); // 75 > 53
+    assert_eq!(problem.check_x(0), None); // x must be positive
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_edge_case_c_less_than_a() {
+    // c < a: no valid x since a*1^2 = a > c
+    let problem = QuadraticDiophantineEquations::new(10, 1, 5);
+    assert_eq!(problem.dims(), vec![0]);
+}
+
+#[test]
+fn test_quadratic_diophantine_equations_paper_example() {
+    // From issue: a=3, b=5, c=53. x=1: y=10, x=4: y=1.
+    let problem = QuadraticDiophantineEquations::new(3, 5, 53);
+    // Verify the claimed solution x=1 (config[0]=0)
+    assert_eq!(problem.evaluate(&[0]), Or(true));
+    // Verify x=4 (config[0]=3) also works
+    assert_eq!(problem.evaluate(&[3]), Or(true));
+
+    let solver = BruteForce::new();
+    let all = solver.find_all_witnesses(&problem);
+    assert_eq!(all.len(), 2);
+}
+
+#[test]
+#[should_panic(expected = "Coefficient a must be positive")]
+fn test_quadratic_diophantine_equations_panics_on_zero_a() {
+    QuadraticDiophantineEquations::new(0, 5, 53);
+}
+
+#[test]
+#[should_panic(expected = "Coefficient b must be positive")]
+fn test_quadratic_diophantine_equations_panics_on_zero_b() {
+    QuadraticDiophantineEquations::new(3, 0, 53);
+}
+
+#[test]
+#[should_panic(expected = "Right-hand side c must be positive")]
+fn test_quadratic_diophantine_equations_panics_on_zero_c() {
+    QuadraticDiophantineEquations::new(3, 5, 0);
+}