[proposal] functions for k-powerfree numbers

[trizen]

Hi. I would like to suggest several new functions for k-powerfree numbers.

Definition:

• a positive integer is k-powerfree iff all the exponents e in the prime-power factorization of n are less than k. (i.e.: A051903(n) < k)

Example:

• k = 2 gives the squarefree numbers (OEIS: A005117)
• k = 3 gives the cubefree numbers (OEIS: A004709)

Proposed functions:

is_powerfree(n, k=2)         true if n is k-powerfree
powerfree_count(n, k=2)      count of k-powerfree numbers <= n

Additional functions:

powerfree_sum(n, k=2)          sum of k-powerfree numbers <= n
powerfree_part(n, k=2)         k-powerfree part of n
powerfree_part_sum(n, k=2)     sum of k-powerfree part of integers in the range 1..n

Formulas:

powerfree_count(n, k)    = Sum_{v=1..floor(n^(1/k))} moebius(v) * floor(n / v^k)
powerfree_sum(n, k)      = Sum_{v=1..floor(n^(1/k))} moebius(v) * v^k * T(floor(n / v^k)), where T(n) = n*(n+1)/2
powerfree_part(n, k)     = multiplificative with powerfree_part(p^e, k) = p^(e mod k)
powerfree_part_sum(n, k) = Sum_{v=1..floor(n^(1/k))} f(v,k) * T(floor(n / v^k)), where f(n,k) = Prod_{p|n} (1 - p^k)

Code example:

use 5.020;
use strict;
use warnings;

use ntheory qw(:all);
use experimental qw(signatures);

sub T ($n) {    # n-th triangular number
    divint(mulint($n, addint($n, 1)), 2);
}

sub is_powerfree ($n, $k = 2) {
    (vecall { $_->[1] < $k } factor_exp($n)) ? 1 : 0;
}

sub powerfree_count ($n, $k = 2) {
    my $count = 0;
    forsquarefree {
        $count += ((scalar(@_) & 1) ? -1 : 1) * divint($n, powint($_, $k));
    } rootint($n, $k);
    return $count;
}

sub powerfree_sum ($n, $k = 2) {
    my $sum = 0;
    forsquarefree {
        $sum = addint($sum, vecprod(((scalar(@_) & 1) ? -1 : 1), powint($_, $k), T(divint($n, powint($_, $k)))));
    } rootint($n, $k);
    return $sum;
}

sub powerfree_part ($n, $k = 2) {
    return 0 if ($n == 0);
    vecprod(map { powint($_->[0], $_->[1] % $k) } factor_exp($n));
}

sub f ($n, $r) {
    vecprod(map { 1 - powint($_->[0], $r) } factor_exp($n));
}

sub powerfree_part_sum ($n, $k = 2) {
    my $sum = 0;
    for (1 .. rootint($n, $k)) {
        $sum = addint($sum, mulint(f($_, $k), T(divint($n, powint($_, $k)))));
    }
    return $sum;
}

say join ', ', map { is_powerfree($_, 2) } 1..30;   # A008966
say join ', ', map { is_powerfree($_, 3) } 1..30;   # A212793

say '';

say join ', ', grep { is_powerfree($_, 2) } 1..30;   # A005117
say join ', ', grep { is_powerfree($_, 3) } 1..30;   # A004709
say join ', ', grep { is_powerfree($_, 4) } 1..30;   # A046100

say '';

say join ', ', map { powerfree_count($_, 2) } 1..30;    # A013928
say join ', ', map { powerfree_count($_, 3) } 1..30;    # A060431

say '';

say join ', ', map { powerfree_sum($_, 2) } 1..30;    # A066779
say join ', ', map { powerfree_sum($_, 3) } 1..30;    # not in the OEIS

say '';

say join ', ', map { powerfree_part($_, 2) } 1..30;    # A007913
say join ', ', map { powerfree_part($_, 3) } 1..30;    # A050985
say join ', ', map { powerfree_part($_, 4) } 1..30;    # A053165

say '';

say join ', ', map { powerfree_part_sum($_, 2) } 1..30;    # A069891
say join ', ', map { powerfree_part_sum($_, 3) } 1..30;    # not in the OEIS

say '';

use Test::More tests => 30;

foreach my $k (1..10) {
    my $n = 500;

    is_deeply(
        [map { powerfree_sum($_, $k) } 1..$n],
        [map { vecsum(grep { is_powerfree($_, $k) } 1..$_) } 1..$n],
    );

    is_deeply(
        [map { powerfree_count($_, $k) } 1..$n],
        [map { scalar grep { is_powerfree($_, $k) } 1..$_ } 1..$n],
    );

    is_deeply(
        [map { powerfree_part_sum($_, $k) } 1..$n],
        [map { vecsum(map { powerfree_part($_, $k) } 1..$_) } 1..$n],
    );
}

[danaj]

I just checked in is_powerfree, powerfree_count, and powerfree_sum.

I will add powerfree_part. I'm not sure on the sum.

[trizen]

Thank you! Awesome work as always.

[danaj]

Added powerfree_part and powerfree_part_sum.