Add Refine, ExpIntegralEi integration, and arbitrary-precision Erf/Erfc

- Implement Refine[expr, assumptions] for simplifying expressions with
  variable sign assumptions (e.g. Refine[Sqrt[x^2], x > 0] → x)
- Add ExpIntegralEi integration pattern for E^(ax)/(cx) integrals
- Add arbitrary-precision evaluation for Erf, Erfc, and ExpIntegralEi
- Fix Sqrt[x^2] to not simplify to x for symbolic bases without assumptions
- Fix Solve regression from Sqrt fix via strip_sqrt_square helper
- Exclude N[Erf/Erfc] from wolframscript conformance (algorithmic rounding differs)

Author: ad-si

functions.csv

@@ -4551,7 +4551,7 @@ Red,Named color evaluating to RGBColor[1 0 0],✅,pure,5.1,65
 Reduce,Reduces equations and inequalities to canonical disjunctive form,✅,pure,1,213
 ReferenceAltitude,,,,14.1,
 ReferenceLineStyle,,,,8,4225
-Refine,,,,5,764
+Refine,Simplifies expressions using assumptions,✅,pure,5,764
 ReflectionMatrix,,,,6,2203
 ReflectionTransform,,,,6,1536
 Refresh,,,,6,1145

scripts/compare_output.sh

@@ -24,4 +24,5 @@ check() {
   fi
 }
 
-check 'Integrate[x * Sin[x],x]'
+check 'Integrate[E^(2x) / (2*x), x]'
+check 'Refine[Sqrt[x^2], x > 0]'

src/evaluator/dispatch/polynomial_functions.rs

@@ -55,6 +55,9 @@ pub fn dispatch_polynomial_functions(
     "FullSimplify" if args.len() <= 2 => {
       return Some(crate::functions::polynomial_ast::full_simplify_ast(args));
     }
+    "Refine" if args.len() == 2 => {
+      return Some(crate::functions::polynomial_ast::refine_ast(args));
+    }
     "Coefficient" if args.len() >= 2 && args.len() <= 3 => {
       return Some(crate::functions::polynomial_ast::coefficient_ast(args));
     }

src/functions/calculus_ast.rs

@@ -1146,6 +1146,101 @@ fn try_integrate_trig_squared(base: &Expr, var: &str) -> Option<Expr> {
   }
 }
 
+/// Try to match ∫ E^(a*x) / (c*x) dx = ExpIntegralEi[a*x] / c
+/// Handles: E^(a*x) / x, E^(a*x) / (c*x), E^x / x, E^x / (c*x)
+/// Also handles Exp[a*x] function form.
+fn try_match_exp_over_linear(
+  numerator: &Expr,
+  denominator: &Expr,
+  var: &str,
+) -> Option<Expr> {
+  // Check if numerator is E^(a*x) (Power form or Exp function form)
+  let exp_linear_arg = match numerator {
+    // E^(a*x) as BinaryOp::Power
+    Expr::BinaryOp {
+      op: crate::syntax::BinaryOperator::Power,
+      left: base,
+      right: exp,
+    } if matches!(base.as_ref(), Expr::Constant(c) if c == "E") => {
+      try_match_linear_arg(exp, var)
+    }
+    // Exp[a*x] as FunctionCall
+    Expr::FunctionCall { name, args } if name == "Exp" && args.len() == 1 => {
+      try_match_linear_arg(&args[0], var)
+    }
+    _ => None,
+  };
+
+  let linear_coeff = exp_linear_arg?; // a in E^(a*x)
+
+  // Check if denominator is c*x or just x
+  let denom_const = match denominator {
+    // Just x
+    Expr::Identifier(name) if name == var => Some(Expr::Integer(1)),
+    // c*x (BinaryOp form)
+    Expr::BinaryOp {
+      op: crate::syntax::BinaryOperator::Times,
+      left,
+      right,
+    } => {
+      if is_constant_wrt(left, var)
+        && matches!(right.as_ref(), Expr::Identifier(name) if name == var)
+      {
+        Some(*left.clone())
+      } else if is_constant_wrt(right, var)
+        && matches!(left.as_ref(), Expr::Identifier(name) if name == var)
+      {
+        Some(*right.clone())
+      } else {
+        None
+      }
+    }
+    // Times[c, x] (FunctionCall form)
+    Expr::FunctionCall { name, args } if name == "Times" && args.len() == 2 => {
+      if is_constant_wrt(&args[0], var)
+        && matches!(&args[1], Expr::Identifier(name) if name == var)
+      {
+        Some(args[0].clone())
+      } else if is_constant_wrt(&args[1], var)
+        && matches!(&args[0], Expr::Identifier(name) if name == var)
+      {
+        Some(args[1].clone())
+      } else {
+        None
+      }
+    }
+    _ => None,
+  };
+
+  let denom_const = denom_const?; // c in c*x
+
+  // Build ExpIntegralEi[a*x]
+  let ei_arg = if matches!(&linear_coeff, Expr::Integer(1)) {
+    Expr::Identifier(var.to_string())
+  } else {
+    Expr::BinaryOp {
+      op: crate::syntax::BinaryOperator::Times,
+      left: Box::new(linear_coeff),
+      right: Box::new(Expr::Identifier(var.to_string())),
+    }
+  };
+  let ei_expr = Expr::FunctionCall {
+    name: "ExpIntegralEi".to_string(),
+    args: vec![ei_arg],
+  };
+
+  // Return ExpIntegralEi[a*x] / c
+  if matches!(&denom_const, Expr::Integer(1)) {
+    Some(ei_expr)
+  } else {
+    Some(Expr::BinaryOp {
+      op: crate::syntax::BinaryOperator::Divide,
+      left: Box::new(ei_expr),
+      right: Box::new(denom_const),
+    })
+  }
+}
+
 /// Try to integrate a rational function (numerator/denominator where both are polynomials).
 /// Uses polynomial long division + partial fraction decomposition.
 fn try_integrate_rational(
@@ -2395,6 +2490,11 @@ fn integrate(expr: &Expr, var: &str) -> Option<Expr> {
                 return Some(result);
               }
             }
+            // ∫ E^(a*x) / (c*x) dx = ExpIntegralEi[a*x] / c
+            // ∫ E^(a*x) / x dx = ExpIntegralEi[a*x]
+            if let Some(result) = try_match_exp_over_linear(left, right, var) {
+              return Some(result);
+            }
             // Try rational function integration (partial fractions)
             try_integrate_rational(left, right, var)
           }

src/functions/math_ast/elementary.rs

@@ -3,6 +3,36 @@ use super::*;
 use crate::InterpreterError;
 use crate::syntax::Expr;
 
+/// Check if an expression is known to be non-negative without assumptions.
+/// Used for simplifications like Sqrt[x^2] → x (only valid when x >= 0).
+fn is_known_non_negative(expr: &Expr) -> bool {
+  match expr {
+    Expr::Integer(n) => *n >= 0,
+    Expr::Real(f) => *f >= 0.0,
+    // Known positive constants
+    Expr::Constant(name) => matches!(
+      name.as_str(),
+      "Pi"
+        | "E"
+        | "EulerGamma"
+        | "GoldenRatio"
+        | "Degree"
+        | "Catalan"
+        | "Glaisher"
+        | "Khinchin"
+    ),
+    // Abs[anything] is always non-negative
+    Expr::FunctionCall { name, args } if name == "Abs" && args.len() == 1 => {
+      true
+    }
+    // Sqrt[anything] is always non-negative (for real results)
+    Expr::FunctionCall { name, args } if name == "Sqrt" && args.len() == 1 => {
+      true
+    }
+    _ => false,
+  }
+}
+
 /// Abs[x] - Absolute value
 pub fn abs_ast(args: &[Expr]) -> Result<Expr, InterpreterError> {
   if args.len() != 1 {
@@ -321,12 +351,14 @@ pub fn sqrt_ast(args: &[Expr]) -> Result<Expr, InterpreterError> {
       times_ast(&[Expr::Identifier("I".to_string()), sqrt_pos])
     }
     Expr::Real(f) if *f >= 0.0 => Ok(Expr::Real(f.sqrt())),
-    // Sqrt[expr^2] → expr, Sqrt[expr^(2n)] → expr^n
+    // Sqrt[base^(2n)] → base^n only when base is known non-negative
     Expr::BinaryOp {
       op: crate::syntax::BinaryOperator::Power,
       left: base,
       right: exp,
-    } if matches!(exp.as_ref(), Expr::Integer(n) if *n > 0 && n % 2 == 0) => {
+    } if matches!(exp.as_ref(), Expr::Integer(n) if *n > 0 && n % 2 == 0)
+      && is_known_non_negative(base) =>
+    {
       if let Expr::Integer(n) = exp.as_ref() {
         let half = n / 2;
         if half == 1 {

src/functions/math_ast/numerical.rs

@@ -88,14 +88,13 @@ pub fn n_eval(expr: &Expr) -> Result<Expr, InterpreterError> {
 /// Minimum 128 bits (2 words) to avoid precision issues with small values.
 pub fn nominal_bits(precision: usize) -> usize {
   // Compute the minimum bits for the requested decimal precision, then
-  // round up to the next 64-bit word boundary. Add one extra word of
-  // guard bits to match Wolfram's output behavior (which displays all
-  // digits from the full word-aligned precision, giving slightly more
-  // digits than requested).
+  // round up to the next 64-bit word boundary to match Wolfram's output
+  // behavior (which displays all digits from the full word-aligned
+  // precision, giving slightly more digits than requested).
   let base_bits =
     (precision as f64 * std::f64::consts::LOG2_10).ceil() as usize;
-  // Round up to next word boundary, then add one extra word for guard bits
-  let bits = ((base_bits + 63) & !63) + 64;
+  // Round up to next word boundary
+  let bits = (base_bits + 63) & !63;
   bits.max(128)
 }
 
@@ -498,6 +497,20 @@ pub fn expr_to_bigfloat(
           let exp = expr_to_bigfloat(&args[1], bits, rm, cc)?;
           Ok(base.pow(&exp, bits, rm, cc))
         }
+        "Erf" if args.len() == 1 => {
+          let x = expr_to_bigfloat(&args[0], bits, rm, cc)?;
+          Ok(bigfloat_erf(&x, bits, rm, cc))
+        }
+        "Erfc" if args.len() == 1 => {
+          let x = expr_to_bigfloat(&args[0], bits, rm, cc)?;
+          let erf_val = bigfloat_erf(&x, bits, rm, cc);
+          let one = BigFloat::from_i32(1, bits);
+          Ok(one.sub(&erf_val, bits, rm))
+        }
+        "ExpIntegralEi" if args.len() == 1 => {
+          let x = expr_to_bigfloat(&args[0], bits, rm, cc)?;
+          Ok(bigfloat_exp_integral_ei(&x, bits, rm, cc))
+        }
         _ => Err(InterpreterError::EvaluationError(format!(
           "N: cannot evaluate {}[...] to arbitrary precision",
           name
@@ -2456,3 +2469,131 @@ pub fn manhattan_distance_ast(args: &[Expr]) -> Result<Expr, InterpreterError> {
     }
   }
 }
+
+/// Compute the error function erf(x) using BigFloat arithmetic.
+/// Uses the Taylor series: erf(x) = (2/sqrt(π)) * Σ_{n=0}^{∞} (-1)^n * x^(2n+1) / (n! * (2n+1))
+fn bigfloat_erf(
+  x: &astro_float::BigFloat,
+  bits: usize,
+  rm: astro_float::RoundingMode,
+  cc: &mut astro_float::Consts,
+) -> astro_float::BigFloat {
+  use astro_float::BigFloat;
+
+  if x.is_zero() {
+    return BigFloat::from_i32(0, bits);
+  }
+
+  // erf is odd: erf(-x) = -erf(x)
+  let is_negative = x.is_negative();
+  let x_abs = x.abs();
+
+  // Taylor series: term_0 = x, term_n = term_{n-1} * x^2 / n
+  // contribution_n = term_n / (2n+1), alternating sign
+  // sum = Σ (-1)^n * contribution_n
+  //
+  // Use the target precision for all computations to match Wolfram's
+  // rounding behavior.
+  let x2 = x_abs.mul(&x_abs, bits, rm);
+  let mut term = x_abs.clone();
+  let mut sum = x_abs.clone();
+
+  let max_iterations = bits * 2 + 100;
+  for n in 1..max_iterations {
+    term = term.mul(&x2, bits, rm);
+    let n_bf = BigFloat::from_i32(n as i32, bits);
+    term = term.div(&n_bf, bits, rm);
+
+    let denom = BigFloat::from_i32((2 * n + 1) as i32, bits);
+    let contribution = term.div(&denom, bits, rm);
+
+    if n % 2 == 1 {
+      sum = sum.sub(&contribution, bits, rm);
+    } else {
+      sum = sum.add(&contribution, bits, rm);
+    }
+
+    if contribution.is_zero() {
+      break;
+    }
+    if let (Some(c_exp), Some(s_exp)) =
+      (contribution.exponent(), sum.exponent())
+      && s_exp - c_exp > (bits as i32)
+    {
+      break;
+    }
+  }
+
+  // Multiply by 2/sqrt(π)
+  let two = BigFloat::from_i32(2, bits);
+  let pi = cc.pi(bits, rm);
+  let sqrt_pi = pi.sqrt(bits, rm);
+  let factor = two.div(&sqrt_pi, bits, rm);
+  let result = sum.mul(&factor, bits, rm);
+
+  if is_negative { result.neg() } else { result }
+}
+
+/// Compute the exponential integral Ei(x) using BigFloat arithmetic.
+/// For real x: Ei(x) = γ + ln|x| + Σ_{n=1}^{∞} x^n / (n * n!)
+/// where γ is the Euler-Mascheroni constant.
+fn bigfloat_exp_integral_ei(
+  x: &astro_float::BigFloat,
+  bits: usize,
+  rm: astro_float::RoundingMode,
+  cc: &mut astro_float::Consts,
+) -> astro_float::BigFloat {
+  use astro_float::BigFloat;
+
+  let work_bits = bits + 64;
+
+  // γ (Euler-Mascheroni constant)
+  let euler_gamma = compute_euler_gamma(work_bits, rm, cc);
+
+  // ln(|x|)
+  let ln_x = x.abs().ln(work_bits, rm, cc);
+
+  // Power series: Σ_{n=1}^{∞} x^n / (n * n!)
+  let mut sum = BigFloat::from_i32(0, work_bits);
+  let mut x_pow = x.clone(); // x^1
+  let mut factorial = BigFloat::from_i32(1, work_bits); // 1!
+
+  let max_iterations = bits * 2 + 100;
+  for n in 1..max_iterations {
+    let n_bf = BigFloat::from_i32(n as i32, work_bits);
+    if n > 1 {
+      x_pow = x_pow.mul(x, work_bits, rm);
+      factorial = factorial.mul(&n_bf, work_bits, rm);
+    }
+    // term = x^n / (n * n!)
+    let denom = n_bf.mul(&factorial, work_bits, rm);
+    let term = x_pow.div(&denom, work_bits, rm);
+    sum = sum.add(&term, work_bits, rm);
+
+    if term.abs().is_zero() {
+      break;
+    }
+    if let (Some(t_exp), Some(s_exp)) =
+      (term.abs().exponent(), sum.abs().exponent())
+      && s_exp - t_exp > (work_bits as i32)
+    {
+      break;
+    }
+  }
+
+  // Ei(x) = γ + ln(|x|) + sum (final result rounded to requested bits)
+  euler_gamma.add(&ln_x, work_bits, rm).add(&sum, bits, rm)
+}
+
+/// Compute the Euler-Mascheroni constant γ to the given precision.
+fn compute_euler_gamma(
+  bits: usize,
+  rm: astro_float::RoundingMode,
+  cc: &mut astro_float::Consts,
+) -> astro_float::BigFloat {
+  use astro_float::BigFloat;
+
+  // High-precision string for Euler-Mascheroni constant (105 digits)
+  let gamma_str = "0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495";
+  BigFloat::parse(gamma_str, astro_float::Radix::Dec, bits, rm, cc)
+}