view release on metacpan or  search on metacpan

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

    }

    print "   ... inverse_totient phase 1 complete ...\n" if $_verbose;

    # 2. Process the divisors in reverse order.
    for my $dinfo (reverse @DIVINFO) {
      my($d,@L) = @$dinfo;
      my %todelete;
      my @T;
      # Multiply through by $pp
      for my $dset (@L) {
        if (defined $r{$dset->[0]}) {
          my($d2,$pp,$F) = @$dset;
          push @T, [$F, [map { Mmulint($pp,$_) } @{$r{$d2}}]];
          $todelete{$d2} = 1 if $needed{$d2} >= $d;
        }
      }
      # Delete intermediate data that isn't needed any more
      delete $r{$_} for keys %todelete;
      # Append the multiplied lists.
      push @{$r{$_->[0]}}, @{$_->[1]} for @T;
    }
    undef %needed;
    print "   ... inverse_totient phase 2 complete ...\n" if $_verbose;

    return (defined $r{$n}) ? @{Mvecsorti($r{$n})} : ();
  }
}

sub _euler_phi_range {
  my($lo, $hi) = @_;
  validate_integer($lo);
  validate_integer($hi);

  my @totients;
  while ($lo < 0 && $lo <= $hi) {
    push @totients, 0;
    $lo++;
  }
  return @totients if $hi < $lo;

  if ($hi > 2**30 || $hi-$lo < 100) {
    ($lo,$hi) = (tobigint($lo),tobigint($hi)) if $hi > 2**49;
    push @totients, euler_phi($lo++)  while $lo <= $hi;
  } else {
    my @tot = (0 .. $hi);
    foreach my $i (2 .. $hi) {
      next unless $tot[$i] == $i;
      $tot[$i] = $i-1;
      foreach my $j (2 .. int($hi / $i)) {
        $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;

  if ($n < 900) {     # Simple linear sum for small values.
    my $sum = 0;
    $sum += $_ for Mtotient(1,$n);
    return $sum;
  }

  my $cbrt = Mrootint($n,3);
  my $csize = Mvecprod(4, $cbrt, $cbrt);
  $csize = 50_000_000 if $csize > 50_000_000;  # Limit memory use to ~2.5GB
  my @sumcache = Mtotient(0,$csize);
  $sumcache[$_] += $sumcache[$_-1] for 2 .. $csize;
  _sumtot($n, \@sumcache, {});
}


sub prime_bigomega {
  my($n) = @_;
  validate_integer_abs($n);
  return scalar(Mfactor($n));
}
sub prime_omega {
  my($n) = @_;
  validate_integer_abs($n);
  return scalar(Mfactor_exp($n));
}

sub moebius {
  return _moebius_range(@_) if scalar @_ > 1;
  my($n) = @_;
  validate_integer_abs($n);
  return ($n == 1) ? 1 : 0  if $n <= 1;
  return 0 if ($n >= 49) && (!($n % 4) || !($n % 9) || !($n % 25) || !($n%49) );
  my @factors = Mfactor($n);
  foreach my $i (1 .. $#factors) {
    return 0 if $factors[$i] == $factors[$i-1];
  }
  return ((scalar @factors) % 2) ? -1 : 1;
}
sub is_square_free {
  return (Mmoebius($_[0]) != 0) ? 1 : 0;
}

sub is_odd {
  my($n) = @_;
  validate_integer($n);
  # Note:  If n is a Math::BigInt (Calc), then performance:
  #    0.25x  $n->is_odd()                           # method is fastest
  #    0.34x  (substr($n,-1,1) =~ tr/13579/13579/)   # Perl look at string
  #    0.46x  is_odd($n)                             # XS looks at the string
  #    1.0x   $n % 2 ? 1 : 0