File Coverage

blib/lib/Math/NumSeq/PowerFlip.pm

Criterion	Covered	Total	%
statement	79	85	92.9
branch	12	18	66.6
condition	6	8	75.0
subroutine	18	18	100.0
pod	2	2	100.0
total	117	131	89.3

# Copyright 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::PowerFlip;
use 5.004;
use strict;

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;
*_to_bigint = \&Math::NumSeq::_to_bigint;

use Math::NumSeq::PrimeFactorCount;;
*_prime_factors = \&Math::NumSeq::PrimeFactorCount::_prime_factors;

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


use constant name => Math::NumSeq::__('Prime Exponent Flip');
use constant description => Math::NumSeq::__('Flip each prime and its exponent, so for example 3^8 -> 8^3');
use constant default_i_start => 1;
use constant characteristic_increasing => 0;
use constant characteristic_non_decreasing => 0;
use constant characteristic_integer => 1;
use constant characteristic_smaller => 1;
use constant values_min => 1; # at i=1


#------------------------------------------------------------------------------
# cf A005117 squarefrees have all exponents 1 so value=1
#    A013929 non-squarefrees have some exponent>1 so value>1
#    A005361 product of exponents of prime factorization

use constant oeis_anum => 'A008477';

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

use constant 1.02 _UV_LIMIT => 31**2; # is value=2**31

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

  if (_is_infinite($i)) {
    return $i;
  }
  $i = abs($i);

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

  if (! @primes) {
    return $i;  # 0,1 unchanged
  }

  my $value = 1;
  my $log = 0;

  for (;;) {
    my $p = shift @primes || last;
    my $count = 1;
    while (@primes && $primes[0] == $p) {
      shift @primes;
      $count++;
    }
    $log += $p*log($count);
    if ($log > 31) {
      $count = _to_bigint($count);
    }
    $value *= $count ** $p;
  }
  return $value;
}

sub pred {
  my ($self, $value) = @_;
  ### PowerFlip pred(): $value

  # ! is_square_free()

  unless ($value >= 0 && $value <= 0xFFFF_FFFF) {
    return undef;
  }
  if ($value != int($value)) {
    return 0;
  }
  $value = "$value"; # numize Math::BigInt for speed

  if ($value < 2) {
    return 1;
  }

  my $limit = sqrt($value) + 1;

  for (my $p = 2; $p <= $limit; $p += 2-($p==2)) {
    next if ($value % $p);
    ### prime factor: $p
    $value /= $p;

    if (($value % $p) == 0) {
      # found a square factor
      return 1;
    }

    $limit = sqrt($value) + 1;
    ### divided out: "$p, new limit $limit"
  }

  ### final: $value
  # $value now either 1 or a prime, no square factor found

  return 0;
}

1;
__END__