[Nickle] nickle: Branch 'master' - 11 commits
Keith Packard
keithp at keithp.com
Thu Oct 9 13:07:34 PDT 2025
builtin-complex.c | 57 ++++++++++
builtin-namespaces.h | 1
builtin-toplevel.c | 38 ++++--
builtin.5c | 2
builtin.c | 2
complex.5c | 176 +++++++++++++++++++++++++++++++
complex.c | 211 +++++++++++++++++++++++++++++++++++++
debian/changelog | 12 ++
file.c | 3
float.c | 4
gamma.5c | 33 +++--
gram.y | 3
lex.l | 1
math-ext.5c | 162 ++++++++++++++++++++++++++++
math.5c | 125 ++++++++++++++--------
meson.build | 9 -
nickle.h | 1
test/cbrttest.5c | 35 ++++++
test/etest.5c | 48 ++++++++
test/exptest.5c | 62 +++++++++++
test/gammatest.5c | 44 +++++++
test/jntest.5c | 24 ++++
test/logtest.5c | 66 +++++++++++
test/math-bits.5c | 86 +++++++++++++++
test/math-tables.sh | 25 ++++
test/math.5c | 286 ---------------------------------------------------
test/meson.build | 10 +
test/sincostest.5c | 29 +++++
test/sqrttest.5c | 34 ++++++
test/tantest.5c | 32 +++++
type.c | 20 +--
value.c | 21 +++
value.h | 44 ++++---
33 files changed, 1319 insertions(+), 387 deletions(-)
New commits:
commit 4d44159c2cef9e20e20fcc7451e31ed307c5a547
Author: Keith Packard <keithp at keithp.com>
Date: Fri Oct 3 13:46:21 2025 -0700
Update to version 2.105
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/debian/changelog b/debian/changelog
index 2312df0..c840dd1 100644
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,15 @@
+nickle (2.105) unstable; urgency=medium
+
+ * Add more math functions
+ * Improve performance of log for large values.
+ * Improve performance of gamma for negative values.
+ * Test gamma using Legendre duplication formula.
+ * Improve accuracy and avoid infinite loop in sqrt and cbrt.
+ * Test sqrt and cbrt.
+ * Add 'complex' type.
+
+ -- Keith Packard <keithp at keithp.com> Sun, 05 Oct 2025 13:10:14 -0700
+
nickle (2.104) unstable; urgency=medium
* Replace #include SYMBOL with conditional includes.
diff --git a/meson.build b/meson.build
index fae4c36..a43ffcf 100644
--- a/meson.build
+++ b/meson.build
@@ -12,10 +12,10 @@ project('nickle', 'c',
],
license : 'BSD',
meson_version : '>= 1.0',
- version: '2.103'
+ version: '2.105'
)
-release_date = '2025-04-28'
+release_date = '2025-10-05'
prefix = get_option('prefix')
include = prefix / get_option('includedir')
commit ed6a50586140a438e87eaf408179abcf90e4cd97
Author: Keith Packard <keithp at keithp.com>
Date: Thu Oct 9 12:48:02 2025 -0700
test: Split math tests into separate files
These can take a while to run, so do them in parallel
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/test/cbrttest.5c b/test/cbrttest.5c
new file mode 100644
index 0000000..35387f4
--- /dev/null
+++ b/test/cbrttest.5c
@@ -0,0 +1,35 @@
+/*
+ * Nickle test suite
+ *
+ * cbrt test
+ */
+library "math-bits.5c"
+
+void check_cbrt()
+{
+ void check_value (real z, int prec) {
+ z = imprecise(z, prec);
+ real s = cbrt(z);
+ real err = imprecise(z, prec+10) - imprecise(s, prec+10)**3;
+ int errdist = exponent(z) - exponent(err);
+ if (err != 0 && errdist < prec - 2) {
+ printf("%g = cbrt(%g): error too large (differ at %d bits instead of %d bits)\n",
+ s, z, errdist, prec);
+ errors++;
+ }
+ }
+
+ printf("Checking cbrt...");
+ for (int prec = 10; prec <= 1000; prec *= 10) {
+ printf(" %d", prec);
+ File::flush(stdout);
+ for (real z = 1e-200; z < 1e200; z *= 1.3) {
+ check_value(z, prec);
+ check_value(-z, prec);
+ }
+ }
+ printf("\n");
+}
+
+check_cbrt();
+exit(errors);
diff --git a/test/etest.5c b/test/etest.5c
new file mode 100644
index 0000000..a226fb0
--- /dev/null
+++ b/test/etest.5c
@@ -0,0 +1,48 @@
+/*
+ * Nickle test suite
+ *
+ * test value of the 'e' constant.
+ */
+library "math-bits.5c"
+
+/* 1000 digits of e */
+string e_digits = ("2." +
+ "71828182845904523536028747135266249775724709369995" +
+ "95749669676277240766303535475945713821785251664274" +
+ "27466391932003059921817413596629043572900334295260" +
+ "59563073813232862794349076323382988075319525101901" +
+ "15738341879307021540891499348841675092447614606680" +
+ "82264800168477411853742345442437107539077744992069" +
+ "55170276183860626133138458300075204493382656029760" +
+ "67371132007093287091274437470472306969772093101416" +
+ "92836819025515108657463772111252389784425056953696" +
+ "77078544996996794686445490598793163688923009879312" +
+ "77361782154249992295763514822082698951936680331825" +
+ "28869398496465105820939239829488793320362509443117" +
+ "30123819706841614039701983767932068328237646480429" +
+ "53118023287825098194558153017567173613320698112509" +
+ "96181881593041690351598888519345807273866738589422" +
+ "87922849989208680582574927961048419844436346324496" +
+ "84875602336248270419786232090021609902353043699418" +
+ "49146314093431738143640546253152096183690888707016" +
+ "76839642437814059271456354906130310720851038375051" +
+ "01157477041718986106873969655212671546889570350353");
+
+void check_e()
+{
+ printf("Checking e\n");
+ File::sscanf(e_digits, "%f", &(real e_good));
+ for (int prec = 10; prec < 3000; prec *= 2) {
+ real e_tmp = Math::e_value(prec);
+ real error = abs (e_tmp - e_good);
+ /* prec-2 bits right of the decimal point */
+ real limit = 2**(-(prec-2));
+ if (error > limit) {
+ printf("e_value(%d) error %g limit %g value %g\n", prec, error, limit, e_tmp);
+ errors++;
+ }
+ }
+}
+
+check_e();
+exit(errors);
diff --git a/test/exptest.5c b/test/exptest.5c
new file mode 100644
index 0000000..b3a8864
--- /dev/null
+++ b/test/exptest.5c
@@ -0,0 +1,62 @@
+/*
+ * Nickle test suite
+ *
+ * exp test
+ */
+library "math-bits.5c"
+
+void check_exp()
+{
+ printf("Checking exp precision\n");
+ for (int prec = 10; prec <= 1000; prec *= 10) {
+ real eval = Math::e_value(1000);
+ for (int x = 0; x < 100; x++) {
+ real expx = exp(imprecise(x, prec));
+ real powx = imprecise(eval ** x, prec);
+ real error = abs((imprecise(expx,prec*10) - imprecise(powx,prec*10)) / imprecise(expx,prec*10));
+ real limit = 2 ** (-(prec - 8));
+
+ if (error > limit) {
+ printf ("exp(i): prec %d error %g limit %g x %.-g val %.-g %.-g\n",
+ prec, error, limit, x, expx, powx);
+ errors++;
+ }
+ }
+ for (int p = 0; p < prec // 3; p++) {
+ real x = imprecise(2 ** -p, prec);
+ if (exp(x) - 1 == 0) {
+ printf("exp(x): prec %d is one for p %d x %.-g\n", prec, p, x);
+ errors++;
+ }
+ }
+ for (real x = imprecise(0.001, prec); x < 20; x *= 1.5) {
+ for (real y = imprecise(0.001, prec); y < 20; y *= 1.5) {
+
+ real expx = exp(x);
+ real expy = exp(y);
+ real expx_expy = exp(x) * exp(y);
+ real expxy = exp(x+y);
+
+ int bitdiff = ceil(abs(log2(abs(x)) - log2(abs(y))));
+
+ real error = abs((expxy - expx_expy) / expxy);
+ real limit = 2 ** (-(prec - 5 - bitdiff));
+
+ if (error > limit) {
+ printf("exp(x+y): prec %d error %g limit %g x %.-g y %.-g\n",
+ prec, error, limit, x, y);
+ errors++;
+ }
+ }
+ }
+ }
+ printf("Checking exp\n");
+ for (int i = 0; i < dim(exp_table); i++) {
+ real in = exp_table[i].in;
+ real exp_value = exp_table[i].exp;
+ check("exp", exp, in, exp_value);
+ }
+}
+
+check_exp();
+exit(errors);
diff --git a/test/gammatest.5c b/test/gammatest.5c
new file mode 100644
index 0000000..0441892
--- /dev/null
+++ b/test/gammatest.5c
@@ -0,0 +1,44 @@
+/*
+ * Nickle test suite
+ *
+ * gamma test
+ */
+library "math-bits.5c"
+
+void check_gamma()
+{
+ printf("Checking gamma using the Legendre duplication formula\n");
+
+ void check_value(real z, int prec) {
+ int iprec = prec + 20 + abs(exponent(imprecise(z)));
+ real zi = imprecise(z, iprec);
+ real left = gamma(zi) * gamma(zi + 1/2);
+ real right = 2**(1-2*zi) * sqrt(pi_value(iprec)) * gamma(zi * 2);
+
+ real error = left - right;
+ if (error != 0 && exponent(error) > exponent(left) - prec) {
+ printf("%g: Γ(z)Γ(z+½) (%g) != 2^(1-2z) √π Γ(2z) (%g): prec %d error %g\n",
+ z, left, right, prec, error);
+ errors++;
+ }
+ }
+
+ for (int prec = 10; prec <= 100; prec *= 10) {
+ for (real z = 1e-38; z < 1e38; z *= 11) {
+ bool invalid = false;
+ check_value(z, prec);
+ try {
+ check_value(-z, prec);
+ } catch invalid_argument(string err, int i, real z) {
+ invalid = true;
+ }
+ if (invalid != (floor(z) == z)) {
+ printf("%g: negative wrong %v\n", z, invalid);
+ errors++;
+ }
+ }
+ }
+}
+
+check_gamma();
+exit(errors);
diff --git a/test/jntest.5c b/test/jntest.5c
new file mode 100644
index 0000000..7c1f0b7
--- /dev/null
+++ b/test/jntest.5c
@@ -0,0 +1,24 @@
+/*
+ * Nickle test suite
+ *
+ * jn test
+ */
+library "math-bits.5c"
+
+void check_jn()
+{
+ real jn_fixed(real n, real x) {
+ return jn(floor(n), x);
+ }
+
+ printf("Checking jn\n");
+ for (int i = 0; i < dim(jn_table); i++) {
+ int n = jn_table[i].n;
+ real in = jn_table[i].in;
+ real jn_value = jn_table[i].jn;
+ check2("jn", jn_fixed, n, in, jn_value);
+ }
+}
+
+check_jn();
+exit(errors);
diff --git a/test/logtest.5c b/test/logtest.5c
new file mode 100644
index 0000000..5af433e
--- /dev/null
+++ b/test/logtest.5c
@@ -0,0 +1,66 @@
+/*
+ * Nickle test suite
+ *
+ * log test
+ */
+library "math-bits.5c"
+
+void check_log()
+{
+ printf("Checking log precision\n");
+ int log_errors = 0;
+ for (int prec = 10; prec <= 1000; prec *= 10) {
+ real eval = Math::e_value(prec);
+ for (real x = imprecise(0.001, prec); x < 20; x *= 1.5) {
+ real lnx = log(x);
+
+ real explnx = exp(lnx);
+
+ real rterror = abs((x - explnx) / x);
+ real rtlimit = 2 ** (-(prec - 6));
+
+ if (rterror > rtlimit) {
+ printf("exp(log(x)) = x: prec %d error %g limit %g x %.-g\n",
+ prec, rterror, rtlimit, x);
+ errors++;
+ }
+
+ for (real y = imprecise(0.001, prec); y < 20; y *= 1.5) {
+ real lny = log(y);
+
+ if (sign (lnx) != sign(lny))
+ continue;
+
+ real lnxy = log(x * y);
+ real lnx_lny = lnx + lny;
+ real error = abs ((lnxy - lnx_lny) / lnxy);
+
+ int plnxy = precision(lnxy);
+ int plnx_lny = precision(lnx_lny);
+
+ if (plnxy != plnx_lny)
+ printf("precision difference %d %d\n",
+ plnxy, plnx_lny);
+
+ int bitdiff = max(2, ceil(abs(log2(abs(lnx)) - log2(abs(lny)))));
+
+ real limit = 2 ** -(prec - bitdiff - 5);
+
+ if (error > limit) {
+ printf("log(x*y) = log(x)+log(y): prec %d error %g limit %g x %.-g y %.-g\n",
+ prec, error, limit, x, y);
+ errors++;
+ }
+ }
+ }
+ }
+ printf("Checking log\n");
+ for (int i = 0; i < dim(log_table); i++) {
+ real in = log_table[i].in;
+ real log_value = log_table[i].log;
+ check("log", log, in, log_value);
+ }
+}
+
+check_log();
+exit(errors);
diff --git a/test/math-bits.5c b/test/math-bits.5c
new file mode 100644
index 0000000..7954955
--- /dev/null
+++ b/test/math-bits.5c
@@ -0,0 +1,86 @@
+/*
+ * Nickle test suite
+ *
+ * Trig operator tests
+ */
+
+load "math-tables.5c"
+
+int test_precision = 333;
+int errors = 0;
+real error_bound = 1e-70;
+
+real checked(real(real) f, real a) {
+ try {
+ real r = f(imprecise (a, test_precision));
+ if (r > 1e10)
+ return 9999;
+ if (r < -1e10)
+ return 9999;
+ return r;
+ } catch invalid_argument(string type, int i, real v) {
+ return 10000;
+ }
+ return 0;
+}
+
+real checked2(real(real, real) f, real a, real b) {
+ try {
+ real r = f(imprecise(a, test_precision), imprecise(b, test_precision));
+
+ /* Just return the same value for large magnitudes
+ * as small variations in input cause the sign to flip
+ */
+
+ if (r > 1e10)
+ return 9999;
+ if (r < -1e10)
+ return 9999;
+ return r;
+ } catch invalid_argument(string type, int i, real v) {
+ return 10000;
+ }
+ return 0;
+}
+
+void check(string op, real(real) f, real input, real correct)
+{
+ real value = checked(f, input);
+
+ if (is_number (value))
+ {
+ real error = abs (value - correct);
+
+ if (abs(correct) >= 1 && error < abs(correct) * error_bound)
+ return;
+ if (abs(correct) < 1 && error < error_bound)
+ return;
+ printf ("error %v value %v\n", error, value);
+ }
+ if (value != correct) {
+ printf ("check failed %s(%v) (was %.-g, should be %.-g)\n",
+ op, input, value, correct);
+ errors++;
+ }
+}
+
+void check2(string op, real(real, real) f, real input1, real input2, real correct)
+{
+ real value = checked2(f, input1, input2);
+
+ if (is_number (value))
+ {
+ real error = abs (value - correct);
+
+ if (abs(correct) >= 1 && error < abs(correct) * error_bound)
+ return;
+ if (abs(correct) < 1 && error < error_bound)
+ return;
+ printf ("error %v value %v\n", error, value);
+ }
+ if (value != correct) {
+ printf ("check failed %s(%.-g, %.-g) (was %.-g, should be %.-g)\n",
+ op, input1, input2, value, correct);
+ errors++;
+ }
+}
diff --git a/test/math.5c b/test/math.5c
deleted file mode 100644
index 0188c05..0000000
--- a/test/math.5c
+++ /dev/null
@@ -1,391 +0,0 @@
-/*
- * Nickle test suite
- *
- * Trig operator tests
- */
-
-load "math-tables.5c"
-
-int test_precision = 333;
-int errors = 0;
-real error_bound = 1e-70;
-
-real checked(real(real) f, real a) {
- try {
- real r = f(imprecise (a, test_precision));
- if (r > 1e10)
- return 9999;
- if (r < -1e10)
- return 9999;
- return r;
- } catch invalid_argument(string type, int i, real v) {
- return 10000;
- }
- return 0;
-}
-
-real checked2(real(real, real) f, real a, real b) {
- try {
- real r = f(imprecise(a, test_precision), imprecise(b, test_precision));
-
- /* Just return the same value for large magnitudes
- * as small variations in input cause the sign to flip
- */
-
- if (r > 1e10)
- return 9999;
- if (r < -1e10)
- return 9999;
- return r;
- } catch invalid_argument(string type, int i, real v) {
- return 10000;
- }
- return 0;
-}
-
-void check(string op, real(real) f, real input, real correct)
-{
- real value = checked(f, input);
-
- if (is_number (value))
- {
- real error = abs (value - correct);
-
- if (abs(correct) >= 1 && error < abs(correct) * error_bound)
- return;
- if (abs(correct) < 1 && error < error_bound)
- return;
- printf ("error %v value %v\n", error, value);
- }
- if (value != correct) {
- printf ("check failed %s(%v) (was %.-g, should be %.-g)\n",
- op, input, value, correct);
- errors++;
- }
-}
-
-void check2(string op, real(real, real) f, real input1, real input2, real correct)
-{
- real value = checked2(f, input1, input2);
-
- if (is_number (value))
- {
- real error = abs (value - correct);
-
- if (abs(correct) >= 1 && error < abs(correct) * error_bound)
- return;
- if (abs(correct) < 1 && error < error_bound)
- return;
- printf ("error %v value %v\n", error, value);
- }
- if (value != correct) {
- printf ("check failed %s(%.-g, %.-g) (was %.-g, should be %.-g)\n",
- op, input1, input2, value, correct);
- errors++;
- }
-}
-
-void check_sqrt()
-{
- void check_value (real z, int prec) {
- z = imprecise(z, prec);
- real s = sqrt(z);
- real err = z - s*s;
- int errdist = exponent(z) - exponent(err);
- if (err != 0 && errdist < prec - 1) {
- printf("%g = sqrt(%g): error too large (differ at %d bits instead of %d bits)\n",
- s, z, errdist, prec);
- errors++;
- }
- }
-
- printf("Checking sqrt...");
- for (int prec = 10; prec <= 1000; prec *= 10) {
- printf(" %d", prec);
- File::flush(stdout);
- for (real z = 1e-200; z < 1e200; z *= 1.3) {
- check_value(z, prec);
- }
- }
- printf("\n");
-}
-
-void check_cbrt()
-{
- void check_value (real z, int prec) {
- z = imprecise(z, prec);
- real s = cbrt(z);
- real err = imprecise(z, prec+10) - imprecise(s, prec+10)**3;
- int errdist = exponent(z) - exponent(err);
- if (err != 0 && errdist < prec - 2) {
- printf("%g = cbrt(%g): error too large (differ at %d bits instead of %d bits)\n",
- s, z, errdist, prec);
- errors++;
- }
- }
-
- printf("Checking cbrt...");
- for (int prec = 10; prec <= 1000; prec *= 10) {
- printf(" %d", prec);
- File::flush(stdout);
- for (real z = 1e-200; z < 1e200; z *= 1.3) {
- check_value(z, prec);
- check_value(-z, prec);
- }
- }
- printf("\n");
-}
-
-void check_sin_cos()
-{
- printf("Checking sin and cos\n");
- for (int i = 0; i < dim(sin_cos_table); i++) {
- real angle = sin_cos_table[i].angle;
- real sin_value = sin_cos_table[i].sin;
- real cos_value = sin_cos_table[i].cos;
- check("sin", sin, angle, sin_value);
- check("cos", cos, angle, cos_value);
- }
- printf("Checking asin and acos\n");
- for (int i = 0; i < dim(asin_acos_table); i++) {
- real ratio = asin_acos_table[i].ratio;
- real asin_value = asin_acos_table[i].asin;
- real acos_value = asin_acos_table[i].acos;
- check("asin", asin, ratio, asin_value);
- check("acos", acos, ratio, acos_value);
- }
-}
-
-void check_tan()
-{
- printf("Checking tan\n");
- for (int i = 0; i < dim(tan_table); i++) {
- real angle = tan_table[i].angle;
- real tan_value = tan_table[i].tan;
- check("tan", tan, angle, tan_value);
- }
- printf("Checking atan\n");
- for (int i = 0; i < dim(atan_table); i++) {
- real ratio = atan_table[i].ratio;
- real atan_value = atan_table[i].atan;
- check("atan", atan, ratio, atan_value);
- }
- printf("Checking atan2\n");
- for (int i = 0; i < dim(atan2_table); i++) {
- real y = atan2_table[i].y;
- real x = atan2_table[i].x;
- real atan2_value = atan2_table[i].atan2;
- check2("atan2", atan2, y, x, atan2_value);
- }
-}
-
-/* 1000 digits of e */
-string e_digits = ("2." +
- "71828182845904523536028747135266249775724709369995" +
- "95749669676277240766303535475945713821785251664274" +
- "27466391932003059921817413596629043572900334295260" +
- "59563073813232862794349076323382988075319525101901" +
- "15738341879307021540891499348841675092447614606680" +
- "82264800168477411853742345442437107539077744992069" +
- "55170276183860626133138458300075204493382656029760" +
- "67371132007093287091274437470472306969772093101416" +
- "92836819025515108657463772111252389784425056953696" +
- "77078544996996794686445490598793163688923009879312" +
- "77361782154249992295763514822082698951936680331825" +
- "28869398496465105820939239829488793320362509443117" +
- "30123819706841614039701983767932068328237646480429" +
- "53118023287825098194558153017567173613320698112509" +
- "96181881593041690351598888519345807273866738589422" +
- "87922849989208680582574927961048419844436346324496" +
- "84875602336248270419786232090021609902353043699418" +
- "49146314093431738143640546253152096183690888707016" +
- "76839642437814059271456354906130310720851038375051" +
- "01157477041718986106873969655212671546889570350353");
-
-void check_e()
-{
- printf("Checking e\n");
- File::sscanf(e_digits, "%f", &(real e_good));
- for (int prec = 10; prec < 3000; prec *= 2) {
- real e_tmp = Math::e_value(prec);
- real error = abs (e_tmp - e_good);
- /* prec-2 bits right of the decimal point */
- real limit = 2**(-(prec-2));
- if (error > limit) {
- printf("e_value(%d) error %g limit %g value %g\n", prec, error, limit, e_tmp);
- errors++;
- }
- }
-}
-
-void check_exp()
-{
- printf("Checking exp precision\n");
- for (int prec = 10; prec <= 1000; prec *= 10) {
- real eval = Math::e_value(1000);
- for (int x = 0; x < 100; x++) {
- real expx = exp(imprecise(x, prec));
- real powx = imprecise(eval ** x, prec);
- real error = abs((imprecise(expx,prec*10) - imprecise(powx,prec*10)) / imprecise(expx,prec*10));
- real limit = 2 ** (-(prec - 8));
-
- if (error > limit) {
- printf ("exp(i): prec %d error %g limit %g x %.-g val %.-g %.-g\n",
- prec, error, limit, x, expx, powx);
- errors++;
- }
- }
- for (int p = 0; p < prec // 3; p++) {
- real x = imprecise(2 ** -p, prec);
- if (exp(x) - 1 == 0) {
- printf("exp(x): prec %d is one for p %d x %.-g\n", prec, p, x);
- errors++;
- }
- }
- for (real x = imprecise(0.001, prec); x < 20; x *= 1.5) {
- for (real y = imprecise(0.001, prec); y < 20; y *= 1.5) {
-
- real expx = exp(x);
- real expy = exp(y);
- real expx_expy = exp(x) * exp(y);
- real expxy = exp(x+y);
-
- int bitdiff = ceil(abs(log2(abs(x)) - log2(abs(y))));
-
- real error = abs((expxy - expx_expy) / expxy);
- real limit = 2 ** (-(prec - 5 - bitdiff));
-
- if (error > limit) {
- printf("exp(x+y): prec %d error %g limit %g x %.-g y %.-g\n",
- prec, error, limit, x, y);
- errors++;
- }
- }
- }
- }
- printf("Checking exp\n");
- for (int i = 0; i < dim(exp_table); i++) {
- real in = exp_table[i].in;
- real exp_value = exp_table[i].exp;
- check("exp", exp, in, exp_value);
- }
-}
-
-void check_log()
-{
- printf("Checking log precision\n");
- int log_errors = 0;
- for (int prec = 10; prec <= 1000; prec *= 10) {
- real eval = Math::e_value(prec);
- for (real x = imprecise(0.001, prec); x < 20; x *= 1.5) {
- real lnx = log(x);
-
- real explnx = exp(lnx);
-
- real rterror = abs((x - explnx) / x);
- real rtlimit = 2 ** (-(prec - 6));
-
- if (rterror > rtlimit) {
- printf("exp(log(x)) = x: prec %d error %g limit %g x %.-g\n",
- prec, rterror, rtlimit, x);
- errors++;
- }
-
- for (real y = imprecise(0.001, prec); y < 20; y *= 1.5) {
- real lny = log(y);
-
- if (sign (lnx) != sign(lny))
- continue;
-
- real lnxy = log(x * y);
- real lnx_lny = lnx + lny;
- real error = abs ((lnxy - lnx_lny) / lnxy);
-
- int plnxy = precision(lnxy);
- int plnx_lny = precision(lnx_lny);
-
- if (plnxy != plnx_lny)
- printf("precision difference %d %d\n",
- plnxy, plnx_lny);
-
- int bitdiff = max(2, ceil(abs(log2(abs(lnx)) - log2(abs(lny)))));
-
- real limit = 2 ** -(prec - bitdiff - 5);
-
- if (error > limit) {
- printf("log(x*y) = log(x)+log(y): prec %d error %g limit %g x %.-g y %.-g\n",
- prec, error, limit, x, y);
- errors++;
- }
- }
- }
- }
- printf("Checking log\n");
- for (int i = 0; i < dim(log_table); i++) {
- real in = log_table[i].in;
- real log_value = log_table[i].log;
- check("log", log, in, log_value);
- }
-}
-
-void check_gamma()
-{
- printf("Checking gamma using the Legendre duplication formula\n");
-
- void check_value(real z, int prec) {
- int iprec = prec + 20 + abs(exponent(imprecise(z)));
- real zi = imprecise(z, iprec);
- real left = gamma(zi) * gamma(zi + 1/2);
- real right = 2**(1-2*zi) * sqrt(pi_value(iprec)) * gamma(zi * 2);
-
- real error = left - right;
- if (error != 0 && exponent(error) > exponent(left) - prec) {
- printf("%g: Γ(z)Γ(z+½) (%g) != 2^(1-2z) √π Γ(2z) (%g): prec %d error %g\n",
- z, left, right, prec, error);
- errors++;
- }
- }
-
- for (int prec = 10; prec <= 100; prec *= 10) {
- for (real z = 1e-38; z < 1e38; z *= 11) {
- bool invalid = false;
- check_value(z, prec);
- try {
- check_value(-z, prec);
- } catch invalid_argument(string err, int i, real z) {
- invalid = true;
- }
- if (invalid != (floor(z) == z)) {
- printf("%g: negative wrong %v\n", z, invalid);
- errors++;
- }
- }
- }
-}
-
-void check_jn()
-{
- real jn_fixed(real n, real x) {
- return jn(floor(n), x);
- }
-
- printf("Checking jn\n");
- for (int i = 0; i < dim(jn_table); i++) {
- int n = jn_table[i].n;
- real in = jn_table[i].in;
- real jn_value = jn_table[i].jn;
- check2("jn", jn_fixed, n, in, jn_value);
- }
-}
-
-check_sqrt();
-check_cbrt();
-check_sin_cos();
-check_tan();
-check_exp();
-check_log();
-check_e();
-check_gamma();
-check_jn();
-
-exit (errors);
diff --git a/test/meson.build b/test/meson.build
index 963c4de..a9f529c 100644
--- a/test/meson.build
+++ b/test/meson.build
@@ -18,7 +18,15 @@ tests = [
'signal.5c',
'round.5c',
'fround.5c',
- 'math.5c',
+ 'sqrttest.5c',
+ 'cbrttest.5c',
+ 'sincostest.5c',
+ 'tantest.5c',
+ 'exptest.5c',
+ 'logtest.5c',
+ 'etest.5c',
+ 'gammatest.5c',
+ 'jntest.5c',
'factorial.5c',
'is_type.5c',
'jsontest.5c',
diff --git a/test/sincostest.5c b/test/sincostest.5c
new file mode 100644
index 0000000..d895ba5
--- /dev/null
+++ b/test/sincostest.5c
@@ -0,0 +1,29 @@
+/*
+ * Nickle test suite
+ *
+ * sin, cos, asin and acos tests test
+ */
+library "math-bits.5c"
+
+void check_sin_cos()
+{
+ printf("Checking sin and cos\n");
+ for (int i = 0; i < dim(sin_cos_table); i++) {
+ real angle = sin_cos_table[i].angle;
+ real sin_value = sin_cos_table[i].sin;
+ real cos_value = sin_cos_table[i].cos;
+ check("sin", sin, angle, sin_value);
+ check("cos", cos, angle, cos_value);
+ }
+ printf("Checking asin and acos\n");
+ for (int i = 0; i < dim(asin_acos_table); i++) {
+ real ratio = asin_acos_table[i].ratio;
+ real asin_value = asin_acos_table[i].asin;
+ real acos_value = asin_acos_table[i].acos;
+ check("asin", asin, ratio, asin_value);
+ check("acos", acos, ratio, acos_value);
+ }
+}
+
+check_sin_cos();
+exit(errors);
diff --git a/test/sqrttest.5c b/test/sqrttest.5c
new file mode 100644
index 0000000..97587d4
--- /dev/null
+++ b/test/sqrttest.5c
@@ -0,0 +1,34 @@
+/*
+ * Nickle test suite
+ *
+ * sqrt test
+ */
+library "math-bits.5c"
+
+void check_sqrt()
+{
+ void check_value (real z, int prec) {
+ z = imprecise(z, prec);
+ real s = sqrt(z);
+ real err = z - s*s;
+ int errdist = exponent(z) - exponent(err);
+ if (err != 0 && errdist < prec - 1) {
+ printf("%g = sqrt(%g): error too large (differ at %d bits instead of %d bits)\n",
+ s, z, errdist, prec);
+ errors++;
+ }
+ }
+
+ printf("Checking sqrt...");
+ for (int prec = 10; prec <= 1000; prec *= 10) {
+ printf(" %d", prec);
+ File::flush(stdout);
+ for (real z = 1e-200; z < 1e200; z *= 1.3) {
+ check_value(z, prec);
+ }
+ }
+ printf("\n");
+}
+
+check_sqrt();
+exit(errors);
diff --git a/test/tantest.5c b/test/tantest.5c
new file mode 100644
index 0000000..10c8906
--- /dev/null
+++ b/test/tantest.5c
@@ -0,0 +1,32 @@
+/*
+ * Nickle test suite
+ *
+ * tan, atan and atan2 test
+ */
+library "math-bits.5c"
+
+void check_tan()
+{
+ printf("Checking tan\n");
+ for (int i = 0; i < dim(tan_table); i++) {
+ real angle = tan_table[i].angle;
+ real tan_value = tan_table[i].tan;
+ check("tan", tan, angle, tan_value);
+ }
+ printf("Checking atan\n");
+ for (int i = 0; i < dim(atan_table); i++) {
+ real ratio = atan_table[i].ratio;
+ real atan_value = atan_table[i].atan;
+ check("atan", atan, ratio, atan_value);
+ }
+ printf("Checking atan2\n");
+ for (int i = 0; i < dim(atan2_table); i++) {
+ real y = atan2_table[i].y;
+ real x = atan2_table[i].x;
+ real atan2_value = atan2_table[i].atan2;
+ check2("atan2", atan2, y, x, atan2_value);
+ }
+}
+
+check_tan();
+exit(errors);
commit efbd076f61d3e6dc42216471963a5980a796b0dd
Author: Keith Packard <keithp at keithp.com>
Date: Mon Oct 6 20:31:43 2025 -0700
test: Test cbrt and sqrt
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/test/math.5c b/test/math.5c
index b4310c3..0188c05 100644
--- a/test/math.5c
+++ b/test/math.5c
@@ -85,6 +85,57 @@ void check2(string op, real(real, real) f, real input1, real input2, real correc
}
}
+void check_sqrt()
+{
+ void check_value (real z, int prec) {
+ z = imprecise(z, prec);
+ real s = sqrt(z);
+ real err = z - s*s;
+ int errdist = exponent(z) - exponent(err);
+ if (err != 0 && errdist < prec - 1) {
+ printf("%g = sqrt(%g): error too large (differ at %d bits instead of %d bits)\n",
+ s, z, errdist, prec);
+ errors++;
+ }
+ }
+
+ printf("Checking sqrt...");
+ for (int prec = 10; prec <= 1000; prec *= 10) {
+ printf(" %d", prec);
+ File::flush(stdout);
+ for (real z = 1e-200; z < 1e200; z *= 1.3) {
+ check_value(z, prec);
+ }
+ }
+ printf("\n");
+}
+
+void check_cbrt()
+{
+ void check_value (real z, int prec) {
+ z = imprecise(z, prec);
+ real s = cbrt(z);
+ real err = imprecise(z, prec+10) - imprecise(s, prec+10)**3;
+ int errdist = exponent(z) - exponent(err);
+ if (err != 0 && errdist < prec - 2) {
+ printf("%g = cbrt(%g): error too large (differ at %d bits instead of %d bits)\n",
+ s, z, errdist, prec);
+ errors++;
+ }
+ }
+
+ printf("Checking cbrt...");
+ for (int prec = 10; prec <= 1000; prec *= 10) {
+ printf(" %d", prec);
+ File::flush(stdout);
+ for (real z = 1e-200; z < 1e200; z *= 1.3) {
+ check_value(z, prec);
+ check_value(-z, prec);
+ }
+ }
+ printf("\n");
+}
+
void check_sin_cos()
{
printf("Checking sin and cos\n");
@@ -327,6 +378,8 @@ void check_jn()
}
}
+check_sqrt();
+check_cbrt();
check_sin_cos();
check_tan();
check_exp();
commit 389785e1c6ba86a3cbe7127c51ab21db27915421
Author: Keith Packard <keithp at keithp.com>
Date: Sun Oct 5 13:07:16 2025 -0700
test: Test jn function
Use bc 'j' function to generate test vectors.
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/test/math-tables.sh b/test/math-tables.sh
index 4cd8bab..936b8b6 100755
--- a/test/math-tables.sh
+++ b/test/math-tables.sh
@@ -19,6 +19,7 @@ make_tan=y
make_atan=y
make_log=y
make_exp=y
+make_jn=y
if [ $# -gt 0 ]; then
make_sin=n
@@ -27,6 +28,7 @@ if [ $# -gt 0 ]; then
make_atan=n
make_log=n
make_exp=n
+ make_jn=n
for i in "$@"; do
case "$i" in
sin)
@@ -47,6 +49,9 @@ if [ $# -gt 0 ]; then
exp)
make_exp=y
;;
+ jn)
+ make_jn=y;
+ ;;
esac
done
fi
@@ -170,3 +175,23 @@ while [ "$r" -le 1000 -a $make_exp = "y" ]; do
r=`expr "$r" + "$inc"`
done
echo "};"
+
+# jn table
+
+echo "typedef struct { int n; real in, jn; } jn_t;"
+echo "jn_t[] jn_table = {"
+
+n="-5"
+r="-1000"
+while [ "$n" -le 5 -a $make_jn = "y" ]; do
+ exp=`echo "j($n, $r / 100)" | bc -l "$dir"/math-funcs.bc | fmt -w 500 | tr -d '\\\n ' `
+ echo " { .n = $n, .in = $r / 100.0,"
+ echo " .jn = $exp },"
+ r=`expr "$r" + "$inc" '*' 4`
+ if [ "$r" -gt 1000 ]; then
+ n=`expr "$n" + 1`
+ r=-1000
+ fi
+done
+echo "};"
+
diff --git a/test/math.5c b/test/math.5c
index 409ee20..b4310c3 100644
--- a/test/math.5c
+++ b/test/math.5c
@@ -305,18 +305,34 @@ void check_gamma()
invalid = true;
}
if (invalid != (floor(z) == z)) {
- printf("%g: negative wrong %b\n", z, invalid);
+ printf("%g: negative wrong %v\n", z, invalid);
errors++;
}
}
}
}
+void check_jn()
+{
+ real jn_fixed(real n, real x) {
+ return jn(floor(n), x);
+ }
+
+ printf("Checking jn\n");
+ for (int i = 0; i < dim(jn_table); i++) {
+ int n = jn_table[i].n;
+ real in = jn_table[i].in;
+ real jn_value = jn_table[i].jn;
+ check2("jn", jn_fixed, n, in, jn_value);
+ }
+}
+
check_sin_cos();
check_tan();
check_exp();
check_log();
check_e();
check_gamma();
+check_jn();
exit (errors);
commit 51a0cb0567cd711d3df56d88273819d2a6b6e1be
Author: Keith Packard <keithp at keithp.com>
Date: Sun Oct 5 10:35:25 2025 -0700
test: Test gamma using Legendre's duplication formula
Ensure the gamma implementation at least satisfies the Legendre
duplication formula for a range of values.
Γ(z)Γ(z+½) = 2^(1-2z) √π Γ(2z)
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/test/math.5c b/test/math.5c
index 9b30768..409ee20 100644
--- a/test/math.5c
+++ b/test/math.5c
@@ -277,10 +277,46 @@ void check_log()
}
}
+void check_gamma()
+{
+ printf("Checking gamma using the Legendre duplication formula\n");
+
+ void check_value(real z, int prec) {
+ int iprec = prec + 20 + abs(exponent(imprecise(z)));
+ real zi = imprecise(z, iprec);
+ real left = gamma(zi) * gamma(zi + 1/2);
+ real right = 2**(1-2*zi) * sqrt(pi_value(iprec)) * gamma(zi * 2);
+
+ real error = left - right;
+ if (error != 0 && exponent(error) > exponent(left) - prec) {
+ printf("%g: Γ(z)Γ(z+½) (%g) != 2^(1-2z) √π Γ(2z) (%g): prec %d error %g\n",
+ z, left, right, prec, error);
+ errors++;
+ }
+ }
+
+ for (int prec = 10; prec <= 100; prec *= 10) {
+ for (real z = 1e-38; z < 1e38; z *= 11) {
+ bool invalid = false;
+ check_value(z, prec);
+ try {
+ check_value(-z, prec);
+ } catch invalid_argument(string err, int i, real z) {
+ invalid = true;
+ }
+ if (invalid != (floor(z) == z)) {
+ printf("%g: negative wrong %b\n", z, invalid);
+ errors++;
+ }
+ }
+ }
+}
+
check_sin_cos();
check_tan();
check_exp();
check_log();
check_e();
+check_gamma();
exit (errors);
commit 85ee6d237ebf6dc9cc337d8df61ed40810830677
Author: Keith Packard <keithp at keithp.com>
Date: Tue Oct 7 13:47:09 2025 -0700
Add 'complex' type.
Provide a complex type to simplify use of complex numbers in nickle
programs.
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/builtin-complex.c b/builtin-complex.c
new file mode 100644
index 0000000..4fa9aa3
--- /dev/null
+++ b/builtin-complex.c
@@ -0,0 +1,57 @@
+/*
+ * Copyright © 2025 Keith Packard and Bart Massey.
+ * All Rights Reserved. See the file COPYING in this directory
+ * for licensing information.
+ */
+
+#include "builtin.h"
+
+NamespacePtr ComplexNamespace;
+
+Value do_Complex_creal (Value z);
+Value do_Complex_cimag (Value z);
+
+void
+import_Complex_namespace()
+{
+ ENTER();
+ static const struct fbuiltin_1 funcs_1[] = {
+ { do_Complex_creal, "creal", "N", "C", "\n"
+ " real creal(complex c)\n"
+ "\n"
+ " Returns the real component of a complex value\n" },
+ { do_Complex_cimag, "cimag", "N", "C", "\n"
+ " real cimag(complex c)\n"
+ "\n"
+ " Returns the imaginary component of a complex value\n" },
+ { 0 }
+ };
+
+ static const struct fbuiltin_2 funcs_2[] = {
+ { NewValueComplex, "cmplx", "C", "NN", "\n"
+ " complex cmplx(real r, real i)\n"
+ "\n"
+ " Returns a new complex value given real and imaginary components\n" },
+ { 0 },
+ };
+
+ ComplexNamespace = BuiltinNamespace(/*parent*/ 0, "Complex")->namespace.namespace;
+
+ BuiltinFuncs1 (&ComplexNamespace, funcs_1);
+ BuiltinFuncs2 (&ComplexNamespace, funcs_2);
+ EXIT();
+}
+
+Value do_Complex_creal (Value z)
+{
+ if (ValueIsComplex(z))
+ return z->complex.r;
+ return z;
+}
+
+Value do_Complex_cimag (Value z)
+{
+ if (ValueIsComplex(z))
+ return z->complex.i;
+ return Zero;
+}
diff --git a/builtin-namespaces.h b/builtin-namespaces.h
index 6f58c4f..e1ef2d2 100644
--- a/builtin-namespaces.h
+++ b/builtin-namespaces.h
@@ -24,3 +24,4 @@ extern void import_Socket_namespace(void);
extern void import_Foreign_namespace(void);
extern void import_PID_namespace(void);
extern void import_Date_namespace(void);
+extern void import_Complex_namespace(void);
diff --git a/builtin-toplevel.c b/builtin-toplevel.c
index 946e234..6e216ab 100644
--- a/builtin-toplevel.c
+++ b/builtin-toplevel.c
@@ -148,8 +148,8 @@ import_Toplevel_namespace(void)
" int numerator (rational r)\n"
"\n"
" Return the numerator of 'r'.\n" },
- { do_precision, "precision", "i", "R", "\n"
- " int precision (real r)\n"
+ { do_precision, "precision", "i", "C", "\n"
+ " int precision (complex r)\n"
"\n"
" Return the number of bits in the\n"
" representation of the mantissa of 'r'.\n" },
@@ -228,9 +228,9 @@ import_Toplevel_namespace(void)
};
static const struct fbuiltin_v funcs_v[] = {
- { do_imprecise, "imprecise", "R", "R.i", "\n"
- " real imprecise (real r)\n"
- " real imprecise (real r, int precision)\n"
+ { do_imprecise, "imprecise", "C", "C.i", "\n"
+ " complex imprecise (complex r)\n"
+ " complex imprecise (complex r, int precision)\n"
"\n"
" Return an imprecise number.\n"
" The precision will be 'precision' if supplied, else 256.\n" },
@@ -378,6 +378,8 @@ do_imprecise (int n, Value *p)
Value float_prec;
if (ValueIsFloat(v))
RETURN(v);
+ if (ValueIsComplex(v) && ValueIsFloat(v->complex.r) && ValueIsFloat(v->complex.i))
+ RETURN(v);
prec = DEFAULT_FLOAT_PREC;
float_prec = lookupVar(0, "float_precision");
if (float_prec)
@@ -387,8 +389,12 @@ do_imprecise (int n, Value *p)
prec = default_prec;
}
}
-
- RETURN (NewValueFloat (v, prec));
+ if (ValueIsComplex(v)) {
+ v = NewValueComplex(NewValueFloat(v->complex.r, prec), NewValueFloat(v->complex.i, prec));
+ } else {
+ v = NewValueFloat (v, prec);
+ }
+ RETURN(v);
}
Value
@@ -558,10 +564,20 @@ do_precision (Value av)
ENTER ();
unsigned prec;
- if (ValueIsFloat(av))
+ switch (ValueTag(av)) {
+ case rep_float:
prec = av->floats.prec;
- else
+ break;
+ case rep_complex:
+ prec = ValueInt(do_precision(av->complex.i));
+ unsigned iprec = ValueInt(do_precision(av->complex.r));
+ if (iprec > prec)
+ prec = iprec;
+ break;
+ default:
prec = 0;
+ break;
+ }
RETURN (NewInt (prec));
}
@@ -722,6 +738,7 @@ do_is_number (Value av)
case rep_integer:
case rep_rational:
case rep_float:
+ case rep_complex:
av = TrueVal;
break;
default:
@@ -964,9 +981,8 @@ Value do_hash_test (Value hv, Value key)
{
return HashTest (hv, key);
}
-
+
Value do_hash_keys (Value hv)
{
return HashKeys (hv);
}
-
diff --git a/builtin.5c b/builtin.5c
index e1d0159..5e153e9 100644
--- a/builtin.5c
+++ b/builtin.5c
@@ -124,6 +124,7 @@ library "math-ext.5c"
import Math;
library "scanf.5c";
library "socket.5c";
+library "complex.5c";
# Now autoload/autoimport the bonus stuff
# (But why? Let the users do it.)
diff --git a/builtin.c b/builtin.c
index aae62d9..fbe310a 100644
--- a/builtin.c
+++ b/builtin.c
@@ -183,6 +183,7 @@ BuiltinType (char *format, Type **type, Bool arg)
case 'p': t = typePoly; break;
case 'n': t = typePrim[rep_float]; break;
case 'N': t = typePrim[rep_float]; break;
+ case 'C': t = typePrim[rep_complex]; break;
case 'E': t = typeFileError; break;
case 'R': t = typePrim[rep_float]; break;
case 'r': t = typePrim[rep_rational]; break;
@@ -348,6 +349,7 @@ BuiltinInit (void)
import_Foreign_namespace ();
import_PID_namespace ();
import_Date_namespace();
+ import_Complex_namespace();
/* Import builtin strings with predefined values */
BuiltinStrings (svars);
diff --git a/complex.5c b/complex.5c
new file mode 100644
index 0000000..b1df539
--- /dev/null
+++ b/complex.5c
@@ -0,0 +1,176 @@
+/*
+ * Copyright © 2025 Keith Packard and Bart Massey.
+ * All Rights Reserved. See the file COPYING in this directory
+ * for licensing information.
+ */
+
+extend namespace Complex {
+ protected complex(real r, real i) new = cmplx;
+
+ protected complex i = new(0,1);
+
+ public complex polar(real r, real t)
+ {
+ return new(r * cos(t), r * sin(t));
+ }
+
+ public complex conj(complex z)
+ {
+ return new(creal(z), -cimag(z));
+ }
+
+ protected real abs(complex z)
+ {
+ return Math::sqrt(creal(z)**2 + cimag(z)**2);
+ }
+
+ public real arg(complex z)
+ {
+ return Math::atan2(cimag(z), creal(z));
+ }
+
+ protected complex log(complex z)
+ {
+ return new(Math::log(Math::sqrt(creal(z)**2 + cimag(z)**2)),
+ Math::atan2(cimag(z), creal(z)));
+ }
+
+ protected complex exp(complex z)
+ {
+ real r = Math::exp(creal(z));
+ return new(r * Math::cos(cimag(z)),
+ r * Math::sin(cimag(z)));
+ }
+
+ protected complex pow(complex x, complex y)
+ {
+ x = imprecise(x);
+ y = imprecise(y);
+ /* add precision to cover the atan2 + sin/cos operations */
+ int prec = min(precision(x), precision(y));
+ int iprec = prec + abs(exponent(creal(x))) + abs(exponent(cimag(x)));
+ complex xi = imprecise(x, iprec);
+ complex yi = imprecise(y, iprec);
+ return imprecise(exp(log(xi) * yi), prec);
+ }
+
+ protected complex sqrt(complex z)
+ {
+ if (cimag(z) == 0)
+ return Math::sqrt(creal(z));
+
+ real mod = abs(z);
+ real sgn = cimag(z) > 0 ? 1 : -1;
+ return new(Math::sqrt((mod + creal(z)) / 2),
+ sgn * Math::sqrt((mod - creal(z)) / 2));
+ }
+
+ protected complex atan(complex z)
+ {
+ z = imprecise(z);
+ int prec = precision(z);
+ int e = exponent(cimag(z));
+ int iprec = prec + 64;
+ if (e < 0)
+ iprec -= 2*e;
+
+ complex zp = imprecise(z, iprec);
+ complex i_over_2 = new(0, 0.5);
+ return imprecise(i_over_2 * log((i+z) / (i-z)), prec);
+ }
+
+ protected complex asin(complex z)
+ {
+ if (Math::abs(creal(z)) == 1 && cimag(z) == 0)
+ raise invalid_argument("complex asin of 1", 0, z);
+
+ z = imprecise(z);
+ int prec = precision(z);
+ int e = min(exponent(creal(z)), exponent(cimag(z)));
+ int iprec = prec + 14;
+ if (e < 0)
+ iprec -= 2*e;
+ complex zp = imprecise(z, iprec);
+ return imprecise(atan(zp / sqrt(1-zp*zp)), prec);
+ }
+
+ protected complex acos(complex z)
+ {
+ if (Math::abs(creal(z)) == 1 && cimag(z) == 0)
+ raise invalid_argument("complex acos of 1", 0, z);
+
+ z = imprecise(z);
+ int prec = precision(z);
+ int iprec = prec + 14;
+ complex asp = asin(imprecise(z, iprec));
+ return imprecise(new(pi_value(iprec)/2 - creal(asp), -cimag(asp)), prec);
+ }
+
+ protected complex cos(complex z)
+ {
+ return new(Math::cos(creal(z)) * Math::cosh(cimag(z)),
+ -Math::sin(creal(z)) * Math::sinh(cimag(z)));
+ }
+
+ protected complex sin(complex z)
+ {
+ return new(Math::sin(creal(z)) * Math::cosh(cimag(z)),
+ Math::cos(creal(z)) * Math::sinh(cimag(z)));
+ }
+
+ protected complex tan(complex z)
+ {
+ return sin(z) / cos(z);
+ }
+
+ protected complex cosh(complex z)
+ {
+ return new(Math::cosh(creal(z)) * Math::cos(cimag(z)),
+ Math::sinh(creal(z)) * Math::sin(cimag(z)));
+ }
+
+ protected complex sinh(complex z)
+ {
+ return new(Math::sinh(creal(z)) * Math::cos(cimag(z)),
+ Math::cosh(creal(z)) * Math::sin(cimag(z)));
+ }
+
+ protected complex tanh(complex z)
+ {
+ return sinh(z) / cosh(z);
+ }
+
+ protected complex atanh(complex z)
+ {
+ z = imprecise(z);
+ int prec = precision(z);
+ int iprec = prec + 512;
+ int e = min(exponent(creal(z)), exponent(cimag(z)));
+ if (e < 0)
+ iprec += e;
+ complex zp = imprecise(z, iprec);
+ return imprecise(0.5 * log((1 + zp)/(1 - zp)), prec);
+ }
+
+ protected complex asinh(complex z)
+ {
+ int prec = precision(z);
+ int iprec = prec + 512;
+ int e = min(exponent(creal(z)), exponent(cimag(z)));
+ if (e < 0)
+ iprec += e;
+ complex zp = imprecise(z, iprec);
+ return imprecise(log(zp + sqrt(zp*zp + 1)), prec);
+ }
+
+ protected complex acosh(complex z)
+ {
+ int prec = precision(z);
+ int iprec = prec + 512;
+ int e = min(exponent(creal(z)), exponent(cimag(z)));
+ if (e < 0)
+ iprec += e;
+ complex zp = imprecise(z, iprec);
+ return imprecise(log(zp + sqrt(zp + 1) * sqrt(zp - 1)), prec);
+ }
+}
diff --git a/complex.c b/complex.c
new file mode 100644
index 0000000..1f11e73
--- /dev/null
+++ b/complex.c
@@ -0,0 +1,211 @@
+/*
+ * Copyright © 2025 Keith Packard and Bart Massey.
+ * All Rights Reserved. See the file COPYING in this directory
+ * for licensing information.
+ */
+
+#include <math.h>
+#include "nickle.h"
+
+static Value
+ComplexPlus (Value av, Value bv, int expandOk)
+{
+ (void) expandOk;
+ ENTER();
+ Value ret = NewValueComplex(Plus(av->complex.r, bv->complex.r),
+ Plus(av->complex.i, bv->complex.i));
+ RETURN(ret);
+}
+
+static Value
+ComplexMinus (Value av, Value bv, int expandOk)
+{
+ (void) expandOk;
+ ENTER();
+ Value ret = NewValueComplex(Minus(av->complex.r, bv->complex.r),
+ Minus(av->complex.i, bv->complex.i));
+ RETURN(ret);
+}
+
+static Value
+ComplexTimes (Value av, Value bv, int expandOk)
+{
+ (void) expandOk;
+ ENTER();
+ Value ret = NewValueComplex(Minus(Times(av->complex.r, bv->complex.r),
+ Times(av->complex.i, bv->complex.i)),
+ Plus(Times(av->complex.r, bv->complex.i),
+ Times(av->complex.i, bv->complex.r)));
+ RETURN(ret);
+}
+
+static Value
+ComplexRecip (Value av)
+{
+ ENTER();
+ Value modsq = Plus(Times(av->complex.r, av->complex.r),
+ Times(av->complex.i, av->complex.i));
+ Value ret = NewValueComplex(Divide(av->complex.r, modsq),
+ Divide(Negate(av->complex.i), modsq));
+ RETURN(ret);
+}
+
+static Value
+ComplexDivide (Value av, Value bv, int expandOk)
+{
+ (void) expandOk;
+ ENTER();
+ Value recip = ComplexRecip(bv);
+ Value ret = ComplexTimes(av, recip, expandOk);
+ RETURN(ret);
+}
+
+static Value
+ComplexLess (Value av, Value bv, int expandOk)
+{
+ (void) expandOk;
+ /* Dictionary ordering is at least consistent */
+ if (True(Less(av->complex.r, bv->complex.r)))
+ return TrueVal;
+ if (True(Equal(av->complex.r, bv->complex.r)))
+ return Less(av->complex.i, bv->complex.i);
+ return FalseVal;
+}
+
+static Value
+ComplexEqual (Value av, Value bv, int expandOk)
+{
+ (void) expandOk;
+ Value ret;
+ ret = Equal(av->complex.r, bv->complex.r);
+ if (True(ret))
+ ret = Equal(av->complex.i, bv->complex.i);
+ return ret;
+}
+
+static Value
+ComplexNegate (Value av, int expandOk)
+{
+ (void) expandOk;
+ ENTER();
+ Value ret = NewValueComplex(Negate(av->complex.r), Negate(av->complex.i));
+ RETURN(ret);
+}
+
+static Value
+ComplexFloor (Value av, int expandOk)
+{
+ (void) expandOk;
+ ENTER();
+ Value ret = NewValueComplex(Floor(av->complex.r), Floor(av->complex.i));
+ RETURN(ret);
+}
+
+static Value
+ComplexCeil (Value av, int expandOk)
+{
+ (void) expandOk;
+ ENTER();
+ Value ret = NewValueComplex(Ceil(av->complex.r), Ceil(av->complex.i));
+ RETURN(ret);
+}
+
+static Value
+ComplexPromote (Value av, Value bv)
+{
+ (void) bv;
+ ENTER ();
+ if (!ValueIsComplex(av))
+ av = NewComplex(av, Zero);
+ RETURN(av);
+}
+
+static Value
+ComplexReduce (Value av)
+{
+ if (Zerop(av->complex.i))
+ return av->complex.r;
+ return av;
+}
+
+static Bool
+ComplexPrint(Value f, Value fv, char format, int base, int width, int prec, int fill)
+{
+ Bool value = format == 'v';
+
+ if (value)
+ FilePuts(f, "cmplx(");
+ if (!Print(f, fv->complex.r, format, base, width, prec, fill))
+ return False;
+ if (value)
+ FilePuts(f, ", ");
+ else
+ FilePuts(f, "+i");
+ if (!Print(f, fv->complex.i, format, base, width, prec, fill))
+ return False;
+ if (value)
+ FilePuts(f, ")");
+ return True;
+}
+
+static HashValue
+ComplexHash (Value av)
+{
+ return ValueInt(ValueHash(av->complex.r)) ^ ValueInt(ValueHash(av->complex.i));
+}
+
+static void
+ComplexMark (void *object)
+{
+ Complex *c = object;
+
+ MemReference(c->r);
+ MemReference(c->i);
+}
+
+ValueRep ComplexRep = {
+ { ComplexMark, 0, "ComplexRep" }, /* base */
+ rep_complex, /* tag */
+ { /* binary */
+ ComplexPlus,
+ ComplexMinus,
+ ComplexTimes,
+ ComplexDivide,
+ NumericDiv,
+ NumericMod,
+ ComplexLess,
+ ComplexEqual,
+ 0,
+ 0,
+ },
+ { /* unary */
+ ComplexNegate,
+ ComplexFloor,
+ ComplexCeil,
+ },
+ ComplexPromote,
+ ComplexReduce,
+ ComplexPrint,
+ 0,
+ ComplexHash,
+};
+
+Value
+NewComplex (Value r, Value i)
+{
+ ENTER();
+ Value ret;
+
+ ret = ALLOCATE (&ComplexRep.data, sizeof (Complex));
+ ret->complex.r = r;
+ ret->complex.i = i;
+ RETURN (ret);
+}
+
+Value
+NewValueComplex (Value r, Value i)
+{
+ if (!Zerop(i))
+ r = NewComplex(r, i);
+ return r;
+}
diff --git a/file.c b/file.c
index 7341263..dc935f0 100644
--- a/file.c
+++ b/file.c
@@ -1443,6 +1443,9 @@ FilePutRep (Value f, Rep tag, Bool minimal)
case rep_float:
FilePuts (f, "real");
break;
+ case rep_complex:
+ FilePuts (f, "complex");
+ break;
case rep_string:
FilePuts (f, "string");
break;
diff --git a/float.c b/float.c
index 4821f99..cd2edce 100644
--- a/float.c
+++ b/float.c
@@ -644,7 +644,7 @@ FloatPromote (Value av, Value bv)
ENTER ();
int prec;
- if (!ValueIsFloat(av))
+ if (!ValueIsFloat(av) && !ValueIsComplex(av))
{
prec = DEFAULT_FLOAT_PREC;
if (bv && ValueIsFloat(bv))
@@ -1289,6 +1289,8 @@ NewValueFloat (Value av, unsigned prec)
av = NewRationalFloat (&av->rational, prec);
break;
case rep_float:
+ if (prec == av->floats.prec)
+ break;
av = NewFloat (av->floats.mant, av->floats.exp, prec);
break;
default:
diff --git a/gram.y b/gram.y
index 694edca..23e3dbd 100644
--- a/gram.y
+++ b/gram.y
@@ -138,7 +138,7 @@ ParseNewSymbol (Publish publish, Class class, Type *type, Atom name);
%token SEMI MOD OC CC DOLLAR DOTDOTDOT ENDFILE
%token <class> GLOBAL AUTO STATIC CONST
-%token <type> POLY INTEGER NATURAL RATIONAL REAL STRING FOREIGN
+%token <type> POLY INTEGER NATURAL RATIONAL REAL COMPLEX STRING FOREIGN
%token <type> FILET MUTEX SEMAPHORE CONTINUATION THREAD VOID BOOL
%token FUNCTION FUNC EXCEPTION RAISE
%token TYPEDEF IMPORT NEW ANONINIT
@@ -1040,6 +1040,7 @@ basetype : POLY
| INTEGER
| RATIONAL
| REAL
+ | COMPLEX
| STRING
| FILET
| MUTEX
diff --git a/lex.l b/lex.l
index 253b556..90a55e6 100644
--- a/lex.l
+++ b/lex.l
@@ -367,6 +367,7 @@ bool { yylval.type = typePrim[rep_bool]; return BOOL; }
int { yylval.type = typePrim[rep_integer]; return INTEGER; }
rational { yylval.type = typePrim[rep_rational]; return RATIONAL; }
real { yylval.type = typePrim[rep_float]; return REAL; }
+complex { yylval.type = typePrim[rep_complex]; return COMPLEX; }
string { yylval.type = typePrim[rep_string]; return STRING; }
file { yylval.type = typePrim[rep_file]; return FILET; }
semaphore { yylval.type = typePrim[rep_semaphore]; return SEMAPHORE; }
diff --git a/math.5c b/math.5c
index 83f24a9..034d058 100644
--- a/math.5c
+++ b/math.5c
@@ -1149,6 +1149,11 @@ extend namespace Math {
c *= i;
return c / k!;
}
+
+ real(real) _abs = abs;
+ public real(real) abs = _abs;
+ int(real) _precision = precision;
+ public int(real) precision = _precision;
}
/* XXX these shouldn't be here, but it was *convenient* */
diff --git a/meson.build b/meson.build
index 35837d5..fae4c36 100644
--- a/meson.build
+++ b/meson.build
@@ -200,7 +200,7 @@ if history_h != ''
endif
nickle_srcs = [
- 'alarm.c', 'array.c', 'atom.c', 'box.c', 'compile.c', 'debug.c',
+ 'alarm.c', 'array.c', 'atom.c', 'box.c', 'compile.c', 'complex.c', 'debug.c',
'divide.c', 'edit.c', 'error.c', 'execute.c', 'expr.c', 'file.c',
'float.c', 'foreign.c', 'frame.c', 'func.c', 'gcd.c', 'hash.c',
'int.c', 'integer.c', 'io.c', 'main.c', 'mem.c', 'natural.c',
@@ -213,6 +213,7 @@ nickle_srcs = [
'builtin-semaphore.c', 'builtin-sockets.c', 'builtin-string.c',
'builtin-thread.c', 'builtin-toplevel.c', 'builtin-pid.c',
'builtin-date.c', 'builtin.c', 'builtin.h', 'builtin-foreign.c',
+ 'builtin-complex.c',
]
nickle_files = [
@@ -221,7 +222,7 @@ nickle_files = [
'printf.5c', 'history.5c', 'ctype.5c', 'string.5c', 'socket.5c',
'file.5c', 'parse-args.5c', 'svg.5c', 'process.5c',
'prime_sieve.5c', 'factorial.5c', 'gamma.5c', 'sort.5c', 'list.5c', 'skiplist.5c',
- 'json.5c', 'cha-cha.5c', 'math-ext.5c',
+ 'json.5c', 'cha-cha.5c', 'math-ext.5c', 'complex.5c'
]
nickle_headers = [
diff --git a/nickle.h b/nickle.h
index f65cb69..00c166b 100644
--- a/nickle.h
+++ b/nickle.h
@@ -695,6 +695,7 @@ extern NamespacePtr GcdNamespace;
#endif
extern NamespacePtr EnvironNamespace;
extern NamespacePtr DateNamespace;
+extern NamespacePtr ComplexNamespace;
void BuiltinInit (void);
diff --git a/type.c b/type.c
index 60c62ab..5cdf179 100644
--- a/type.c
+++ b/type.c
@@ -653,12 +653,12 @@ TypeBinaryGroup (Type *left, Type *right)
if (TypePoly (left))
{
if (TypePoly (right) || TypeNumeric (right))
- return typePrim[rep_float];
+ return typePrim[rep_complex];
}
else if (TypePoly (right))
{
if (TypeNumeric (left))
- return typePrim[rep_float];
+ return typePrim[rep_complex];
}
else if (TypeNumeric (left) && TypeNumeric (right))
{
@@ -702,12 +702,12 @@ TypeBinaryField (Type *left, Type *right)
if (TypePoly (left))
{
if (TypePoly (right) || TypeNumeric (right))
- return typePrim[rep_float];
+ return typePrim[rep_complex];
}
else if (TypePoly (right))
{
if (TypeNumeric (left))
- return typePrim[rep_float];
+ return typePrim[rep_complex];
}
else if (TypeNumeric (left) && TypeNumeric (right))
{
@@ -728,9 +728,9 @@ static Type *
TypeBinaryDiv (Type *left, Type *right)
{
if (TypePoly (left))
- left = typePrim[rep_float];
+ left = typePrim[rep_complex];
if (TypePoly (right))
- right = typePrim[rep_float];
+ right = typePrim[rep_complex];
if (TypeNumeric (left) && TypeNumeric (right))
{
return typePrim[rep_integer];
@@ -746,14 +746,14 @@ static Type *
TypeBinaryPow (Type *left, Type *right)
{
if (TypePoly (left))
- left = typePrim[rep_float];
+ left = typePrim[rep_complex];
if (TypePoly (right))
- right = typePrim[rep_float];
+ right = typePrim[rep_complex];
if (TypeNumeric (left) && TypeNumeric (right))
{
if (TypeIntegral (right))
return left;
- return typePrim[rep_float];
+ return typePrim[rep_complex];
}
return 0;
}
@@ -821,7 +821,7 @@ static Type *
TypeUnaryGroup (Type *type)
{
if (TypePoly (type))
- return typePrim[rep_float];
+ return typePrim[rep_complex];
else if (TypeNumeric (type))
return type;
return 0;
diff --git a/value.c b/value.c
index 441efd0..a7074d7 100644
--- a/value.c
+++ b/value.c
@@ -19,6 +19,23 @@ volatile Bool signaling;
#ifndef Numericp
Bool
Numericp (Rep t)
+{
+ switch (t) {
+ case rep_int:
+ case rep_integer:
+ case rep_rational:
+ case rep_float:
+ case rep_complex:
+ return True;
+ default:;
+ }
+ return False;
+}
+#endif
+
+#ifndef Scalarp
+Bool
+Scalarp (Rep t)
{
switch (t) {
case rep_int:
@@ -58,6 +75,8 @@ Zerop (Value av)
return av->rational.num->length == 0;
case rep_float:
return av->floats.mant->mag->length == 0;
+ case rep_complex:
+ return Zerop(av->complex.r) && Zerop(av->complex.i);
default:;
}
return False;
@@ -75,6 +94,8 @@ Negativep (Value av)
return av->rational.sign == Negative;
case rep_float:
return av->floats.mant->sign == Negative;
+ case rep_complex:
+ return Negativep(av->complex.r);
default:;
}
return False;
diff --git a/value.h b/value.h
index 8536239..6e8291f 100644
--- a/value.h
+++ b/value.h
@@ -279,34 +279,36 @@ typedef enum _rep {
rep_integer = 1,
rep_rational = 2,
rep_float = 3,
- rep_string = 4,
- rep_file = 5,
- rep_thread = 6,
- rep_semaphore = 7,
- rep_continuation = 8,
- rep_bool = 9,
- rep_foreign = 10,
- rep_void = 11,
+ rep_complex = 4,
+ rep_string = 5,
+ rep_file = 6,
+ rep_thread = 7,
+ rep_semaphore = 8,
+ rep_continuation = 9,
+ rep_bool = 10,
+ rep_foreign = 11,
+ rep_void = 12,
/* composite types */
- rep_ref = 12,
- rep_func = 13,
+ rep_ref = 13,
+ rep_func = 14,
/* mutable type */
- rep_array = 14,
- rep_struct = 15,
- rep_union = 16,
- rep_hash = 17
+ rep_array = 15,
+ rep_struct = 16,
+ rep_union = 17,
+ rep_hash = 18
} Rep;
/* because rep_undef is -1, using (unsigned) makes these a single compare */
-#define Numericp(t) ((unsigned) (t) <= (unsigned) rep_float)
+#define Numericp(t) ((unsigned) (t) <= (unsigned) rep_complex)
+#define Scalarp(t) ((unsigned) (t) <= (unsigned) rep_float)
#define Integralp(t) ((unsigned) (t) <= (unsigned) rep_integer)
#define Mutablep(t) ((t) >= rep_array)
-extern ValueRep IntRep, IntegerRep, RationalRep, FloatRep;
+extern ValueRep IntRep, IntegerRep, RationalRep, FloatRep, ComplexRep;
extern ValueRep StringRep, ArrayRep, FileRep;
extern ValueRep RefRep, StructRep, UnionRep, HashRep;
extern ValueRep FuncRep, ThreadRep;
@@ -369,6 +371,7 @@ static inline ValueRep *_ValueRep(Value v);
#define ValueIsInteger(v) (ValueRep(v) == &IntegerRep)
#define ValueIsRational(v) (ValueRep(v) == &RationalRep)
#define ValueIsFloat(v) (ValueRep(v) == &FloatRep)
+#define ValueIsComplex(v) (ValueRep(v) == &ComplexRep)
#define ValueIsString(v) (ValueRep(v) == &StringRep)
#define ValueIsArray(v) (ValueRep(v) == &ArrayRep)
#define ValueIsFile(v) (ValueRep(v) == &FileRep)
@@ -596,6 +599,11 @@ typedef struct _float {
unsigned prec;
} Float;
+typedef struct _complex {
+ BaseValue base;
+ Value r, i;
+} Complex;
+
typedef struct _string {
BaseValue base;
size_t length;
@@ -887,6 +895,7 @@ typedef union _value {
Integer integer;
Rational rational;
Float floats;
+ Complex complex;
String string;
Array array;
File file;
@@ -1114,6 +1123,8 @@ Value NewNaturalFloat (Sign sign, Natural *n, unsigned prec);
Value NewRationalFloat (Rational *r, unsigned prec);
Value NewValueFloat (Value av, unsigned prec);
Value NewDoubleFloat (double d);
+Value NewComplex (Value r, Value i);
+Value NewValueComplex (Value r, Value i);
Value NewContinuation (ContinuationPtr continuation, InstPtr pc);
Value NewForeign (const char *id, void *data, void (*mark)(void *data), void (*free)(void *data));
@@ -1182,6 +1193,7 @@ Value UnaryOperate (Value v, UnaryOp operator);
Value NumericDiv (Value av, Value bv, int expandOk);
Value NumericMod (Value av, Value bv, int expandOk);
+
# define OK_TRUNC 1
extern Value Blank, Elementless, Void, TrueVal, FalseVal;
commit fc75fd6b1fe3410491ebf7ef9ab93155ebd02803
Author: Keith Packard <keithp at keithp.com>
Date: Tue Oct 7 10:32:04 2025 -0700
Check trig operands for lack of precision
Trig functions reduce their argment to the range -π < x <= π. If the
operand is large and doesn't have sufficient precision to represent a
specific value in that range, raise an invalid_argument instead of
mapping it to zero.
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/math.5c b/math.5c
index e87acb7..83f24a9 100644
--- a/math.5c
+++ b/math.5c
@@ -691,6 +691,8 @@ extend namespace Math {
aa = imprecise (aa);
my_pi = pi_value (precision (aa));
+ if (aa + my_pi == aa)
+ raise invalid_argument ("argument not precise enough", 0, aa);
aa %= 2 * my_pi;
if (aa > my_pi)
aa -= 2 * my_pi;
commit 24397499036b9d855d4b1358cd462ed9386ee526
Author: Keith Packard <keithp at keithp.com>
Date: Mon Oct 6 19:34:32 2025 -0700
Improve accuracy of sqrt and cbrt and avoid infinite loops
Iterate cbrt and sqrt until they either find a fixed point or
oscillate between two values. This avoids the imprecise loop count
estimate from sqrt and the infinite loop in cbrt.
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/math.5c b/math.5c
index 27362df..e87acb7 100644
--- a/math.5c
+++ b/math.5c
@@ -21,20 +21,33 @@ extend namespace Math {
{
real err;
real prev, cur;
- int n, iter;
v = imprecise (v);
- err = 1/(2**(precision (v)+3));
- prev = imprecise (1 / (2**(floor (exponent(v)/2))), precision(v));
- iter = precision (v) + 15;
- while (iter-- > 0)
- {
- cur = 0.5 * prev * (3 - v * prev**2);
- if (abs (cur - prev) < err)
- break;
+ int prec = precision(v);
+ int iprec = prec + 5 * bit_width(prec) + 10;
+ real vi = imprecise(v, iprec);
+ /* Initial estimate is 1/2 ** log2(v)/2 */
+ cur = imprecise(0.5 * (2**(exponent(v) ⫽ 2)), iprec);
+ int iter = iprec + 100;
+ do {
+ /*
+ * Newton's method iteration
+ *
+ * f(y) = y² - x
+ * f'(y) = 2y
+ *
+ * y' = y - f(y)/f'(y)
+ * = y - (y² - x) / 2y
+ * = y - y/2 + x/2y
+ * = y/2 + x/2y
+ */
prev = cur;
- }
- return abs (1/cur);
+ cur = prev/2 + vi/(2*prev);
+ /* bail if we take too long */
+ if (--iter <= 0)
+ break;
+ } while (imprecise(cur, prec+2) != imprecise(prev, prec+2));
+ return abs(imprecise(cur, prec));
}
if (v == 0)
@@ -52,9 +65,7 @@ extend namespace Math {
num = floor (num_s + 0.5);
den = floor (den_s + 0.5);
if (num * num == numerator (v) && den * den == denominator (v))
- {
return num/den;
- }
}
return real_sqrt (v);
}
@@ -71,20 +82,40 @@ extend namespace Math {
int s = sign (v);
v = imprecise (abs (v));
- int result_bits = precision (v);
- int intermediate_bits = result_bits + 10;
- v = imprecise (v, intermediate_bits);
- real err = imprecise (1/(2**(result_bits+3)), intermediate_bits);
-
- cur = imprecise (1 / (0.75 * 2**(floor (exponent(v)/3))),
- intermediate_bits);
+ int prec = precision (v);
+ /*
+ * Estimate about log2(prec) iterations, add 5 bits
+ * of precision for each one (for the 5 operations)
+ */
+ int iprec = prec + 5 * bit_width(prec) + 10;
+ v = imprecise (v, iprec);
+ /* Initial estimate is 3/4 ** log2(v)/3 */
+ cur = imprecise (0.75 * 2**(exponent(v) ⫽ 3), iprec);
+ int iter = iprec + 100;
do {
prev = cur;
- cur = 1/3 * (2 * prev + v / prev**2);
- } while (abs (cur - prev) > err);
- return s * imprecise (abs (cur), result_bits);
+ /*
+ * Newton's method iteration
+ *
+ * f(y) = y³ - x
+ * f'(y) = 3y²
+ *
+ * y' = y - f(y)/f'(y)
+ * = y - (y³ - x) / 3y²
+ * = y - y/3 + x/3y²
+ * = 2y/3 + x/3y²
+ */
+ cur = 0.{6} * prev + 0.{3} * v / (prev*prev);
+ /* bail if we take too long */
+ if (--iter <= 0)
+ break;
+ } while (imprecise(cur, prec+2) != imprecise(prev, prec+2));
+ return s * imprecise (abs (cur), prec);
}
+ if (v == 0)
+ return 0;
+
if (is_rational (v))
{
int num, den;
@@ -96,11 +127,8 @@ extend namespace Math {
den_s = real_cbrt (imprecise (den, bit_width(den) + 128));
num = floor (num_s + 0.5);
den = floor (den_s + 0.5);
- /* printf ("num %g den %g\n", num, den); */
if (num ** 3 == numerator (v) && den ** 3 == denominator (v))
- {
return num/den;
- }
}
return real_cbrt (v);
}
@@ -584,7 +612,7 @@ extend namespace Math {
low = 1;
while (high - low > 1)
{
- p = (high + low) // 2;
+ p = (high + low) ⫽ 2;
if (avail_prec (v, p) > prec)
high = p;
else
@@ -609,7 +637,7 @@ extend namespace Math {
p = 3;
q = 5;
mden = imprecise (den, digits * 4) ** 4;
- l = loops (1 / den, digits) // 2;
+ l = loops (1 / den, digits) ⫽ 2;
/*
* Need at least digits + log10(loops) for all intermediate
* computations
@@ -1083,7 +1111,7 @@ extend namespace Math {
/* binary search phase */
int ll = 0;
while (ul > ll + 1) {
- int step = (ul - ll) // 2;
+ int step = (ul - ll) ⫽ 2;
if (mask(b, ul - step)) {
ul -= step;
continue;
@@ -1106,7 +1134,7 @@ extend namespace Math {
/* This keeps the size of the intermediate terms
smaller below, and also makes the bounds check
slightly easier. */
- if (k > (n + 1) // 2)
+ if (k > (n + 1) ⫽ 2)
k = n - k;
if (n < 0 || k < 0)
return 0;
commit a6c231f54c6adc08bf35707eb237f0ed6cbd3554
Author: Keith Packard <keithp at keithp.com>
Date: Sun Oct 5 10:34:21 2025 -0700
Use Euler's reflection formula for gamma < -0.5
This avoids a long iteration for large negative values.
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/gamma.5c b/gamma.5c
index 533488b..17dc9cb 100644
--- a/gamma.5c
+++ b/gamma.5c
@@ -77,24 +77,27 @@ extend namespace Math {
return (floor(n)-1)!;
}
+ n = imprecise(n);
+ int bits = precision(n);
+
/*
- * Use Γ(z) = Γ(z+1)/z
- * to solve for negative arguments
+ * Use Euler's reflection formula to solve for arguments < -½
+ *
+ * Γ(z) = π/(sin(πz)Γ(1-z))
*/
- if (n < 0) {
- int i = -floor(n) + 1;
- real g = gamma(n + i);
-
- while (i-- > 0) {
- g = g / n;
- n++;
- }
- return g;
+ if (n < -0.5) {
+ int ibits = bits + 15;
+ real n_ibits = imprecise(n, ibits);
+ real pi_ibits = pi_value(ibits);
+ real sin_n_pi = sin(pi_ibits * n_ibits);
+ real inv_gamma = gamma(1 - n_ibits);
+ /* convert divide-by-zero into invalid_argument */
+ real div = sin_n_pi * inv_gamma;
+ if (div == 0)
+ raise invalid_argument("gamma of non-positive integer", 0, n);
+ return imprecise(pi_ibits / div, bits);
}
- n = imprecise(n);
- int bits = precision(n);
-
/*
* Smaller numbers take a lot more coefficients to
* get the desired number of bits. Make the value
@@ -102,7 +105,7 @@ extend namespace Math {
* to make the result converge faster.
*/
if (n < 10000) {
- int new_bits = bits + 14 - floor(log2(n));
+ int new_bits = bits + 20 - exponent(n);
real n_new = imprecise(n, new_bits) + 10000;
real g = gamma(n_new);
for (int i = 0; i < 10000; i++) {
commit c5d7fbcb2b70e9e556dae4e09e2c45944e4bc802
Author: Keith Packard <keithp at keithp.com>
Date: Sat Oct 4 16:51:21 2025 -0700
Use exponent/mantissa for log range reduction
This performs all of the log computations on the range of 1-2 instead of
an arbitrary range and makes log of large values much faster.
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/math.5c b/math.5c
index ec86bfa..27362df 100644
--- a/math.5c
+++ b/math.5c
@@ -441,6 +441,16 @@ extend namespace Math {
int prec = precision(a);
+ int e = 0;
+
+ /*
+ * Bring all values within the range of 1 <= a <= 2.
+ */
+ if (a > 2) {
+ e = exponent(a) - 1;
+ a = mantissa(a) * 2;
+ }
+
/*
* estimate = bsd_log (a). This gives 53 bits
*/
@@ -468,16 +478,6 @@ extend namespace Math {
v = imprecise (v, iprec);
a = imprecise (a, iprec);
- /* reduce to range 1-2 to improve convergence speed */
- real epow = imprecise(floor (v) - 2, iprec);
-
-
- /*
- * compute log(a) = log(a/(exp(epow))) + epow;
- */
- v = v - epow;
- a /= exp(epow);
-
/*
* Newton's method
*
@@ -497,9 +497,15 @@ extend namespace Math {
while (maxiter-- > 0)
v = (v - one) + a / exp(v);
+ }
- /* Add the range reduction value back in */
- v = v + epow;
+ /* Mix in the log of the exponent */
+ if (e != 0) {
+ int eprec = prec + 16;
+ static real log2 = 0;
+ if (log2 == 0 || precision(log2) < eprec)
+ log2 = log(imprecise(2, eprec));
+ v += imprecise(e, eprec) * log2;
}
return imprecise (v, prec);
commit 3f08c19478415605d56a3c19393aa4074233f9d7
Author: Keith Packard <keithp at keithp.com>
Date: Fri Oct 3 13:44:30 2025 -0700
Add a bunch more math functions
Hyperbolic trig, erf and jn
Signed-off-by: Keith Packard <keithp at keithp.com>
diff --git a/builtin.5c b/builtin.5c
index 318c2d8..e1d0159 100644
--- a/builtin.5c
+++ b/builtin.5c
@@ -120,6 +120,7 @@ library "math.5c";
library "prime_sieve.5c"
library "factorial.5c"
library "gamma.5c"
+library "math-ext.5c"
import Math;
library "scanf.5c";
library "socket.5c";
diff --git a/math-ext.5c b/math-ext.5c
new file mode 100644
index 0000000..24b7e19
--- /dev/null
+++ b/math-ext.5c
@@ -0,0 +1,162 @@
+/*
+ * Copyright © 2025 Keith Packard and Bart Massey.
+ * All Rights Reserved. See the file COPYING in this directory
+ * for licensing information.
+ */
+
+extend namespace Math {
+
+ public real acosh(x)
+ /*
+ * Returns inverse hyperbolic cosine of 'x'.
+ */
+ {
+ if (x < 1)
+ raise invalid_argument("acosh arg < 1", 0, x);
+ return log(x + sqrt(x*x-1));
+ }
+
+ public real asinh(x)
+ /*
+ * Returns inverse hyperbolic sine of 'x'.
+ */
+ {
+ if (x == 0) return 0;
+ real sign = 1;
+ if (x < 0) {
+ sign = -1;
+ x = -x;
+ }
+ return sign * log(x + sqrt(x*x+1));
+ }
+
+ public real atanh(x)
+ /*
+ * Returns inverse hyperbolic tangent of 'x'.
+ */
+ {
+ if (abs(x) > 1)
+ raise invalid_argument("atanh arg > 1", 0, x);
+ return 0.5 * log((1 + x) / (1 - x));
+ }
+
+ public real cosh(x)
+ /*
+ * Returns hyperbolic cosine of 'x'.
+ */
+ {
+ return (exp(x) + exp(-x)) / 2;
+ }
+
+ public real sinh(x)
+ /*
+ * Returns hyperbolic sine of 'x'.
+ */
+ {
+ return (exp(x) - exp(-x)) / 2;
+ }
+
+ public real tanh(x)
+ /*
+ * Returns hyperbolic tangent of 'x'.
+ */
+ {
+ return sinh(x) / cosh(x);
+ }
+
+ real _erf(real x, real off)
+ {
+ if (x == off)
+ return 0;
+ x = imprecise(x);
+ int bits = precision(x);
+ int obits = bits + 512;
+ real factor = 2 / sqrt(pi_value(obits));
+
+ x = imprecise(x, obits);
+ off = imprecise(off, obits) / factor;
+ real val = x - off;
+
+ for (int n = 1; ; n++) {
+ int f = 2 * n + 1;
+ real a = imprecise(((-1)**n * x**f) / (n! * f), obits);
+ val += a;
+ if (exponent(val) - exponent(a) > obits)
+ break;
+ }
+ return imprecise(val * factor, bits);
+ }
+
+ public real erf(real x)
+ /*
+ * Returns Gauss error function of 'x'.
+ */
+ {
+ return _erf(x, 0);
+ }
+
+ public real erfc(real x)
+ /*
+ * Returns Gauss complementary error function of 'x'.
+ */
+ {
+ return -_erf(x, 1);
+ }
+
+ public real jn(int n, real x)
+ /*
+ * Returns Bessel function of the first kind, order 'n' of 'x'
+ */
+ {
+ if (n < 0) {
+ real v = jn(-n, x);
+ if (n % 2 == 1)
+ v = -v;
+ return v;
+ }
+
+ x = imprecise(x);
+ int bits = precision(x);
+ int obits = bits + 64;
+
+ x = imprecise(x, obits);
+ real val = imprecise(0, obits);
+
+ int s = 1;
+
+ int mf = 1;
+ int mnf = n!;
+ real x_2 = x/2;
+ real x_2_pow = x_2 ** n;
+ real x_2_2 = x_2 ** 2;
+
+ for (int m = 0; ;) {
+ real a = (s / (mf * mnf)) * x_2_pow;
+ val += a;
+ if (exponent(val) - exponent(a) > obits)
+ break;
+ s = s > 0 ? -1 : 1;
+ m++;
+ mf *= m;
+ mnf *= (m+n);
+ x_2_pow *= x_2_2;
+ }
+ return imprecise(val, bits);
+ }
+
+ public real j0(real x)
+ /*
+ * Returns Bessel function of the first kind, order '0' of 'x'
+ */
+ {
+ return jn(0, x);
+ }
+
+ public real j1(real x)
+ /*
+ * Returns Bessel function of the first kind, order '1' of 'x'
+ */
+ {
+ return jn(1, x);
+ }
+}
diff --git a/meson.build b/meson.build
index 3606e73..35837d5 100644
--- a/meson.build
+++ b/meson.build
@@ -221,7 +221,7 @@ nickle_files = [
'printf.5c', 'history.5c', 'ctype.5c', 'string.5c', 'socket.5c',
'file.5c', 'parse-args.5c', 'svg.5c', 'process.5c',
'prime_sieve.5c', 'factorial.5c', 'gamma.5c', 'sort.5c', 'list.5c', 'skiplist.5c',
- 'json.5c', 'cha-cha.5c',
+ 'json.5c', 'cha-cha.5c', 'math-ext.5c',
]
nickle_headers = [
More information about the Nickle
mailing list