diff --git a/.gitignore b/.gitignore index b04a8c8..abe52eb 100644 --- a/.gitignore +++ b/.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 diff --git a/Gemfile.lock b/Gemfile.lock index b81f8fd..07886ce 100644 --- a/Gemfile.lock +++ b/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) diff --git a/Rakefile b/Rakefile index f0b2328..7144d90 100644 --- a/Rakefile +++ b/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] diff --git a/ephem.gemspec b/ephem.gemspec index d4a59bd..612d014 100644 --- a/ephem.gemspec +++ b/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 diff --git a/ext/ephem/chebyshev/chebyshev.c b/ext/ephem/chebyshev/chebyshev.c new file mode 100644 index 0000000..d555cc6 --- /dev/null +++ b/ext/ephem/chebyshev/chebyshev.c @@ -0,0 +1,143 @@ +#include + +/* 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); +} diff --git a/ext/ephem/chebyshev/extconf.rb b/ext/ephem/chebyshev/extconf.rb new file mode 100644 index 0000000..89e9658 --- /dev/null +++ b/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") diff --git a/lib/ephem/computation/chebyshev_polynomial.rb b/lib/ephem/computation/chebyshev_polynomial.rb index c6973fa..2875224 100644 --- a/lib/ephem/computation/chebyshev_polynomial.rb +++ b/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