- Added the Number bphi(n) method.

The bi-unitary analog of Euler's totient function of `n`. (OEIS: A116550)
trizen · Oct 22, 2023 · a121ecd · a121ecd
1 parent e662b70
commit a121ecd
Show file tree

Hide file tree

Showing 4 changed files with 45 additions and 1 deletion.
diff --git a/NUMBER_THEORY_TUTORIAL.md b/NUMBER_THEORY_TUTORIAL.md
@@ -663,6 +663,7 @@ say 20.by { .is_proth_prime(3) }   # numbers n such that 2^n * 3 + 1 is prime
 say map(1..50, { .psi })
 say map(1..50, { .phi })
 say map(1..50, { .iphi })
+say map(1..50, { .bphi })
 say map(1..50, { .uphi })
 say map(1..50, { .nuphi })
 say map(1..50, { .sigma })

diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -24081,7 +24081,7 @@ package Sidef::Types::Number::Number {
                     Math::GMPz::Rmpz_ui_pow_ui($pp, $p, $e);
                 }
                 else {
-                    Math::GMPz::Rmpz_set_str($pp, $p, 10);
+                    Math::GMPz::Rmpz_set_str($pp, "$p", 10);
                     Math::GMPz::Rmpz_pow_ui($pp, $pp, $e);
                 }
             }
@@ -26235,6 +26235,40 @@ package Sidef::Types::Number::Number {
         bless \_binsplit(\@terms, \&__mul__);
     }
 
+    sub bphi {    # OEIS: A116550 -- bi-unitary analog of Euler's totient function of n.
+        my ($n) = @_;
+
+        # Formula:
+        #   a(n) = Sum{d|n, d is unitary} (-1)^omega(d) * Sum{k|d, k is squarefree, k <= n/d} mu(k) * floor(n/(d*k)).
+
+        my @final_terms;
+
+        foreach my $v (@{$n->udivisors}) {
+            my $sign = (-1)**(scalar _factor_exp($$v));
+            my $x    = $n->idiv($v);
+            my $nod  = "$$x";
+
+            if ($nod eq '1') {
+                push @final_terms, $sign;
+                next;
+            }
+
+            my @terms;
+            foreach my $u (@{$v->squarefree_divisors}) {
+                my $k   = $$u;
+                my $div = Math::Prime::Util::GMP::divint($nod, $k);
+                last if ($div eq '0');
+                my $mu = ((HAS_PRIME_UTIL and !ref($k)) ? Math::Prime::Util::moebius($k) : Math::Prime::Util::GMP::moebius($k));
+                push @terms, Math::Prime::Util::GMP::mulint($mu, $div);
+            }
+
+            push @final_terms, Math::Prime::Util::GMP::mulint($sign, ((scalar(@terms) == 1) ? $terms[0] : Math::Prime::Util::GMP::vecsum(@terms)));
+        }
+
+        @final_terms || return ZERO;
+        _set_int(Math::Prime::Util::GMP::vecsum(@final_terms));
+    }
+
     sub prime_power_sigma {
         my ($n, $k) = @_;
 

diff --git a/lib/Sidef/Types/Number/Number.pod b/lib/Sidef/Types/Number/Number.pod
@@ -1225,6 +1225,14 @@ Returns C<nil> if C<n> cannot be truncated to an integer or if C<k> is negative.
 
 =cut
 
+=head2 bphi
+
+    bphi(n)
+
+The bi-unitary analog of Euler's totient function of C<n>. (OEIS: A116550)
+
+=cut
+
 =head2 bsearch
 
     bsearch(n, {...})

diff --git a/scripts/Tests/number_methods.sf b/scripts/Tests/number_methods.sf
@@ -1373,5 +1373,6 @@ for n in ([@|(0..20), 120, 480, 5040, (2**64 + 1)**3, 2**64 + 1, 2**127 - 1, (2*
 }
 
 assert_eq(map(1..70, {.nuphi}), %n[1,2,3,2,5,6,7,4,6,10,11,6,13,14,15,8,17,12,19,10,21,22,23,12,20,26,18,14,29,30,31,16,33,34,35,12,37,38,39,20,41,42,43,22,30,46,47,24,42,40,51,26,53,36,55,28,57,58,59,30,61,62,42,32,65,66,67,34,69,70])
+assert_eq(map(1..74, {.bphi}), %n[1,1,2,3,4,3,6,7,8,6,10,8,12,9,9,15,16,12,18,14,14,15,22,17,24,18,26,21,28,15,30,31,23,24,25,29,36,27,28,31,40,21,42,35,34,33,46,36,48,36,37,42,52,39,42,46,42,42,58,34,60,45,51,63,50,35,66,56,51,38,70,62,72,54])
 
 say "** Tests passed!"