Add C extension for Chebyshev polynomial evaluation

.gitignore

@@ -9,3 +9,13 @@
 
 # rspec failure tracking
 .rspec_status
+
+# C extension build artifacts
+/lib/ephem/*.so
+/lib/ephem/*.bundle
+/lib/ephem/*.dll
+/ext/ephem/chebyshev/Makefile
+/ext/ephem/chebyshev/*.o
+/ext/ephem/chebyshev/*.bundle
+/ext/ephem/chebyshev/*.so
+/ext/ephem/chebyshev/*.dSYM

Gemfile.lock

@@ -38,6 +38,8 @@ GEM
     racc (1.8.1)
     rainbow (3.1.1)
     rake (13.3.1)
+    rake-compiler (1.3.1)
+      rake
     rdoc (6.13.1)
       psych (>= 4.0.0)
     regexp_parser (2.11.3)
@@ -104,6 +106,7 @@ DEPENDENCIES
   irb (~> 1.15)
   parallel (~> 1.26)
   rake (~> 13.0)
+  rake-compiler (~> 1.2)
   rspec (~> 3.13)
   standard (~> 1.43)
 

Rakefile

@@ -9,4 +9,17 @@ RSpec::Core::RakeTask.new(:spec)
 
 require "standard/rake"
 
+begin
+  require "rake/extensiontask"
+
+  Rake::ExtensionTask.new("chebyshev") do |ext|
+    ext.lib_dir = "lib/ephem"
+    ext.ext_dir = "ext/ephem/chebyshev"
+  end
+
+  task spec: :compile
+rescue LoadError
+  # rake-compiler not available; skip extension compilation tasks.
+end
+
 task default: %i[spec standard]

ephem.gemspec

@@ -40,6 +40,7 @@ Gem::Specification.new do |spec|
   spec.bindir = "bin"
   spec.executables = ["ruby-ephem"]
   spec.require_paths = ["lib"]
+  spec.extensions = ["ext/ephem/chebyshev/extconf.rb"]
 
   spec.add_dependency "minitar", "~> 0.12"
   spec.add_dependency "numo-narray", "~> 0.9.2.1"
@@ -51,5 +52,6 @@ Gem::Specification.new do |spec|
   spec.add_development_dependency "rake", "~> 13.0"
   spec.add_development_dependency "rspec", "~> 3.13"
   spec.add_development_dependency "benchmark-ips", "~> 2.14"
+  spec.add_development_dependency "rake-compiler", "~> 1.2"
   spec.add_development_dependency "standard", "~> 1.43"
 end

ext/ephem/chebyshev/chebyshev.c

@@ -0,0 +1,143 @@
+#include <ruby.h>
+
+/* Matches Ephem::Core::Constants::Time::SECONDS_PER_DAY */
+#define SECONDS_PER_DAY 86400.0
+
+/*
+ * Ephem::Computation::ChebyshevPolynomial.evaluate(coeffs, t)
+ *
+ * Evaluates a 3D Chebyshev polynomial at a given normalized time
+ * using the Clenshaw recurrence algorithm.
+ *
+ * coeffs: Array of Arrays, shape [n_terms][3]
+ * t:      Float in [-1, 1]
+ *
+ * Returns: Array of 3 Floats [x, y, z]
+ */
+static VALUE
+chebyshev_evaluate(VALUE self, VALUE rb_coeffs, VALUE rb_t)
+{
+    long n, k;
+    double t, t2;
+    double b1x, b1y, b1z;
+    double b2x, b2y, b2z;
+    double c0, c1, c2;
+    double tx, ty, tz;
+    VALUE c_ary;
+
+    Check_Type(rb_coeffs, T_ARRAY);
+    t = NUM2DBL(rb_t);
+
+    n = RARRAY_LEN(rb_coeffs);
+
+    b1x = b1y = b1z = 0.0;
+    b2x = b2y = b2z = 0.0;
+
+    t2 = 2.0 * t;
+
+    for (k = n - 1; k > 0; k--) {
+        c_ary = rb_ary_entry(rb_coeffs, k);
+        Check_Type(c_ary, T_ARRAY);
+        c0 = NUM2DBL(rb_ary_entry(c_ary, 0));
+        c1 = NUM2DBL(rb_ary_entry(c_ary, 1));
+        c2 = NUM2DBL(rb_ary_entry(c_ary, 2));
+
+        tx = t2 * b1x - b2x + c0;
+        ty = t2 * b1y - b2y + c1;
+        tz = t2 * b1z - b2z + c2;
+
+        b2x = b1x; b2y = b1y; b2z = b1z;
+        b1x = tx;  b1y = ty;  b1z = tz;
+    }
+
+    c_ary = rb_ary_entry(rb_coeffs, 0);
+    Check_Type(c_ary, T_ARRAY);
+    c0 = NUM2DBL(rb_ary_entry(c_ary, 0));
+    c1 = NUM2DBL(rb_ary_entry(c_ary, 1));
+    c2 = NUM2DBL(rb_ary_entry(c_ary, 2));
+
+    return rb_ary_new_from_args(3,
+        DBL2NUM(t * b1x - b2x + c0),
+        DBL2NUM(t * b1y - b2y + c1),
+        DBL2NUM(t * b1z - b2z + c2));
+}
+
+/*
+ * Ephem::Computation::ChebyshevPolynomial.evaluate_derivative(coeffs, t, radius)
+ *
+ * Evaluates the time derivative of a 3D Chebyshev polynomial
+ * using the Clenshaw recurrence algorithm.
+ *
+ * coeffs: Array of Arrays, shape [n_terms][3]
+ * t:      Float in [-1, 1]
+ * radius: Float (half-interval in days)
+ *
+ * Returns: Array of 3 Floats [vx, vy, vz] in units per second
+ */
+static VALUE
+chebyshev_evaluate_derivative(VALUE self, VALUE rb_coeffs, VALUE rb_t,
+                              VALUE rb_radius)
+{
+    long n, k;
+    double t, t2, radius, scale, k2;
+    double d1x, d1y, d1z;
+    double d2x, d2y, d2z;
+    double c0, c1, c2;
+    double tx, ty, tz;
+    VALUE c_ary;
+
+    Check_Type(rb_coeffs, T_ARRAY);
+    t = NUM2DBL(rb_t);
+    radius = NUM2DBL(rb_radius);
+
+    n = RARRAY_LEN(rb_coeffs);
+
+    if (n < 2) {
+        return rb_ary_new_from_args(3,
+            DBL2NUM(0.0), DBL2NUM(0.0), DBL2NUM(0.0));
+    }
+
+    d1x = d1y = d1z = 0.0;
+    d2x = d2y = d2z = 0.0;
+
+    t2 = 2.0 * t;
+
+    for (k = n - 1; k > 0; k--) {
+        c_ary = rb_ary_entry(rb_coeffs, k);
+        Check_Type(c_ary, T_ARRAY);
+        c0 = NUM2DBL(rb_ary_entry(c_ary, 0));
+        c1 = NUM2DBL(rb_ary_entry(c_ary, 1));
+        c2 = NUM2DBL(rb_ary_entry(c_ary, 2));
+
+        k2 = 2.0 * (double)k;
+        tx = t2 * d1x - d2x + k2 * c0;
+        ty = t2 * d1y - d2y + k2 * c1;
+        tz = t2 * d1z - d2z + k2 * c2;
+
+        d2x = d1x; d2y = d1y; d2z = d1z;
+        d1x = tx;  d1y = ty;  d1z = tz;
+    }
+
+    scale = SECONDS_PER_DAY / (2.0 * radius);
+
+    return rb_ary_new_from_args(3,
+        DBL2NUM(d1x * scale),
+        DBL2NUM(d1y * scale),
+        DBL2NUM(d1z * scale));
+}
+
+void
+Init_chebyshev(void)
+{
+    VALUE mEphem, mComputation, mChebyshevPolynomial;
+
+    mEphem = rb_define_module("Ephem");
+    mComputation = rb_define_module_under(mEphem, "Computation");
+    mChebyshevPolynomial = rb_define_module_under(mComputation,
+                                                  "ChebyshevPolynomial");
+
+    rb_define_module_function(mChebyshevPolynomial, "evaluate",
+                              chebyshev_evaluate, 2);
+    rb_define_module_function(mChebyshevPolynomial, "evaluate_derivative",
+                              chebyshev_evaluate_derivative, 3);
+}

ext/ephem/chebyshev/extconf.rb

@@ -0,0 +1,11 @@
+# frozen_string_literal: true
+
+require "mkmf"
+
+# Disable FMA (fused multiply-add) contraction so that the C extension
+# produces bit-identical results to the pure-Ruby implementation.
+# Without this, ARM64 compilers emit fmadd instructions that skip
+# intermediate rounding, causing ULP-level differences.
+append_cflags("-ffp-contract=off")
+
+create_makefile("ephem/chebyshev")

lib/ephem/computation/chebyshev_polynomial.rb

@@ -88,3 +88,12 @@ def self.evaluate_derivative(coeffs, t, radius)
     end
   end
 end
+
+# Attempt to load the C extension, which redefines the methods above.
+# If compilation is not available (e.g., JRuby, TruffleRuby, or missing
+# compiler), the pure-Ruby implementation remains in place.
+begin
+  require "ephem/chebyshev"
+rescue LoadError
+  # C extension not available; pure-Ruby fallback is already loaded above.
+end