Add squarefree_kernel

Changes

@@ -88,6 +88,7 @@ Revision history for Perl module Math::Prime::Util
     - powerfree_sum(n[,k])            sum of k-powerfree numbers <= n
     - powerfree_part(n[,k])           remove excess powers from n
     - powerfree_part_sum(n[,k])       sum of k-powerfree parts for 1 to n
+    - squarefree_kernel(n)            integer radical of |n|
     - perfect_power_count([lo,] hi)   count of perfect powers
     - smooth_count(n,k)               count of k-smooth numbers <= n
     - rough_count(n,k)                count of k-rough numbers <= n

XS.xs

@@ -1905,6 +1905,19 @@ void is_square(IN SV* svn)
     }
     return; /* skip implicit PUTBACK */
 
+void squarefree_kernel(IN SV* svn)
+  PREINIT:
+    int status;
+    UV n;
+  PPCODE:
+    status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
+    if (status == -1)
+      XSRETURN_IV( -(IV)squarefree_kernel(-(IV)n) );
+    if (status == 1)
+      XSRETURN_UV( squarefree_kernel(n) );
+    _vcallsub_with_pp("squarefree_kernel");
+    return;
+
 void is_powerfree(IN SV* svn, IN int k = 2)
   ALIAS:
     powerfree_sum = 1
@@ -3099,7 +3112,6 @@ void qnr(IN SV* svn)
     objectify_result(aTHX_ svn, ST(0));
     return;
 
-
 void
 is_smooth(IN SV* svn, IN SV* svk)
   ALIAS:

lib/Math/Prime/Util.pm

@@ -77,7 +77,7 @@ our @EXPORT_OK =
       perfect_power_count_lower perfect_power_count_upper
       nth_perfect_power nth_perfect_power_approx
       nth_perfect_power_lower nth_perfect_power_upper
-      nth_powerfree powerfree_count powerfree_sum
+      nth_powerfree powerfree_count powerfree_sum squarefree_kernel
       powerfree_part powerfree_part_sum
       smooth_count rough_count powersum
       lucky_count lucky_count_approx lucky_count_lower lucky_count_upper
@@ -2938,6 +2938,16 @@ but substantially faster.
 With C<k == 2> this produces the sequence
 L<OEIS A069891|http://oeis.org/A069891>.
 
+=head2 squarefree_kernel
+
+Given an integer C<n>, returns the squarefree kernel of C<n>.  This is
+also known as the integer radical.  It is the largest squarefree divisor
+of C<n>, which is also the product of the distinct primes dividing C<n>.
+
+We choose to accept negative inputs, with the result matching the input sign.
+
+This is the L<OEIS series A007947|http://oeis.org/A007947>.
+
 =head2 sqrtint
 
 Given a non-negative integer input C<n>, returns the integer square root.

lib/Math/Prime/Util/PP.pm

@@ -2103,6 +2103,13 @@ sub powerfree_part_sum {
   $sum;
 }
 
+sub squarefree_kernel {
+  my($n) = @_;
+  validate_integer($n);
+  return -1 * Mlcm(Mfactor(-$n)) if $n < 0;
+  Mlcm(Mfactor($n));
+}
+
 sub is_perfect_power {
   my($n) = @_;
   validate_integer($n);

lib/Math/Prime/Util/PPFE.pm

@@ -176,6 +176,7 @@ sub entropy_bytes {
 *powerfree_sum = \&Math::Prime::Util::PP::powerfree_sum;
 *powerfree_part = \&Math::Prime::Util::PP::powerfree_part;
 *powerfree_part_sum = \&Math::Prime::Util::PP::powerfree_part_sum;
+*squarefree_kernel = \&Math::Prime::Util::PP::squarefree_kernel;
 # TODO: Should this do validation here?
 *powersum = \&Math::Prime::Util::PP::powersum;
 *next_chen_prime = \&Math::Prime::Util::PP::next_chen_prime;

lib/ntheory.pm

@@ -388,6 +388,7 @@ Tags:
   powerfree_sum(n[,k])                sum of k-powerfree numbers <= n
   powerfree_part(n[,k])               remove excess powers so n is k-free
   powerfree_part_sum(n[,k])           sum of k-powerfree parts for 1 to n
+  squarefree_kernel(n)                integer radical of |n|
   powersum(n,k)                       sum of kth powers from 1 to n
 
 =head2 NON-INTEGER MATH

powerfree.c

@@ -268,3 +268,16 @@ UV nth_powerfree(UV n, uint32_t k)
   }
   return qk;
 }
+
+/******************************************************************************/
+
+UV squarefree_kernel(UV n)
+{
+  UV P, fac[MPU_MAX_FACTORS+1];
+  int i, nfactors;
+
+  nfactors = factor_exp(n, fac, 0);
+  for (P = 1, i = 0; i < nfactors; i++)
+    P *= fac[i];
+  return P;
+}

powerfree.h

@@ -12,4 +12,6 @@ extern UV powerfree_part_sum(UV n, uint32_t k);
 
 extern UV nth_powerfree(UV n, uint32_t k);
 
+extern UV  squarefree_kernel(UV n);
+
 #endif

t/02-can.t

@@ -71,7 +71,7 @@ my @functions =  qw(
       perfect_power_count_lower perfect_power_count_upper
       nth_perfect_power nth_perfect_power_approx
       nth_perfect_power_lower nth_perfect_power_upper
-      nth_powerfree powerfree_count powerfree_sum
+      nth_powerfree powerfree_count powerfree_sum squarefree_kernel
       powerfree_part powerfree_part_sum
       smooth_count rough_count powersum
       lucky_count lucky_count_approx lucky_count_lower lucky_count_upper

t/26-powerfree.t

@@ -5,18 +5,32 @@ use warnings;
 use Test::More;
 use Math::Prime::Util qw/is_powerfree powerfree_count powerfree_sum
                          powerfree_part powerfree_part_sum nth_powerfree
-                         is_square_free vecsum vecmax factor_exp/;
+                         squarefree_kernel is_square_free
+                         vecsum vecmax factor_exp/;
 
 my @simple = (0 .. 16,
               758096738,434420340,870589313,695486396,602721315,418431087,
               752518565,723570005,506916483,617459403);
 
+my @T = (             # powerfree part, squarefree_kernel
+ [0, [0, 0]],
+ [1, [1, 1]],
+ [2, [2, 2]],
+ [ 48, [3, 6]],
+ [-48, [-3, -6]],
+ [2*2*3*5, [3*5, 2*3*5]],
+ [2*3*3*3*5*5*7, [2*3*7, 2*3*5*7]],
+ ["54713282649239", [5471, 547116413]],
+ ["4000000000000027", ["4162330905307", "129032258064517"]],
+);
+
 plan tests => 3     # simple is square free
             + 11*2  # powerfree_count, powerfree_sum
             + 6+2   # ""
             + 7     # nth_powerfree
             + 2+8   # powerfree_part
-            + 8*2;  # powerfree_part_sum
+            + 8*2   # powerfree_part_sum
+            + 2;    # powerfree_part and squarefree_kernel
 
 ##### is_powerfree
 
@@ -98,6 +112,11 @@ for my $k (0 .. 7) {
   is( powerfree_part_sum(654321,$k), $pfpst[$k], "powerfree_part_sum(654321,$k) = $pfpst[$k]" );
 }
 
+##### powerfree_part and squarefree_kernel
+
+is_deeply( [map { powerfree_part($_->[0]) } @T], [map { $_->[1]->[0] } @T], "powerfree_part" );
+is_deeply( [map { squarefree_kernel($_->[0]) } @T], [map { $_->[1]->[1] } @T], "squarefree_kernel" );
+
 ##### subs
 
 sub ipf {

t/92-release-pod-coverage.t

@@ -104,7 +104,7 @@ sub mpu_public_regex {
       perfect_power_count_lower perfect_power_count_upper
       nth_perfect_power nth_perfect_power_approx
       nth_perfect_power_lower nth_perfect_power_upper
-      nth_powerfree powerfree_count powerfree_sum
+      nth_powerfree powerfree_count powerfree_sum squarefree_kernel
       powerfree_part powerfree_part_sum
       smooth_count rough_count powersum
       lucky_count lucky_count_approx lucky_count_lower lucky_count_upper