comb.pl 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152
  1. #!/usr/bin/perl
  2. use 5.014;
  3. use ntheory qw(forcomb carmichael_lambda divisors is_prime binomial);
  4. use Math::Prime::Util::GMP qw(is_carmichael vecprod);
  5. #~ 1.8173 record: 64075459460541239985
  6. #~ 1.8274 record: 14370921672011352376545
  7. #~ 1.8670 record: 100724996094964745183793345
  8. #~ 1.8937 record: 5981802520294676451270038145
  9. #~ 1.9103 record: 141934960805533535384100259905
  10. #~ 1.9103 record: 2609668563371076823446624111609234945
  11. #~ 1.9176 record: 871005264581548962932418919531990803260747265
  12. #~ 1.9296 record: 74875990602425226938138664024259158432269224416705
  13. #~ 1.9332 record: 1018329989576989704159754753835889677652405111555820545
  14. #~ 1.9483 record: 165502573617193579886078524884554371013787876466850820614145
  15. #~ 1.9498 record: 14701083488299057530174696885922722686956286485260401577115414785
  16. #~ 1.9540 record: 321030150905393790929751720043602006651739765349595158023943724894346115208705
  17. #my @p = (3, 5, 17, 23, 29, 83, 89, 353, 449);
  18. #my @p = (3, 5, 17, 23, 29);
  19. #my @p = (3, 5, 17, 23, 29);
  20. #my @p = (3, 5, 17, 23, 29, 83, 89);
  21. #my @p = (3, 5, 17, 23, 29, 83, 89, 113, 197, 353, 449, 617);
  22. #my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 353, 449);
  23. #my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617);
  24. #my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113);
  25. #my @p = (3, 5, 17, 23, 29, 83, 89, 113, 197, 353, 449, 617, 3137);
  26. #my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 4019, 656657);
  27. #my @p = (3, 5, 17, 23, 29, 53, 83, 89);
  28. #my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113);
  29. #my @p = (3, 5, 17, 23, 29, 43, 53, 89, 113, 127);
  30. #my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 617);
  31. #my @p = (3, 5, 17, 23, 29);
  32. #my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617);
  33. #my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617);
  34. #my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 1409, 3137, 10193, 16073, 23297, 88397, 896897, 1500929, 18386369);
  35. my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 617); # seems promising
  36. #my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 617);
  37. my $P = vecprod(@p);
  38. #my $C = "37839385943068863406967633413004957540054532539686888463944906014566240419460804270776358938980032660929917901837033235462145";
  39. #my $C = "772459017179480479061611372132330246001039753130436193419524315193543873326133868681083905";
  40. my $C = "321030150905393790929751720043602006651739765349595158023943724894346115208705";
  41. #my $C = "463149379463251167706230703520995680069543921166639491398457345";
  42. my $L = carmichael_lambda($C);
  43. #$L = 4352299952*128;
  44. #$L = vecprod(2, 2, 2, 2, 2, 2, 2, 2, 7, 7, 7, 11, 13, 41);
  45. #$L *= 2;
  46. #$L *= 7;
  47. #$L = 294181888;
  48. say "lambda = $L";
  49. my @divisors = divisors($L);
  50. @divisors = grep { $_ > 1 and is_prime($_+1) } @divisors;
  51. @divisors = map{ $_ + 1 } @divisors;
  52. @divisors = grep { "@p" !~ /\b$_\b/ } @divisors;
  53. #@divisors = (53, 257, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 8009, 9857, 10193, 13313, 16073, 23297, 50177, 50513, 68993, 88397, 93809, 202049, 275969, 375233, 656657, 896897, 1500929, 3232769, 18386369, 22629377);
  54. say "Primes = {", join(', ', @divisors), '}';
  55. my $len = scalar(@divisors);
  56. printf("binomial(%s, %s) = %s\n", $len, $len>>1, binomial($len, $len>>1));
  57. foreach my $k(1..scalar(@divisors)) {
  58. say "Testing: $k";
  59. forcomb {
  60. if (is_carmichael(vecprod(@divisors[@_],$P))) {
  61. say vecprod(@divisors[@_], $P);
  62. }
  63. } $len, $k;
  64. }
  65. __END__
  66. 166320 -> 3649
  67. 15120 -> 2972
  68. 110880 -> 2925
  69. 55440 -> 2828
  70. 277200 -> 2775
  71. 25200 -> 2685
  72. 196560 -> 2586
  73. 332640 -> 2572
  74. 75600 -> 2483
  75. 65520 -> 2432
  76. 131040 -> 2274
  77. 720720 -> 2217
  78. 327600 -> 2185
  79. 100800 -> 2114
  80. 45360 -> 1964
  81. 221760 -> 1957
  82. 30240 -> 1940
  83. 831600 -> 1937
  84. 257040 -> 1917
  85. 60480 -> 1899
  86. 7201568 -> 1
  87. 4352299952 -> 1
  88. 5789168 -> 1
  89. 2618528 -> 1
  90. 2700544 -> 1
  91. 285824 -> 1
  92. 1350272 -> 1
  93. 2520 -> 226
  94. 3960 -> 141
  95. 3360 -> 128
  96. 3600 -> 127
  97. 2160 -> 119
  98. 4680 -> 114
  99. 6120 -> 106
  100. 6336 -> 95
  101. 5940 -> 86
  102. 4032 -> 82
  103. 6048 -> 67
  104. 1680 -> 65
  105. 9000 -> 63
  106. 9072 -> 61
  107. 4200 -> 60
  108. 3024 -> 53
  109. 4320 -> 52
  110. 1260 -> 52
  111. 4800 -> 52
  112. 9720 -> 48
  113. 1093200 -> 4
  114. 1012480 -> 3
  115. 1568400 -> 3
  116. 19153848 -> 2
  117. 1086840 -> 2
  118. 6984864 -> 2
  119. 2495532 -> 2
  120. 1320840 -> 2
  121. 1649640 -> 2
  122. 1461564 -> 2
  123. 1754970 -> 2
  124. 1840230 -> 2
  125. 5145300 -> 2
  126. 4156860 -> 2
  127. 2716452 -> 2
  128. 3715860 -> 2
  129. 1402092 -> 2
  130. 6531960 -> 2
  131. 1236300 -> 2
  132. 1626576 -> 2