Math-Prime-Util

 view release on metacpan or  search on metacpan

lib/Math/Prime/Util.pm  view on Meta::CPAN

  # Exponential of Mangoldt function
  say "lamba(49) = ", log(exp_mangoldt(49));
  # Some more number theoretical functions
  say liouville(4292384);
  say chebyshev_psi(234984);
  say chebyshev_theta(92384234);
  say partitions(1000);
  # Show all prime partitions of 25
  forpart { say "@_" unless scalar grep { !is_prime($_) } @_ } 25;
  # List all 3-way combinations of an array
  my @cdata = qw/apple bread curry donut eagle/;
  forcomb { say "@cdata[@_]" } @cdata, 3;
  # or all permutations
  forperm { say "@cdata[@_]" } @cdata;

  # divisor sum
  my $sigma  = divisor_sum( $n );       # sum of divisors
  my $sigma0 = divisor_sum( $n, 0 );    # count of divisors
  my $sigmak = divisor_sum( $n, $k );
  my $sigmaf = divisor_sum( $n, sub { log($_[0]) } ); # arbitrary func

  # primorial n#, primorial p(n)#, and lcm
  say "The product of primes below 47 is ",     primorial(47);
  say "The product of the first 47 primes is ", pn_primorial(47);

lib/Math/Prime/Util/PP.pm  view on Meta::CPAN

        $tot[$i*$j] -= $tot[$i*$j]/$i;
      }
    }
    splice(@tot, 0, $lo) if $lo > 0;
    push @totients, @tot;
  }
  @totients;
}

sub _sumtot {
  my($n, $cdata, $ecache) = @_;
  return $cdata->[$n] if $n <= 0+$#$cdata;
  return $ecache->{$n} if defined $ecache->{$n};

  my $sum = Mmulint($n, $n+1) >> 1;
  my $s = sqrtint($n);
  my $lim = Mdivint($n, $s+1);

  my($x, $nextx) = ($n, Mdivint($n,2));
  $sum -= Mmulint($x - $nextx, $cdata->[1]);
  for my $k (2 .. $lim) {
    ($x,$nextx) = ($nextx, Mdivint($n,$k+1));
    $sum -= ($x <= 0+$#$cdata) ? $cdata->[$x] : _sumtot($x, $cdata, $ecache);
    $sum -= Mmulint($x - $nextx,
                    ($k <= $#$cdata) ? $cdata->[$k] : _sumtot($k, $cdata, $ecache));
  }
  if ($s > $lim) {
    ($x,$nextx) = ($nextx, Mdivint($n,$s+1));
    $sum -= Mmulint($x - $nextx,
                    ($s <= 0+$#$cdata) ? $cdata->[$s] : _sumtot($s, $cdata, $ecache));
  }
  $ecache->{$n} = $sum;
  $sum;
}

sub sumtotient {
  my($n) = @_;
  validate_integer_nonneg($n);
  return $n if $n <= 2;

totients.c  view on Meta::CPAN

/* Recursive method using a cache.  */

typedef struct {
  UV    hsize;
  UV   *nhash;  /* n value */
  UV   *shash;  /* sum for n */
} sumt_hash_t;
#define _CACHED_SUMT(x) \
  (((x)<csize)  ?  cdata[x]  :  _sumt((x), cdata, csize, thash))

static UV _sumt(UV n, const UV *cdata, UV csize, sumt_hash_t thash) {
  UV sum, s, k, lim, hn;
  uint32_t const probes = 8;
  uint32_t hashk;

  if (n < csize) return cdata[n];

  hn = n % thash.hsize;
  for (hashk = 0; hashk <= probes; hashk++) {
    if (thash.nhash[hn] == n)  return thash.shash[hn];
    if (thash.nhash[hn] == 0)  break;
    hn = (hn+1 < thash.hsize) ? hn+1 : 0;
  }

  s = isqrt(n);
  lim = n/(s+1);

totients.c  view on Meta::CPAN

  Safefree(thash.nhash);
  Safefree(thash.shash);
  Safefree(sumcache);
  return sum;
}



#if HAVE_SUMTOTIENT_128
#define _CACHED_SUMT128(x) \
  (((x)<csize)  ?  (uint128_t)cdata[x]  :  _sumt128((x), cdata, csize, thash))
typedef struct {
  UV         hsize;
  UV        *nhash;  /* n value */
  uint128_t *shash;  /* sum for n */
} sumt_hash_128_t;
static uint128_t _sumt128(UV n, const UV *cdata, UV csize, sumt_hash_128_t thash) {
  uint128_t sum;
  UV s, k, lim, hn;
  uint32_t const probes = 16;
  uint32_t const hinc = 1 + ((n >> 8) & 15);  /* mitigate clustering */
  uint32_t hashk;

  if (n < csize) return cdata[n];

  hn = n % thash.hsize;
  for (hashk = 0; hashk <= probes; hashk++) {
    if (thash.nhash[hn] == n)  return thash.shash[hn];
    if (thash.nhash[hn] == 0)  break;
    hn = (hn+hinc < thash.hsize) ? hn+hinc : hn+hinc-thash.hsize;
  }

  s = isqrt(n);
  lim = n/(s+1);