Complex number support

examples/complex_impedance.nbt

@@ -0,0 +1,114 @@
+# Complex Impedance in AC Circuits
+#
+# In alternating-current (AC) circuits, resistors, capacitors, and inductors
+# each oppose current flow differently. Their combined effect is described by
+# a complex-valued impedance Z, measured in ohms:
+#
+#   Resistor:  Z_R = R            (purely real)
+#   Capacitor: Z_C = 1 / (i·ω·C) (purely imaginary, negative)
+#   Inductor:  Z_L = i·ω·L       (purely imaginary, positive)
+#
+# where ω = 2π·f is the angular frequency. Complex arithmetic lets us combine
+# these naturally — no separate magnitude/phase bookkeeping needed.
+#
+# This example analyses a series RLC band-pass filter and an RC low-pass filter.
+
+# --- Component values ---
+
+let R₁ = 100 Ω
+let C₁ = 10 nF
+let L₁ = 1 mH
+
+# Resonant frequency: f₀ = 1 / (2π √(LC))
+let f_resonant = 1 / (2π × sqrt(L₁ × C₁))
+
+print("=========================================")
+print("  Series RLC Band-Pass Filter Analysis")
+print("=========================================")
+print()
+print("Components: R = {R₁}, C = {C₁}, L = {L₁}")
+print("Resonant frequency f₀ = {f_resonant -> kHz}")
+print()
+
+# --- Impedance as a function of frequency ---
+# Z = R + iωL + 1/(iωC)
+
+fn Z_total(f: Frequency) -> ElectricResistance =
+    R₁ + i × (2π × f) × L₁ + 1 / (i × (2π × f) × C₁)
+
+# --- Below resonance (capacitive regime) ---
+
+let f_low = f_resonant / 10
+let Z_low = Z_total(f_low)
+
+print("  f = f₀/10 = {f_low -> kHz}  (below resonance)")
+print("    Z   = {Z_low}")
+print("    |Z| = {abs(Z_low)}")
+print("    arg = {arg(Z_low / Ω) -> °}")
+print()
+
+# --- At resonance (purely resistive) ---
+
+let Z_res = Z_total(f_resonant)
+
+print("  f = f₀   = {f_resonant -> kHz}  (at resonance)")
+print("    Z   = {Z_res}")
+print("    |Z| = {abs(Z_res)}")
+print("    arg = {arg(Z_res / Ω) -> °}")
+print()
+
+# At resonance the reactive parts cancel, leaving only the resistance R.
+assert_eq(re(Z_res), R₁, 1 mΩ)
+assert_eq(abs(im(Z_res)), 0 Ω, 1 mΩ)
+
+# --- Above resonance (inductive regime) ---
+
+let f_high = f_resonant × 10
+let Z_high = Z_total(f_high)
+
+print("  f = 10·f₀ = {f_high -> kHz}  (above resonance)")
+print("    Z   = {Z_high}")
+print("    |Z| = {abs(Z_high)}")
+print("    arg = {arg(Z_high / Ω) -> °}")
+print()
+
+# --- RC Low-Pass Filter ---
+#
+# A resistor and capacitor in series form a low-pass filter. The transfer
+# function (voltage divider) is:
+#
+#   H(f) = Z_C / (Z_R + Z_C) = 1 / (1 + i·2π·f·R·C)
+#
+# The gain |H(f)| rolls off above the cutoff frequency f_c = 1/(2π·RC).
+
+print("=========================================")
+print("  RC Low-Pass Filter")
+print("=========================================")
+print()
+
+let R_lp = 1 kΩ
+let C_lp = 100 nF
+let f_cutoff = 1 / (2π × R_lp × C_lp)
+
+print("Components: R = {R_lp}, C = {C_lp}")
+print("Cutoff frequency f_c = {f_cutoff -> kHz}")
+print()
+
+# Complex transfer function
+fn transfer(f: Frequency) -> Scalar = 1 / (1 + i × 2π × f × R_lp × C_lp)
+
+# Gain = |H(f)|
+fn gain(f: Frequency) -> Scalar = abs(transfer(f))
+
+# Phase = arg(transfer(f))
+fn phase(f: Frequency) -> Scalar = arg(transfer(f))
+
+print("  f = f_c/10 :  gain = {gain(f_cutoff / 10)},  phase = {phase(f_cutoff / 10) -> °}")
+print("  f = f_c    :  gain = {gain(f_cutoff)},  phase = {phase(f_cutoff) -> °}")
+print("  f = 10·f_c :  gain = {gain(f_cutoff × 10)},  phase = {phase(f_cutoff × 10) -> °}")
+print()
+
+# At the cutoff frequency the gain is 1/√2 ≈ 0.707 (the −3 dB point)
+# and the phase is −45°.
+assert_eq(gain(f_cutoff), 1 / sqrt(2), 0.001)
+assert_eq(phase(f_cutoff), -45°, 0.01°)

numbat/modules/core/functions.nbt

@@ -5,6 +5,12 @@ use core::scalar
 @example("id(8 kg)")
 fn id<A>(x: A) -> A = x
 
+@name("Is real")
+@description("Returns true if the input is a real number (has no imaginary part).")
+@example("is_real(3)")
+@example("is_real(i)")
+fn is_real<T: Dim>(n: T) -> Bool
+
 @name("Absolute value")
 @description("Return the absolute value $|x|$ of the input. This works for quantities, too: `abs(-5 m) = 5 m`.")
 @url("https://doc.rust-lang.org/std/primitive.f64.html#method.abs")
@@ -21,7 +27,11 @@ fn sqrt<D: Dim>(x: D^2) -> D = x^(1/2)
 @description("Return the cube root $\\sqrt[3]{{x}}$ of the input: `cbrt(8 m^3) = 2 m`.")
 @url("https://en.wikipedia.org/wiki/Cube_root")
 @example("cbrt(8 L) -> cm")
-fn cbrt<D: Dim>(x: D^3) -> D = if x > 0 then x^(1/3) else - (-x)^(1/3)
+fn cbrt<D: Dim>(x: D^3) -> D =
+    if is_real(x) then
+        (if x >= 0 then x^(1/3) else -((-x)^(1/3)))
+    else
+        x^(1/3)
 
 @name("Square function")
 @description("Return the square of the input, $x^2$: `sqr(5 m) = 25 m^2`.")
@@ -105,6 +115,29 @@ fn mod<T: Dim>(a: T, b: T) -> T
 @example("parse(\"0xFF\") -> bin")
 fn parse<T: Dim>(input: String) -> T
 
+@name("Real part")
+@description("Return the real part of the input. For quantities with units, the unit is preserved.")
+@example("re(3 + 2i)")
+@example("re(5 m)")
+fn re<T: Dim>(z: T) -> T
+
+@name("Imaginary part")
+@description("Return the imaginary part of the input. For quantities with units, the unit is preserved.")
+@example("im(3 + 2i)")
+@example("im(5 m)")
+fn im<T: Dim>(z: T) -> T
+
+@name("Complex conjugate")
+@description("Return the complex conjugate of the input. For quantities with units, the unit is preserved.")
+@example("conj(3 + 2i)")
+fn conj<T: Dim>(z: T) -> T
+
+@name("Phase angle")
+@description("Return the phase angle (argument) of the complex number, in radians.")
+@example("arg(i)")
+@example("arg(-1)")
+fn arg(z: Scalar) -> Scalar
+
 @name("Command-line arguments")
 @description("Returns the command-line arguments passed to the script. The first argument is the name of the script itself.")
 @example("let xs = tail(args())", "Get a list of all arguments except the script name.")

numbat/modules/core/numbers.nbt

@@ -31,4 +31,4 @@ fn is_nonzero<D: Dim>(value: D) -> Bool = !is_zero(value)
 @description("Returns true if the input is an integer.")
 @example("is_integer(3)")
 @example("is_integer(pi)")
-fn is_integer(x: Scalar) -> Bool = is_zero(fract(x))
+fn is_integer(x: Scalar) -> Bool = is_real(x) && is_zero(fract(x))

numbat/modules/math/constants.nbt

@@ -1,4 +1,5 @@
 use core::scalar
+use core::functions
 
 ### Mathematical
 
@@ -21,6 +22,10 @@ let e = 2.71828182845904523536028747135266250
 @aliases(golden_ratio)
 let φ = 1.61803398874989484820458683436563811
 
+@name("Imaginary unit")
+@url("https://en.wikipedia.org/wiki/Imaginary_unit")
+let i = sqrt(-1)
+
 ### Named numbers
 
 #### Large numbers

numbat/src/ffi/datetime.rs

@@ -97,7 +97,10 @@ pub fn from_unixtime_us(
     mut args: Args,
     _return_type: &TypeScheme,
 ) -> Result<Value, Box<RuntimeErrorKind>> {
-    let us = quantity_arg!(args).unsafe_value().to_f64() as i64;
+    let q = quantity_arg!(args);
+    let us = q.unsafe_value().try_as_real().ok_or_else(|| {
+        RuntimeErrorKind::ExpectedRealNumberInFunction("from_unixtime".into())
+    })? as i64;
 
     let dt = Timestamp::from_microsecond(us)
         .map_err(|_| RuntimeErrorKind::DateTimeOutOfRange)?
@@ -112,7 +115,10 @@ fn calendar_add(
     to_span: fn(i64) -> std::result::Result<Span, jiff::Error>,
 ) -> Result<Value, Box<RuntimeErrorKind>> {
     let dt = datetime_arg!(args);
-    let n = quantity_arg!(args).unsafe_value().to_f64();
+    let q = quantity_arg!(args);
+    let n = q.unsafe_value().try_as_real().ok_or_else(|| {
+        RuntimeErrorKind::ExpectedRealNumberInFunction(format!("calendar_add ({unit_name})"))
+    })?;
 
     if n.fract() != 0.0 {
         return Err(Box::new(RuntimeErrorKind::UserError(format!(

numbat/src/ffi/functions.rs

@@ -83,6 +83,12 @@ pub(crate) fn functions() -> &'static HashMap<&'static str, ForeignFunction> {
 
         insert_function!(random, 0..=0);
 
+        insert_function!("re", re_fn, 1..=1);
+        insert_function!("im", im_fn, 1..=1);
+        insert_function!("conj", conj, 1..=1);
+        insert_function!("arg", arg_fn, 1..=1);
+        insert_function!(is_real, 1..=1);
+
         // Lists
         insert_function!(len, 1..=1);
         insert_function!(head, 1..=1);
@@ -166,7 +172,12 @@ fn value_of(
 ) -> Result<Value, Box<RuntimeErrorKind>> {
     let quantity = quantity_arg!(args);
 
-    return_scalar!(quantity.unsafe_value().to_f64())
+    let n = quantity.unsafe_value().try_as_real().ok_or_else(|| {
+        Box::new(RuntimeErrorKind::ExpectedRealNumberInFunction(
+            "value_of".into(),
+        ))
+    })?;
+    return_scalar!(n)
 }
 
 fn base_unit_of(

numbat/src/ffi/macros.rs

@@ -56,6 +56,16 @@ macro_rules! return_quantity {
 }
 pub(crate) use return_quantity;
 
+macro_rules! return_complex_scalar {
+    ($re:expr, $im:expr) => {
+        Ok(Value::Quantity(Quantity::new(
+            crate::number::Number::new($re, $im),
+            crate::unit::Unit::scalar(),
+        )))
+    };
+}
+pub(crate) use return_complex_scalar;
+
 macro_rules! return_boolean {
     ($value:expr) => {
         Ok(Value::Boolean($value))