File Coverage

blib/lib/Math/NumSeq/Totient.pm
Criterion Covered Total %
statement 84 96 87.5
branch 22 32 68.7
condition n/a
subroutine 18 19 94.7
pod 2 2 100.0
total 126 149 84.5


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014 Kevin Ryde
2              
3             # This file is part of Math-NumSeq.
4             #
5             # Math-NumSeq is free software; you can redistribute it and/or modify
6             # it under the terms of the GNU General Public License as published by the
7             # Free Software Foundation; either version 3, or (at your option) any later
8             # version.
9             #
10             # Math-NumSeq is distributed in the hope that it will be useful, but
11             # WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12             # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13             # for more details.
14             #
15             # You should have received a copy of the GNU General Public License along
16             # with Math-NumSeq. If not, see .
17              
18              
19             package Math::NumSeq::Totient;
20 3     3   6313 use 5.004;
  3         7  
21 3     3   9 use strict;
  3         4  
  3         60  
22 3     3   785 use Math::Prime::XS 0.23 'is_prime'; # version 0.23 fix for 1928099
  3         22504  
  3         143  
23 3     3   754 use Math::Factor::XS 'factors';
  3         17138  
  3         136  
24              
25 3     3   16 use vars '$VERSION', '@ISA';
  3         5  
  3         127  
26             $VERSION = 72;
27              
28 3     3   370 use Math::NumSeq;
  3         4  
  3         59  
29 3     3   722 use Math::NumSeq::Base::IterateIth;
  3         5  
  3         109  
30             @ISA = ('Math::NumSeq::Base::IterateIth',
31             'Math::NumSeq');
32             *_is_infinite = \&Math::NumSeq::_is_infinite;
33              
34 3     3   696 use Math::NumSeq::PrimeFactorCount;;
  3         5  
  3         104  
35             *_prime_factors = \&Math::NumSeq::PrimeFactorCount::_prime_factors;
36              
37             # uncomment this to run the ### lines
38             #use Smart::Comments;
39              
40              
41             # use constant name => Math::NumSeq::__('Totient');
42 3     3   12 use constant description => Math::NumSeq::__('Euler totient function, the count of how many numbers coprime to N.');
  3         4  
  3         8  
43 3     3   11 use constant characteristic_count => 1;
  3         3  
  3         116  
44 3     3   10 use constant characteristic_smaller => 1;
  3         3  
  3         104  
45 3     3   9 use constant characteristic_increasing => 0;
  3         3  
  3         111  
46 3     3   15 use constant values_min => 1;
  3         3  
  3         104  
47 3     3   9 use constant default_i_start => 1;
  3         3  
  3         103  
48              
49             #------------------------------------------------------------------------------
50             # cf A007617 non-totients, all odds, plus evens per A005277
51             # A005277 even non-totients, n s.t. n==phi(something) no solution
52             # A058980 non-totients 0mod4
53             # A056595 - sum non-square divisors
54             # A007614 totients ascending, with multiplicity
55              
56             # Dressler (1970) N(x) = num phi(n)<=x, then N(x)/x -> A
57             # A = zeta(2)*zeta(3)/zeta(6) = product primes 1+1/(p*(p-1))
58             #
59             # 2p is a non-totient if 2p+1 composite (p not an S-G prime)
60             # 4p is a non-totient iff 2p+1 and 4p+1 both composite
61             # if n non-totient and 2n+1 composite then 2n also non-totient
62             #
63 3     3   10 use constant oeis_anum => 'A000010';
  3         3  
  3         949  
64              
65             sub ith {
66 0     0 1 0 my ($self, $i) = @_;
67 0         0 return _totient($i);
68             }
69             sub _totient {
70 1135     1135   948 my ($n) = @_;
71             ### _totient(): $n
72              
73 1135 50       1275 if (_is_infinite($n)) {
74 0         0 return $n;
75             }
76 1135 100       1350 if ($n == 0) {
77 98         127 return 0;
78             }
79 1037         1216 my ($good, @primes) = _prime_factors($n);
80 1037 50       1247 return undef unless $good;
81              
82 1037         602 my $prev = 0;
83 1037         589 my $ret = 1;
84 1037         806 foreach my $p (@primes) {
85 1467 100       1340 if ($p == $prev) {
86 362         269 $ret *= $p;
87             } else {
88 1105         707 $ret *= $p - 1;
89 1105         900 $prev = $p;
90             }
91             }
92 1037         1415 return $ret;
93             }
94              
95             # totient(x)=p^a.q^b.r^c=n
96             # seek a prime w for x with w-1 dividing in n
97             # combinations of the primes of n to make w-1
98             #
99             # factor 2*f of n, arising from prime 2*f+1
100             #
101             # 8 arises from totient(15=3*5) = (3-1)*(5-1)=2*4
102             # 484=2*2*11*11 2*11=23 prime
103             #
104             sub pred {
105 2     2 1 25237 my ($self, $value) = @_;
106             ### Totient pred(): $value
107              
108 2 50       5 if ($value <= 1) {
109 0         0 return ($value == 1); # $value==0 no, $value==1 yes
110             }
111 2 50       92 if ($value % 2) {
112             ### no because odd ...
113 0         0 return 0;
114             }
115 2 50       127 unless ($value <= 0xFFFF_FFFF) {
116 0         0 return undef;
117             }
118 2         68 $value = "$value"; # numize any Math::BigInt for factors()
119 2 50       23 if (_pred_f($value,$value)) {
120 2         8 return 1;
121             }
122 0         0 return 0;
123             }
124             sub _pred_f {
125 16     16   13 my ($n, $prev_factor) = @_;
126             ### _pred_f(): "n=$n prev=$prev_factor"
127              
128 16 100       22 if ($n & 1) {
129             ### no odd ...
130 6         10 return 0;
131             }
132              
133 10         5 $n >>= 1;
134             ### halved: $n
135 10 50       13 if ($n == 1) {
136 0         0 return 1; # totient(3)=2 occurs
137             }
138              
139 10         23 foreach my $f (1, factors($n)) {
140             ### at: "n=$n f=$f"
141 14 100       18 if ($f >= $prev_factor) {
142             ### f too big, chop search ...
143 6         13 return 0;
144             }
145 8         8 my $p = 2*$f+1;
146             ### $p
147              
148 8         6 my $r = $n / $f;
149             ### divide out: "f=$f to r=$r"
150              
151 8 50       15 unless (is_prime($p)) {
152             ### no, not prime ...
153 0         0 next;
154             }
155              
156 8         6 for (;;) {
157 14 100       17 if (_pred_f ($r, $f)) { # recurse
158 2         4 return 1;
159             }
160 12 100       17 if ($r % $p) {
161 4         5 last;
162             }
163 8 100       10 if ($r == $p) {
164 2         4 return 1;
165             }
166 6         4 $r /= $p;
167             ### divide out prime: "p=$p to r=$r"
168             }
169             }
170              
171             ### whole: "n=$n p=".(2*$n+1)
172 0 0         if ($n >= $prev_factor) {
173             ### f too big, chop search ...
174 0           return 0;
175             }
176 0           return is_prime(2*$n+1);
177             }
178              
179              
180             # sub _totient {
181             # my ($x) = @_;
182             # my $count = (($x >= 1) # y=1 always
183             # + ($x > 2 && ($x&1)) # y=2 if $x odd
184             # + ($x > 3 && ($x % 3) != 0) # y=3
185             # + ($x > 4 && ($x&1)) # y=4 if $x odd
186             # );
187             # for (my $y = 5; $y < $x; $y++) {
188             # $count += _coprime($x,$y);
189             # }
190             # return $count;
191             # }
192             # sub _coprime { # for x
193             # my ($x, $y) = @_;
194             # #### _coprime(): "$x,$y"
195             # if ($y > $x) {
196             # return 0;
197             # }
198             # for (;;) {
199             # if ($y <= 1) {
200             # return ($y == 1);
201             # }
202             # ($x,$y) = ($y, $x % $y);
203             # }
204             # }
205              
206              
207             1;
208             __END__