File Coverage

blib/lib/Math/NumSeq/Totient.pm

Criterion	Covered	Total	%
statement	85	97	87.6
branch	22	32	68.7
condition			n/a
subroutine	18	19	94.7
pod	2	2	100.0
total	127	150	84.6

# Copyright 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::Totient;
use 5.004;
use strict;
use Math::Prime::XS 0.23 'is_prime'; # version 0.23 fix for 1928099
use Math::Factor::XS 'factors';

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

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


# use constant name => Math::NumSeq::__('Totient');
use constant description => Math::NumSeq::__('Euler totient function, the count of how many numbers coprime to N.');
use constant characteristic_count => 1;
use constant characteristic_smaller => 1;
use constant characteristic_increasing => 0;
use constant values_min => 1;
use constant default_i_start => 1;

#------------------------------------------------------------------------------
# cf A007617 non-totients, all odds, plus evens per A005277
#    A005277 even non-totients, n s.t. n==phi(something) no solution
#    A058980 non-totients 0mod4
#    A056595 - sum non-square divisors
#    A007614 totients ascending, with multiplicity

# Dressler (1970) N(x) = num phi(n)<=x, then N(x)/x -> A
# A = zeta(2)*zeta(3)/zeta(6) = product primes 1+1/(p*(p-1))
#
# 2p is a non-totient if 2p+1 composite (p not an S-G prime)
# 4p is a non-totient iff 2p+1 and 4p+1 both composite
# if n non-totient and 2n+1 composite then 2n also non-totient
#
use constant oeis_anum => 'A000010';

sub ith {
  my ($self, $i) = @_;
  return _totient($i);
}
sub _totient {
  my ($n) = @_;
  ### _totient(): $n

  if (_is_infinite($n)) {
    return $n;
  }
  if ($n == 0) {
    return 0;
  }
  my ($good, @primes) = _prime_factors($n);
  return undef unless $good;

  my $prev = 0;
  my $ret = 1;
  foreach my $p (@primes) {
    if ($p == $prev) {
      $ret *= $p;
    } else {
      $ret *= $p - 1;
      $prev = $p;
    }
  }
  return $ret;
}

# totient(x)=p^a.q^b.r^c=n
# seek a prime w for x with w-1 dividing in n
# combinations of the primes of n to make w-1
#
# factor 2*f of n, arising from prime 2*f+1
#
# 8 arises from totient(15=3*5) = (3-1)*(5-1)=2*4
# 484=2*2*11*11   2*11=23 prime
#
sub pred {
  my ($self, $value) = @_;
  ### Totient pred(): $value

  if ($value <= 1) {
    return ($value == 1); # $value==0 no, $value==1 yes
  }
  if ($value % 2) {
    ### no because odd ...
    return 0;
  }
  unless ($value <= 0xFFFF_FFFF) {
    return undef;
  }
  $value = "$value"; # numize any Math::BigInt for factors()
  if (_pred_f($value,$value)) {
    return 1;
  }
  return 0;
}
sub _pred_f {
  my ($n, $prev_factor) = @_;
  ### _pred_f(): "n=$n  prev=$prev_factor"

  if ($n & 1) {
    ### no odd ...
    return 0;
  }

  $n >>= 1;
  ### halved: $n
  if ($n == 1) {
    return 1;  # totient(3)=2 occurs
  }

  foreach my $f (1, factors($n)) {
    ### at: "n=$n f=$f"
    if ($f >= $prev_factor) {
      ### f too big, chop search ...
      return 0;
    }
    my $p = 2*$f+1;
    ### $p

    my $r = $n / $f;
    ### divide out: "f=$f to r=$r"

    unless (is_prime($p)) {
      ### no, not prime ...
      next;
    }

    for (;;) {
      if (_pred_f ($r, $f)) {  # recurse
        return 1;
      }
      if ($r % $p) {
        last;
      }
      if ($r == $p) {
        return 1;
      }
      $r /= $p;
      ### divide out prime: "p=$p to r=$r"
    }
  }

  ### whole: "n=$n  p=".(2*$n+1)
  if ($n >= $prev_factor) {
    ### f too big, chop search ...
    return 0;
  }
  return is_prime(2*$n+1);
}


# sub _totient {
#   my ($x) = @_;
#   my $count = (($x >= 1)                    # y=1 always
#                + ($x > 2 && ($x&1))         # y=2 if $x odd
#                + ($x > 3 && ($x % 3) != 0)  # y=3
#                + ($x > 4 && ($x&1))         # y=4 if $x odd
#               );
#   for (my $y = 5; $y < $x; $y++) {
#     $count += _coprime($x,$y);
#   }
#   return $count;
# }
# sub _coprime {        # for x
#   my ($x, $y) = @_;
#   #### _coprime(): "$x,$y"
#   if ($y > $x) {
#     return 0;
#   }
#   for (;;) {
#     if ($y <= 1) {
#       return ($y == 1);
#     }
#     ($x,$y) = ($y, $x % $y);
#   }
# }


1;
__END__