Fix parsing, Apply, CosIntegral, and playground splitting (#60, #61, #63, #64)

- Fix #60: Allow Part extraction in increment/decrement operators
  (++x[[i]], --x[[i]], x[[i]]++, x[[i]]--) by updating grammar rules
  and adding Part evaluation support in the increment/decrement handler.

- Fix #61: Implement CosIntegral and SinIntegral with special values
  (0, Infinity, -Infinity) and numeric evaluation via power series and
  asymptotic expansion.

- Fix #63: Make Apply work with any expression type (BinaryOp, UnaryOp,
  Rule, Association) by adding expr_children() helper that extracts
  children in Wolfram canonical form.

- Fix #64: Fix split_into_statements to recognize /; (Condition) and
  operator continuations at line endings, matching insert_statement_separators
  behavior for consistent playground/file parsing.

Author: ad-si

functions.csv

@@ -1128,7 +1128,7 @@ Cosh,Returns the hyperbolic cosine,✅,pure,1,259
 CoshIntegral,,,,3,1347
 CosineDistance,,,,6,2911
 CosineWindow,,,,9,3862
-CosIntegral,,,,2,1178
+CosIntegral,Cosine integral Ci(z),✅,pure,2,1178
 Cot,Returns the cotangent,✅,pure,1,358
 CotDegrees,,,,14.1,
 Coth,Returns the hyperbolic cotangent,✅,pure,1,592
@@ -4994,7 +4994,7 @@ SingularValueList,,,,5,1892
 SingularValuePlot,,,,8,3755
 Sinh,Returns the hyperbolic sine,✅,pure,1,278
 SinhIntegral,,,,3,1309
-SinIntegral,,,,2,1148
+SinIntegral,Sine integral Si(z),✅,pure,2,1148
 SixJSymbol,,,,2,1217
 Skeleton,Represents omitted elements in output display as <<n>>,✅,none,1,4479
 SkeletonTransform,,,,8,3522

src/evaluator/core_eval.rs

@@ -1050,37 +1050,59 @@ pub fn evaluate_expr_to_expr_inner(
         || name == "PreIncrement"
         || name == "PreDecrement")
         && args.len() == 1
-        && let Expr::Identifier(var_name) = &args[0]
       {
-        let current = ENV.with(|e| e.borrow().get(var_name).cloned());
-        let current_val = match current {
-          Some(StoredValue::ExprVal(e)) => e,
-          Some(StoredValue::Raw(s)) => {
-            crate::syntax::string_to_expr(&s).unwrap_or(Expr::Integer(0))
+        if let Expr::Identifier(var_name) = &args[0] {
+          let current = ENV.with(|e| e.borrow().get(var_name).cloned());
+          let current_val = match current {
+            Some(StoredValue::ExprVal(e)) => e,
+            Some(StoredValue::Raw(s)) => {
+              crate::syntax::string_to_expr(&s).unwrap_or(Expr::Integer(0))
+            }
+            _ => Expr::Integer(0),
+          };
+          let delta = if name == "Increment" || name == "PreIncrement" {
+            Expr::Integer(1)
+          } else {
+            Expr::Integer(-1)
+          };
+          let new_val = evaluate_expr_to_expr(&Expr::BinaryOp {
+            op: crate::syntax::BinaryOperator::Plus,
+            left: Box::new(current_val.clone()),
+            right: Box::new(delta),
+          })?;
+          ENV.with(|e| {
+            e.borrow_mut().insert(
+              var_name.clone(),
+              StoredValue::Raw(crate::syntax::expr_to_string(&new_val)),
+            );
+          });
+          // Post-increment/decrement returns old value; pre returns new value
+          if name == "PreIncrement" || name == "PreDecrement" {
+            return Ok(new_val);
           }
-          _ => Expr::Integer(0),
-        };
-        let delta = if name == "Increment" || name == "PreIncrement" {
-          Expr::Integer(1)
-        } else {
-          Expr::Integer(-1)
-        };
-        let new_val = evaluate_expr_to_expr(&Expr::BinaryOp {
-          op: crate::syntax::BinaryOperator::Plus,
-          left: Box::new(current_val.clone()),
-          right: Box::new(delta),
-        })?;
-        ENV.with(|e| {
-          e.borrow_mut().insert(
-            var_name.clone(),
-            StoredValue::Raw(crate::syntax::expr_to_string(&new_val)),
-          );
-        });
-        // Post-increment/decrement returns old value; pre returns new value
-        if name == "PreIncrement" || name == "PreDecrement" {
-          return Ok(new_val);
+          return Ok(current_val);
+        }
+        // Handle Part expressions: ++x[[i]], --x[[i]], x[[i]]++, x[[i]]--
+        if let Expr::Part { .. } = &args[0] {
+          // Get the current value at the part position
+          let current_val = evaluate_expr_to_expr(&args[0])?;
+          let delta = if name == "Increment" || name == "PreIncrement" {
+            Expr::Integer(1)
+          } else {
+            Expr::Integer(-1)
+          };
+          let new_val = evaluate_expr_to_expr(&Expr::BinaryOp {
+            op: crate::syntax::BinaryOperator::Plus,
+            left: Box::new(current_val.clone()),
+            right: Box::new(delta),
+          })?;
+          // Use Set to assign the new value back to the part
+          crate::evaluator::assignment::set_ast(&args[0], &new_val)?;
+          if name == "PreIncrement" || name == "PreDecrement" {
+            return Ok(new_val);
+          }
+          return Ok(current_val);
         }
-        return Ok(current_val);
       }
       // Special handling for Unset - x =. (removes definition)
       if name == "Unset" && args.len() == 1 {

src/evaluator/dispatch/math_functions.rs

@@ -409,6 +409,12 @@ pub fn dispatch_math_functions(
     "ExpIntegralE" if args.len() == 2 => {
       return Some(crate::functions::math_ast::exp_integral_e_ast(args));
     }
+    "CosIntegral" if args.len() == 1 => {
+      return Some(crate::functions::math_ast::cos_integral_ast(args));
+    }
+    "SinIntegral" if args.len() == 1 => {
+      return Some(crate::functions::math_ast::sin_integral_ast(args));
+    }
     "BesselI" if args.len() == 2 => {
       return Some(crate::functions::math_ast::bessel_i_ast(args));
     }

src/functions/list_helpers_ast/functional.rs

@@ -298,12 +298,69 @@ pub fn nest_while_list_ast(
   Ok(Expr::List(results))
 }
 
+/// Extract the children of any expression in Wolfram canonical form.
+/// Returns None for atomic expressions (Integer, Real, String, Symbol).
+fn expr_children(expr: &Expr) -> Option<Vec<Expr>> {
+  match expr {
+    Expr::List(items) => Some(items.clone()),
+    Expr::FunctionCall { args, .. } => Some(args.clone()),
+    Expr::BinaryOp { op, left, right } => {
+      use crate::syntax::BinaryOperator;
+      match op {
+        // a - b → Plus[a, Times[-1, b]]
+        BinaryOperator::Minus => Some(vec![
+          left.as_ref().clone(),
+          Expr::FunctionCall {
+            name: "Times".to_string(),
+            args: vec![Expr::Integer(-1), right.as_ref().clone()],
+          },
+        ]),
+        // 1/b → Power[b, -1]; a/b → Times[a, Power[b, -1]]
+        BinaryOperator::Divide => {
+          if matches!(left.as_ref(), Expr::Integer(1)) {
+            Some(vec![right.as_ref().clone(), Expr::Integer(-1)])
+          } else {
+            Some(vec![
+              left.as_ref().clone(),
+              Expr::FunctionCall {
+                name: "Power".to_string(),
+                args: vec![right.as_ref().clone(), Expr::Integer(-1)],
+              },
+            ])
+          }
+        }
+        // Plus, Times, Power, And, Or, etc.
+        _ => Some(vec![left.as_ref().clone(), right.as_ref().clone()]),
+      }
+    }
+    Expr::UnaryOp { operand, .. } => Some(vec![operand.as_ref().clone()]),
+    Expr::Rule {
+      pattern,
+      replacement,
+    }
+    | Expr::RuleDelayed {
+      pattern,
+      replacement,
+    } => Some(vec![pattern.as_ref().clone(), replacement.as_ref().clone()]),
+    Expr::Association(pairs) => Some(
+      pairs
+        .iter()
+        .map(|(k, v)| Expr::Rule {
+          pattern: Box::new(k.clone()),
+          replacement: Box::new(v.clone()),
+        })
+        .collect(),
+    ),
+    // Atomic expressions have no children
+    _ => None,
+  }
+}
+
 /// Apply[f, list] - applies f to the elements of list (f @@ list)
 pub fn apply_ast(func: &Expr, list: &Expr) -> Result<Expr, InterpreterError> {
-  let items = match list {
-    Expr::List(items) => items.clone(),
-    Expr::FunctionCall { args, .. } => args.clone(),
-    _ => {
+  let items = match expr_children(list) {
+    Some(items) => items,
+    None => {
       return Err(InterpreterError::EvaluationError(
         "Apply expects a list or expression as second argument".into(),
       ));

src/functions/math_ast/special_functions.rs

@@ -2467,6 +2467,187 @@ pub fn exp_integral_ei_numeric(x: f64) -> f64 {
   }
 }
 
+/// CosIntegral[z] - Cosine integral Ci(z)
+/// Ci(z) = γ + ln(z) + ∫₀ᶻ (cos(t)-1)/t dt
+pub fn cos_integral_ast(args: &[Expr]) -> Result<Expr, InterpreterError> {
+  if args.len() != 1 {
+    return Err(InterpreterError::EvaluationError(
+      "CosIntegral expects exactly 1 argument".into(),
+    ));
+  }
+
+  match &args[0] {
+    // CosIntegral[0] = -Infinity
+    Expr::Integer(0) => Ok(Expr::FunctionCall {
+      name: "Times".to_string(),
+      args: vec![Expr::Integer(-1), Expr::Identifier("Infinity".to_string())],
+    }),
+    // CosIntegral[Infinity] = 0
+    Expr::Identifier(s) if s == "Infinity" => Ok(Expr::Integer(0)),
+    // Numeric evaluation
+    Expr::Real(x) => Ok(Expr::Real(cos_integral_numeric(*x))),
+    // Check for -Infinity
+    other => {
+      if is_neg_infinity(other) {
+        // CosIntegral[-Infinity] = I*Pi
+        return Ok(Expr::FunctionCall {
+          name: "Times".to_string(),
+          args: vec![
+            Expr::Identifier("I".to_string()),
+            Expr::Constant("Pi".to_string()),
+          ],
+        });
+      }
+      // Unevaluated
+      Ok(Expr::FunctionCall {
+        name: "CosIntegral".to_string(),
+        args: args.to_vec(),
+      })
+    }
+  }
+}
+
+/// Compute Ci(z) numerically for real z
+/// Ci(z) = γ + ln|z| + ∫₀ᶻ (cos(t)-1)/t dt
+/// = γ + ln|z| + Σ_{n=1}^∞ (-1)^n z^(2n) / (2n · (2n)!)
+fn cos_integral_numeric(z: f64) -> f64 {
+  let euler_gamma = 0.5772156649015329;
+
+  if z.abs() < 40.0 {
+    // Power series: Ci(z) = γ + ln|z| + Σ (-1)^n z^(2n) / (2n · (2n)!)
+    let mut sum = euler_gamma + z.abs().ln();
+    let z2 = z * z;
+    let mut term = 1.0; // Will build up z^(2n) / (2n)!
+    for n in 1..200 {
+      let n2 = (2 * n) as f64;
+      // term *= z^2 / ((2n-1) * 2n)
+      term *= z2 / ((n2 - 1.0) * n2);
+      let contrib = if n % 2 == 1 { -term } else { term };
+      sum += contrib / n2;
+      if (contrib / n2).abs() < 1e-16 * sum.abs().max(1e-300) {
+        break;
+      }
+    }
+    sum
+  } else {
+    // Asymptotic expansion for large |z|
+    // Ci(z) ≈ sin(z)/z * f(z) - cos(z)/z * g(z)
+    // where f(z) = Σ (-1)^n (2n)! / z^(2n)
+    //       g(z) = Σ (-1)^n (2n+1)! / z^(2n+1)
+    let mut f = 0.0;
+    let mut g = 0.0;
+    let mut f_term = 1.0;
+    let mut g_term = 1.0 / z;
+    for n in 0..100 {
+      f += if n % 2 == 0 { f_term } else { -f_term };
+      g += if n % 2 == 0 { g_term } else { -g_term };
+      let n2 = (2 * n + 1) as f64;
+      let n2p = (2 * n + 2) as f64;
+      f_term *= n2 * n2p / (z * z);
+      g_term *= n2p * (n2p + 1.0) / (z * z);
+      if f_term.abs() < 1e-16 && g_term.abs() < 1e-16 {
+        break;
+      }
+      // Divergent series: stop when terms start growing
+      if n > 0
+        && (f_term.abs() > (n2 * n2p / (z * z) * f_term).abs()
+          || g_term.abs() > 1e10)
+      {
+        break;
+      }
+    }
+    f / z * z.sin() - g * z.cos()
+  }
+}
+
+/// SinIntegral[z] - Sine integral Si(z)
+/// Si(z) = ∫₀ᶻ sin(t)/t dt
+pub fn sin_integral_ast(args: &[Expr]) -> Result<Expr, InterpreterError> {
+  if args.len() != 1 {
+    return Err(InterpreterError::EvaluationError(
+      "SinIntegral expects exactly 1 argument".into(),
+    ));
+  }
+
+  match &args[0] {
+    // SinIntegral[0] = 0
+    Expr::Integer(0) => Ok(Expr::Integer(0)),
+    // SinIntegral[Infinity] = Pi/2
+    Expr::Identifier(s) if s == "Infinity" => Ok(Expr::FunctionCall {
+      name: "Times".to_string(),
+      args: vec![
+        crate::functions::math_ast::make_rational(1, 2),
+        Expr::Constant("Pi".to_string()),
+      ],
+    }),
+    // Numeric evaluation
+    Expr::Real(x) => Ok(Expr::Real(sin_integral_numeric(*x))),
+    // Check for -Infinity
+    other => {
+      if is_neg_infinity(other) {
+        // SinIntegral[-Infinity] = -Pi/2
+        return Ok(Expr::FunctionCall {
+          name: "Times".to_string(),
+          args: vec![
+            crate::functions::math_ast::make_rational(-1, 2),
+            Expr::Constant("Pi".to_string()),
+          ],
+        });
+      }
+      // Unevaluated
+      Ok(Expr::FunctionCall {
+        name: "SinIntegral".to_string(),
+        args: args.to_vec(),
+      })
+    }
+  }
+}
+
+/// Compute Si(z) numerically for real z
+/// Si(z) = Σ_{n=0}^∞ (-1)^n z^(2n+1) / ((2n+1) · (2n+1)!)
+fn sin_integral_numeric(z: f64) -> f64 {
+  if z.abs() < 40.0 {
+    let z2 = z * z;
+    let mut sum = z;
+    let mut term = z;
+    for n in 1..200 {
+      let n2 = (2 * n) as f64;
+      let n2p1 = n2 + 1.0;
+      term *= -z2 / (n2 * n2p1);
+      sum += term / n2p1;
+      if (term / n2p1).abs() < 1e-16 * sum.abs().max(1e-300) {
+        break;
+      }
+    }
+    sum
+  } else {
+    // Asymptotic expansion for large |z|
+    // Si(z) ≈ π/2 - cos(z)/z * f(z) - sin(z)/z * g(z)
+    let sign = if z > 0.0 { 1.0 } else { -1.0 };
+    let az = z.abs();
+    let mut f = 0.0;
+    let mut g = 0.0;
+    let mut f_term = 1.0;
+    let mut g_term = 1.0 / az;
+    for n in 0..100 {
+      f += if n % 2 == 0 { f_term } else { -f_term };
+      g += if n % 2 == 0 { g_term } else { -g_term };
+      let n2 = (2 * n + 1) as f64;
+      let n2p = (2 * n + 2) as f64;
+      let old_f = f_term;
+      f_term *= n2 * n2p / (az * az);
+      g_term *= n2p * (n2p + 1.0) / (az * az);
+      if f_term.abs() < 1e-16 && g_term.abs() < 1e-16 {
+        break;
+      }
+      if f_term.abs() > old_f.abs() {
+        break;
+      }
+    }
+    sign * (std::f64::consts::FRAC_PI_2 - f / az * az.cos() - g * az.sin())
+  }
+}
+
 /// ExpIntegralE[n, z] - Generalized exponential integral E_n(z)
 pub fn exp_integral_e_ast(args: &[Expr]) -> Result<Expr, InterpreterError> {
   if args.len() != 2 {