123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566 |
- #!/usr/bin/perl
- # Carmichael numbers divisible by the sum of their prime factors, sopfr (A001414).
- # https://oeis.org/A309003
- use 5.020;
- use strict;
- use warnings;
- use experimental qw(signatures);
- use Math::GMPz;
- use ntheory qw(:all);
- use Math::AnyNum qw(is_smooth);
- sub odd_part ($n) {
- $n >> valuation($n, 2);
- }
- my %seen;
- my @nums;
- while (<>) {
- next if /^\h*#/;
- /\S/ or next;
- my $n = (split(' ', $_))[-1];
- $n || next;
- length($n) > 15 and next;
- if ($n < 3240392401) {
- next;
- }
- next if $seen{$n}++;
- # is_smooth($n, 1e7) || next;
- is_carmichael($n) || next;
- say "Candidate: $n";
- if ($n > ~0) {
- $n = Math::GMPz->new("$n");
- }
- push @nums, $n;
- }
- @nums = sort { $a <=> $b } @nums;
- say "There are ", scalar(@nums), " total numbers";
- foreach my $n (@nums) {
- #say "Testing: $n";
- #if (odd_part($n - 1) == odd_part(euler_phi($n))) {
- # die "Found counter-example: $n";
- #}
- if ($n % vecsum(factor($n)) == 0) {
- print $n, ", ";
- }
- }
- __END__
- 3240392401, 13577445505, 14446721521, 84127131361, 203340265921, 241420757761, 334797586201, 381334973041, 461912170321, 1838314142785, 3636869821201, 10285271821441, 17624045440981, 18773053896961, 20137015596061, 24811804945201, 26863480687681, 35598629998801
|