mirror of
https://github.com/storytold/gridcraft.git
synced 2026-10-27 12:16:54 -06:00
refactor(functions): move engineering.rs items to engineering/
This commit is contained in:
aiongg
@@ -197,283 +197,6 @@ fn bitrshift(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
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<Cx> {
|
||||
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<Cx> {
|
||||
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<Cx> {
|
||||
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<f64> {
|
||||
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::<f64>().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<char>)> {
|
||||
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<char>)> {
|
||||
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<char>)> {
|
||||
cx_of(&a.get(i).map(scalar).unwrap_or(Value::Empty))
|
||||
}
|
||||
|
||||
fn merge_suffix(acc: Option<char>, s: Option<char>) -> R<Option<char>> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
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<Cx>) -> R<Value> {
|
||||
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<f64>) -> R<Value> {
|
||||
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<Value> {
|
||||
@@ -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<Value> {
|
||||
im_real_fn(a, |z| Ok(z.abs()))
|
||||
}
|
||||
fn imreal(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
im_real_fn(a, |z| Ok(z.re))
|
||||
}
|
||||
fn imaginary(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
im_real_fn(a, |z| Ok(z.im))
|
||||
}
|
||||
fn imargument(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
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<Cx>) -> R<Value> {
|
||||
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<Value> {
|
||||
im_binary(a, Cx::div)
|
||||
}
|
||||
fn imsub(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
im_binary(a, |x, y| Ok(x.sub(y)))
|
||||
}
|
||||
fn impower(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
im_fold(a, Cx::new(0.0, 0.0), Cx::add)
|
||||
}
|
||||
fn improduct(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
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<f64> {
|
||||
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<f64> {
|
||||
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::<f64>();
|
||||
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<f64> {
|
||||
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<f64> {
|
||||
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<f64> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
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<FnSpec> {
|
||||
vec![
|
||||
f!("BIN2DEC", 1, 1, Engineering, S, "BIN2DEC(number)", "Converts a binary number to decimal.", bin2dec),
|
||||
|
||||
@@ -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<f64> {
|
||||
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<f64> {
|
||||
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::<f64>();
|
||||
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<f64> {
|
||||
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<f64> {
|
||||
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<f64> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
let (x, n) = bessel_args(a)?;
|
||||
if x <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(bessel_k(n, x).ok_or(CellError::Num)?)
|
||||
}
|
||||
@@ -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<Cx> {
|
||||
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<Cx> {
|
||||
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<Cx> {
|
||||
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<f64> {
|
||||
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::<f64>().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<char>)> {
|
||||
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<char>)> {
|
||||
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<char>)> {
|
||||
cx_of(&a.get(i).map(scalar).unwrap_or(Value::Empty))
|
||||
}
|
||||
|
||||
fn merge_suffix(acc: Option<char>, s: Option<char>) -> R<Option<char>> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
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<Cx>) -> R<Value> {
|
||||
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<f64>) -> R<Value> {
|
||||
let (z, _) = cx_arg(a, 0)?;
|
||||
num_val(f(z)?)
|
||||
}
|
||||
|
||||
pub(super) fn imabs(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
im_real_fn(a, |z| Ok(z.abs()))
|
||||
}
|
||||
|
||||
pub(super) fn imreal(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
im_real_fn(a, |z| Ok(z.re))
|
||||
}
|
||||
|
||||
pub(super) fn imaginary(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
im_real_fn(a, |z| Ok(z.im))
|
||||
}
|
||||
|
||||
pub(super) fn imargument(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
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<Cx>) -> R<Value> {
|
||||
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<Value> {
|
||||
im_binary(a, Cx::div)
|
||||
}
|
||||
|
||||
pub(super) fn imsub(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
im_binary(a, |x, y| Ok(x.sub(y)))
|
||||
}
|
||||
|
||||
pub(super) fn impower(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
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<Value> {
|
||||
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<Value> {
|
||||
im_fold(a, Cx::new(0.0, 0.0), Cx::add)
|
||||
}
|
||||
|
||||
pub(super) fn improduct(a: &[Arg], _c: &mut dyn Ctx) -> R<Value> {
|
||||
im_fold(a, ONE, Cx::mul)
|
||||
}
|
||||
@@ -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<Value> {
|
||||
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))
|
||||
}
|
||||
Reference in New Issue
Block a user