prog.pl 1.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566
  1. #!/usr/bin/perl
  2. # Carmichael numbers divisible by the sum of their prime factors, sopfr (A001414).
  3. # https://oeis.org/A309003
  4. use 5.020;
  5. use strict;
  6. use warnings;
  7. use experimental qw(signatures);
  8. use Math::GMPz;
  9. use ntheory qw(:all);
  10. use Math::AnyNum qw(is_smooth);
  11. sub odd_part ($n) {
  12. $n >> valuation($n, 2);
  13. }
  14. my %seen;
  15. my @nums;
  16. while (<>) {
  17. next if /^\h*#/;
  18. /\S/ or next;
  19. my $n = (split(' ', $_))[-1];
  20. $n || next;
  21. length($n) > 15 and next;
  22. if ($n < 3240392401) {
  23. next;
  24. }
  25. next if $seen{$n}++;
  26. # is_smooth($n, 1e7) || next;
  27. is_carmichael($n) || next;
  28. say "Candidate: $n";
  29. if ($n > ~0) {
  30. $n = Math::GMPz->new("$n");
  31. }
  32. push @nums, $n;
  33. }
  34. @nums = sort { $a <=> $b } @nums;
  35. say "There are ", scalar(@nums), " total numbers";
  36. foreach my $n (@nums) {
  37. #say "Testing: $n";
  38. #if (odd_part($n - 1) == odd_part(euler_phi($n))) {
  39. # die "Found counter-example: $n";
  40. #}
  41. if ($n % vecsum(factor($n)) == 0) {
  42. print $n, ", ";
  43. }
  44. }
  45. __END__
  46. 3240392401, 13577445505, 14446721521, 84127131361, 203340265921, 241420757761, 334797586201, 381334973041, 461912170321, 1838314142785, 3636869821201, 10285271821441, 17624045440981, 18773053896961, 20137015596061, 24811804945201, 26863480687681, 35598629998801