diff --git a/crates/functions/src/engineering.rs b/crates/functions/src/engineering.rs index 86ab7d7c7..ed2572f90 100644 --- a/crates/functions/src/engineering.rs +++ b/crates/functions/src/engineering.rs @@ -197,283 +197,6 @@ fn bitrshift(a: &[Arg], _c: &mut dyn Ctx) -> R { shift(a, false) } -// --------------------------------------------------------------------------------------------- -// CONVERT - -#[derive(Clone, Copy, PartialEq, Eq, Debug)] -enum Dim { - Mass, - Length, - Time, - Pressure, - Force, - Energy, - Power, - Magnetism, - Temperature, - Volume, - Area, - Information, - Speed, -} - -/// How a unit accepts prefixes: none, SI (raised to `power` for areas/volumes) or SI + binary. -#[derive(Clone, Copy, PartialEq, Eq)] -enum Pre { - No, - Si(i32), - Bin, -} - -struct Unit { - names: &'static [&'static str], - dim: Dim, - factor: f64, - pre: Pre, -} - -const fn u(names: &'static [&'static str], dim: Dim, factor: f64, pre: Pre) -> Unit { - Unit { names, dim, factor, pre } -} - -const IN: f64 = 0.0254; -const FT: f64 = 0.3048; -const YD: f64 = 0.9144; -const MI: f64 = 1609.344; -const NMI: f64 = 1852.0; -const LY: f64 = 9_460_730_472_580_800.0; -const PICA_PT: f64 = 0.0254 / 72.0; -const PICA: f64 = 0.0254 / 6.0; -const GAL: f64 = 0.003_785_411_784; - -const UNITS: &[Unit] = &[ - // Mass (gram) - u(&["g"], Dim::Mass, 1.0, Pre::Si(1)), - u(&["sg"], Dim::Mass, 14_593.902_937_206_4, Pre::No), - u(&["lbm"], Dim::Mass, 453.592_37, Pre::No), - u(&["u"], Dim::Mass, 1.660_539_066_6e-24, Pre::Si(1)), - u(&["ozm"], Dim::Mass, 28.349_523_125, Pre::No), - u(&["grain"], Dim::Mass, 0.064_798_91, Pre::No), - u(&["cwt", "shweight"], Dim::Mass, 45_359.237, Pre::No), - u(&["uk_cwt", "lcwt", "hweight"], Dim::Mass, 50_802.345_44, Pre::No), - u(&["stone"], Dim::Mass, 6_350.293_18, Pre::No), - u(&["ton"], Dim::Mass, 907_184.74, Pre::No), - u(&["uk_ton", "LTON", "brton"], Dim::Mass, 1_016_046.908_8, Pre::No), - // Length (metre) - u(&["m"], Dim::Length, 1.0, Pre::Si(1)), - u(&["mi"], Dim::Length, MI, Pre::No), - u(&["Nmi"], Dim::Length, NMI, Pre::No), - u(&["in"], Dim::Length, IN, Pre::No), - u(&["ft"], Dim::Length, FT, Pre::No), - u(&["yd"], Dim::Length, YD, Pre::No), - u(&["ang"], Dim::Length, 1e-10, Pre::Si(1)), - u(&["ell"], Dim::Length, 1.143, Pre::No), - u(&["ly"], Dim::Length, LY, Pre::Si(1)), - u(&["parsec", "pc"], Dim::Length, 30_856_775_814_913_673.0, Pre::Si(1)), - u(&["Picapt", "Pica"], Dim::Length, PICA_PT, Pre::No), - u(&["pica"], Dim::Length, PICA, Pre::No), - u(&["survey_mi"], Dim::Length, 1_609.347_218_694_437, Pre::No), - // Time (second) - u(&["yr"], Dim::Time, 31_557_600.0, Pre::No), - u(&["day", "d"], Dim::Time, 86_400.0, Pre::No), - u(&["hr"], Dim::Time, 3_600.0, Pre::No), - u(&["mn", "min"], Dim::Time, 60.0, Pre::No), - u(&["sec", "s"], Dim::Time, 1.0, Pre::Si(1)), - // Pressure (pascal) - u(&["Pa", "p"], Dim::Pressure, 1.0, Pre::Si(1)), - u(&["atm", "at"], Dim::Pressure, 101_325.0, Pre::Si(1)), - u(&["mmHg"], Dim::Pressure, 133.322_387_415, Pre::Si(1)), - u(&["psi"], Dim::Pressure, 6_894.757_293_168_36, Pre::No), - u(&["Torr"], Dim::Pressure, 101_325.0 / 760.0, Pre::No), - // Force (newton) - u(&["N"], Dim::Force, 1.0, Pre::Si(1)), - u(&["dyn", "dy"], Dim::Force, 1e-5, Pre::Si(1)), - u(&["lbf"], Dim::Force, 4.448_221_615_260_5, Pre::No), - u(&["pond"], Dim::Force, 0.009_806_65, Pre::Si(1)), - // Energy (joule) - u(&["J"], Dim::Energy, 1.0, Pre::Si(1)), - u(&["e"], Dim::Energy, 1e-7, Pre::Si(1)), - u(&["c"], Dim::Energy, 4.184, Pre::Si(1)), - u(&["cal"], Dim::Energy, 4.1868, Pre::Si(1)), - u(&["eV", "ev"], Dim::Energy, 1.602_176_634e-19, Pre::Si(1)), - u(&["HPh", "hh"], Dim::Energy, 2_684_519.537_696_17, Pre::No), - u(&["Wh", "wh"], Dim::Energy, 3_600.0, Pre::Si(1)), - u(&["flb"], Dim::Energy, 1.355_817_948_331_4, Pre::No), - u(&["BTU", "btu"], Dim::Energy, 1_055.055_852_62, Pre::No), - // Power (watt) - u(&["W", "w"], Dim::Power, 1.0, Pre::Si(1)), - u(&["HP", "h"], Dim::Power, 745.699_871_582_270_2, Pre::No), - u(&["PS"], Dim::Power, 735.498_75, Pre::No), - // Magnetism (tesla) - u(&["T"], Dim::Magnetism, 1.0, Pre::Si(1)), - u(&["ga"], Dim::Magnetism, 1e-4, Pre::Si(1)), - // Temperature (handled specially; factor unused) - u(&["C", "cel"], Dim::Temperature, 1.0, Pre::No), - u(&["F", "fah"], Dim::Temperature, 1.0, Pre::No), - u(&["K", "kel"], Dim::Temperature, 1.0, Pre::Si(1)), - u(&["Rank"], Dim::Temperature, 1.0, Pre::No), - u(&["Reau"], Dim::Temperature, 1.0, Pre::No), - // Volume (cubic metre) - u(&["l", "L", "lt"], Dim::Volume, 0.001, Pre::Si(1)), - u(&["tsp"], Dim::Volume, GAL / 768.0, Pre::No), - u(&["tspm"], Dim::Volume, 5e-6, Pre::No), - u(&["tbs"], Dim::Volume, GAL / 256.0, Pre::No), - u(&["oz"], Dim::Volume, GAL / 128.0, Pre::No), - u(&["cup"], Dim::Volume, GAL / 16.0, Pre::No), - u(&["pt", "us_pt"], Dim::Volume, GAL / 8.0, Pre::No), - u(&["uk_pt"], Dim::Volume, 0.000_568_261_25, Pre::No), - u(&["qt"], Dim::Volume, GAL / 4.0, Pre::No), - u(&["uk_qt"], Dim::Volume, 0.001_136_522_5, Pre::No), - u(&["gal"], Dim::Volume, GAL, Pre::No), - u(&["uk_gal"], Dim::Volume, 0.004_546_09, Pre::No), - u(&["m3"], Dim::Volume, 1.0, Pre::Si(3)), - u(&["mi3"], Dim::Volume, MI * MI * MI, Pre::No), - u(&["Nmi3"], Dim::Volume, NMI * NMI * NMI, Pre::No), - u(&["yd3"], Dim::Volume, YD * YD * YD, Pre::No), - u(&["ft3"], Dim::Volume, FT * FT * FT, Pre::No), - u(&["in3"], Dim::Volume, IN * IN * IN, Pre::No), - u(&["ang3"], Dim::Volume, 1e-30, Pre::Si(3)), - u(&["ly3"], Dim::Volume, LY * LY * LY, Pre::No), - u(&["Picapt3", "Pica3"], Dim::Volume, PICA_PT * PICA_PT * PICA_PT, Pre::No), - u(&["barrel"], Dim::Volume, 0.158_987_294_928, Pre::No), - u(&["bushel"], Dim::Volume, 0.035_239_070_166_88, Pre::No), - u(&["MTON"], Dim::Volume, 1.132_673_863_68, Pre::No), - u(&["GRT", "regton"], Dim::Volume, 2.831_684_659_2, Pre::No), - // Area (square metre) - u(&["m2"], Dim::Area, 1.0, Pre::Si(2)), - u(&["ha"], Dim::Area, 10_000.0, Pre::No), - u(&["ar"], Dim::Area, 100.0, Pre::Si(1)), - u(&["ft2"], Dim::Area, FT * FT, Pre::No), - u(&["in2"], Dim::Area, IN * IN, Pre::No), - u(&["yd2"], Dim::Area, YD * YD, Pre::No), - u(&["mi2"], Dim::Area, MI * MI, Pre::No), - u(&["Nmi2"], Dim::Area, NMI * NMI, Pre::No), - u(&["ang2"], Dim::Area, 1e-20, Pre::Si(2)), - u(&["ly2"], Dim::Area, LY * LY, Pre::No), - u(&["Picapt2", "Pica2"], Dim::Area, PICA_PT * PICA_PT, Pre::No), - u(&["uk_acre"], Dim::Area, 4_046.856_422_4, Pre::No), - u(&["us_acre"], Dim::Area, 4_046.872_609_874_252, Pre::No), - u(&["Morgen"], Dim::Area, 2_500.0, Pre::No), - // Information (bit) - u(&["bit"], Dim::Information, 1.0, Pre::Bin), - u(&["byte"], Dim::Information, 8.0, Pre::Bin), - // Speed (metres per second) - u(&["m/s", "m/sec"], Dim::Speed, 1.0, Pre::Si(1)), - u(&["m/h", "m/hr"], Dim::Speed, 1.0 / 3600.0, Pre::Si(1)), - u(&["mph"], Dim::Speed, MI / 3600.0, Pre::No), - u(&["kn"], Dim::Speed, NMI / 3600.0, Pre::No), - u(&["admkn"], Dim::Speed, 6080.0 * FT / 3600.0, Pre::No), -]; - -const SI_PREFIXES: &[(&str, f64)] = &[ - ("da", 1e1), - ("Y", 1e24), - ("Z", 1e21), - ("E", 1e18), - ("P", 1e15), - ("T", 1e12), - ("G", 1e9), - ("M", 1e6), - ("k", 1e3), - ("h", 1e2), - ("e", 1e1), - ("d", 1e-1), - ("c", 1e-2), - ("m", 1e-3), - ("u", 1e-6), - ("n", 1e-9), - ("p", 1e-12), - ("f", 1e-15), - ("a", 1e-18), - ("z", 1e-21), - ("y", 1e-24), -]; - -const BIN_PREFIXES: &[(&str, f64)] = &[ - ("ki", 1024.0), - ("Mi", 1_048_576.0), - ("Gi", 1_073_741_824.0), - ("Ti", 1_099_511_627_776.0), - ("Pi", 1_125_899_906_842_624.0), - ("Ei", 1_152_921_504_606_846_976.0), - ("Zi", 1_180_591_620_717_411_303_424.0), - ("Yi", 1_208_925_819_614_629_174_706_176.0), -]; - -fn find_unit(name: &str) -> Option<&'static Unit> { - UNITS.iter().find(|u| u.names.contains(&name)) -} - -/// Resolves a unit name to (unit, prefix multiplier). Exact names win over prefixed ones. -fn resolve(name: &str) -> Option<(&'static Unit, f64)> { - let name = name.replace("^2", "2").replace("^3", "3"); - if let Some(u) = find_unit(&name) { - return Some((u, 1.0)); - } - for (p, f) in BIN_PREFIXES { - if let Some(rest) = name.strip_prefix(p) - && let Some(u) = find_unit(rest) - && u.pre == Pre::Bin - { - return Some((u, *f)); - } - } - for (p, f) in SI_PREFIXES { - if let Some(rest) = name.strip_prefix(p) - && let Some(u) = find_unit(rest) - { - match u.pre { - Pre::Si(pow) => return Some((u, f.powi(pow))), - Pre::Bin => return Some((u, *f)), - Pre::No => {} - } - } - } - None -} - -fn to_kelvin(name: &str, x: f64) -> f64 { - match name { - "C" | "cel" => x + 273.15, - "F" | "fah" => (x - 32.0) * 5.0 / 9.0 + 273.15, - "Rank" => x * 5.0 / 9.0, - "Reau" => x * 1.25 + 273.15, - _ => x, - } -} - -fn from_kelvin(name: &str, k: f64) -> f64 { - match name { - "C" | "cel" => k - 273.15, - "F" | "fah" => (k - 273.15) * 9.0 / 5.0 + 32.0, - "Rank" => k * 9.0 / 5.0, - "Reau" => (k - 273.15) / 1.25, - _ => k, - } -} - -fn convert(a: &[Arg], _c: &mut dyn Ctx) -> R { - let x = num(a, 0)?; - let from = text(a, 1)?; - let to = text(a, 2)?; - let (fu, fp) = resolve(&from).ok_or(CellError::NA)?; - let (tu, tp) = resolve(&to).ok_or(CellError::NA)?; - if fu.dim != tu.dim { - return Err(CellError::NA); - } - if fu.dim == Dim::Temperature { - let fname = fu.names.first().copied().unwrap_or("K"); - let tname = tu.names.first().copied().unwrap_or("K"); - if fname == tname && fp == tp { - return num_val(x); - } - let k = to_kelvin(fname, x * fp); - return num_val(from_kelvin(tname, k) / tp); - } - let base = x * fu.factor * fp; - num_val(base / (tu.factor * tp)) -} - // --------------------------------------------------------------------------------------------- // DELTA, GESTEP, ERF @@ -502,201 +225,6 @@ fn erfc_fn(a: &[Arg], _c: &mut dyn Ctx) -> R { num_val(erfc(num(a, 0)?)) } -// --------------------------------------------------------------------------------------------- -// Complex numbers - -#[derive(Clone, Copy, Debug, PartialEq)] -struct Cx { - re: f64, - im: f64, -} - -impl Cx { - fn new(re: f64, im: f64) -> Cx { - Cx { re, im } - } - fn add(self, o: Cx) -> Cx { - Cx::new(self.re + o.re, self.im + o.im) - } - fn sub(self, o: Cx) -> Cx { - Cx::new(self.re - o.re, self.im - o.im) - } - fn mul(self, o: Cx) -> Cx { - Cx::new(self.re * o.re - self.im * o.im, self.re * o.im + self.im * o.re) - } - fn div(self, o: Cx) -> R { - let d = o.re * o.re + o.im * o.im; - if d == 0.0 { - return Err(CellError::Num); - } - Ok(Cx::new((self.re * o.re + self.im * o.im) / d, (self.im * o.re - self.re * o.im) / d)) - } - fn abs(self) -> f64 { - self.re.hypot(self.im) - } - fn arg(self) -> f64 { - self.im.atan2(self.re) - } - fn exp(self) -> Cx { - let e = self.re.exp(); - Cx::new(e * self.im.cos(), e * self.im.sin()) - } - fn ln(self) -> R { - if self.re == 0.0 && self.im == 0.0 { - return Err(CellError::Num); - } - Ok(Cx::new(self.abs().ln(), self.arg())) - } - fn sin(self) -> Cx { - Cx::new(self.re.sin() * self.im.cosh(), self.re.cos() * self.im.sinh()) - } - fn cos(self) -> Cx { - Cx::new(self.re.cos() * self.im.cosh(), -self.re.sin() * self.im.sinh()) - } - fn sinh(self) -> Cx { - Cx::new(self.re.sinh() * self.im.cos(), self.re.cosh() * self.im.sin()) - } - fn cosh(self) -> Cx { - Cx::new(self.re.cosh() * self.im.cos(), self.re.sinh() * self.im.sin()) - } - fn powf(self, n: f64) -> R { - if self.re == 0.0 && self.im == 0.0 { - return if n > 0.0 { Ok(Cx::new(0.0, 0.0)) } else { Err(CellError::Num) }; - } - let r = self.abs().powf(n); - let t = self.arg() * n; - Ok(Cx::new(r * t.cos(), r * t.sin())) - } -} - -const ONE: Cx = Cx { re: 1.0, im: 0.0 }; - -/// Parses a real number written in a complex string (no spaces, optional exponent). -fn parse_real(s: &str) -> Option { - if s.is_empty() || s.contains(char::is_whitespace) { - return None; - } - let ok = s.chars().all(|c| c.is_ascii_digit() || matches!(c, '.' | 'e' | 'E' | '+' | '-')); - if !ok { - return None; - } - s.parse::().ok().filter(|v| v.is_finite()) -} - -/// Parses "a+bi" forms. Returns the number and its suffix ('i', 'j' or none). -fn parse_complex(s: &str) -> R<(Cx, Option)> { - if s.is_empty() { - return Ok((Cx::new(0.0, 0.0), None)); - } - let suffix = match s.chars().last() { - Some(c @ ('i' | 'j')) => Some(c), - _ => None, - }; - let Some(sfx) = suffix else { - return parse_real(s).map(|r| (Cx::new(r, 0.0), None)).ok_or(CellError::Num); - }; - let body = &s[..s.len() - 1]; - // Split at the last sign that is not at the start and not part of an exponent. - let bytes = body.as_bytes(); - let mut split = None; - for i in (1..bytes.len()).rev() { - let c = bytes.get(i).copied().unwrap_or(0); - let prev = bytes.get(i - 1).copied().unwrap_or(0); - if (c == b'+' || c == b'-') && prev != b'e' && prev != b'E' { - split = Some(i); - break; - } - } - let (re_s, im_s) = match split { - Some(i) => (&body[..i], &body[i..]), - None => ("", body), - }; - let re = if re_s.is_empty() { 0.0 } else { parse_real(re_s).ok_or(CellError::Num)? }; - let im = match im_s { - "" | "+" => 1.0, - "-" => -1.0, - t => parse_real(t).ok_or(CellError::Num)?, - }; - Ok((Cx::new(re, im), Some(sfx))) -} - -fn cx_of(v: &Value) -> R<(Cx, Option)> { - match v { - Value::Number(n) => Ok((Cx::new(*n, 0.0), None)), - Value::Empty => Ok((Cx::new(0.0, 0.0), None)), - Value::Bool(_) => Err(CellError::Value), - Value::Error(e) => Err(*e), - Value::Text(t) => parse_complex(t), - Value::Array(a) => cx_of(a.data.first().unwrap_or(&Value::Empty)), - } -} - -fn cx_arg(a: &[Arg], i: usize) -> R<(Cx, Option)> { - cx_of(&a.get(i).map(scalar).unwrap_or(Value::Empty)) -} - -fn merge_suffix(acc: Option, s: Option) -> R> { - match (acc, s) { - (Some(x), Some(y)) if x != y => Err(CellError::Value), - (Some(x), _) => Ok(Some(x)), - (None, y) => Ok(y), - } -} - -fn fmt_complex(z: Cx, suffix: char) -> R { - if !z.re.is_finite() || !z.im.is_finite() { - return Err(CellError::Num); - } - let re = clean(z.re); - let im = clean(z.im); - let mut s = String::new(); - if im == 0.0 { - s.push_str(&number_to_text(re)); - return text_val(s); - } - if re != 0.0 { - s.push_str(&number_to_text(re)); - if im > 0.0 { - s.push('+'); - } - } - if im == -1.0 { - s.push('-'); - } else if im != 1.0 { - s.push_str(&number_to_text(im)); - } - s.push(suffix); - text_val(s) -} - -/// Rounds to 15 significant digits so values like 0.999999999999999999 print as 1. -fn clean(x: f64) -> f64 { - let s = format!("{:.14e}", x); - s.parse().unwrap_or(x) -} - -fn complex(a: &[Arg], _c: &mut dyn Ctx) -> R { - let re = num(a, 0)?; - let im = num(a, 1)?; - let sfx = if has(a, 2) { text(a, 2)? } else { "i".into() }; - let sfx = match sfx.as_str() { - "i" | "" => 'i', - "j" => 'j', - _ => return Err(CellError::Value), - }; - fmt_complex(Cx::new(re, im), sfx) -} - -fn im_unary(a: &[Arg], f: fn(Cx) -> R) -> R { - let (z, s) = cx_arg(a, 0)?; - fmt_complex(f(z)?, s.unwrap_or('i')) -} - -fn im_real_fn(a: &[Arg], f: fn(Cx) -> R) -> R { - let (z, _) = cx_arg(a, 0)?; - num_val(f(z)?) -} - macro_rules! im1 { ($name:ident, $f:expr) => { fn $name(a: &[Arg], _c: &mut dyn Ctx) -> R { @@ -705,6 +233,16 @@ macro_rules! im1 { }; } +mod bessel; +mod complex; +mod units; + +use bessel::{besseli, besselj, besselk, bessely}; +#[cfg(test)] +use complex::parse_complex; +use complex::{Cx, ONE, complex, im_unary, imabs, imaginary, imargument, imdiv, impower, improduct, imreal, imsub, imsum}; +use units::convert; + im1!(imconjugate, |z: Cx| Ok(Cx::new(z.re, -z.im))); im1!(imexp, |z: Cx| Ok(z.exp())); im1!(imln, |z: Cx| z.ln()); @@ -727,310 +265,6 @@ im1!(imcosh, |z: Cx| Ok(z.cosh())); im1!(imsech, |z: Cx| ONE.div(z.cosh())); im1!(imcsch, |z: Cx| ONE.div(z.sinh())); -fn imabs(a: &[Arg], _c: &mut dyn Ctx) -> R { - im_real_fn(a, |z| Ok(z.abs())) -} -fn imreal(a: &[Arg], _c: &mut dyn Ctx) -> R { - im_real_fn(a, |z| Ok(z.re)) -} -fn imaginary(a: &[Arg], _c: &mut dyn Ctx) -> R { - im_real_fn(a, |z| Ok(z.im)) -} -fn imargument(a: &[Arg], _c: &mut dyn Ctx) -> R { - im_real_fn(a, |z| if z.re == 0.0 && z.im == 0.0 { Err(CellError::Div0) } else { Ok(z.arg()) }) -} - -fn im_binary(a: &[Arg], f: fn(Cx, Cx) -> R) -> R { - let (x, s1) = cx_arg(a, 0)?; - let (y, s2) = cx_arg(a, 1)?; - let s = merge_suffix(s1, s2)?; - fmt_complex(f(x, y)?, s.unwrap_or('i')) -} -fn imdiv(a: &[Arg], _c: &mut dyn Ctx) -> R { - im_binary(a, Cx::div) -} -fn imsub(a: &[Arg], _c: &mut dyn Ctx) -> R { - im_binary(a, |x, y| Ok(x.sub(y))) -} -fn impower(a: &[Arg], _c: &mut dyn Ctx) -> R { - let (z, s) = cx_arg(a, 0)?; - let n = num(a, 1)?; - // Integer powers by repeated multiplication are exact for small Gaussian integers. - let r = if n == n.trunc() && n.abs() <= 64.0 && !(z.re == 0.0 && z.im == 0.0) { - let mut acc = ONE; - for _ in 0..(n.abs() as u32) { - acc = acc.mul(z); - } - if n < 0.0 { ONE.div(acc)? } else { acc } - } else { - z.powf(n)? - }; - fmt_complex(r, s.unwrap_or('i')) -} - -fn im_fold(a: &[Arg], init: Cx, f: fn(Cx, Cx) -> Cx) -> R { - let mut acc = init; - let mut sfx = None; - for arg in a { - let direct = !matches!(arg.value, Value::Array(_)); - for v in flatten(std::slice::from_ref(arg)) { - if matches!(v, Value::Empty) && !direct { - continue; - } - let (z, s) = cx_of(&v)?; - sfx = merge_suffix(sfx, s)?; - acc = f(acc, z); - } - } - fmt_complex(acc, sfx.unwrap_or('i')) -} -fn imsum(a: &[Arg], _c: &mut dyn Ctx) -> R { - im_fold(a, Cx::new(0.0, 0.0), Cx::add) -} -fn improduct(a: &[Arg], _c: &mut dyn Ctx) -> R { - im_fold(a, ONE, Cx::mul) -} - -// --------------------------------------------------------------------------------------------- -// Bessel functions - -const MAX_BESSEL_POINTS: f64 = 4_000_000.0; - -/// J_n(x) for integer n ≥ 0 via Bessel's integral, J_n(x) = (1/2π)∫₀^{2π} cos(nτ − x sin τ) dτ, -/// evaluated with the trapezoid rule, which converges exponentially for periodic integrands. -fn bessel_j(n: u32, x: f64) -> Option { - let pts = (2.0 * (x.abs() + n as f64) + 64.0).ceil(); - if pts > MAX_BESSEL_POINTS { - return if x.abs() > 1e6 { hankel(n, x.abs()).map(|(j, _)| if x < 0.0 && n % 2 == 1 { -j } else { j }) } else { None }; - } - let m = pts as usize; - let h = 2.0 * PI / m as f64; - let nf = n as f64; - let mut sum = 0.0; - for k in 0..m { - let t = k as f64 * h; - sum += (nf * t - x * t.sin()).cos(); - } - Some(sum / m as f64) -} - -/// Hankel asymptotic expansion: (J_n(x), Y_n(x)) for large x. -fn hankel(n: u32, x: f64) -> Option<(f64, f64)> { - let mu = 4.0 * (n as f64) * (n as f64); - let mut p = 0.0; - let mut q = 0.0; - let mut term: f64 = 1.0; - let mut last = f64::INFINITY; - for k in 0..200 { - if term.abs() > last { - break; - } - let t = if k % 4 < 2 { term } else { -term }; - if k % 2 == 0 { - p += t; - } else { - q += t; - } - last = term.abs(); - if last < 1e-17 { - break; - } - let kk = (k + 1) as f64; - let odd = 2.0 * kk - 1.0; - term *= (mu - odd * odd) / (kk * 8.0 * x); - } - let chi = x - (n as f64 / 2.0 + 0.25) * PI; - let f = (2.0 / (PI * x)).sqrt(); - let j = f * (p * chi.cos() - q * chi.sin()); - let y = f * (p * chi.sin() + q * chi.cos()); - if j.is_finite() && y.is_finite() { Some((j, y)) } else { None } -} - -const EULER: f64 = std::f64::consts::EULER_GAMMA; - -/// J_0..J_m(x) for 0 < x ≤ ~100 by Miller's backward recurrence, normalised with -/// J_0 + 2ΣJ_{2k} = 1. -fn bessel_j_all(x: f64) -> Vec { - let top = x.ceil() as usize + 30 + (40.0 * x.max(1.0)).sqrt() as usize; - let m = top + top % 2; - let mut j = vec![0.0f64; m + 2]; - if let Some(v) = j.get_mut(m) { - *v = 1e-30; - } - for k in (1..=m).rev() { - let jk = j.get(k).copied().unwrap_or(0.0); - let jk1 = j.get(k + 1).copied().unwrap_or(0.0); - let mut prev = 2.0 * k as f64 / x * jk - jk1; - if prev.abs() > 1e250 { - for v in j.iter_mut() { - *v *= 1e-250; - } - prev *= 1e-250; - } - if let Some(v) = j.get_mut(k - 1) { - *v = prev; - } - } - let norm: f64 = j.first().copied().unwrap_or(0.0) + 2.0 * j.iter().skip(2).step_by(2).sum::(); - if norm != 0.0 { - for v in j.iter_mut() { - *v /= norm; - } - } - j -} - -/// Y_0 and Y_1 by Neumann series over J_k (x ≤ 25) or the asymptotic expansion. -fn bessel_y01(x: f64) -> Option<(f64, f64)> { - if x > 25.0 { - let (_, y0) = hankel(0, x)?; - let (_, y1) = hankel(1, x)?; - return Some((y0, y1)); - } - let j = bessel_j_all(x); - let jj = |k: usize| j.get(k).copied().unwrap_or(0.0); - let l = (x / 2.0).ln() + EULER; - // Y0 = (2/π)(L J0) − (4/π) Σ (−1)^k J_{2k}/k - // Y1 = −Y0' = −(2/π)(J0/x − L J1) + (2/π) Σ (−1)^k (J_{2k−1} − J_{2k+1})/k - let mut s0 = 0.0; - let mut s1 = 0.0; - let mut k = 1usize; - while 2 * k + 1 < j.len() { - let sign = if k.is_multiple_of(2) { 1.0 } else { -1.0 }; - s0 += sign * jj(2 * k) / k as f64; - s1 += sign * (jj(2 * k - 1) - jj(2 * k + 1)) / k as f64; - k += 1; - } - let y0 = (2.0 / PI) * l * jj(0) - (4.0 / PI) * s0; - let y1 = -(2.0 / PI) * (jj(0) / x - l * jj(1)) + (2.0 / PI) * s1; - if y0.is_finite() && y1.is_finite() { Some((y0, y1)) } else { None } -} - -fn bessel_y(n: u32, x: f64) -> Option { - let (y0, y1) = bessel_y01(x)?; - if n == 0 { - return Some(y0); - } - let (mut a, mut b) = (y0, y1); - for k in 1..n { - let c = 2.0 * k as f64 / x * b - a; - a = b; - b = c; - if !b.is_finite() { - return None; - } - } - Some(b) -} - -/// I_n(x) by its power series (all terms positive), computed with a log-scaled first term. -fn bessel_i(n: u32, x: f64) -> Option { - if x == 0.0 { - return Some(if n == 0 { 1.0 } else { 0.0 }); - } - let ax = x.abs(); - let nf = n as f64; - let ln_t0 = nf * (ax / 2.0).ln() - crate::special::ln_gamma(nf + 1.0); - let q = ax * ax / 4.0; - let mut term = 1.0; - let mut sum = 1.0; - let mut k = 0.0; - while k < 1e6 { - k += 1.0; - term *= q / (k * (k + nf)); - sum += term; - if term < 1e-17 * sum { - break; - } - if !sum.is_finite() { - return None; - } - } - let r = (ln_t0 + sum.ln()).exp(); - let r = if x < 0.0 && n % 2 == 1 { -r } else { r }; - if r.is_finite() { Some(r) } else { None } -} - -/// K_n(x) for x > 0: K_0 and K_1 from ∫₀^∞ e^{−x cosh t} cosh(νt) dt (trapezoid rule), then -/// upward recurrence. -fn bessel_k(n: u32, x: f64) -> Option { - let k_int = |nu: f64| -> f64 { - let h: f64 = 0.02; - let mut sum = 0.5 * (-x).exp(); - let mut t: f64 = h; - while t < 800.0 { - let v = (-x * t.cosh() + (nu * t).ln_cosh_safe()).exp(); - sum += v; - if v < 1e-18 * sum && x * t.cosh() > 50.0 { - break; - } - t += h; - } - sum * h - }; - let k0 = k_int(0.0); - if n == 0 { - return Some(k0); - } - let k1 = k_int(1.0); - let (mut a, mut b) = (k0, k1); - for k in 1..n { - let c = a + 2.0 * k as f64 / x * b; - a = b; - b = c; - if !b.is_finite() { - return None; - } - } - if b.is_finite() { Some(b) } else { None } -} - -trait LnCosh { - fn ln_cosh_safe(self) -> f64; -} -impl LnCosh for f64 { - /// ln(cosh(x)) without overflow. - fn ln_cosh_safe(self) -> f64 { - let a = self.abs(); - a + (-2.0 * a).exp().ln_1p() - std::f64::consts::LN_2 - } -} - -fn bessel_args(a: &[Arg]) -> R<(f64, u32)> { - let x = num(a, 0)?; - let n = num(a, 1)?.trunc(); - if n < 0.0 { - return Err(CellError::Num); - } - if n > 100_000.0 { - return Err(CellError::Num); - } - Ok((x, n as u32)) -} - -fn besselj(a: &[Arg], _c: &mut dyn Ctx) -> R { - let (x, n) = bessel_args(a)?; - num_val(bessel_j(n, x).ok_or(CellError::Num)?) -} -fn bessely(a: &[Arg], _c: &mut dyn Ctx) -> R { - let (x, n) = bessel_args(a)?; - if x <= 0.0 { - return Err(CellError::Num); - } - num_val(bessel_y(n, x).ok_or(CellError::Num)?) -} -fn besseli(a: &[Arg], _c: &mut dyn Ctx) -> R { - let (x, n) = bessel_args(a)?; - num_val(bessel_i(n, x).ok_or(CellError::Num)?) -} -fn besselk(a: &[Arg], _c: &mut dyn Ctx) -> R { - let (x, n) = bessel_args(a)?; - if x <= 0.0 { - return Err(CellError::Num); - } - num_val(bessel_k(n, x).ok_or(CellError::Num)?) -} - pub(crate) fn specs() -> Vec { vec![ f!("BIN2DEC", 1, 1, Engineering, S, "BIN2DEC(number)", "Converts a binary number to decimal.", bin2dec), diff --git a/crates/functions/src/engineering/bessel.rs b/crates/functions/src/engineering/bessel.rs new file mode 100644 index 000000000..34be81783 --- /dev/null +++ b/crates/functions/src/engineering/bessel.rs @@ -0,0 +1,245 @@ +use super::*; + +// --------------------------------------------------------------------------------------------- +// Bessel functions + +const MAX_BESSEL_POINTS: f64 = 4_000_000.0; + +/// J_n(x) for integer n ≥ 0 via Bessel's integral, J_n(x) = (1/2π)∫₀^{2π} cos(nτ − x sin τ) dτ, +/// evaluated with the trapezoid rule, which converges exponentially for periodic integrands. +fn bessel_j(n: u32, x: f64) -> Option { + let pts = (2.0 * (x.abs() + n as f64) + 64.0).ceil(); + if pts > MAX_BESSEL_POINTS { + return if x.abs() > 1e6 { hankel(n, x.abs()).map(|(j, _)| if x < 0.0 && n % 2 == 1 { -j } else { j }) } else { None }; + } + let m = pts as usize; + let h = 2.0 * PI / m as f64; + let nf = n as f64; + let mut sum = 0.0; + for k in 0..m { + let t = k as f64 * h; + sum += (nf * t - x * t.sin()).cos(); + } + Some(sum / m as f64) +} + +/// Hankel asymptotic expansion: (J_n(x), Y_n(x)) for large x. +fn hankel(n: u32, x: f64) -> Option<(f64, f64)> { + let mu = 4.0 * (n as f64) * (n as f64); + let mut p = 0.0; + let mut q = 0.0; + let mut term: f64 = 1.0; + let mut last = f64::INFINITY; + for k in 0..200 { + if term.abs() > last { + break; + } + let t = if k % 4 < 2 { term } else { -term }; + if k % 2 == 0 { + p += t; + } else { + q += t; + } + last = term.abs(); + if last < 1e-17 { + break; + } + let kk = (k + 1) as f64; + let odd = 2.0 * kk - 1.0; + term *= (mu - odd * odd) / (kk * 8.0 * x); + } + let chi = x - (n as f64 / 2.0 + 0.25) * PI; + let f = (2.0 / (PI * x)).sqrt(); + let j = f * (p * chi.cos() - q * chi.sin()); + let y = f * (p * chi.sin() + q * chi.cos()); + if j.is_finite() && y.is_finite() { Some((j, y)) } else { None } +} + +const EULER: f64 = std::f64::consts::EULER_GAMMA; + +/// J_0..J_m(x) for 0 < x ≤ ~100 by Miller's backward recurrence, normalised with +/// J_0 + 2ΣJ_{2k} = 1. +fn bessel_j_all(x: f64) -> Vec { + let top = x.ceil() as usize + 30 + (40.0 * x.max(1.0)).sqrt() as usize; + let m = top + top % 2; + let mut j = vec![0.0f64; m + 2]; + if let Some(v) = j.get_mut(m) { + *v = 1e-30; + } + for k in (1..=m).rev() { + let jk = j.get(k).copied().unwrap_or(0.0); + let jk1 = j.get(k + 1).copied().unwrap_or(0.0); + let mut prev = 2.0 * k as f64 / x * jk - jk1; + if prev.abs() > 1e250 { + for v in j.iter_mut() { + *v *= 1e-250; + } + prev *= 1e-250; + } + if let Some(v) = j.get_mut(k - 1) { + *v = prev; + } + } + let norm: f64 = j.first().copied().unwrap_or(0.0) + 2.0 * j.iter().skip(2).step_by(2).sum::(); + if norm != 0.0 { + for v in j.iter_mut() { + *v /= norm; + } + } + j +} + +/// Y_0 and Y_1 by Neumann series over J_k (x ≤ 25) or the asymptotic expansion. +fn bessel_y01(x: f64) -> Option<(f64, f64)> { + if x > 25.0 { + let (_, y0) = hankel(0, x)?; + let (_, y1) = hankel(1, x)?; + return Some((y0, y1)); + } + let j = bessel_j_all(x); + let jj = |k: usize| j.get(k).copied().unwrap_or(0.0); + let l = (x / 2.0).ln() + EULER; + // Y0 = (2/π)(L J0) − (4/π) Σ (−1)^k J_{2k}/k + // Y1 = −Y0' = −(2/π)(J0/x − L J1) + (2/π) Σ (−1)^k (J_{2k−1} − J_{2k+1})/k + let mut s0 = 0.0; + let mut s1 = 0.0; + let mut k = 1usize; + while 2 * k + 1 < j.len() { + let sign = if k.is_multiple_of(2) { 1.0 } else { -1.0 }; + s0 += sign * jj(2 * k) / k as f64; + s1 += sign * (jj(2 * k - 1) - jj(2 * k + 1)) / k as f64; + k += 1; + } + let y0 = (2.0 / PI) * l * jj(0) - (4.0 / PI) * s0; + let y1 = -(2.0 / PI) * (jj(0) / x - l * jj(1)) + (2.0 / PI) * s1; + if y0.is_finite() && y1.is_finite() { Some((y0, y1)) } else { None } +} + +fn bessel_y(n: u32, x: f64) -> Option { + let (y0, y1) = bessel_y01(x)?; + if n == 0 { + return Some(y0); + } + let (mut a, mut b) = (y0, y1); + for k in 1..n { + let c = 2.0 * k as f64 / x * b - a; + a = b; + b = c; + if !b.is_finite() { + return None; + } + } + Some(b) +} + +/// I_n(x) by its power series (all terms positive), computed with a log-scaled first term. +fn bessel_i(n: u32, x: f64) -> Option { + if x == 0.0 { + return Some(if n == 0 { 1.0 } else { 0.0 }); + } + let ax = x.abs(); + let nf = n as f64; + let ln_t0 = nf * (ax / 2.0).ln() - crate::special::ln_gamma(nf + 1.0); + let q = ax * ax / 4.0; + let mut term = 1.0; + let mut sum = 1.0; + let mut k = 0.0; + while k < 1e6 { + k += 1.0; + term *= q / (k * (k + nf)); + sum += term; + if term < 1e-17 * sum { + break; + } + if !sum.is_finite() { + return None; + } + } + let r = (ln_t0 + sum.ln()).exp(); + let r = if x < 0.0 && n % 2 == 1 { -r } else { r }; + if r.is_finite() { Some(r) } else { None } +} + +/// K_n(x) for x > 0: K_0 and K_1 from ∫₀^∞ e^{−x cosh t} cosh(νt) dt (trapezoid rule), then +/// upward recurrence. +fn bessel_k(n: u32, x: f64) -> Option { + let k_int = |nu: f64| -> f64 { + let h: f64 = 0.02; + let mut sum = 0.5 * (-x).exp(); + let mut t: f64 = h; + while t < 800.0 { + let v = (-x * t.cosh() + (nu * t).ln_cosh_safe()).exp(); + sum += v; + if v < 1e-18 * sum && x * t.cosh() > 50.0 { + break; + } + t += h; + } + sum * h + }; + let k0 = k_int(0.0); + if n == 0 { + return Some(k0); + } + let k1 = k_int(1.0); + let (mut a, mut b) = (k0, k1); + for k in 1..n { + let c = a + 2.0 * k as f64 / x * b; + a = b; + b = c; + if !b.is_finite() { + return None; + } + } + if b.is_finite() { Some(b) } else { None } +} + +trait LnCosh { + fn ln_cosh_safe(self) -> f64; +} + +impl LnCosh for f64 { + /// ln(cosh(x)) without overflow. + fn ln_cosh_safe(self) -> f64 { + let a = self.abs(); + a + (-2.0 * a).exp().ln_1p() - std::f64::consts::LN_2 + } +} + +fn bessel_args(a: &[Arg]) -> R<(f64, u32)> { + let x = num(a, 0)?; + let n = num(a, 1)?.trunc(); + if n < 0.0 { + return Err(CellError::Num); + } + if n > 100_000.0 { + return Err(CellError::Num); + } + Ok((x, n as u32)) +} + +pub(super) fn besselj(a: &[Arg], _c: &mut dyn Ctx) -> R { + let (x, n) = bessel_args(a)?; + num_val(bessel_j(n, x).ok_or(CellError::Num)?) +} + +pub(super) fn bessely(a: &[Arg], _c: &mut dyn Ctx) -> R { + let (x, n) = bessel_args(a)?; + if x <= 0.0 { + return Err(CellError::Num); + } + num_val(bessel_y(n, x).ok_or(CellError::Num)?) +} + +pub(super) fn besseli(a: &[Arg], _c: &mut dyn Ctx) -> R { + let (x, n) = bessel_args(a)?; + num_val(bessel_i(n, x).ok_or(CellError::Num)?) +} + +pub(super) fn besselk(a: &[Arg], _c: &mut dyn Ctx) -> R { + let (x, n) = bessel_args(a)?; + if x <= 0.0 { + return Err(CellError::Num); + } + num_val(bessel_k(n, x).ok_or(CellError::Num)?) +} diff --git a/crates/functions/src/engineering/complex.rs b/crates/functions/src/engineering/complex.rs new file mode 100644 index 000000000..86ab6b535 --- /dev/null +++ b/crates/functions/src/engineering/complex.rs @@ -0,0 +1,268 @@ +use super::*; + +// --------------------------------------------------------------------------------------------- +// Complex numbers + +#[derive(Clone, Copy, Debug, PartialEq)] +pub(super) struct Cx { + pub(super) re: f64, + pub(super) im: f64, +} + +impl Cx { + pub(super) fn new(re: f64, im: f64) -> Cx { + Cx { re, im } + } + fn add(self, o: Cx) -> Cx { + Cx::new(self.re + o.re, self.im + o.im) + } + fn sub(self, o: Cx) -> Cx { + Cx::new(self.re - o.re, self.im - o.im) + } + fn mul(self, o: Cx) -> Cx { + Cx::new(self.re * o.re - self.im * o.im, self.re * o.im + self.im * o.re) + } + pub(super) fn div(self, o: Cx) -> R { + let d = o.re * o.re + o.im * o.im; + if d == 0.0 { + return Err(CellError::Num); + } + Ok(Cx::new((self.re * o.re + self.im * o.im) / d, (self.im * o.re - self.re * o.im) / d)) + } + pub(super) fn abs(self) -> f64 { + self.re.hypot(self.im) + } + fn arg(self) -> f64 { + self.im.atan2(self.re) + } + pub(super) fn exp(self) -> Cx { + let e = self.re.exp(); + Cx::new(e * self.im.cos(), e * self.im.sin()) + } + pub(super) fn ln(self) -> R { + if self.re == 0.0 && self.im == 0.0 { + return Err(CellError::Num); + } + Ok(Cx::new(self.abs().ln(), self.arg())) + } + pub(super) fn sin(self) -> Cx { + Cx::new(self.re.sin() * self.im.cosh(), self.re.cos() * self.im.sinh()) + } + pub(super) fn cos(self) -> Cx { + Cx::new(self.re.cos() * self.im.cosh(), -self.re.sin() * self.im.sinh()) + } + pub(super) fn sinh(self) -> Cx { + Cx::new(self.re.sinh() * self.im.cos(), self.re.cosh() * self.im.sin()) + } + pub(super) fn cosh(self) -> Cx { + Cx::new(self.re.cosh() * self.im.cos(), self.re.sinh() * self.im.sin()) + } + fn powf(self, n: f64) -> R { + if self.re == 0.0 && self.im == 0.0 { + return if n > 0.0 { Ok(Cx::new(0.0, 0.0)) } else { Err(CellError::Num) }; + } + let r = self.abs().powf(n); + let t = self.arg() * n; + Ok(Cx::new(r * t.cos(), r * t.sin())) + } +} + +pub(super) const ONE: Cx = Cx { re: 1.0, im: 0.0 }; + +/// Parses a real number written in a complex string (no spaces, optional exponent). +fn parse_real(s: &str) -> Option { + if s.is_empty() || s.contains(char::is_whitespace) { + return None; + } + let ok = s.chars().all(|c| c.is_ascii_digit() || matches!(c, '.' | 'e' | 'E' | '+' | '-')); + if !ok { + return None; + } + s.parse::().ok().filter(|v| v.is_finite()) +} + +/// Parses "a+bi" forms. Returns the number and its suffix ('i', 'j' or none). +pub(super) fn parse_complex(s: &str) -> R<(Cx, Option)> { + if s.is_empty() { + return Ok((Cx::new(0.0, 0.0), None)); + } + let suffix = match s.chars().last() { + Some(c @ ('i' | 'j')) => Some(c), + _ => None, + }; + let Some(sfx) = suffix else { + return parse_real(s).map(|r| (Cx::new(r, 0.0), None)).ok_or(CellError::Num); + }; + let body = &s[..s.len() - 1]; + // Split at the last sign that is not at the start and not part of an exponent. + let bytes = body.as_bytes(); + let mut split = None; + for i in (1..bytes.len()).rev() { + let c = bytes.get(i).copied().unwrap_or(0); + let prev = bytes.get(i - 1).copied().unwrap_or(0); + if (c == b'+' || c == b'-') && prev != b'e' && prev != b'E' { + split = Some(i); + break; + } + } + let (re_s, im_s) = match split { + Some(i) => (&body[..i], &body[i..]), + None => ("", body), + }; + let re = if re_s.is_empty() { 0.0 } else { parse_real(re_s).ok_or(CellError::Num)? }; + let im = match im_s { + "" | "+" => 1.0, + "-" => -1.0, + t => parse_real(t).ok_or(CellError::Num)?, + }; + Ok((Cx::new(re, im), Some(sfx))) +} + +fn cx_of(v: &Value) -> R<(Cx, Option)> { + match v { + Value::Number(n) => Ok((Cx::new(*n, 0.0), None)), + Value::Empty => Ok((Cx::new(0.0, 0.0), None)), + Value::Bool(_) => Err(CellError::Value), + Value::Error(e) => Err(*e), + Value::Text(t) => parse_complex(t), + Value::Array(a) => cx_of(a.data.first().unwrap_or(&Value::Empty)), + } +} + +fn cx_arg(a: &[Arg], i: usize) -> R<(Cx, Option)> { + cx_of(&a.get(i).map(scalar).unwrap_or(Value::Empty)) +} + +fn merge_suffix(acc: Option, s: Option) -> R> { + match (acc, s) { + (Some(x), Some(y)) if x != y => Err(CellError::Value), + (Some(x), _) => Ok(Some(x)), + (None, y) => Ok(y), + } +} + +fn fmt_complex(z: Cx, suffix: char) -> R { + if !z.re.is_finite() || !z.im.is_finite() { + return Err(CellError::Num); + } + let re = clean(z.re); + let im = clean(z.im); + let mut s = String::new(); + if im == 0.0 { + s.push_str(&number_to_text(re)); + return text_val(s); + } + if re != 0.0 { + s.push_str(&number_to_text(re)); + if im > 0.0 { + s.push('+'); + } + } + if im == -1.0 { + s.push('-'); + } else if im != 1.0 { + s.push_str(&number_to_text(im)); + } + s.push(suffix); + text_val(s) +} + +/// Rounds to 15 significant digits so values like 0.999999999999999999 print as 1. +fn clean(x: f64) -> f64 { + let s = format!("{:.14e}", x); + s.parse().unwrap_or(x) +} + +pub(super) fn complex(a: &[Arg], _c: &mut dyn Ctx) -> R { + let re = num(a, 0)?; + let im = num(a, 1)?; + let sfx = if has(a, 2) { text(a, 2)? } else { "i".into() }; + let sfx = match sfx.as_str() { + "i" | "" => 'i', + "j" => 'j', + _ => return Err(CellError::Value), + }; + fmt_complex(Cx::new(re, im), sfx) +} + +pub(super) fn im_unary(a: &[Arg], f: fn(Cx) -> R) -> R { + let (z, s) = cx_arg(a, 0)?; + fmt_complex(f(z)?, s.unwrap_or('i')) +} + +fn im_real_fn(a: &[Arg], f: fn(Cx) -> R) -> R { + let (z, _) = cx_arg(a, 0)?; + num_val(f(z)?) +} + +pub(super) fn imabs(a: &[Arg], _c: &mut dyn Ctx) -> R { + im_real_fn(a, |z| Ok(z.abs())) +} + +pub(super) fn imreal(a: &[Arg], _c: &mut dyn Ctx) -> R { + im_real_fn(a, |z| Ok(z.re)) +} + +pub(super) fn imaginary(a: &[Arg], _c: &mut dyn Ctx) -> R { + im_real_fn(a, |z| Ok(z.im)) +} + +pub(super) fn imargument(a: &[Arg], _c: &mut dyn Ctx) -> R { + im_real_fn(a, |z| if z.re == 0.0 && z.im == 0.0 { Err(CellError::Div0) } else { Ok(z.arg()) }) +} + +fn im_binary(a: &[Arg], f: fn(Cx, Cx) -> R) -> R { + let (x, s1) = cx_arg(a, 0)?; + let (y, s2) = cx_arg(a, 1)?; + let s = merge_suffix(s1, s2)?; + fmt_complex(f(x, y)?, s.unwrap_or('i')) +} + +pub(super) fn imdiv(a: &[Arg], _c: &mut dyn Ctx) -> R { + im_binary(a, Cx::div) +} + +pub(super) fn imsub(a: &[Arg], _c: &mut dyn Ctx) -> R { + im_binary(a, |x, y| Ok(x.sub(y))) +} + +pub(super) fn impower(a: &[Arg], _c: &mut dyn Ctx) -> R { + let (z, s) = cx_arg(a, 0)?; + let n = num(a, 1)?; + // Integer powers by repeated multiplication are exact for small Gaussian integers. + let r = if n == n.trunc() && n.abs() <= 64.0 && !(z.re == 0.0 && z.im == 0.0) { + let mut acc = ONE; + for _ in 0..(n.abs() as u32) { + acc = acc.mul(z); + } + if n < 0.0 { ONE.div(acc)? } else { acc } + } else { + z.powf(n)? + }; + fmt_complex(r, s.unwrap_or('i')) +} + +fn im_fold(a: &[Arg], init: Cx, f: fn(Cx, Cx) -> Cx) -> R { + let mut acc = init; + let mut sfx = None; + for arg in a { + let direct = !matches!(arg.value, Value::Array(_)); + for v in flatten(std::slice::from_ref(arg)) { + if matches!(v, Value::Empty) && !direct { + continue; + } + let (z, s) = cx_of(&v)?; + sfx = merge_suffix(sfx, s)?; + acc = f(acc, z); + } + } + fmt_complex(acc, sfx.unwrap_or('i')) +} + +pub(super) fn imsum(a: &[Arg], _c: &mut dyn Ctx) -> R { + im_fold(a, Cx::new(0.0, 0.0), Cx::add) +} + +pub(super) fn improduct(a: &[Arg], _c: &mut dyn Ctx) -> R { + im_fold(a, ONE, Cx::mul) +} diff --git a/crates/functions/src/engineering/units.rs b/crates/functions/src/engineering/units.rs new file mode 100644 index 000000000..e54abd87a --- /dev/null +++ b/crates/functions/src/engineering/units.rs @@ -0,0 +1,286 @@ +use super::*; + +// --------------------------------------------------------------------------------------------- +// CONVERT + +#[derive(Clone, Copy, PartialEq, Eq, Debug)] +enum Dim { + Mass, + Length, + Time, + Pressure, + Force, + Energy, + Power, + Magnetism, + Temperature, + Volume, + Area, + Information, + Speed, +} + +/// How a unit accepts prefixes: none, SI (raised to `power` for areas/volumes) or SI + binary. +#[derive(Clone, Copy, PartialEq, Eq)] +enum Pre { + No, + Si(i32), + Bin, +} + +struct Unit { + names: &'static [&'static str], + dim: Dim, + factor: f64, + pre: Pre, +} + +const fn u(names: &'static [&'static str], dim: Dim, factor: f64, pre: Pre) -> Unit { + Unit { names, dim, factor, pre } +} + +const IN: f64 = 0.0254; + +const FT: f64 = 0.3048; + +const YD: f64 = 0.9144; + +const MI: f64 = 1609.344; + +const NMI: f64 = 1852.0; + +const LY: f64 = 9_460_730_472_580_800.0; + +const PICA_PT: f64 = 0.0254 / 72.0; + +const PICA: f64 = 0.0254 / 6.0; + +const GAL: f64 = 0.003_785_411_784; + +const UNITS: &[Unit] = &[ + // Mass (gram) + u(&["g"], Dim::Mass, 1.0, Pre::Si(1)), + u(&["sg"], Dim::Mass, 14_593.902_937_206_4, Pre::No), + u(&["lbm"], Dim::Mass, 453.592_37, Pre::No), + u(&["u"], Dim::Mass, 1.660_539_066_6e-24, Pre::Si(1)), + u(&["ozm"], Dim::Mass, 28.349_523_125, Pre::No), + u(&["grain"], Dim::Mass, 0.064_798_91, Pre::No), + u(&["cwt", "shweight"], Dim::Mass, 45_359.237, Pre::No), + u(&["uk_cwt", "lcwt", "hweight"], Dim::Mass, 50_802.345_44, Pre::No), + u(&["stone"], Dim::Mass, 6_350.293_18, Pre::No), + u(&["ton"], Dim::Mass, 907_184.74, Pre::No), + u(&["uk_ton", "LTON", "brton"], Dim::Mass, 1_016_046.908_8, Pre::No), + // Length (metre) + u(&["m"], Dim::Length, 1.0, Pre::Si(1)), + u(&["mi"], Dim::Length, MI, Pre::No), + u(&["Nmi"], Dim::Length, NMI, Pre::No), + u(&["in"], Dim::Length, IN, Pre::No), + u(&["ft"], Dim::Length, FT, Pre::No), + u(&["yd"], Dim::Length, YD, Pre::No), + u(&["ang"], Dim::Length, 1e-10, Pre::Si(1)), + u(&["ell"], Dim::Length, 1.143, Pre::No), + u(&["ly"], Dim::Length, LY, Pre::Si(1)), + u(&["parsec", "pc"], Dim::Length, 30_856_775_814_913_673.0, Pre::Si(1)), + u(&["Picapt", "Pica"], Dim::Length, PICA_PT, Pre::No), + u(&["pica"], Dim::Length, PICA, Pre::No), + u(&["survey_mi"], Dim::Length, 1_609.347_218_694_437, Pre::No), + // Time (second) + u(&["yr"], Dim::Time, 31_557_600.0, Pre::No), + u(&["day", "d"], Dim::Time, 86_400.0, Pre::No), + u(&["hr"], Dim::Time, 3_600.0, Pre::No), + u(&["mn", "min"], Dim::Time, 60.0, Pre::No), + u(&["sec", "s"], Dim::Time, 1.0, Pre::Si(1)), + // Pressure (pascal) + u(&["Pa", "p"], Dim::Pressure, 1.0, Pre::Si(1)), + u(&["atm", "at"], Dim::Pressure, 101_325.0, Pre::Si(1)), + u(&["mmHg"], Dim::Pressure, 133.322_387_415, Pre::Si(1)), + u(&["psi"], Dim::Pressure, 6_894.757_293_168_36, Pre::No), + u(&["Torr"], Dim::Pressure, 101_325.0 / 760.0, Pre::No), + // Force (newton) + u(&["N"], Dim::Force, 1.0, Pre::Si(1)), + u(&["dyn", "dy"], Dim::Force, 1e-5, Pre::Si(1)), + u(&["lbf"], Dim::Force, 4.448_221_615_260_5, Pre::No), + u(&["pond"], Dim::Force, 0.009_806_65, Pre::Si(1)), + // Energy (joule) + u(&["J"], Dim::Energy, 1.0, Pre::Si(1)), + u(&["e"], Dim::Energy, 1e-7, Pre::Si(1)), + u(&["c"], Dim::Energy, 4.184, Pre::Si(1)), + u(&["cal"], Dim::Energy, 4.1868, Pre::Si(1)), + u(&["eV", "ev"], Dim::Energy, 1.602_176_634e-19, Pre::Si(1)), + u(&["HPh", "hh"], Dim::Energy, 2_684_519.537_696_17, Pre::No), + u(&["Wh", "wh"], Dim::Energy, 3_600.0, Pre::Si(1)), + u(&["flb"], Dim::Energy, 1.355_817_948_331_4, Pre::No), + u(&["BTU", "btu"], Dim::Energy, 1_055.055_852_62, Pre::No), + // Power (watt) + u(&["W", "w"], Dim::Power, 1.0, Pre::Si(1)), + u(&["HP", "h"], Dim::Power, 745.699_871_582_270_2, Pre::No), + u(&["PS"], Dim::Power, 735.498_75, Pre::No), + // Magnetism (tesla) + u(&["T"], Dim::Magnetism, 1.0, Pre::Si(1)), + u(&["ga"], Dim::Magnetism, 1e-4, Pre::Si(1)), + // Temperature (handled specially; factor unused) + u(&["C", "cel"], Dim::Temperature, 1.0, Pre::No), + u(&["F", "fah"], Dim::Temperature, 1.0, Pre::No), + u(&["K", "kel"], Dim::Temperature, 1.0, Pre::Si(1)), + u(&["Rank"], Dim::Temperature, 1.0, Pre::No), + u(&["Reau"], Dim::Temperature, 1.0, Pre::No), + // Volume (cubic metre) + u(&["l", "L", "lt"], Dim::Volume, 0.001, Pre::Si(1)), + u(&["tsp"], Dim::Volume, GAL / 768.0, Pre::No), + u(&["tspm"], Dim::Volume, 5e-6, Pre::No), + u(&["tbs"], Dim::Volume, GAL / 256.0, Pre::No), + u(&["oz"], Dim::Volume, GAL / 128.0, Pre::No), + u(&["cup"], Dim::Volume, GAL / 16.0, Pre::No), + u(&["pt", "us_pt"], Dim::Volume, GAL / 8.0, Pre::No), + u(&["uk_pt"], Dim::Volume, 0.000_568_261_25, Pre::No), + u(&["qt"], Dim::Volume, GAL / 4.0, Pre::No), + u(&["uk_qt"], Dim::Volume, 0.001_136_522_5, Pre::No), + u(&["gal"], Dim::Volume, GAL, Pre::No), + u(&["uk_gal"], Dim::Volume, 0.004_546_09, Pre::No), + u(&["m3"], Dim::Volume, 1.0, Pre::Si(3)), + u(&["mi3"], Dim::Volume, MI * MI * MI, Pre::No), + u(&["Nmi3"], Dim::Volume, NMI * NMI * NMI, Pre::No), + u(&["yd3"], Dim::Volume, YD * YD * YD, Pre::No), + u(&["ft3"], Dim::Volume, FT * FT * FT, Pre::No), + u(&["in3"], Dim::Volume, IN * IN * IN, Pre::No), + u(&["ang3"], Dim::Volume, 1e-30, Pre::Si(3)), + u(&["ly3"], Dim::Volume, LY * LY * LY, Pre::No), + u(&["Picapt3", "Pica3"], Dim::Volume, PICA_PT * PICA_PT * PICA_PT, Pre::No), + u(&["barrel"], Dim::Volume, 0.158_987_294_928, Pre::No), + u(&["bushel"], Dim::Volume, 0.035_239_070_166_88, Pre::No), + u(&["MTON"], Dim::Volume, 1.132_673_863_68, Pre::No), + u(&["GRT", "regton"], Dim::Volume, 2.831_684_659_2, Pre::No), + // Area (square metre) + u(&["m2"], Dim::Area, 1.0, Pre::Si(2)), + u(&["ha"], Dim::Area, 10_000.0, Pre::No), + u(&["ar"], Dim::Area, 100.0, Pre::Si(1)), + u(&["ft2"], Dim::Area, FT * FT, Pre::No), + u(&["in2"], Dim::Area, IN * IN, Pre::No), + u(&["yd2"], Dim::Area, YD * YD, Pre::No), + u(&["mi2"], Dim::Area, MI * MI, Pre::No), + u(&["Nmi2"], Dim::Area, NMI * NMI, Pre::No), + u(&["ang2"], Dim::Area, 1e-20, Pre::Si(2)), + u(&["ly2"], Dim::Area, LY * LY, Pre::No), + u(&["Picapt2", "Pica2"], Dim::Area, PICA_PT * PICA_PT, Pre::No), + u(&["uk_acre"], Dim::Area, 4_046.856_422_4, Pre::No), + u(&["us_acre"], Dim::Area, 4_046.872_609_874_252, Pre::No), + u(&["Morgen"], Dim::Area, 2_500.0, Pre::No), + // Information (bit) + u(&["bit"], Dim::Information, 1.0, Pre::Bin), + u(&["byte"], Dim::Information, 8.0, Pre::Bin), + // Speed (metres per second) + u(&["m/s", "m/sec"], Dim::Speed, 1.0, Pre::Si(1)), + u(&["m/h", "m/hr"], Dim::Speed, 1.0 / 3600.0, Pre::Si(1)), + u(&["mph"], Dim::Speed, MI / 3600.0, Pre::No), + u(&["kn"], Dim::Speed, NMI / 3600.0, Pre::No), + u(&["admkn"], Dim::Speed, 6080.0 * FT / 3600.0, Pre::No), +]; + +const SI_PREFIXES: &[(&str, f64)] = &[ + ("da", 1e1), + ("Y", 1e24), + ("Z", 1e21), + ("E", 1e18), + ("P", 1e15), + ("T", 1e12), + ("G", 1e9), + ("M", 1e6), + ("k", 1e3), + ("h", 1e2), + ("e", 1e1), + ("d", 1e-1), + ("c", 1e-2), + ("m", 1e-3), + ("u", 1e-6), + ("n", 1e-9), + ("p", 1e-12), + ("f", 1e-15), + ("a", 1e-18), + ("z", 1e-21), + ("y", 1e-24), +]; + +const BIN_PREFIXES: &[(&str, f64)] = &[ + ("ki", 1024.0), + ("Mi", 1_048_576.0), + ("Gi", 1_073_741_824.0), + ("Ti", 1_099_511_627_776.0), + ("Pi", 1_125_899_906_842_624.0), + ("Ei", 1_152_921_504_606_846_976.0), + ("Zi", 1_180_591_620_717_411_303_424.0), + ("Yi", 1_208_925_819_614_629_174_706_176.0), +]; + +fn find_unit(name: &str) -> Option<&'static Unit> { + UNITS.iter().find(|u| u.names.contains(&name)) +} + +/// Resolves a unit name to (unit, prefix multiplier). Exact names win over prefixed ones. +fn resolve(name: &str) -> Option<(&'static Unit, f64)> { + let name = name.replace("^2", "2").replace("^3", "3"); + if let Some(u) = find_unit(&name) { + return Some((u, 1.0)); + } + for (p, f) in BIN_PREFIXES { + if let Some(rest) = name.strip_prefix(p) + && let Some(u) = find_unit(rest) + && u.pre == Pre::Bin + { + return Some((u, *f)); + } + } + for (p, f) in SI_PREFIXES { + if let Some(rest) = name.strip_prefix(p) + && let Some(u) = find_unit(rest) + { + match u.pre { + Pre::Si(pow) => return Some((u, f.powi(pow))), + Pre::Bin => return Some((u, *f)), + Pre::No => {} + } + } + } + None +} + +fn to_kelvin(name: &str, x: f64) -> f64 { + match name { + "C" | "cel" => x + 273.15, + "F" | "fah" => (x - 32.0) * 5.0 / 9.0 + 273.15, + "Rank" => x * 5.0 / 9.0, + "Reau" => x * 1.25 + 273.15, + _ => x, + } +} + +fn from_kelvin(name: &str, k: f64) -> f64 { + match name { + "C" | "cel" => k - 273.15, + "F" | "fah" => (k - 273.15) * 9.0 / 5.0 + 32.0, + "Rank" => k * 9.0 / 5.0, + "Reau" => (k - 273.15) / 1.25, + _ => k, + } +} + +pub(super) fn convert(a: &[Arg], _c: &mut dyn Ctx) -> R { + let x = num(a, 0)?; + let from = text(a, 1)?; + let to = text(a, 2)?; + let (fu, fp) = resolve(&from).ok_or(CellError::NA)?; + let (tu, tp) = resolve(&to).ok_or(CellError::NA)?; + if fu.dim != tu.dim { + return Err(CellError::NA); + } + if fu.dim == Dim::Temperature { + let fname = fu.names.first().copied().unwrap_or("K"); + let tname = tu.names.first().copied().unwrap_or("K"); + if fname == tname && fp == tp { + return num_val(x); + } + let k = to_kelvin(fname, x * fp); + return num_val(from_kelvin(tname, k) / tp); + } + let base = x * fu.factor * fp; + num_val(base / (tu.factor * tp)) +}