File Coverage

blib/lib/Math/NumSeq/PrimeFactorCount.pm

Criterion	Covered	Total	%
statement	88	112	78.5
branch	42	58	72.4
condition	10	17	58.8
subroutine	23	24	95.8
pod	5	5	100.0
total	168	216	77.7

# Copyright 2010, 2011, 2012, 2013, 2014 Kevin Ryde

# This file is part of Math-NumSeq.
#
# Math-NumSeq is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3, or (at your option) any later
# version.
#
# Math-NumSeq is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License along
# with Math-NumSeq.  If not, see .


# values_type => 'mod2'


package Math::NumSeq::PrimeFactorCount;
use 5.004;
use strict;
use List::Util 'min', 'max';

use vars '$VERSION','@ISA';
$VERSION = 71;

use Math::NumSeq;
use Math::NumSeq::Base::IterateIth;
@ISA = ('Math::NumSeq::Base::IterateIth',
        'Math::NumSeq');
*_is_infinite = \&Math::NumSeq::_is_infinite;

use Math::Prime::XS 'is_prime';
use Math::Factor::XS 'prime_factors';

use Math::NumSeq::Fibonacci;
*_blog2_estimate = \&Math::NumSeq::Fibonacci::_blog2_estimate;

# uncomment this to run the ### lines
#use Smart::Comments;


# cf. Untouchables, not sum of proper divisors of any other integer
# p*q sum S=1+p+q
# so sums up to hi need factorize to (hi^2)/4
#

use constant values_min => 0;
use constant i_start => 1;

sub values_max {
  my ($self) = @_;
  if ($self->{'values_type'} eq 'mod2') {
    return 1;
  } else {
    return undef;
  }
}
use constant characteristic_count => 1;
use constant characteristic_smaller => 1;
use constant characteristic_increasing => 0;

use constant parameter_info_array =>
  [
   { name    => 'prime_type',
     display => Math::NumSeq::__('Prime Type'),
     type    => 'enum',
     default => 'all',
     choices => ['all','odd','4k+1','4k+3',
                 'twin','SG','safe'],
     choices_display => [Math::NumSeq::__('All'),
                         Math::NumSeq::__('Odd'),
                         # TRANSLATORS: "4k+1" meaning numbers 1,5,9,13 etc, probably no need to translate except into another script if Latin letter "k" won't be recognised
                         Math::NumSeq::__('4k+1'),
                         Math::NumSeq::__('4k+3'),
                         Math::NumSeq::__('Twin'),
                         Math::NumSeq::__('SG'),
                         Math::NumSeq::__('Safe'),
                        ],
     description => Math::NumSeq::__('The type of primes to count.
twin=P where P+2 or P-2 also prime.
SG=Sophie Germain P where 2P+1 also prime.
safe=P where (P-1)/2 also prime (the "other" of the SGs).'),
   },
   { name    => 'multiplicity',
     display => Math::NumSeq::__('Multiplicity'),
     type    => 'enum',
     default => 'repeated',
     choices => ['repeated','distinct'],
     choices_display => [Math::NumSeq::__('Repeated'),
                         Math::NumSeq::__('Distinct'),
                        ],
     description => Math::NumSeq::__('Whether to include repeated prime factors, or only distinct prime factors.'),
   },

   # not documented yet
   { name    => 'values_type',
     share_key => 'values_type_cm2',
     display => Math::NumSeq::__('Values Type'),
     type    => 'enum',
     default => 'count',
     choices => ['count','mod2'],
     choices_display => [Math::NumSeq::__('Count'),
                         Math::NumSeq::__('Mod2'),
                        ],
     # description => Math::NumSeq::__('...'),
   },
  ];

sub description {
  my ($self) = @_;
  if (ref $self) {
    return ($self->{'multiplicity'} eq 'repeated'
            ? Math::NumSeq::__('Count of prime factors, including repetitions.')
            : Math::NumSeq::__('Count of distinct prime factors.'))
      . ($self->{'prime_type'} eq 'odd' ? "\nOdd primes only."
         : $self->{'prime_type'} eq '4k+1' ? "\nPrimes of form 4k+1 only."
         : $self->{'prime_type'} eq '4k+3' ? "\nPrimes of form 4k+3 only."
         : $self->{'prime_type'} eq 'twin' ? "\nTwin primes only."
         : $self->{'prime_type'} eq 'SG' ? "\nSophie Germain primes only (2P+1 also prime)."
         : $self->{'prime_type'} eq 'SG' ? "\nSafe primes only ((P-1)/2 also prime)."
         : "");
  } else {
    # class method
    return Math::NumSeq::__('Count of prime factors.');
  }
}

#------------------------------------------------------------------------------
#
# count 1-bits in exponents of primes
# A000028,A000379 seqs
#    A133008  characteristic
#    A131181,A026416  same, but 1 in "B" class
#    A064547  count 1 bits in prime exponents
#    A066724  so a(i)*a(j) not in seq
#    A026477  so a(i)*a(j)*a(k) not in seq
#    A050376  prime^(2^k)
#    A084400  smallest not dividing product a(1)..a(n-1), is prime^(2^k)

my %oeis_anum = (repeated => { all    => 'A001222',
                               odd    => 'A087436',
                               '4k+1' => 'A083025',
                               '4k+3' => 'A065339',
                             },
                 distinct => { all    => 'A001221',
                               odd    => 'A005087',
                               '4k+1' => 'A005089',
                               '4k+3' => 'A005091',
                               SG     => 'A156542',
                             },
                );
# OEIS-Catalogue: A001222
# OEIS-Catalogue: A087436 prime_type=odd
# OEIS-Catalogue: A083025 prime_type=4k+1
# OEIS-Catalogue: A065339 prime_type=4k+3

# OEIS-Catalogue: A001221 multiplicity=distinct
# OEIS-Catalogue: A005087 multiplicity=distinct prime_type=odd
# OEIS-Catalogue: A005089 multiplicity=distinct prime_type=4k+1
# OEIS-Catalogue: A005091 multiplicity=distinct prime_type=4k+3
# OEIS-Catalogue: A156542 multiplicity=distinct prime_type=SG

sub oeis_anum {
  my ($self) = @_;
  return $oeis_anum{$self->{'multiplicity'}}->{$self->{'prime_type'}};
}

#------------------------------------------------------------------------------

# prime_factors() is about 5x faster
#
sub ith {
  my ($self, $i) = @_;
  $i = abs($i);

  my ($good, @primes) = _prime_factors($i);
  return undef unless $good;

  my $multiplicity = ($self->{'multiplicity'} ne 'distinct');
  my $prime_type = $self->{'prime_type'};
  my $count = 0;

  while (@primes) {
    my $p = shift @primes;
    my $c = 1;
    while (@primes && $primes[0] == $p) {
      shift @primes;
      $c += $multiplicity;
    }

    if ($prime_type eq 'odd') {
      next unless $p & 1;
    } elsif ($prime_type eq '4k+1') {
      next unless ($p&3)==1;
    } elsif ($prime_type eq '4k+3') {
      next unless ($p&3)==3;
    } elsif ($prime_type eq 'twin') {
      next unless _is_twin_prime($p);
    } elsif ($prime_type eq 'SG') {
      next unless _is_SG_prime($p);
    } elsif ($prime_type eq 'safe') {
      next unless _is_safe_prime($p);

    # } elsif ($prime_type eq 'twin_first') {
    #   next unless is_prime($p+2);
    # } elsif ($prime_type eq 'twin_second') {
    #   next unless is_prime($p-2);
    }
    $count += $c;
  }

  if ($self->{'values_type'} eq 'mod2') {
    $count %= 2;
  }
  return $count;
}

# Return ($good, $prime,$prime,$prime,...).
# $good is true if a full factorization is found.
# $good is false if cannot factorize because $n is too big or infinite.
#
# If $n==0 or $n==1 then there are no prime factors and the return is
# $good=1 and an empty list of primes.
#
sub _prime_factors {
  my ($n) = @_;
  ### _prime_factors(): $n

  unless ($n >= 0) {
    return 0;
  }
  if (_is_infinite($n)) {
    return 0;
  }

  if ($n <= 0xFFFF_FFFF) {
    return (1, prime_factors($n));
  }

  my @ret;
  until ($n % 2) {
    ### div2: $n
    $n /= 2;
    push @ret, 2;
  }

  # Stop at when prime $p reaches $limit and when no prime factor has been
  # found for the last 20 attempted $p.  Stopping only after a run of no
  # factors found allows big primorials 2*3*5*7*13*... to be divided out.
  # If the divisions are making progress reducing $i then continue.
  #
  # Would like $p and $gap to count primes, not just odd numbers.  Perhaps
  # a table of small primes.  The first gap of 36 odds between primes
  # occurs at prime=31469.  cf A000230 smallest prime p for gap 2n.

  my $limit = 10_000 / (_blog2_estimate($n) || 1);
  my $gap = 0;
  for (my $p = 3; $gap < 36 || $p <= $limit ; $p += 2) {
    if ($n % $p) {
      $gap++;
    } else {
      do  {
        ### prime: $p
        $n /= $p;
        push @ret, $p;
      } until ($n % $p);

      if ($n <= 1) {
        ### all factors found ...
        return (1, @ret);
      }
      if ($n < 0xFFFF_FFFF) {
        ### remaining factors by XS ...
        return (1, @ret, prime_factors($n));
      }
      $gap = 0;
    }
  }
  return 0;  # factors too big
}

sub _is_twin_prime {
  my ($n) = @_;
  ### assert: $n >= 2
  ### assert: is_prime($n)
  return (is_prime($n+2) || is_prime($n-2));
}
sub _is_SG_prime {
  my ($n) = @_;
  ### assert: is_prime($n)
  return is_prime(2*$n+1);
}
sub _is_safe_prime {
  my ($n) = @_;
  ### assert: is_prime($n)
  return (($n&1) && is_prime(($n-1)/2));
}

sub pred {
  my ($self, $value) = @_;
  return ($value >= 0 && $value == int($value));
}

1;
__END__