feat: symbolic amplitude and probability formatting

5 files changed, 162 insertions, 141 deletions

diff --git a/include/maths/format.h b/include/maths/format.h
new file mode 100644
index 0000000..65c1680
--- /dev/null
+++ b/include/maths/format.h
@@ -0,0 +1,8 @@
+#pragma once
+
+#include <stddef.h>
+
+#include "maths/complex.h"
+
+char *psi_format_amplitude(struct PsiComplex c, char *out, size_t cap);
+char *psi_format_probability(double p, char *out, size_t cap);
diff --git a/include/psi.h b/include/psi.h
index 42797e3..45cf8a1 100644
--- a/include/psi.h
+++ b/include/psi.h
@@ -9,6 +9,7 @@ const char *psi_version(void);
 #include "maths/complex.h"
 #include "maths/vector.h"
 #include "maths/matrix.h"
+#include "maths/format.h"
 
 #include "core/quantum_components.h"
 #include "core/gates.h"
diff --git a/src/maths/format.c b/src/maths/format.c
new file mode 100644
index 0000000..ad80ff0
--- /dev/null
+++ b/src/maths/format.c
@@ -0,0 +1,147 @@
+#include "maths/format.h"
+
+#include <math.h>
+#include <stdbool.h>
+#include <stdio.h>
+
+static const double EPSILON = 1e-10;
+static const double SQRT_2 = 1.4142135623730951;
+static const double INV_SQRT_2 = 0.7071067811865475;
+static const double INV_SQRT_8 = 0.3535533905932738;
+static const double INV_SQRT_32 = 0.1767766952966369;
+
+static bool approx_eq(double a, double b)
+{
+	return fabs(a - b) < EPSILON;
+}
+
+static bool real_symbolic(double v, char *out, size_t cap)
+{
+	double abs_v = fabs(v);
+	const char *sign = v < 0.0 ? "-" : "";
+
+	if (approx_eq(abs_v, 0.0))
+	{
+		snprintf(out, cap, "%s", "0");
+		return true;
+	}
+
+	const char *sym = NULL;
+	if (approx_eq(abs_v, 1.0))
+		sym = "1";
+	else if (approx_eq(abs_v, 0.5))
+		sym = "½";
+	else if (approx_eq(abs_v, 0.25))
+		sym = "¼";
+	else if (approx_eq(abs_v, 0.75))
+		sym = "¾";
+	else if (approx_eq(abs_v, 0.125))
+		sym = "⅛";
+	else if (approx_eq(abs_v, SQRT_2))
+		sym = "√2";
+	else if (approx_eq(abs_v, INV_SQRT_2))
+		sym = "¹⁄√2";
+	else if (approx_eq(abs_v, INV_SQRT_8))
+		sym = "¹⁄√8";
+	else if (approx_eq(abs_v, INV_SQRT_32))
+		sym = "¹⁄√32";
+	else if (approx_eq(abs_v, 2.0))
+		sym = "2";
+	else if (approx_eq(abs_v, 1.0 / 3.0))
+		sym = "⅓";
+	else if (approx_eq(abs_v, 2.0 / 3.0))
+		sym = "⅔";
+	else
+		return false;
+
+	snprintf(out, cap, "%s%s", sign, sym);
+	return true;
+}
+
+char *psi_format_amplitude(struct PsiComplex c, char *out, size_t cap)
+{
+	double re = c.real;
+	double im = c.imaginary;
+
+	bool re_zero = approx_eq(fabs(re), 0.0);
+	bool im_zero = approx_eq(fabs(im), 0.0);
+
+	if (re_zero && im_zero)
+	{
+		snprintf(out, cap, "%s", "0");
+		return out;
+	}
+
+	if (im_zero)
+	{
+		if (!real_symbolic(re, out, cap))
+			snprintf(out, cap, "%.4f", re);
+
+		return out;
+	}
+
+	if (re_zero)
+	{
+		if (approx_eq(fabs(im), 1.0))
+		{
+			snprintf(out, cap, "%s", im > 0.0 ? "i" : "-i");
+			return out;
+		}
+
+		char sym[32];
+		if (real_symbolic(im, sym, sizeof sym))
+			snprintf(out, cap, "%si", sym);
+		else
+			snprintf(out, cap, "%.4fi", im);
+
+		return out;
+	}
+
+	char re_str[32];
+	if (!real_symbolic(re, re_str, sizeof re_str))
+		snprintf(re_str, sizeof re_str, "%.4f", re);
+
+	char im_str[32];
+	if (approx_eq(fabs(im), 1.0))
+	{
+		snprintf(im_str, sizeof im_str, "%s", im > 0.0 ? "+i" : "-i");
+	}
+	else
+	{
+		const char *sign = im > 0.0 ? "+" : "-";
+		char sym[32];
+		if (real_symbolic(fabs(im), sym, sizeof sym))
+			snprintf(im_str, sizeof im_str, "%s%si", sign, sym);
+		else
+			snprintf(im_str, sizeof im_str, "%s%.4fi", sign, fabs(im));
+	}
+
+	snprintf(out, cap, "%s%s", re_str, im_str);
+	return out;
+}
+
+char *psi_format_probability(double p, char *out, size_t cap)
+{
+	if (approx_eq(p, 0.0))
+		snprintf(out, cap, "%s", "0");
+	else if (approx_eq(p, 1.0))
+		snprintf(out, cap, "%s", "1");
+	else if (approx_eq(p, 0.5))
+		snprintf(out, cap, "%s", "½");
+	else if (approx_eq(p, 0.25))
+		snprintf(out, cap, "%s", "¼");
+	else if (approx_eq(p, 0.75))
+		snprintf(out, cap, "%s", "¾");
+	else if (approx_eq(p, 0.125))
+		snprintf(out, cap, "%s", "⅛");
+	else if (approx_eq(p, 0.0625))
+		snprintf(out, cap, "%s", "¹⁄₁₆");
+	else if (approx_eq(p, 1.0 / 3.0))
+		snprintf(out, cap, "%s", "⅓");
+	else if (approx_eq(p, 2.0 / 3.0))
+		snprintf(out, cap, "%s", "⅔");
+	else
+		snprintf(out, cap, "%.4f", p);
+
+	return out;
+}
diff --git a/src/maths/format.rs b/src/maths/format.rs
deleted file mode 100644
index b6b0191..0000000
--- a/src/maths/format.rs
+++ /dev/null
@@ -1,141 +0,0 @@
-use crate::Complex;
-
-const EPSILON: f64 = 1e-10;
-const SQRT_2: f64 = std::f64::consts::SQRT_2;
-const INV_SQRT_2: f64 = 0.7071067811865475;
-const INV_SQRT_8: f64 = 0.3535533905932738;
-const INV_SQRT_32: f64 = 0.1767766952966369;
-
-fn approx_eq(a: f64, b: f64) -> bool {
-    (a - b).abs() < EPSILON
-}
-
-fn format_real_symbolic(v: f64) -> Option<String> {
-    let abs_v = v.abs();
-    let sign = if v < 0.0 { "-" } else { "" };
-
-    if approx_eq(abs_v, 0.0) {
-        return Some("0".to_string());
-    }
-    if approx_eq(abs_v, 1.0) {
-        return Some(format!("{}1", sign));
-    }
-    if approx_eq(abs_v, 0.5) {
-        return Some(format!("{}½", sign));
-    }
-    if approx_eq(abs_v, 0.25) {
-        return Some(format!("{}¼", sign));
-    }
-    if approx_eq(abs_v, 0.75) {
-        return Some(format!("{}¾", sign));
-    }
-    if approx_eq(abs_v, 0.125) {
-        return Some(format!("{}⅛", sign));
-    }
-    if approx_eq(abs_v, SQRT_2) {
-        return Some(format!("{}√2", sign));
-    }
-    if approx_eq(abs_v, INV_SQRT_2) {
-        return Some(format!("{}¹⁄√2", sign));
-    }
-    if approx_eq(abs_v, INV_SQRT_8) {
-        return Some(format!("{}¹⁄√8", sign));
-    }
-    if approx_eq(abs_v, INV_SQRT_32) {
-        return Some(format!("{}¹⁄√32", sign));
-    }
-    if approx_eq(abs_v, 2.0) {
-        return Some(format!("{}2", sign));
-    }
-    if approx_eq(abs_v, 1.0 / 3.0) {
-        return Some(format!("{}⅓", sign));
-    }
-    if approx_eq(abs_v, 2.0 / 3.0) {
-        return Some(format!("{}⅔", sign));
-    }
-
-    None
-}
-
-pub fn format_amplitude(c: &Complex<f64>) -> String {
-    let re = c.real;
-    let im = c.imaginary;
-
-    let re_zero = approx_eq(re.abs(), 0.0);
-    let im_zero = approx_eq(im.abs(), 0.0);
-
-    if re_zero && im_zero {
-        return "0".to_string();
-    }
-
-    if im_zero {
-        if let Some(s) = format_real_symbolic(re) {
-            return s;
-        }
-        return format!("{:.4}", re);
-    }
-
-    if re_zero {
-        if approx_eq(im.abs(), 1.0) {
-            return if im > 0.0 {
-                "i".to_string()
-            } else {
-                "-i".to_string()
-            };
-        }
-        if let Some(s) = format_real_symbolic(im) {
-            return format!("{}i", s);
-        }
-        return format!("{:.4}i", im);
-    }
-
-    let re_str = format_real_symbolic(re).unwrap_or_else(|| format!("{:.4}", re));
-    let im_str = if approx_eq(im.abs(), 1.0) {
-        if im > 0.0 {
-            "+i".to_string()
-        } else {
-            "-i".to_string()
-        }
-    } else {
-        let im_sym = format_real_symbolic(im.abs());
-        let sign = if im > 0.0 { "+" } else { "-" };
-        match im_sym {
-            Some(s) => format!("{}{}i", sign, s.trim_start_matches('-')),
-            None => format!("{}{:.4}i", sign, im.abs()),
-        }
-    };
-
-    format!("{}{}", re_str, im_str)
-}
-
-pub fn format_probability(p: f64) -> String {
-    if approx_eq(p, 0.0) {
-        return "0".to_string();
-    }
-    if approx_eq(p, 1.0) {
-        return "1".to_string();
-    }
-    if approx_eq(p, 0.5) {
-        return "½".to_string();
-    }
-    if approx_eq(p, 0.25) {
-        return "¼".to_string();
-    }
-    if approx_eq(p, 0.75) {
-        return "¾".to_string();
-    }
-    if approx_eq(p, 0.125) {
-        return "⅛".to_string();
-    }
-    if approx_eq(p, 0.0625) {
-        return "¹⁄₁₆".to_string();
-    }
-    if approx_eq(p, 1.0 / 3.0) {
-        return "⅓".to_string();
-    }
-    if approx_eq(p, 2.0 / 3.0) {
-        return "⅔".to_string();
-    }
-
-    format!("{:.4}", p)
-}
diff --git a/tester/main.c b/tester/main.c
index c8077e9..346b16a 100644
--- a/tester/main.c
+++ b/tester/main.c
@@ -208,6 +208,12 @@ int main(void)
 	psi_circuit_probabilities(&ghz, probs);
 	printf("p(000)  = %g, p(111) = %g\n", probs[0], probs[7]);
 
+	char abuf[64];
+	char pbuf[64];
+	printf("amp     = |000> %s, p|111> %s\n",
+		psi_format_amplitude(ghz_state->data[0], abuf, sizeof abuf),
+		psi_format_probability(probs[7], pbuf, sizeof pbuf));
+
 	psi_free_quantum_circuit(&ghz);
 }