- Added the Number faulhaber_range(a, b, k) method.

Sum of powers: `a^k + (a+1)^k + (a+2)^k + ... + b^k`, using Faulhaber's summation formula. The value for `k` must be a non-negative native integer. Example: faulhaber_range(50, 100, 2) # 50^2 + 51^2 + ... + 100^2 = 297925
trizen · May 16, 2021 · 99e7ca4 · 99e7ca4
1 parent 0f7b113
commit 99e7ca4
Show file tree

Hide file tree

Showing 3 changed files with 37 additions and 2 deletions.
diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -9235,6 +9235,13 @@ package Sidef::Types::Number::Number {
         bless \__div__(__sub__(__pow__($r, __add__($n, $ONE)), $ONE), __sub__($r, $ONE));
     }
 
+    sub faulhaber_range {
+        my ($from, $to, $k) = @_;
+        _valid(\$to, \$k);
+        return ZERO if $to->lt($from);
+        return $to->faulhaber_sum($k)->sub($from->dec->faulhaber_sum($k));
+    }
+
     sub faulhaber_sum {
         my ($n, $p) = @_;
 
@@ -17279,7 +17286,7 @@ package Sidef::Types::Number::Number {
         }
 
         my $n = _any2mpz($$from) // return ZERO;
-        $k = _any2ui($$k)  // return ZERO;
+        $k = _any2ui($$k) // return ZERO;
 
         if ($k < 2 or Math::GMPz::Rmpz_sgn($n) <= 0) {
             return ZERO;
@@ -17367,7 +17374,7 @@ package Sidef::Types::Number::Number {
         }
 
         my $n = _any2mpz($$from) // return ZERO;
-        $k = _any2ui($$k)  // return ZERO;
+        $k = _any2ui($$k) // return ZERO;
 
         if (Math::GMPz::Rmpz_sgn($n) <= 0) {
             return ZERO;

diff --git a/lib/Sidef/Types/Number/Number.pod b/lib/Sidef/Types/Number/Number.pod
@@ -2415,6 +2415,20 @@ Defined in terms of the Bernoulli polynomials, as:
 
 =cut
 
+=head2 faulhaber_range
+
+    faulhaber_range(a, b, k)
+
+Sum of powers: C<a^k + (a+1)^k + (a+2)^k + ... + b^k>, using Faulhaber's summation formula.
+
+The value for C<k> must be a non-negative native integer.
+
+Example:
+
+    faulhaber_range(50, 100, 2)   # 50^2 + 51^2 + ... + 100^2 = 297925
+
+=cut
+
 =head2 fermat_factor
 
     n.fermat_factor(k=1e4)

diff --git a/scripts/Tests/quaternion_integers.sf b/scripts/Tests/quaternion_integers.sf
@@ -116,6 +116,20 @@ do {
     assert_eq(20.of { b(_).abs**2 -> round }, %n[1, 1, 2, 5, 9, 20, 41, 85, 178, 369, 769, 1600, 3329, 6929, 14418, 30005, 62441, 129940, 270409, 562725])
 }
 
+do {
+    func b(n) is cached {
+        return 1 if (n <= 1)
+        b(n-1) + b(n-2)*Quaternion(0, 1/Quadratic(0, 1, 3), 1/Quadratic(0, 1, 3), 1/Quadratic(0, 1, 3))
+    }
+
+    # OEIS: A105309
+    assert_eq(20.of { b(_).norm }, %n[1, 1, 2, 5, 9, 20, 41, 85, 178, 369, 769, 1600, 3329, 6929, 14418, 30005, 62441, 129940, 270409, 562725])
+
+    assert_eq(b(151).norm, 579074356283752148309391378541196518313641586965)
+    assert_eq(b(240).norm, 12282939451223766602706295263727875255095274997784797700575322884638891494401)
+    assert_eq(b(300).norm, 153413999717187674323079946300077491667777701038030587210963604285451850164508905918344669013249)
+}
+
 assert_eq(  # OEIS: A213421
     with(Quaternion(2, 1, 1, 1)) {|q| 20.of {|n| q**n -> a } },
     %n[1, 2, 1, -10, -47, -118, -143, 254, 2017, 6290, 11041, 134, -76751, -307942, -694511, -622450, 2371777, 13844258, 38774593, 58188566]