- Extended Number harmonic and harmreal method to accept an optio…

…nal argument, returning the k-th Harmonic number. Example: say harmonic(100, 5) # 100-th Harmonic number 5-th order
trizen · Jun 6, 2021 · 3ee4ba4 · 3ee4ba4
1 parent b3f38b2
commit 3ee4ba4
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 4 deletions.
diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -4655,7 +4655,17 @@ package Sidef::Types::Number::Number {
     *bernoulli_log = \&lnbernreal;
 
     sub harmfrac {
-        my ($n) = @_;
+        my ($n, $k) = @_;
+
+        # Formula in terms of the Harmonic numbers, due to Conway and Guy (1996),
+        # for computing the Harmonic numbers of the k-th order.
+        if (defined($k)) {
+
+            my $km1   = $k->dec;
+            my $npkm1 = $n->add($km1);
+
+            return $npkm1->binomial($km1)->mul($npkm1->harmfrac->sub($km1->harmfrac));
+        }
 
         $n = _any2ui($$n) // goto &nan;
         $n || return ZERO();
@@ -4673,12 +4683,22 @@ package Sidef::Types::Number::Number {
     *harmonic_number = \&harmfrac;
 
     sub harmreal {
-        my ($x) = @_;
+        my ($n, $k) = @_;
 
-        $x = _any2mpfr($$x);
+        # Formula in terms of the Harmonic numbers, due to Conway and Guy (1996),
+        # for computing the Harmonic numbers of the k-th order.
+        if (defined($k)) {
+
+            my $km1   = $k->dec;
+            my $npkm1 = $n->add($km1);
+
+            return $npkm1->binomial($km1)->mul($npkm1->harmreal->sub($km1->harmreal));
+        }
+
+        $n = _any2mpfr($$n);
 
         my $r = Math::MPFR::Rmpfr_init2(CORE::int($PREC));
-        Math::MPFR::Rmpfr_add_ui($r, $x, 1, $ROUND);
+        Math::MPFR::Rmpfr_add_ui($r, $n, 1, $ROUND);
         Math::MPFR::Rmpfr_digamma($r, $r, $ROUND);
 
         my $t = Math::MPFR::Rmpfr_init2(CORE::int($PREC));

diff --git a/lib/Sidef/Types/Number/Number.pod b/lib/Sidef/Types/Number/Number.pod
@@ -2585,23 +2585,29 @@ Returns the Hamming distance (number of bit-positions where the bits differ) bet
 =head2 harm
 
     harmonic(n)
+    harmonic(n, k)
 
 Returns the n-th Harmonic number C<H_n>. The harmonic numbers are the sum of reciprocals of the first C<n> natural numbers: C<1 + 1/2 + 1/3 + ... + 1/n>.
 
+When an additional argument is given, it returns the n-th Harmonic number of the k-th order.
+
 Aliases: I<harmfrac>, I<harmonic>, I<harmonic_number>
 
 =cut
 
 =head2 harmreal
 
     harmreal(n)
+    harmreal(n, k)
 
 Returns the n-th Harmonic number C<H_n> as a floating-point value, defined as:
 
     harmreal(n) = digamma(n+1) + γ
 
 where C<γ> is the Euler-Mascheroni constant.
 
+When an additional argument is given, it returns the n-th Harmonic number of the k-th order.
+
 =cut
 
 =head2 hermite_H