File Coverage

blib/lib/Math/NumSeq/MobiusFunction.pm

Criterion	Covered	Total	%
statement	68	69	98.5
branch	16	18	88.8
condition	5	8	62.5
subroutine	16	16	100.0
pod	2	2	100.0
total	107	113	94.6

# 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 .

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

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

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

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

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


# use constant name => Math::NumSeq::__('Mobius Function');
use constant description => Math::NumSeq::__('The Mobius function, being 1 for an even number of prime factors, -1 for an odd number, or 0 if any repeated factors (ie. not square-free).');
use constant characteristic_increasing => 0;
use constant characteristic_integer => 1;
use constant values_min => -1;
use constant values_max => 1;
use constant default_i_start => 1;

#------------------------------------------------------------------------------
# cf A030059 the -1 positions, being odd number of distinct primes
#    A030229 the 1 positions, being even number of distinct primes
#    A013929 the 0 positions, being square factor, ie. the non-square-frees
#    A005117 square free numbers, mobius -1 or +1
#
use constant oeis_anum => 'A008683'; # mobius -1,0,1 starting i=1


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

sub ith {
  my ($self, $i) = @_;
  ### MobiusFunction ith(): $i

  my $ret = 0;

  if (_is_infinite($i) || $i < 0) {
    return undef;
  }

  if (($i % 2) == 0) {
    $i /= 2;
    if (($i % 2) == 0) {
      return 0;  # square factor
    }
    $ret = 1;
  }

  if ($i <= 0xFFFF_FFFF) {
    $i = "$i"; # numize Math::BigInt for speed
  }

  my $sqrt = int(sqrt($i));
  my $limit = min ($sqrt,
                   10_000 / (_blog2_estimate($i) || 1));
  ### $sqrt
  ### $limit

  for (my $p = 3; $p <= $limit; $p += 2) {
    if (($i % $p) == 0) {
      $i /= $p;
      if (($i % $p) == 0) {
        ### square factor, zero ...
        return 0;
      }
      $ret ^= 1;
      $sqrt = int(sqrt($i));  # new smaller limit
      $limit = min ($sqrt, $limit);
      ### factor: "$p new ret $ret new limit $limit"
    }
  }
  if ($sqrt > $limit) {
    ### not all factors found up to limit ...
    # ENHANCE-ME: prime_factors() here if <2^32
    return undef;
  }
  if ($i != 1) {
    $ret ^= 1;
  }
  return ($ret ? -1 : 1);
}

sub pred {
  my ($self, $value) = @_;
  return ($value == 0 || $value == 1 || $value == -1);
}

1;
__END__