use super::{CodegenOptions, Precision};
use crate::base::arena::Arena;
use crate::base::errors::SymplexError;
use crate::base::node::{ExprId, ExprNode};
use num_traits::ToPrimitive;
use rustc_hash::FxHashMap;
use std::collections::BTreeSet;
pub(crate) fn to_c_fn(
arena: &mut Arena,
expr: ExprId,
name: &str,
args: &[&str],
) -> Result<String, SymplexError> {
to_c_fn_with_options(arena, expr, name, args, &CodegenOptions::default())
}
pub(crate) fn to_c_fn_with_options(
arena: &mut Arena,
expr: ExprId,
name: &str,
args: &[&str],
options: &CodegenOptions,
) -> Result<String, SymplexError> {
let ty = c_type(options.precision);
let (bindings, final_expr) = if options.cse {
let r = crate::output::cse::cse(arena, expr);
(r.bindings, r.expr)
} else {
(Vec::new(), expr)
};
let mut cse_slots: FxHashMap<ExprId, usize> = FxHashMap::default();
for (i, (name_id, _)) in bindings.iter().enumerate() {
cse_slots.insert(*name_id, i);
}
let mut body: Vec<String> = Vec::new();
for (i, (_, value)) in bindings.iter().enumerate() {
let code = emit_expr(arena, *value, args, &cse_slots, options)?;
body.push(format!(
" const {ty} t{i} = {};",
strip_outer_parens(&code)
));
}
let result = emit_expr(arena, final_expr, args, &cse_slots, options)?;
body.push(format!(" return {};", strip_outer_parens(&result)));
let params = if args.is_empty() {
"void".to_string()
} else {
args.iter()
.map(|a| format!("{ty} {a}"))
.collect::<Vec<_>>()
.join(", ")
};
let storage = if options.inline { "static inline " } else { "" };
let mut out = String::new();
out.push_str("/* Generated by symplex. */\n#include <math.h>\n");
if options.checked_domain {
out.push_str("#include <assert.h>\n");
}
let body_text = body.join("\n");
if options.emit_runtime {
let used = used_c_helpers(&body_text);
if !used.is_empty() {
out.push('\n');
out.push_str(&helper_source(&used));
}
}
out.push('\n');
out.push_str(&format!("{storage}{ty} {name}({params}) {{\n"));
out.push_str(&body_text);
out.push_str("\n}\n");
Ok(out)
}
pub(crate) fn c_runtime_source() -> String {
let all: BTreeSet<String> = C_HELPERS.iter().map(|h| h.name.to_string()).collect();
let mut out = String::from("/* symplex C runtime. */\n#include <math.h>\n\n");
out.push_str(&helper_source(&all));
out
}
fn c_type(p: Precision) -> &'static str {
match p {
Precision::F64 => "double",
Precision::F32 => "float",
}
}
#[derive(Clone, Copy, PartialEq, Eq)]
enum Kind {
Value,
ValueStripNeg,
Bool,
}
enum Plan {
Unary(&'static str),
Binary(&'static str),
Neg,
Add {
signs: Vec<bool>,
fma: Vec<bool>,
},
ExpM1 {
signs: Vec<bool>,
},
Mul {
n_num: usize,
n_den: usize,
neg: bool,
},
FmaSplit {
n_rest: usize,
},
PowInt(i64),
RealRoot {
odd_numer: bool,
exp: String,
},
PowF,
MinMax(&'static str, usize),
Piecewise(usize),
Rel(&'static str),
AndOr(&'static str, usize),
Not,
BoolToNum,
NumToBool,
RtUnary(&'static str),
RtBinary(&'static str),
RtOrdered(&'static str, i32),
}
struct Frame {
id: ExprId,
kind: Kind,
plan: Option<Plan>,
}
struct CEmitter<'a> {
arena: &'a Arena,
args: &'a [&'a str],
cse_slots: &'a FxHashMap<ExprId, usize>,
options: &'a CodegenOptions,
suffix: &'static str,
work: Vec<Frame>,
values: Vec<String>,
}
fn emit_expr(
arena: &Arena,
root: ExprId,
args: &[&str],
cse_slots: &FxHashMap<ExprId, usize>,
options: &CodegenOptions,
) -> Result<String, SymplexError> {
let mut em = CEmitter {
arena,
args,
cse_slots,
options,
suffix: match options.precision {
Precision::F64 => "",
Precision::F32 => "f",
},
work: Vec::new(),
values: Vec::new(),
};
em.work.push(Frame {
id: root,
kind: Kind::Value,
plan: None,
});
while let Some(frame) = em.work.pop() {
match frame.plan {
None => em.visit(frame.id, frame.kind)?,
Some(plan) => em.combine(plan)?,
}
}
em.values
.pop()
.ok_or_else(|| SymplexError::ComputationFailed {
operation: "to_c_fn",
reason: "internal error: empty value stack".to_string(),
})
}
impl<'a> CEmitter<'a> {
fn lit(&self, v: f64) -> String {
format_c_float(v, self.suffix)
}
fn mf(&self, name: &str) -> String {
format!("{name}{}", self.suffix)
}
fn pop_n(&mut self, n: usize) -> Vec<String> {
let at = self.values.len().saturating_sub(n);
self.values.split_off(at)
}
fn schedule(&mut self, id: ExprId, plan: Plan, children: &[(ExprId, Kind)]) {
self.work.push(Frame {
id,
kind: Kind::Value,
plan: Some(plan),
});
for &(c, k) in children.iter().rev() {
self.work.push(Frame {
id: c,
kind: k,
plan: None,
});
}
}
fn unsupported(&self, what: &str) -> SymplexError {
SymplexError::NotImplemented(format!("cannot generate C code for `{what}`"))
}
fn visit(&mut self, id: ExprId, kind: Kind) -> Result<(), SymplexError> {
let arena = self.arena;
let node = arena.node(id).clone();
let is_bool_node = matches!(
node,
ExprNode::BoolTrue
| ExprNode::BoolFalse
| ExprNode::Gt(_, _)
| ExprNode::Ge(_, _)
| ExprNode::Eq_(_, _)
| ExprNode::Ne(_, _)
| ExprNode::And(_)
| ExprNode::Or(_)
| ExprNode::Not(_)
);
if kind == Kind::Bool && !is_bool_node {
self.schedule(id, Plan::NumToBool, &[(id, Kind::Value)]);
return Ok(());
}
if kind != Kind::Bool && is_bool_node {
self.schedule(id, Plan::BoolToNum, &[(id, Kind::Bool)]);
return Ok(());
}
let v = Kind::Value;
match node {
ExprNode::Num(nid) => {
let r = arena.num(nid);
let n = r.numer().to_f64().unwrap_or(f64::NAN);
let d = r.denom().to_f64().unwrap_or(f64::NAN);
let s = self.lit(n / d);
self.values.push(s);
}
ExprNode::Symbol(sid) => {
if let Some(&slot) = self.cse_slots.get(&id) {
self.values.push(format!("t{slot}"));
} else {
let name = arena.symbol_name(sid);
if self.args.contains(&name) {
self.values.push(name.to_string());
} else {
return Err(SymplexError::FreeSymbol {
name: name.to_string(),
});
}
}
}
ExprNode::Pi => {
let s = self.lit(std::f64::consts::PI);
self.values.push(s);
}
ExprNode::E => {
let s = self.lit(std::f64::consts::E);
self.values.push(s);
}
ExprNode::EulerGamma => {
let s = self.lit(crate::output::codegen::numeric_rt::EULER_GAMMA_F64);
self.values.push(s);
}
ExprNode::Catalan => {
let s = self.lit(crate::output::codegen::numeric_rt::CATALAN_F64);
self.values.push(s);
}
ExprNode::GoldenRatio => {
let s = self.lit(crate::output::codegen::numeric_rt::GOLDEN_RATIO_F64);
self.values.push(s);
}
ExprNode::PhysicalConstant(_, value_id) => self.work.push(Frame {
id: value_id,
kind: v,
plan: None,
}),
ExprNode::Infinity => self.values.push("INFINITY".to_string()),
ExprNode::NegInfinity => self.values.push("(-INFINITY)".to_string()),
ExprNode::NaN | ExprNode::ComplexInfinity => self.values.push("NAN".to_string()),
ExprNode::ImaginaryUnit => return Err(self.unsupported("ImaginaryUnit")),
ExprNode::Re(_) => return Err(self.unsupported("re")),
ExprNode::Im(_) => return Err(self.unsupported("im")),
ExprNode::Conjugate(_) => return Err(self.unsupported("conjugate")),
ExprNode::Arg(_) => return Err(self.unsupported("arg")),
ExprNode::Si(_) => return Err(self.unsupported("Si")),
ExprNode::Ci(_) => return Err(self.unsupported("Ci")),
ExprNode::Ei(_) => return Err(self.unsupported("Ei")),
ExprNode::Li(_) => return Err(self.unsupported("li")),
ExprNode::Zeta(_) => return Err(self.unsupported("zeta")),
ExprNode::Polygamma(_, _) => return Err(self.unsupported("polygamma")),
ExprNode::KroneckerDelta(_, _) => return Err(self.unsupported("KroneckerDelta")),
ExprNode::Add(children) => {
if children.is_empty() {
let s = self.lit(0.0);
self.values.push(s);
return Ok(());
}
let exp_idx = children
.iter()
.position(|&c| matches!(arena.node(c), ExprNode::Exp(_)));
let neg_one_idx = children.iter().position(|&c| is_const(arena, c, -1.0));
if let (Some(ei), Some(ni)) = (exp_idx, neg_one_idx)
&& ei != ni
&& let ExprNode::Exp(inner) = arena.node(children[ei]).clone()
{
let mut kids = vec![(inner, v)];
let mut signs = Vec::new();
for (i, &c) in children.iter().enumerate() {
if i != ei && i != ni {
let (neg, child) = split_sign(arena, c);
signs.push(neg);
kids.push(child);
}
}
self.schedule(id, Plan::ExpM1 { signs }, &kids);
return Ok(());
}
let mut kids: Vec<(ExprId, Kind)> = Vec::new();
let mut signs = Vec::new();
let mut fma = Vec::new();
for &c in children.iter() {
let (neg, child) = split_sign(arena, c);
let is_fma = self.options.use_mul_add
&& !neg
&& matches!(arena.node(c), ExprNode::Mul(f) if f.len() >= 2)
&& !is_neg_one_mul(arena, c);
signs.push(neg);
fma.push(is_fma);
kids.push(child);
}
let any_fma = fma.iter().filter(|&&b| b).count() >= 1 && children.len() >= 2;
if !any_fma {
fma.iter_mut().for_each(|b| *b = false);
}
self.work.push(Frame {
id,
kind: v,
plan: Some(Plan::Add {
signs,
fma: fma.clone(),
}),
});
for (i, &(c, k)) in kids.iter().enumerate().rev() {
if fma[i] {
let n_rest = arena.node(c).children().len() - 1;
let factors: Vec<ExprId> = arena.node(c).children().to_vec();
self.work.push(Frame {
id: c,
kind: v,
plan: Some(Plan::FmaSplit { n_rest }),
});
for &f in factors.iter().rev() {
self.work.push(Frame {
id: f,
kind: v,
plan: None,
});
}
} else {
self.work.push(Frame {
id: c,
kind: k,
plan: None,
});
}
}
}
ExprNode::Mul(children) => {
if children.is_empty() {
let s = self.lit(1.0);
self.values.push(s);
return Ok(());
}
let mut num: Vec<(ExprId, Kind)> = Vec::new();
let mut den: Vec<(ExprId, Kind)> = Vec::new();
let mut neg = false;
for (i, &c) in children.iter().enumerate() {
if i == 0 && is_const(arena, c, -1.0) && children.len() > 1 {
neg = kind != Kind::ValueStripNeg;
continue;
}
if let ExprNode::Pow(b, e) = arena.node(c)
&& is_const(arena, *e, -1.0)
{
den.push((*b, v));
continue;
}
num.push((c, v));
}
let mut kids = num.clone();
kids.extend(den.iter().copied());
self.schedule(
id,
Plan::Mul {
n_num: num.len(),
n_den: den.len(),
neg,
},
&kids,
);
}
ExprNode::Pow(base, exp) => {
if let Some(r) = arena.as_num(exp) {
if r.is_integer()
&& let Some(n) = r.numer().to_i64()
{
self.schedule(id, Plan::PowInt(n), &[(base, v)]);
return Ok(());
}
let one = num_bigint::BigInt::from(1);
let two = num_bigint::BigInt::from(2);
if *r.numer() == one && *r.denom() == two {
self.schedule(id, Plan::Unary("sqrt"), &[(base, v)]);
return Ok(());
}
if *r.numer() == one && *r.denom() == num_bigint::BigInt::from(3) {
self.schedule(id, Plan::Unary("cbrt"), &[(base, v)]);
return Ok(());
}
if r.denom() % &two != num_bigint::BigInt::from(0) {
let odd_numer = r.numer() % &two != num_bigint::BigInt::from(0);
let e = r.numer().to_f64().unwrap_or(f64::NAN)
/ r.denom().to_f64().unwrap_or(f64::NAN);
let exp = self.lit(e);
self.schedule(id, Plan::RealRoot { odd_numer, exp }, &[(base, v)]);
return Ok(());
}
}
self.schedule(id, Plan::PowF, &[(base, v), (exp, v)]);
}
ExprNode::Neg(x) => self.schedule(id, Plan::Neg, &[(x, v)]),
ExprNode::Floor(x) => self.schedule(id, Plan::Unary("floor"), &[(x, v)]),
ExprNode::Ceiling(x) => self.schedule(id, Plan::Unary("ceil"), &[(x, v)]),
ExprNode::Sin(x) => self.schedule(id, Plan::Unary("sin"), &[(x, v)]),
ExprNode::Cos(x) => self.schedule(id, Plan::Unary("cos"), &[(x, v)]),
ExprNode::Tan(x) => self.schedule(id, Plan::Unary("tan"), &[(x, v)]),
ExprNode::Exp(x) => self.schedule(id, Plan::Unary("exp"), &[(x, v)]),
ExprNode::Ln(x) => {
if let ExprNode::Add(ch) = arena.node(x)
&& ch.len() == 2
{
if is_const(arena, ch[0], 1.0) {
self.schedule(id, Plan::Unary("log1p"), &[(ch[1], v)]);
return Ok(());
}
if is_const(arena, ch[1], 1.0) {
self.schedule(id, Plan::Unary("log1p"), &[(ch[0], v)]);
return Ok(());
}
}
self.schedule(id, Plan::Unary("log"), &[(x, v)]);
}
ExprNode::Abs(x) => self.schedule(id, Plan::Unary("fabs"), &[(x, v)]),
ExprNode::Asin(x) => self.schedule(id, Plan::Unary("asin"), &[(x, v)]),
ExprNode::Acos(x) => self.schedule(id, Plan::Unary("acos"), &[(x, v)]),
ExprNode::Atan(x) => self.schedule(id, Plan::Unary("atan"), &[(x, v)]),
ExprNode::Atan2(y, x) => self.schedule(id, Plan::Binary("atan2"), &[(y, v), (x, v)]),
ExprNode::Sinh(x) => self.schedule(id, Plan::Unary("sinh"), &[(x, v)]),
ExprNode::Cosh(x) => self.schedule(id, Plan::Unary("cosh"), &[(x, v)]),
ExprNode::Tanh(x) => self.schedule(id, Plan::Unary("tanh"), &[(x, v)]),
ExprNode::Asinh(x) => self.schedule(id, Plan::Unary("asinh"), &[(x, v)]),
ExprNode::Acosh(x) => self.schedule(id, Plan::Unary("acosh"), &[(x, v)]),
ExprNode::Atanh(x) => self.schedule(id, Plan::Unary("atanh"), &[(x, v)]),
ExprNode::Sign(x) => self.schedule(id, Plan::RtUnary("sign"), &[(x, v)]),
ExprNode::Heaviside(x) => self.schedule(id, Plan::RtUnary("heaviside"), &[(x, v)]),
ExprNode::DiracDelta(_) => {
let s = self.lit(0.0);
self.values.push(s);
}
ExprNode::Min(ch) => {
if ch.is_empty() {
self.values.push("INFINITY".to_string());
return Ok(());
}
let kids: Vec<(ExprId, Kind)> = ch.iter().map(|&c| (c, v)).collect();
self.schedule(id, Plan::MinMax("fmin", kids.len()), &kids);
}
ExprNode::Max(ch) => {
if ch.is_empty() {
self.values.push("(-INFINITY)".to_string());
return Ok(());
}
let kids: Vec<(ExprId, Kind)> = ch.iter().map(|&c| (c, v)).collect();
self.schedule(id, Plan::MinMax("fmax", kids.len()), &kids);
}
ExprNode::Gamma(x) => self.schedule(id, Plan::Unary("tgamma"), &[(x, v)]),
ExprNode::LogGamma(x) => self.schedule(id, Plan::Unary("lgamma"), &[(x, v)]),
ExprNode::Erf(x) => self.schedule(id, Plan::Unary("erf"), &[(x, v)]),
ExprNode::Erfc(x) => self.schedule(id, Plan::Unary("erfc"), &[(x, v)]),
ExprNode::Digamma(x) => self.schedule(id, Plan::RtUnary("digamma"), &[(x, v)]),
ExprNode::LambertW(x) => self.schedule(id, Plan::RtUnary("lambert_w0"), &[(x, v)]),
ExprNode::Factorial(x) => self.schedule(id, Plan::RtUnary("factorial"), &[(x, v)]),
ExprNode::Beta(a, b) => self.schedule(id, Plan::RtBinary("beta"), &[(a, v), (b, v)]),
ExprNode::Binomial(n, k) => {
self.schedule(id, Plan::RtBinary("binomial"), &[(n, v), (k, v)])
}
ExprNode::Piecewise(branches) => {
if branches.is_empty() {
self.values.push("NAN".to_string());
return Ok(());
}
let mut kids = Vec::with_capacity(branches.len() * 2);
for &(val, cond) in branches.iter() {
kids.push((cond, Kind::Bool));
kids.push((val, v));
}
self.schedule(id, Plan::Piecewise(branches.len()), &kids);
}
ExprNode::BoolTrue => self.values.push("1".to_string()),
ExprNode::BoolFalse => self.values.push("0".to_string()),
ExprNode::Gt(a, b) => self.schedule(id, Plan::Rel(">"), &[(a, v), (b, v)]),
ExprNode::Ge(a, b) => self.schedule(id, Plan::Rel(">="), &[(a, v), (b, v)]),
ExprNode::Eq_(a, b) => self.schedule(id, Plan::Rel("=="), &[(a, v), (b, v)]),
ExprNode::Ne(a, b) => self.schedule(id, Plan::Rel("!="), &[(a, v), (b, v)]),
ExprNode::And(ch) => {
if ch.is_empty() {
self.values.push("1".to_string());
return Ok(());
}
let kids: Vec<(ExprId, Kind)> = ch.iter().map(|&c| (c, Kind::Bool)).collect();
self.schedule(id, Plan::AndOr("&&", kids.len()), &kids);
}
ExprNode::Or(ch) => {
if ch.is_empty() {
self.values.push("0".to_string());
return Ok(());
}
let kids: Vec<(ExprId, Kind)> = ch.iter().map(|&c| (c, Kind::Bool)).collect();
self.schedule(id, Plan::AndOr("||", kids.len()), &kids);
}
ExprNode::Not(x) => self.schedule(id, Plan::Not, &[(x, Kind::Bool)]),
ExprNode::Apply(sid, apply_args) => {
let fname = arena.symbol_name(sid).to_string();
self.visit_apply(id, &fname, &apply_args)?;
}
ExprNode::Derivative(_, _) => return Err(self.unsupported("Derivative")),
ExprNode::Integral(_, _) => return Err(self.unsupported("Integral")),
ExprNode::DefiniteIntegral(_, _, _, _) => {
return Err(self.unsupported("DefiniteIntegral"));
}
ExprNode::Sum(_, _, _, _) => return Err(self.unsupported("Sum")),
ExprNode::Product_(_, _, _, _) => return Err(self.unsupported("Product")),
ExprNode::Limit(_, _, _) => return Err(self.unsupported("Limit")),
ExprNode::Series(_, _, _, _) => return Err(self.unsupported("Series")),
ExprNode::LaplaceTransform(_, _, _) => {
return Err(self.unsupported("LaplaceTransform"));
}
ExprNode::InverseLaplaceTransform(_, _, _) => {
return Err(self.unsupported("InverseLaplaceTransform"));
}
ExprNode::Residue(_, _, _) => return Err(self.unsupported("Residue")),
ExprNode::RootOf(_, _) => return Err(self.unsupported("RootOf")),
ExprNode::RootSum(_, _, _) => return Err(self.unsupported("RootSum")),
ExprNode::DSolve(_, _, _) => return Err(self.unsupported("DSolve")),
ExprNode::ConditionSet(_, _) => return Err(self.unsupported("ConditionSet")),
ExprNode::EmptySet
| ExprNode::UniversalSet
| ExprNode::Interval(_, _, _)
| ExprNode::FiniteSet(_)
| ExprNode::SetUnion(_)
| ExprNode::SetIntersection(_)
| ExprNode::SetComplement(_, _) => return Err(self.unsupported("set expression")),
}
Ok(())
}
fn visit_apply(
&mut self,
id: ExprId,
fname: &str,
args: &[ExprId],
) -> Result<(), SymplexError> {
use crate::base::arena as names;
let v = Kind::Value;
let arity_err = |n: usize| {
SymplexError::NotImplemented(format!(
"cannot generate C code for `{fname}` with {} argument(s) (expected {n})",
args.len()
))
};
let ordered = match fname {
n if n == names::FN_BESSELJ => Some("bessel_j"),
n if n == names::FN_BESSELY => Some("bessel_y"),
n if n == names::FN_BESSELI => Some("bessel_i"),
n if n == names::FN_BESSELK => Some("bessel_k"),
n if n == names::FN_LEGENDRE => Some("legendre_p"),
n if n == names::FN_CHEBYSHEV_T => Some("chebyshev_t"),
n if n == names::FN_CHEBYSHEV_U => Some("chebyshev_u"),
n if n == names::FN_HERMITE => Some("hermite_h"),
n if n == names::FN_LAGUERRE => Some("laguerre_l"),
_ => None,
};
if let Some(helper) = ordered {
if args.len() != 2 {
return Err(arity_err(2));
}
let order = const_order(self.arena, args[0], fname)?;
self.schedule(id, Plan::RtOrdered(helper, order), &[(args[1], v)]);
return Ok(());
}
let unary = match fname {
n if n == names::FN_FIBONACCI => Some("fibonacci"),
n if n == names::FN_LUCAS => Some("lucas"),
n if n == names::FN_HARMONIC => Some("harmonic"),
n if n == names::FN_FACTORIAL2 => Some("factorial2"),
_ => None,
};
if let Some(helper) = unary {
if args.len() != 1 {
return Err(arity_err(1));
}
self.schedule(id, Plan::RtUnary(helper), &[(args[0], v)]);
return Ok(());
}
let binary = match fname {
n if n == names::FN_RISING_FACTORIAL => Some("rising_factorial"),
n if n == names::FN_FALLING_FACTORIAL => Some("falling_factorial"),
_ => None,
};
if let Some(helper) = binary {
if args.len() != 2 {
return Err(arity_err(2));
}
self.schedule(id, Plan::RtBinary(helper), &[(args[0], v), (args[1], v)]);
return Ok(());
}
Err(SymplexError::NotImplemented(format!(
"cannot generate C code for user-defined Apply node `{fname}`"
)))
}
fn domain_condition(&self, func: &str, a: &str) -> Option<String> {
let z = self.lit(0.0);
let one = self.lit(1.0);
Some(match func {
"log" => format!("{a} > {z}"),
"sqrt" => format!("{a} >= {z}"),
"asin" | "acos" => format!("{}({a}) <= {one}", self.mf("fabs")),
"acosh" => format!("{a} >= {one}"),
"atanh" => format!("{}({a}) < {one}", self.mf("fabs")),
"log1p" => format!("{a} > -{one}"),
"tgamma" | "lgamma" => format!("!({a} <= {z} && {}({a}) == {a})", self.mf("floor")),
_ => return None,
})
}
fn with_domain_check(&self, call: String, func: &str, a: &str) -> String {
if !self.options.checked_domain {
return call;
}
match self.domain_condition(func, a) {
Some(cond) => format!("(assert({cond}), {call})"),
None => call,
}
}
fn rt_call(&self, helper: &str, order: Option<i32>, args: &[String]) -> String {
let mut parts: Vec<String> = Vec::new();
if let Some(n) = order {
parts.push(format!("{n}"));
}
match self.options.precision {
Precision::F64 => {
parts.extend(args.iter().cloned());
format!("symplex_{helper}({})", parts.join(", "))
}
Precision::F32 => {
parts.extend(args.iter().map(|a| format!("(double){a}")));
format!("(float)symplex_{helper}({})", parts.join(", "))
}
}
}
fn combine(&mut self, plan: Plan) -> Result<(), SymplexError> {
let out = match plan {
Plan::Unary(f) => {
let a = self.pop_n(1).remove(0);
let call = format!("{}({a})", self.mf(f));
self.with_domain_check(call, f, &a)
}
Plan::Binary(f) => {
let ab = self.pop_n(2);
format!("{}({}, {})", self.mf(f), ab[0], ab[1])
}
Plan::Neg => {
let a = self.pop_n(1).remove(0);
format!("(-{a})")
}
Plan::Add { signs, fma } => {
let slots: usize = fma.iter().map(|&b| if b { 2 } else { 1 }).sum();
let vals = self.pop_n(slots);
let mut parts: Vec<String> = Vec::new(); let mut fma_terms: Vec<(String, String)> = Vec::new();
let mut vi = 0;
for (i, &neg) in signs.iter().enumerate() {
if fma[i] {
fma_terms.push((vals[vi].clone(), vals[vi + 1].clone()));
vi += 2;
} else {
parts.push(if neg {
format!("-{}", vals[vi])
} else {
vals[vi].clone()
});
vi += 1;
}
}
let mut acc = if parts.is_empty() {
let (a, b) = fma_terms.remove(0);
format!("({a} * {b})")
} else if parts.len() == 1 {
parts.remove(0)
} else {
join_terms(&parts)
};
for (a, b) in fma_terms {
acc = format!("{}({a}, {b}, {acc})", self.mf("fma"));
}
acc
}
Plan::ExpM1 { signs } => {
let vals = self.pop_n(signs.len() + 1);
let mut parts = vec![format!("{}({})", self.mf("expm1"), vals[0])];
for (i, &neg) in signs.iter().enumerate() {
parts.push(if neg {
format!("-{}", vals[i + 1])
} else {
vals[i + 1].clone()
});
}
if parts.len() == 1 {
parts.remove(0)
} else {
join_terms(&parts)
}
}
Plan::Mul { n_num, n_den, neg } => {
let vals = self.pop_n(n_num + n_den);
let num = if n_num == 0 {
self.lit(1.0)
} else if n_num == 1 {
vals[0].clone()
} else {
format!("({})", vals[..n_num].join(" * "))
};
let body = if n_den == 0 {
num
} else if n_den == 1 {
format!("({num} / {})", vals[n_num])
} else {
format!("({num} / ({}))", vals[n_num..].join(" * "))
};
if neg { format!("(-{body})") } else { body }
}
Plan::FmaSplit { n_rest } => {
let vals = self.pop_n(n_rest + 1);
let first = vals[0].clone();
let rest = if n_rest == 1 {
vals[1].clone()
} else {
format!("({})", vals[1..].join(" * "))
};
self.values.push(first);
rest
}
Plan::PowInt(n) => {
let b = self.pop_n(1).remove(0);
let simple = is_simple_operand(&b);
match n {
0 => self.lit(1.0),
1 => b,
-1 => format!("({} / {b})", self.lit(1.0)),
2..=4 if simple => {
let reps = vec![b.as_str(); n as usize];
format!("({})", reps.join(" * "))
}
-4..=-2 if simple => {
let reps = vec![b.as_str(); (-n) as usize];
format!("({} / ({}))", self.lit(1.0), reps.join(" * "))
}
_ => format!("{}({b}, {})", self.mf("pow"), self.lit(n as f64)),
}
}
Plan::RealRoot { odd_numer, exp } => {
let b = self.pop_n(1).remove(0);
let mag = format!("{}({}({b}), {exp})", self.mf("pow"), self.mf("fabs"));
if odd_numer {
format!("{}({mag}, {b})", self.mf("copysign"))
} else {
mag
}
}
Plan::PowF => {
let be = self.pop_n(2);
format!("{}({}, {})", self.mf("pow"), be[0], be[1])
}
Plan::MinMax(f, n) => {
let vals = self.pop_n(n);
let mut acc = vals[0].clone();
for v in &vals[1..] {
acc = format!("{}({acc}, {v})", self.mf(f));
}
acc
}
Plan::Piecewise(n) => {
let vals = self.pop_n(2 * n);
let mut acc = "NAN".to_string();
for i in (0..n).rev() {
acc = format!("({} ? {} : {acc})", vals[2 * i], vals[2 * i + 1]);
}
acc
}
Plan::Rel(op) => {
let ab = self.pop_n(2);
format!("({} {op} {})", ab[0], ab[1])
}
Plan::AndOr(op, n) => {
let vals = self.pop_n(n);
format!("({})", vals.join(&format!(" {op} ")))
}
Plan::Not => {
let a = self.pop_n(1).remove(0);
format!("(!{a})")
}
Plan::BoolToNum => {
let c = self.pop_n(1).remove(0);
format!("({c} ? {} : {})", self.lit(1.0), self.lit(0.0))
}
Plan::NumToBool => {
let a = self.pop_n(1).remove(0);
format!("({a} != {})", self.lit(0.0))
}
Plan::RtUnary(h) => {
let a = self.pop_n(1).remove(0);
let call = self.rt_call(h, None, std::slice::from_ref(&a));
if self.options.checked_domain && h == "lambert_w0" {
format!("(assert({a} >= -0.36787944117144233), {call})")
} else {
call
}
}
Plan::RtBinary(h) => {
let ab = self.pop_n(2);
self.rt_call(h, None, &ab)
}
Plan::RtOrdered(h, n) => {
let a = self.pop_n(1).remove(0);
let call = self.rt_call(h, Some(n), std::slice::from_ref(&a));
if self.options.checked_domain && (h == "bessel_y" || h == "bessel_k") {
format!("(assert({a} > {}), {call})", self.lit(0.0))
} else {
call
}
}
};
self.values.push(out);
Ok(())
}
}
fn join_terms(parts: &[String]) -> String {
let mut s = String::from("(");
for (i, p) in parts.iter().enumerate() {
if i == 0 {
s.push_str(p);
} else if let Some(rest) = p.strip_prefix('-') {
s.push_str(" - ");
s.push_str(rest);
} else {
s.push_str(" + ");
s.push_str(p);
}
}
s.push(')');
s
}
fn split_sign(arena: &Arena, c: ExprId) -> (bool, (ExprId, Kind)) {
if let ExprNode::Neg(inner) = arena.node(c) {
return (true, (*inner, Kind::Value));
}
if is_neg_one_mul(arena, c) {
return (true, (c, Kind::ValueStripNeg));
}
(false, (c, Kind::Value))
}
fn is_neg_one_mul(arena: &Arena, id: ExprId) -> bool {
if let ExprNode::Mul(children) = arena.node(id)
&& let Some(&first) = children.first()
{
return children.len() > 1 && is_const(arena, first, -1.0);
}
false
}
fn is_const(arena: &Arena, id: ExprId, value: f64) -> bool {
if let Some(r) = arena.as_num(id) {
let n = r.numer().to_f64().unwrap_or(f64::NAN);
let d = r.denom().to_f64().unwrap_or(f64::NAN);
return n / d == value;
}
if value < 0.0
&& let ExprNode::Neg(inner) = arena.node(id)
{
return is_const(arena, *inner, -value);
}
false
}
fn const_order(arena: &Arena, id: ExprId, what: &str) -> Result<i32, SymplexError> {
if let Some(r) = arena.as_num(id)
&& r.is_integer()
&& let Some(n) = r.numer().to_i32()
{
return Ok(n);
}
if let ExprNode::Neg(inner) = arena.node(id)
&& let Some(r) = arena.as_num(*inner)
&& r.is_integer()
&& let Some(n) = r.numer().to_i32()
{
return Ok(-n);
}
Err(SymplexError::NotImplemented(format!(
"{what} requires a constant integer order/degree, got `{}`",
arena.display(id)
)))
}
fn is_simple_operand(code: &str) -> bool {
!code.is_empty()
&& code
.chars()
.all(|c| c.is_alphanumeric() || c == '_' || c == '.' || c == '-' || c == '+')
&& !code.contains('(')
}
fn format_c_float(v: f64, suffix: &str) -> String {
if v.is_nan() {
return "NAN".to_string();
}
if v.is_infinite() {
return if v > 0.0 {
"INFINITY".to_string()
} else {
"(-INFINITY)".to_string()
};
}
let mut s = format!("{v:?}");
if !s.contains('.') && !s.contains('e') && !s.contains('E') {
s.push_str(".0");
}
if v < 0.0 {
format!("({s}{suffix})")
} else {
format!("{s}{suffix}")
}
}
fn strip_outer_parens(s: &str) -> &str {
let bytes = s.as_bytes();
if bytes.first() != Some(&b'(') || bytes.last() != Some(&b')') {
return s;
}
let inner = &s[1..s.len() - 1];
let mut depth: i32 = 0;
for c in inner.chars() {
match c {
'(' => depth += 1,
')' => {
depth -= 1;
if depth < 0 {
return s;
}
}
_ => {}
}
}
if depth == 0 { inner } else { s }
}
struct CHelper {
name: &'static str,
deps: &'static [&'static str],
src: &'static str,
}
fn used_c_helpers(code: &str) -> BTreeSet<String> {
let mut out = BTreeSet::new();
let mut rest = code;
while let Some(pos) = rest.find("symplex_") {
let after = &rest[pos + "symplex_".len()..];
let ident: String = after
.chars()
.take_while(|c| c.is_alphanumeric() || *c == '_')
.collect();
if !ident.is_empty() && after[ident.len()..].starts_with('(') {
out.insert(ident);
}
rest = after;
}
out
}
fn helper_source(names: &BTreeSet<String>) -> String {
let mut needed: BTreeSet<&str> = BTreeSet::new();
let mut stack: Vec<&str> = names.iter().map(String::as_str).collect();
while let Some(n) = stack.pop() {
if let Some(h) = C_HELPERS.iter().find(|h| h.name == n)
&& needed.insert(h.name)
{
stack.extend(h.deps.iter().copied());
}
}
let mut out = String::new();
for h in C_HELPERS {
if needed.contains(h.name) {
out.push_str(h.src.trim_matches('\n'));
out.push_str("\n\n");
}
}
out
}
const C_HELPERS: &[CHelper] = &[
CHelper {
name: "util",
deps: &[],
src: r#"
static inline int symplex_is_int(double x) { return isfinite(x) && x == floor(x); }
static inline int symplex_is_even(double n) { double h = n * 0.5; return h == floor(h); }
static inline double symplex_sin_pi(double x) {
double n = floor(x + 0.5); double r = x - n; double s = sin(3.141592653589793 * r);
return symplex_is_even(n) ? s : -s;
}
static inline double symplex_cos_pi(double x) {
double n = floor(x + 0.5); double r = x - n; double c = cos(3.141592653589793 * r);
return symplex_is_even(n) ? c : -c;
}
static inline double symplex_gamma_sign(double x) { return (x > 0.0 || symplex_is_even(floor(x))) ? 1.0 : -1.0; }
static inline int symplex_is_gamma_pole(double x) { return x <= 0.0 && symplex_is_int(x); }
"#,
},
CHelper {
name: "sign",
deps: &[],
src: r#"
static inline double symplex_sign(double x) { return x > 0.0 ? 1.0 : (x < 0.0 ? -1.0 : 0.0); }
"#,
},
CHelper {
name: "heaviside",
deps: &[],
src: r#"
static inline double symplex_heaviside(double x) { return x > 0.0 ? 1.0 : (x < 0.0 ? 0.0 : 0.5); }
"#,
},
CHelper {
name: "factorial",
deps: &[],
src: r#"
static inline double symplex_factorial(double x) { return tgamma(x + 1.0); }
"#,
},
CHelper {
name: "digamma",
deps: &["util"],
src: r#"
/* Digamma: reflection for x < 0, recurrence to x >= 10, asymptotic series through x^-14. */
static inline double symplex_digamma(double x) {
if (isnan(x) || x == -INFINITY) return NAN;
if (x == INFINITY) return INFINITY;
if (symplex_is_gamma_pole(x)) return NAN;
double result = 0.0;
if (x < 0.0) { result -= 3.141592653589793 * symplex_cos_pi(x) / symplex_sin_pi(x); x = 1.0 - x; }
while (x < 10.0) { result -= 1.0 / x; x += 1.0; }
double inv = 1.0 / x, inv2 = inv * inv;
double series = inv2 * (1.0/12.0 - inv2 * (1.0/120.0 - inv2 * (1.0/252.0 - inv2 * (1.0/240.0
- inv2 * (1.0/132.0 - inv2 * (691.0/32760.0 - inv2 * (1.0/12.0)))))));
return result + log(x) - 0.5 * inv - series;
}
"#,
},
CHelper {
name: "lambert_w0",
deps: &[],
src: r#"
/* Principal branch of the Lambert W function (branch-point series + Halley). */
static inline double symplex_lambert_w0(double x) {
const double neg_inv_e = -0.36787944117144233;
if (isnan(x) || x < neg_inv_e) return NAN;
if (x == INFINITY) return INFINITY;
if (x == 0.0) return x;
if (x == neg_inv_e) return -1.0;
double w;
if (x < -0.25) {
double p = sqrt(2.0 * (2.718281828459045 * x + 1.0));
w = -1.0 + p * (1.0 - p * (1.0/3.0 - p * (11.0/72.0 - p * (43.0/540.0 - p * (769.0/17280.0
- p * (221.0/8505.0 - p * (680863.0/43545600.0 - p * (1963.0/204120.0))))))));
if (p < 0.02) return w;
} else if (x < 2.718281828459045) {
w = log(1.0 + x);
} else {
double l1 = log(x), l2 = log(l1);
w = l1 - l2 + l2 / l1;
}
for (int iter = 0; iter < 40; iter++) {
double ew = exp(w), f = w * ew - x, wp1 = w + 1.0;
double denom = ew * wp1 - (w + 2.0) * f / (2.0 * wp1);
double dw = f / denom;
w -= dw;
if (fabs(dw) <= 2.0 * 2.220446049250313e-16 * (fabs(w) + 1e-300)) break;
}
return w;
}
"#,
},
CHelper {
name: "beta",
deps: &["util"],
src: r#"
/* Beta function B(a,b) = G(a)G(b)/G(a+b). */
static inline double symplex_beta(double a, double b) {
if (isnan(a) || isnan(b)) return NAN;
double s = a + b;
if (symplex_is_gamma_pole(a) || symplex_is_gamma_pole(b)) return NAN;
if (symplex_is_gamma_pole(s)) return 0.0;
if (a > 0.0 && b > 0.0 && s < 171.0) return tgamma(a) * tgamma(b) / tgamma(s);
double sign = symplex_gamma_sign(a) * symplex_gamma_sign(b) * symplex_gamma_sign(s);
return sign * exp(lgamma(a) + lgamma(b) - lgamma(s));
}
"#,
},
CHelper {
name: "binomial",
deps: &["util"],
src: r#"
/* Generalised binomial coefficient. */
static inline double symplex_binomial(double n, double k) {
if (isnan(n) || isnan(k)) return NAN;
if (symplex_is_int(k)) {
if (k < 0.0) return 0.0;
double kk = k;
if (symplex_is_int(n) && n >= 0.0) { if (k > n) return 0.0; if (n - k < k) kk = n - k; }
if (kk <= 2000.0) {
double acc = 1.0;
for (double i = 1.0; i <= kk; i += 1.0) acc = acc * (n - kk + i) / i;
return acc;
}
}
double a = n + 1.0, b = k + 1.0, c = n - k + 1.0;
if (symplex_is_gamma_pole(b) || symplex_is_gamma_pole(c)) return symplex_is_gamma_pole(a) ? NAN : 0.0;
if (symplex_is_gamma_pole(a)) return NAN;
if (a > 0.0 && b > 0.0 && c > 0.0 && a < 171.0) return tgamma(a) / (tgamma(b) * tgamma(c));
double sign = symplex_gamma_sign(a) * symplex_gamma_sign(b) * symplex_gamma_sign(c);
return sign * exp(lgamma(a) - lgamma(b) - lgamma(c));
}
"#,
},
CHelper {
name: "bessel_core",
deps: &[],
src: r#"
/* Power series for J_n(x), n >= 0, x >= 0 (monotone terms for x <= 2 sqrt(n+1)). */
static inline double symplex_bessel_j_series(int n, double x) {
double half = 0.5 * x, term = 1.0;
for (int i = 1; i <= n; i++) term *= half / (double)i;
if (term == 0.0) return 0.0;
double h2 = half * half, sum = term;
for (double k = 1.0; k < 500.0; k += 1.0) {
term *= -h2 / (k * ((double)n + k));
sum += term;
if (fabs(term) <= 2.220446049250313e-16 * fabs(sum)) break;
}
return sum;
}
/* Miller's backward recurrence: J_n plus Y_0, Y_1 via Neumann series. */
static inline void symplex_bessel_miller(int n, double x, double *jn_out, double *y0_out, double *y1_out) {
int nmax = ((double)n > x) ? n : (int)x;
int m = nmax + 20 + (int)sqrt(40.0 * (double)nmax);
if (m % 2 == 1) m += 1;
double two_over_x = 2.0 / x, bjp = 0.0, bj = 1.0, sum = 2.0;
double mk = (double)(m / 2);
double s0 = ((m / 2) % 2 == 0) ? 1.0 / mk : -1.0 / mk, s1 = 0.0, last_odd = 0.0;
double jn = (n == m) ? 1.0 : 0.0, j0 = 0.0, j1 = 0.0;
for (int j = m; j >= 1; j--) {
double bjm = ((double)j * two_over_x) * bj - bjp;
bjp = bj; bj = bjm;
if (fabs(bj) > 1e250) {
const double sc = 1e-250;
bj *= sc; bjp *= sc; sum *= sc; s0 *= sc; s1 *= sc; last_odd *= sc; jn *= sc; j0 *= sc; j1 *= sc;
}
int order = j - 1;
if (order == n) jn = bj;
if (order % 2 == 0) {
if (order == 0) { sum += bj; j0 = bj; }
else {
sum += 2.0 * bj;
double k = (double)(order / 2);
if ((order / 2) % 2 == 0) s0 += bj / k; else s0 -= bj / k;
}
} else {
int kk = (order + 1) / 2; double k = (double)kk;
if (kk % 2 == 0) s1 += (bj - last_odd) / k; else s1 -= (bj - last_odd) / k;
last_odd = bj;
if (order == 1) j1 = bj;
}
}
double inv = 1.0 / sum;
jn *= inv; j0 *= inv; j1 *= inv; s0 *= inv; s1 *= inv;
double lg = log(0.5 * x) + 0.5772156649015329, two_over_pi = 2.0 / 3.141592653589793;
*jn_out = jn;
*y0_out = two_over_pi * (lg * j0 - 2.0 * s0);
*y1_out = two_over_pi * (lg * j1 - j0 / x + s1);
}
/* Hankel asymptotic expansion for x >= 25; returns 0 when it does not converge. */
static inline int symplex_bessel_hankel(int n, double x, double *j_out, double *y_out) {
double mu = 4.0 * (double)n * (double)n, p = 1.0, q = 0.0, term = 1.0, prev = INFINITY;
int converged = 0;
for (int k = 1; k < 80; k++) {
double kf = (double)k;
term *= (mu - (2.0 * kf - 1.0) * (2.0 * kf - 1.0)) / (kf * 8.0 * x);
if (fabs(term) >= prev) break;
prev = fabs(term);
if (k % 2 == 1) { if ((k / 2) % 2 == 0) q += term; else q -= term; }
else { if ((k / 2) % 2 == 0) p += term; else p -= term; }
if (prev < 1e-17) { converged = 1; break; }
}
if (!converged && prev > 1e-15) return 0;
double chi = x - (0.5 * (double)n + 0.25) * 3.141592653589793;
double c = cos(chi), s = sin(chi), pref = sqrt(2.0 / (3.141592653589793 * x));
*j_out = pref * (p * c - q * s);
*y_out = pref * (p * s + q * c);
return 1;
}
/* Power series for Y_0, Y_1 with 0 < x <= 2. */
static inline void symplex_bessel_y01_series(double x, double *y0_out, double *y1_out) {
double half = 0.5 * x, h2 = half * half, lg = log(half);
double j0 = symplex_bessel_j_series(0, x), j1 = symplex_bessel_j_series(1, x);
double term = 1.0, hk = 0.0, s0 = 0.0, term1 = 1.0, s1 = 0.0;
const double eg = 0.5772156649015329, eps = 2.220446049250313e-16;
for (double k = 0.0; k < 200.0;) {
double hk1 = hk + 1.0 / (k + 1.0);
double t1 = term1 * (hk + hk1 - 2.0 * eg);
s1 += t1;
double kn = k + 1.0;
term *= h2 / (kn * kn);
hk = hk1;
double t0 = term * hk;
s0 += t0;
term1 *= -h2 / (kn * (kn + 1.0));
if (fabs(t0) <= eps * fabs(s0) && fabs(term1) <= eps * fabs(s1)) break;
term = -term;
k = kn;
}
double two_over_pi = 2.0 / 3.141592653589793;
*y0_out = two_over_pi * ((lg + eg) * j0 + s0);
*y1_out = -two_over_pi / x + two_over_pi * lg * j1 - half * s1 / 3.141592653589793;
}
"#,
},
CHelper {
name: "bessel_j",
deps: &["bessel_core"],
src: r#"
/* Bessel function of the first kind, integer order. */
static inline double symplex_bessel_j(int n, double x) {
if (isnan(x)) return NAN;
double sign = 1.0;
if (n < 0) { if (n % 2 != 0) sign = -sign; n = -n; }
double ax = x;
if (x < 0.0) { if (n % 2 != 0) sign = -sign; ax = -x; }
if (ax == 0.0) return n == 0 ? 1.0 : 0.0;
if (isinf(ax)) return 0.0;
double v, y;
if (ax <= 2.0 * sqrt((double)n + 1.0)) v = symplex_bessel_j_series(n, ax);
else if (ax < 25.0) { double y0, y1; symplex_bessel_miller(n, ax, &v, &y0, &y1); }
else if (!symplex_bessel_hankel(n, ax, &v, &y)) { double y0, y1; symplex_bessel_miller(n, ax, &v, &y0, &y1); }
return sign * v;
}
"#,
},
CHelper {
name: "bessel_y",
deps: &["bessel_core"],
src: r#"
/* Bessel function of the second kind, integer order, x > 0. */
static inline double symplex_bessel_y(int n, double x) {
if (isnan(x) || x < 0.0) return NAN;
if (x == 0.0) return -INFINITY;
if (isinf(x)) return 0.0;
double sign = 1.0;
if (n < 0) { if (n % 2 != 0) sign = -sign; n = -n; }
double y0, y1, jd;
if (x <= 2.0) symplex_bessel_y01_series(x, &y0, &y1);
else if (x < 25.0) symplex_bessel_miller(0, x, &jd, &y0, &y1);
else {
double ja, jb;
if (!symplex_bessel_hankel(0, x, &ja, &y0) || !symplex_bessel_hankel(1, x, &jb, &y1))
symplex_bessel_miller(0, x, &jd, &y0, &y1);
}
if (n == 0) return sign * y0;
double ym = y0, y = y1;
for (int k = 1; k < n; k++) { double yp = (2.0 * (double)k / x) * y - ym; ym = y; y = yp; }
return sign * y;
}
"#,
},
CHelper {
name: "bessel_i",
deps: &[],
src: r#"
/* Modified Bessel function of the first kind, integer order. */
static inline double symplex_bessel_i(int n, double x) {
if (isnan(x)) return NAN;
if (n < 0) n = -n;
double sign = 1.0, ax = x;
if (x < 0.0) { if (n % 2 != 0) sign = -1.0; ax = -x; }
if (ax == 0.0) return n == 0 ? 1.0 : 0.0;
if (isinf(ax)) return sign * INFINITY;
if (ax > 30.0) {
double mu = 4.0 * (double)n * (double)n, sum = 1.0, term = 1.0, prev = INFINITY;
int ok = 0;
for (int k = 1; k < 80; k++) {
double kf = (double)k;
term *= -(mu - (2.0 * kf - 1.0) * (2.0 * kf - 1.0)) / (kf * 8.0 * ax);
if (fabs(term) >= prev) break;
prev = fabs(term);
sum += term;
if (prev < 1e-17 * fabs(sum)) { ok = 1; break; }
}
if (ok) return sign * exp(ax - 0.5 * log(2.0 * 3.141592653589793 * ax)) * sum;
}
double half = 0.5 * ax, term = 1.0;
for (int i = 1; i <= n; i++) term *= half / (double)i;
if (term == 0.0) return 0.0;
double h2 = half * half, sum = term;
for (double k = 1.0; k < 2000.0; k += 1.0) {
term *= h2 / (k * ((double)n + k));
sum += term;
if (term <= 2.220446049250313e-16 * sum) break;
}
return sign * sum;
}
"#,
},
CHelper {
name: "bessel_k",
deps: &[],
src: r#"
static inline double symplex_bessel_k_trapezoid(double nf, double x, double h) {
double r = nf / x, t_peak = log(r + sqrt(r * r + 1.0));
double sum = 0.5 * exp(-x);
for (int j = 1; j < 100000; j++) {
double t = h * (double)j;
double ch = 0.5 * (exp(t) + exp(-t));
double a = -x * ch;
double f = 0.5 * (exp(a + nf * t) + exp(a - nf * t));
sum += f;
if (t > t_peak + 1.0 && f <= 1e-17 * sum) break;
}
return h * sum;
}
/* Modified Bessel function of the second kind, integer order, x > 0. */
static inline double symplex_bessel_k(int n, double x) {
if (isnan(x) || x < 0.0) return NAN;
if (x == 0.0) return INFINITY;
if (isinf(x)) return 0.0;
if (n < 0) n = -n;
double nf = (double)n;
if (x >= 20.0) {
double mu = 4.0 * nf * nf, sum = 1.0, term = 1.0, prev = INFINITY;
int ok = 0;
for (int k = 1; k < 80; k++) {
double kf = (double)k;
term *= (mu - (2.0 * kf - 1.0) * (2.0 * kf - 1.0)) / (kf * 8.0 * x);
if (fabs(term) >= prev) break;
prev = fabs(term);
sum += term;
if (prev < 1e-17 * fabs(sum)) { ok = 1; break; }
}
if (ok) return sqrt(3.141592653589793 / (2.0 * x)) * exp(-x) * sum;
}
double scale = sqrt(x * x + nf * nf);
double h = 3.141592653589793 * 3.141592653589793 / (scale + 45.0);
double prev = symplex_bessel_k_trapezoid(nf, x, h);
for (int refinements = 0; refinements < 8; refinements++) {
h *= 0.5;
double cur = symplex_bessel_k_trapezoid(nf, x, h);
if (fabs(cur - prev) <= 4.0 * 2.220446049250313e-16 * fabs(cur)) return cur;
prev = cur;
}
return prev;
}
"#,
},
CHelper {
name: "legendre_p",
deps: &[],
src: r#"
static inline double symplex_legendre_p(int n, double x) {
if (n < 0) n = -n - 1;
if (n == 0) return 1.0;
double pm = 1.0, p = x;
for (int k = 1; k < n; k++) { double kf = (double)k; double pn = ((2.0 * kf + 1.0) * x * p - kf * pm) / (kf + 1.0); pm = p; p = pn; }
return p;
}
"#,
},
CHelper {
name: "chebyshev_t",
deps: &[],
src: r#"
static inline double symplex_chebyshev_t(int n, double x) {
if (n < 0) n = -n;
if (n == 0) return 1.0;
double tm = 1.0, t = x;
for (int k = 1; k < n; k++) { double tn = 2.0 * x * t - tm; tm = t; t = tn; }
return t;
}
"#,
},
CHelper {
name: "chebyshev_u",
deps: &[],
src: r#"
static inline double symplex_chebyshev_u(int n, double x) {
if (n == -1) return 0.0;
double sign = 1.0;
if (n < 0) { n = -n - 2; sign = -1.0; }
if (n == 0) return sign;
double um = 1.0, u = 2.0 * x;
for (int k = 1; k < n; k++) { double un = 2.0 * x * u - um; um = u; u = un; }
return sign * u;
}
"#,
},
CHelper {
name: "hermite_h",
deps: &[],
src: r#"
static inline double symplex_hermite_h(int n, double x) {
if (n < 0) return NAN;
if (n == 0) return 1.0;
double hm = 1.0, h = 2.0 * x;
for (int k = 1; k < n; k++) { double hn = 2.0 * x * h - 2.0 * (double)k * hm; hm = h; h = hn; }
return h;
}
"#,
},
CHelper {
name: "laguerre_l",
deps: &[],
src: r#"
static inline double symplex_laguerre_l(int n, double x) {
if (n < 0) return NAN;
if (n == 0) return 1.0;
double lm = 1.0, l = 1.0 - x;
for (int k = 1; k < n; k++) { double kf = (double)k; double ln = ((2.0 * kf + 1.0 - x) * l - kf * lm) / (kf + 1.0); lm = l; l = ln; }
return l;
}
"#,
},
CHelper {
name: "fib_pair",
deps: &[],
src: r#"
/* (F_m, F_{m+1}) by fast doubling in double precision (exact for m <= 78). */
static inline void symplex_fib_pair(unsigned long m, double *f, double *f1) {
double a = 0.0, b = 1.0;
for (int bit = 31; bit >= 0; bit--) {
double d = a * (2.0 * b - a), e = a * a + b * b;
a = d; b = e;
if ((m >> bit) & 1UL) { double t = a + b; a = b; b = t; }
}
*f = a; *f1 = b;
}
"#,
},
CHelper {
name: "fibonacci",
deps: &["fib_pair", "util"],
src: r#"
static inline double symplex_fibonacci(double n) {
if (n == INFINITY) return INFINITY;
if (!symplex_is_int(n)) return NAN;
double m = fabs(n);
if (m > 4000000000.0) return (n < 0.0 && symplex_is_even(m)) ? -INFINITY : INFINITY;
double f, f1; symplex_fib_pair((unsigned long)m, &f, &f1);
return (n < 0.0 && symplex_is_even(m)) ? -f : f;
}
"#,
},
CHelper {
name: "lucas",
deps: &["fib_pair", "util"],
src: r#"
static inline double symplex_lucas(double n) {
if (n == INFINITY) return INFINITY;
if (!symplex_is_int(n)) return NAN;
double m = fabs(n);
if (m > 4000000000.0) return (n < 0.0 && !symplex_is_even(m)) ? -INFINITY : INFINITY;
double f, f1; symplex_fib_pair((unsigned long)m, &f, &f1);
double v = 2.0 * f1 - f;
return (n < 0.0 && !symplex_is_even(m)) ? -v : v;
}
"#,
},
CHelper {
name: "harmonic",
deps: &["digamma", "util"],
src: r#"
static inline double symplex_harmonic(double n) {
if (isnan(n)) return NAN;
if (n == INFINITY) return INFINITY;
if (symplex_is_int(n)) {
if (n < 0.0) return NAN;
if (n <= 100.0) { double sum = 0.0; for (double k = n; k >= 1.0; k -= 1.0) sum += 1.0 / k; return sum; }
}
return symplex_digamma(n + 1.0) + 0.5772156649015329;
}
"#,
},
CHelper {
name: "factorial2",
deps: &["util"],
src: r#"
static inline double symplex_factorial2(double n) {
if (!symplex_is_int(n)) return NAN;
if (n >= -1.0) { double acc = 1.0; for (double k = n; k > 1.0; k -= 2.0) acc *= k; return acc; }
if (symplex_is_even(n)) return NAN;
double acc = 1.0;
for (double k = n + 2.0; k <= -1.0; k += 2.0) acc *= k;
return 1.0 / acc;
}
"#,
},
CHelper {
name: "rising_factorial",
deps: &["util"],
src: r#"
static inline double symplex_rising_factorial(double x, double n) {
if (isnan(x) || isnan(n)) return NAN;
if (symplex_is_int(n) && fabs(n) <= 1000.0) {
double acc = 1.0;
if (n >= 0.0) { for (double i = 0.0; i < n; i += 1.0) acc *= x + i; return acc; }
for (double i = 1.0; i <= -n; i += 1.0) acc *= x - i;
return 1.0 / acc;
}
double top = x + n;
if (symplex_is_gamma_pole(x)) return symplex_is_gamma_pole(top) ? NAN : 0.0;
if (symplex_is_gamma_pole(top)) return NAN;
return symplex_gamma_sign(top) * symplex_gamma_sign(x) * exp(lgamma(top) - lgamma(x));
}
"#,
},
CHelper {
name: "falling_factorial",
deps: &["rising_factorial"],
src: r#"
static inline double symplex_falling_factorial(double x, double n) { return symplex_rising_factorial(x - n + 1.0, n); }
"#,
},
];
#[cfg(test)]
mod tests {
use super::*;
fn gen_c(arena: &mut Arena, e: ExprId, args: &[&str]) -> String {
to_c_fn(arena, e, "f", args).unwrap()
}
#[test]
fn basic_arithmetic_and_cse() {
let mut a = Arena::new();
let x = a.symbol("x");
let s = a.sin(x);
let two = a.int(2);
let s2 = a.pow(s, two);
let e = a.add(&[s2, s]);
let code = gen_c(&mut a, e, &["x"]);
assert!(code.contains("#include <math.h>"));
assert!(code.contains("double f(double x) {"), "{code}");
assert!(code.contains("const double t0 = sin(x);"), "{code}");
assert!(code.contains("return"), "{code}");
assert!(
code.contains("t0 * t0"),
"small integer powers are multiplied out:\n{code}"
);
}
#[test]
fn division_and_powers() {
let mut a = Arena::new();
let x = a.symbol("x");
let y = a.symbol("y");
let m1 = a.int(-1);
let inv_y = a.pow(y, m1);
let e = a.mul(&[x, inv_y]);
let code = gen_c(&mut a, e, &["x", "y"]);
assert!(code.contains("return x / y;"), "{code}");
let seven = a.int(7);
let x7 = a.pow(x, seven);
let code = gen_c(&mut a, x7, &["x"]);
assert!(code.contains("pow(x, 7.0)"), "{code}");
let half = a.rational(1, 2);
let sq = a.pow(x, half);
let code = gen_c(&mut a, sq, &["x"]);
assert!(code.contains("sqrt(x)"), "{code}");
let two_fifths = a.rational(2, 5);
let r = a.pow(x, two_fifths);
let code = gen_c(&mut a, r, &["x"]);
assert!(
code.contains("pow(fabs(x), 0.4)"),
"even numerator → no sign:\n{code}"
);
let three_fifths = a.rational(3, 5);
let r = a.pow(x, three_fifths);
let code = gen_c(&mut a, r, &["x"]);
assert!(
code.contains("copysign(pow(fabs(x), 0.6), x)"),
"odd numerator → sign:\n{code}"
);
}
#[test]
fn fma_and_options() {
let mut a = Arena::new();
let x = a.symbol("x");
let y = a.symbol("y");
let xy = a.mul(&[x, y]);
let s = a.sin(x);
let e = a.add(&[xy, s]);
let code = gen_c(&mut a, e, &["x", "y"]);
assert!(code.contains("fma(x, y, sin(x))"), "{code}");
let opts = CodegenOptions {
use_mul_add: false,
precision: Precision::F32,
inline: true,
checked_domain: true,
..Default::default()
};
let ln = a.ln(x);
let e2 = a.add(&[xy, ln]);
let code = to_c_fn_with_options(&mut a, e2, "g", &["x", "y"], &opts).unwrap();
assert!(
code.contains("static inline float g(float x, float y) {"),
"{code}"
);
assert!(!code.contains("fma"), "{code}");
assert!(code.contains("#include <assert.h>"));
assert!(code.contains("(assert(x > 0.0f), logf(x))"), "{code}");
}
#[test]
fn special_functions_and_helpers() {
let mut a = Arena::new();
let x = a.symbol("x");
let g = a.gamma(x);
let code = gen_c(&mut a, g, &["x"]);
assert!(code.contains("tgamma(x)"), "{code}");
assert!(
!code.contains("static inline"),
"no helpers needed for tgamma:\n{code}"
);
let w = a.lambertw(x);
let code = gen_c(&mut a, w, &["x"]);
assert!(
code.contains("static inline double symplex_lambert_w0(double x)"),
"{code}"
);
assert!(code.contains("return symplex_lambert_w0(x);"), "{code}");
let two = a.int(2);
let k = a.besselk(two, x);
let code = gen_c(&mut a, k, &["x"]);
assert!(code.contains("symplex_bessel_k(2, x)"), "{code}");
assert!(code.contains("symplex_bessel_k_trapezoid"), "{code}");
assert!(
!code.contains("symplex_bessel_miller"),
"K does not need J core:\n{code}"
);
let h = a.harmonic(x);
let code = gen_c(&mut a, h, &["x"]);
let i_util = code.find("symplex_is_int(double").unwrap();
let i_dig = code.find("symplex_digamma(double").unwrap();
let i_harm = code.find("symplex_harmonic(double").unwrap();
assert!(i_util < i_dig && i_dig < i_harm, "{code}");
}
#[test]
fn piecewise_and_bools() {
let mut a = Arena::new();
let x = a.symbol("x");
let zero = a.int(0);
let gt = a.gt(x, zero);
let neg = a.neg(x);
let t = a.bool_true();
let pw = a.piecewise(&[(x, gt), (neg, t)]);
let code = gen_c(&mut a, pw, &["x"]);
assert!(
code.contains("return (x > 0.0) ? x : (1 ? (-x) : NAN);"),
"{code}"
);
let code = gen_c(&mut a, gt, &["x"]);
assert!(code.contains("return (x > 0.0) ? 1.0 : 0.0;"), "{code}");
}
#[test]
fn errors() {
let mut a = Arena::new();
let x = a.symbol("x");
let y = a.symbol("y");
let e = a.add(&[x, y]);
assert!(matches!(
to_c_fn(&mut a, e, "f", &["x"]),
Err(SymplexError::FreeSymbol { .. })
));
let i = a.i_unit;
assert!(matches!(
to_c_fn(&mut a, i, "f", &[]),
Err(SymplexError::NotImplemented(_))
));
let n = a.symbol("n");
let j = a.besselj(n, x);
assert!(matches!(
to_c_fn(&mut a, j, "f", &["x", "n"]),
Err(SymplexError::NotImplemented(_))
));
}
#[test]
fn full_runtime_is_well_formed() {
let rt = c_runtime_source();
for h in C_HELPERS {
assert!(rt.contains(h.src.trim_matches('\n')), "{}", h.name);
}
assert_eq!(rt.matches('{').count(), rt.matches('}').count());
assert_eq!(rt.matches('(').count(), rt.matches(')').count());
}
}