Circular primes: Difference between revisions
m (added an Oxford comma.) |
|||
Line 7: | Line 7: | ||
'''1193''' is a circular prime, since '''1931''', '''9311''' and '''3119''' are all also prime. |
'''1193''' is a circular prime, since '''1931''', '''9311''' and '''3119''' are all also prime. |
||
Note that a number which is a cyclic permutation of a smaller circular prime is not considered to be itself a circular prime. So '''13''' is a circular prime but '''31''' is not. |
Note that a number which is a cyclic permutation of a smaller circular prime is not considered to be itself a circular prime. So '''13''' is a circular prime, but '''31''' is not. |
||
Revision as of 19:51, 24 May 2020
You are encouraged to solve this task according to the task description, using any language you may know.
- Definitions
A circular prime is a prime number with the property that the number generated at each intermediate step when cyclically permuting its (base 10) digits will also be prime.
For example: 1193 is a circular prime, since 1931, 9311 and 3119 are all also prime.
Note that a number which is a cyclic permutation of a smaller circular prime is not considered to be itself a circular prime. So 13 is a circular prime, but 31 is not.
A repunit (denoted by R) is a number whose base 10 representation contains only the digit 1.
For example: R(2) = 11 and R(5) = 11111 are repunits.
- Task
- Find the first 19 circular primes.
- If your language has access to arbitrary precision integer arithmetic, given that they are all repunits, find the next 4 circular primes.
- (Stretch) Determine which of the following repunits are probably circular primes: R(5003), R(9887), R(15073), R(25031), R(35317) and R(49081). The larger ones may take a long time to process so just do as many as you reasonably can.
- See also
- Wikipedia article - Circular primes.
- Wikipedia article - Repunit.
- OEIS sequence A016114 - Circular primes.
C
<lang c>#include <stdbool.h>
- include <stdint.h>
- include <stdio.h>
- include <stdlib.h>
- include <string.h>
- include <gmp.h>
bool is_prime(uint32_t n) {
if (n == 2) return true; if (n < 2 || n % 2 == 0) return false; for (uint32_t p = 3; p * p <= n; p += 2) { if (n % p == 0) return false; } return true;
}
// e.g. returns 2341 if n = 1234 uint32_t cycle(uint32_t n) {
uint32_t m = n, p = 1; while (m >= 10) { p *= 10; m /= 10; } return m + 10 * (n % p);
}
bool is_circular_prime(uint32_t p) {
if (!is_prime(p)) return false; uint32_t p2 = cycle(p); while (p2 != p) { if (p2 < p || !is_prime(p2)) return false; p2 = cycle(p2); } return true;
}
void test_repunit(uint32_t digits) {
char* str = malloc(digits + 1); if (str == 0) { fprintf(stderr, "Out of memory\n"); exit(1); } memset(str, '1', digits); str[digits] = 0; mpz_t bignum; mpz_init_set_str(bignum, str, 10); free(str); if (mpz_probab_prime_p(bignum, 10)) printf("R(%u) is probably prime.\n", digits); else printf("R(%u) is not prime.\n", digits); mpz_clear(bignum);
}
int main() {
uint32_t p = 2; printf("First 19 circular primes:\n"); for (int count = 0; count < 19; ++p) { if (is_circular_prime(p)) { if (count > 0) printf(", "); printf("%u", p); ++count; } } printf("\n"); printf("Next 4 circular primes:\n"); uint32_t repunit = 1, digits = 1; for (; repunit < p; ++digits) repunit = 10 * repunit + 1; mpz_t bignum; mpz_init_set_ui(bignum, repunit); for (int count = 0; count < 4; ) { if (mpz_probab_prime_p(bignum, 15)) { if (count > 0) printf(", "); printf("R(%u)", digits); ++count; } ++digits; mpz_mul_ui(bignum, bignum, 10); mpz_add_ui(bignum, bignum, 1); } mpz_clear(bignum); printf("\n"); test_repunit(5003); test_repunit(9887); test_repunit(15073); test_repunit(25031); test_repunit(35317); test_repunit(49081); return 0;
}</lang>
- Output:
With GMP 6.2.0, execution time on my system is about 13 minutes (3.2 GHz Quad-Core Intel Core i5, macOS 10.15.4).
First 19 circular primes: 2, 3, 5, 7, 11, 13, 17, 37, 79, 113, 197, 199, 337, 1193, 3779, 11939, 19937, 193939, 199933 Next 4 circular primes: R(19), R(23), R(317), R(1031) R(5003) is not prime. R(9887) is not prime. R(15073) is not prime. R(25031) is not prime. R(35317) is not prime. R(49081) is probably prime.
C++
<lang cpp>#include <cstdint>
- include <algorithm>
- include <iostream>
- include <sstream>
- include <gmpxx.h>
typedef mpz_class integer;
bool is_prime(const integer& n, int reps = 50) {
return mpz_probab_prime_p(n.get_mpz_t(), reps);
}
std::string to_string(const integer& n) {
std::ostringstream out; out << n; return out.str();
}
bool is_circular_prime(const integer& p) {
if (!is_prime(p)) return false; std::string str(to_string(p)); for (size_t i = 0, n = str.size(); i + 1 < n; ++i) { std::rotate(str.begin(), str.begin() + 1, str.end()); integer p2(str, 10); if (p2 < p || !is_prime(p2)) return false; } return true;
}
integer next_repunit(const integer& n) {
integer p = 1; while (p < n) p = 10 * p + 1; return p;
}
integer repunit(int digits) {
std::string str(digits, '1'); integer p(str); return p;
}
void test_repunit(int digits) {
if (is_prime(repunit(digits), 10)) std::cout << "R(" << digits << ") is probably prime\n"; else std::cout << "R(" << digits << ") is not prime\n";
}
int main() {
integer p = 2; std::cout << "First 19 circular primes:\n"; for (int count = 0; count < 19; ++p) { if (is_circular_prime(p)) { if (count > 0) std::cout << ", "; std::cout << p; ++count; } } std::cout << '\n'; std::cout << "Next 4 circular primes:\n"; p = next_repunit(p); std::string str(to_string(p)); int digits = str.size(); for (int count = 0; count < 4; ) { if (is_prime(p, 15)) { if (count > 0) std::cout << ", "; std::cout << "R(" << digits << ")"; ++count; } p = repunit(++digits); } std::cout << '\n'; test_repunit(5003); test_repunit(9887); test_repunit(15073); test_repunit(25031); test_repunit(35317); test_repunit(49081); return 0;
}</lang>
- Output:
First 19 circular primes: 2, 3, 5, 7, 11, 13, 17, 37, 79, 113, 197, 199, 337, 1193, 3779, 11939, 19937, 193939, 199933 Next 4 circular primes: R(19), R(23), R(317), R(1031) R(5003) is not prime R(9887) is not prime R(15073) is not prime R(25031) is not prime R(35317) is not prime R(49081) is probably prime
D
<lang D>import std.bigint; import std.stdio;
immutable PRIMES = [
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997
];
bool isPrime(BigInt n) {
if (n < 2) { return false; }
foreach (p; PRIMES) { if (n == p) { return true; } if (n % p == 0) { return false; } if (p * p > n) { return true; } }
for (auto m = BigInt(PRIMES[$ - 1]); m * m <= n ; m += 2) { if (n % m == 0) { return false; } }
return true;
}
// e.g. returns 2341 if n = 1234 BigInt cycle(BigInt n) {
BigInt m = n; BigInt p = 1; while (m >= 10) { p *= 10; m /= 10; } return m + 10 * (n % p);
}
bool isCircularPrime(BigInt p) {
if (!isPrime(p)) { return false; } for (auto p2 = cycle(p); p2 != p; p2 = cycle(p2)) { if (p2 < p || !isPrime(p2)) { return false; } } return true;
}
BigInt repUnit(int len) {
BigInt n = 0; while (len > 0) { n = 10 * n + 1; len--; } return n;
}
void main() {
writeln("First 19 circular primes:"); int count = 0; foreach (p; PRIMES) { if (isCircularPrime(BigInt(p))) { if (count > 0) { write(", "); } write(p); count++; } } for (auto p = BigInt(PRIMES[$ - 1]) + 2; count < 19; p += 2) { if (isCircularPrime(BigInt(p))) { if (count > 0) { write(", "); } write(p); count++; } } writeln;
}</lang>
- Output:
First 19 circular primes: 2, 3, 5, 7, 11, 13, 17, 37, 79, 113, 197, 199, 337, 1193, 3779, 11939, 19937, 193939, 199933
Factor
Unfortunately Factor's miller-rabin test or bignums aren't quite up to the task of finding the next four circular prime repunits in a reasonable time. It takes ~90 seconds to check R(7)-R(1031).
<lang factor>USING: combinators.short-circuit formatting io kernel lists lists.lazy math math.combinatorics math.functions math.parser math.primes sequences sequences.extras ;
! Create an ordered infinite lazy list of circular prime ! "candidates" -- the numbers 2, 3, 5 followed by numbers ! composed of only the digits 1, 3, 7, and 9.
- candidates ( -- list )
L{ "2" "3" "5" "7" } 2 lfrom [ "1379" swap selections >list ] lmap-lazy lconcat lappend ;
- circular-prime? ( str -- ? )
all-rotations { [ [ infimum ] [ first = ] bi ] [ [ string>number prime? ] all? ] } 1&& ;
- circular-primes ( -- list )
candidates [ circular-prime? ] lfilter ;
- prime-repunits ( -- list )
7 lfrom [ 10^ 1 - 9 / prime? ] lfilter ;
"The first 19 circular primes are:" print 19 circular-primes ltake [ write bl ] leach nl nl
"The next 4 circular primes, in repunit format, are:" print 4 prime-repunits ltake [ "R(%d) " printf ] leach nl</lang>
- Output:
The first 19 circular primes are: 2 3 5 7 11 13 17 37 79 113 197 199 337 1193 3779 11939 19937 193939 199933 The next 4 circular primes, in repunit format, are: R(19) R(23) R(317) R(1031)
Go
<lang go>package main
import (
"fmt" big "github.com/ncw/gmp" "strings"
)
// OK for 'small' numbers. func isPrime(n int) bool {
switch { case n < 2: return false case n%2 == 0: return n == 2 case n%3 == 0: return n == 3 default: d := 5 for d*d <= n { if n%d == 0 { return false } d += 2 if n%d == 0 { return false } d += 4 } return true }
}
func repunit(n int) *big.Int {
ones := strings.Repeat("1", n) b, _ := new(big.Int).SetString(ones, 10) return b
}
var circs = []int{}
// binary search is overkill for a small number of elements func alreadyFound(n int) bool {
for _, i := range circs { if i == n { return true } } return false
}
func isCircular(n int) bool {
nn := n pow := 1 // will eventually contain 10 ^ d where d is number of digits in n for nn > 0 { pow *= 10 nn /= 10 } nn = n for { nn *= 10 f := nn / pow // first digit nn += f * (1 - pow) if alreadyFound(nn) { return false } if nn == n { break } if !isPrime(nn) { return false } } return true
}
func main() {
fmt.Println("The first 19 circular primes are:") digits := [4]int{1, 3, 7, 9} q := []int{1, 2, 3, 5, 7, 9} // queue the numbers to be examined fq := []int{1, 2, 3, 5, 7, 9} // also queue the corresponding first digits count := 0 for { f := q[0] // peek first element fd := fq[0] // peek first digit if isPrime(f) && isCircular(f) { circs = append(circs, f) count++ if count == 19 { break } } copy(q, q[1:]) // pop first element q = q[:len(q)-1] // reduce length by 1 copy(fq, fq[1:]) // ditto for first digit queue fq = fq[:len(fq)-1] if f == 2 || f == 5 { // if digits > 1 can't contain a 2 or 5 continue } // add numbers with one more digit to queue // only numbers whose last digit >= first digit need be added for _, d := range digits { if d >= fd { q = append(q, f*10+d) fq = append(fq, fd) } } } fmt.Println(circs) fmt.Println("\nThe next 4 circular primes, in repunit format, are:") count = 0 var rus []string for i := 7; count < 4; i++ { if repunit(i).ProbablyPrime(10) { count++ rus = append(rus, fmt.Sprintf("R(%d)", i)) } } fmt.Println(rus) fmt.Println("\nThe following repunits are probably circular primes:") for _, i := range []int{5003, 9887, 15073, 25031, 35317, 49081} { fmt.Printf("R(%-5d) : %t\n", i, repunit(i).ProbablyPrime(10)) }
}</lang>
- Output:
The first 19 circular primes are: [2 3 5 7 11 13 17 37 79 113 197 199 337 1193 3779 11939 19937 193939 199933] The next 4 circular primes, in repunit format, are: [R(19) R(23) R(317) R(1031)] The following repunits are probably circular primes: R(5003 ) : false R(9887 ) : false R(15073) : false R(25031) : false R(35317) : false R(49081) : true
Haskell
Uses arithmoi Library: http://hackage.haskell.org/package/arithmoi-0.10.0.0 <lang haskell>import Math.NumberTheory.Primes (Prime, unPrime, nextPrime) import Math.NumberTheory.Primes.Testing (isPrime, millerRabinV) import Text.Printf (printf)
rotated :: [Integer] -> [Integer] rotated xs
| any (< head xs) xs = [] | otherwise = map asNum $ take (pred $ length xs) $ rotate xs where rotate [] = [] rotate (d:ds) = ds <> [d] : rotate (ds <> [d])
asNum :: [Integer] -> Integer asNum [] = 0 asNum n@(d:ds)
| all (==1) n = read $ concatMap show n | otherwise = (d * (10 ^ length ds)) + asNum ds
digits :: Integer -> [Integer] digits 0 = [] digits n = digits d <> [r]
where (d, r) = n `quotRem` 10
isCircular :: Bool -> Integer -> Bool isCircular repunit n
| repunit = millerRabinV 0 n | n < 10 = True | even n = False | null rotations = False | any (<n) rotations = False | otherwise = all isPrime rotations where rotations = rotated $ digits n
repunits :: [Integer] repunits = go 2
where go n = asNum (replicate n 1) : go (succ n)
asRepunit :: Int -> Integer asRepunit n = asNum $ replicate n 1
main :: IO () main = do
printf "The first 19 circular primes are:\n%s\n\n" $ circular primes printf "The next 4 circular primes, in repunit format are:\n" mapM_ (printf "R(%d) ") $ reps repunits printf "\n\nThe following repunits are probably circular primes:\n" mapM_ (uncurry (printf "R(%d) : %s\n") . checkReps) [5003, 9887, 15073, 25031, 35317, 49081] where primes = map unPrime [nextPrime 1..] circular = show . take 19 . filter (isCircular False) reps = map (sum . digits). tail . take 5 . filter (isCircular True) checkReps = (,) <$> id <*> show . isCircular True . asRepunit</lang>
- Output:
The first 19 circular primes are: [2,3,5,7,11,13,17,37,79,113,197,199,337,1193,3779,11939,19937,193939,199933] The next 4 circular primes, in repunit format are: R(19) R(23) R(317) R(1031) The following repunits are probably circular primes: R(5003) : False R(9887) : False R(15073) : False R(25031) : False R(35317) : False R(49081) : True ./circular_primes 277.56s user 1.82s system 97% cpu 4:47.91 total
J
<lang J>
R=: [: ". 'x' ,~ #&'1' assert 11111111111111111111111111111111x -: R 32
Filter=: (#~`)(`:6)
rotations=: (|."0 1~ i.@#)&.(10&#.inv) assert 123 231 312 -: rotations 123
primes_less_than=: i.&.:(p:inv) assert 2 3 5 7 11 -: primes_less_than 12
NB. circular y --> y is the order of magnitude.
circular=: monad define
P25=: ([: -. (0 e. 1 3 7 9 e.~ 10 #.inv ])&>)Filter primes_less_than 10^y NB. Q25 are primes with 1 3 7 9 digits P=: 2 5 , P25 en=: # P group=: en # 0 next=: 1 for_i. i. # group do. if. 0 = i { group do. NB. if untested j =: P i. rotations i { P NB. j are the indexes of the rotated numbers in the list of primes if. en e. j do. NB. if any are unfound j=: j -. en NB. prepare to mark them all as searched, and failed. g=: _1 else. g=: next NB. mark the set as found in a new group. Because we can. next=: >: next end. group=: g j} group NB. apply the tested mark end. end. group </. P
) </lang> J lends itself to investigation. Demonstrate and play with the definitions.
circular 3 NB. the values in the long list have a composite under rotation ┌─┬─┬─┬─┬──┬─────┬─────┬──────────────────────────────────────────────────────────────────────────┬─────┬─────┬───────────┬───────────┬───────────┬───────────┐ │2│5│3│7│11│13 31│17 71│19 137 139 173 179 191 193 313 317 331 379 397 739 773 797 911 937 977 997│37 73│79 97│113 131 311│197 719 971│199 919 991│337 373 733│ └─┴─┴─┴─┴──┴─────┴─────┴──────────────────────────────────────────────────────────────────────────┴─────┴─────┴───────────┴───────────┴───────────┴───────────┘ circular 2 NB. VV ┌─┬─┬─┬─┬──┬─────┬─────┬──┬─────┬─────┐ │2│5│3│7│11│13 31│17 71│19│37 73│79 97│ └─┴─┴─┴─┴──┴─────┴─────┴──┴─────┴─────┘ q: 91 NB. factor the lone entry 7 13 RC=: circular 8 {: RC NB. tail ┌─────────────────────────────────────────┐ │199933 319993 331999 933199 993319 999331│ └─────────────────────────────────────────┘ (< >./)@:(#&>) Filter circular 3 NB. remove the box containing most items ┌─┬─┬─┬─┬──┬─────┬─────┬─────┬─────┬───────────┬───────────┬───────────┬───────────┐ │2│5│3│7│11│13 31│17 71│37 73│79 97│113 131 311│197 719 971│199 919 991│337 373 733│ └─┴─┴─┴─┴──┴─────┴─────┴─────┴─────┴───────────┴───────────┴───────────┴───────────┘ ] CPS=: {.&> (< >./)@:(#&>) Filter RC NB. first 19 circular primes 2 5 3 7 11 13 17 37 79 113 197 199 337 1193 3779 11939 19937 193939 199933 # CPS NB. yes, 19 found. 19
Brief investigation into repunits.
Note 'The current Miller-Rabin test implemented in c is insufficient for this task' R=: ([: ". 'x' ,~ #&'1')&> (;q:@R)&> 500 |limit error | (;q:@R)&>500 ) boxdraw_j_ 0 Filter=: (#~`)(`:6) R=: ([: ". 'x' ,~ #&'1')&> (; q:@R)&> (0 ~: 3&|)Filter >: i. 26 NB. factor some repunits ┌──┬─────────────────────────────────┐ │1 │ │ ├──┼─────────────────────────────────┤ │2 │11 │ ├──┼─────────────────────────────────┤ │4 │11 101 │ ├──┼─────────────────────────────────┤ │5 │41 271 │ ├──┼─────────────────────────────────┤ │7 │239 4649 │ ├──┼─────────────────────────────────┤ │8 │11 73 101 137 │ ├──┼─────────────────────────────────┤ │10│11 41 271 9091 │ ├──┼─────────────────────────────────┤ │11│21649 513239 │ ├──┼─────────────────────────────────┤ │13│53 79 265371653 │ ├──┼─────────────────────────────────┤ │14│11 239 4649 909091 │ ├──┼─────────────────────────────────┤ │16│11 17 73 101 137 5882353 │ ├──┼─────────────────────────────────┤ │17│2071723 5363222357 │ ├──┼─────────────────────────────────┤ │19│1111111111111111111 │ ├──┼─────────────────────────────────┤ │20│11 41 101 271 3541 9091 27961 │ ├──┼─────────────────────────────────┤ │22│11 11 23 4093 8779 21649 513239 │ ├──┼─────────────────────────────────┤ │23│11111111111111111111111 │ ├──┼─────────────────────────────────┤ │25│41 271 21401 25601 182521213001 │ ├──┼─────────────────────────────────┤ │26│11 53 79 859 265371653 1058313049│ └──┴─────────────────────────────────┘ NB. R(2) R(19), R(23) are probably prime.
Julia
Note that the evalpoly function used in this program was added in Julia 1.4 <lang julia>using Lazy, Primes
function iscircularprime(n)
!isprime(n) && return false dig = digits(n) return all(i -> (m = evalpoly(10, circshift(dig, i))) >= n && isprime(m), 1:length(dig)-1)
end
filtcircular(n, rang) = Int.(collect(take(n, filter(iscircularprime, rang)))) isprimerepunit(n) = isprime(evalpoly(BigInt(10), ones(Int, n))) filtrep(n, rang) = collect(take(n, filter(isprimerepunit, rang)))
println("The first 19 circular primes are:\n", filtcircular(19, Lazy.range(2))) print("\nThe next 4 circular primes, in repunit format, are: ",
mapreduce(n -> "R($n) ", *, filtrep(4, Lazy.range(6))))
println("\n\nChecking larger repunits:") for i in [5003, 9887, 15073, 25031, 35317, 49081]
println("R($i) is ", isprimerepunit(i) ? "prime." : "not prime.")
end
</lang>
- Output:
The first 19 circular primes are: [2, 3, 5, 7, 11, 13, 17, 37, 79, 113, 197, 199, 337, 1193, 3779, 11939, 19937, 193939, 199933] The next 4 circular primes, in repunit format, are: R(19) R(23) R(317) R(1031) Checking larger repunits: R(5003) is not prime. R(9887) is not prime. R(15073) is not prime. R(25031) is not prime. R(35317) is not prime. R(49081) is prime.
Perl
<lang perl>use strict; use warnings; use feature 'say'; use List::Util 'min'; use ntheory 'is_prime';
sub rotate { my($i,@a) = @_; join , @a[$i .. @a-1, 0 .. $i-1] }
sub isCircular {
my ($n) = @_; return 0 unless is_prime($n); my @circular = split //, $n; return 0 if min(@circular) < $circular[0]; for (1 .. scalar @circular) { my $r = join , rotate($_,@circular); return 0 unless is_prime($r) and $r >= $n; } 1
}
say "The first 19 circular primes are:"; for ( my $i = 1, my $count = 0; $count < 19; $i++ ) {
++$count and print "$i " if isCircular($i);
}
say "\n\nThe next 4 circular primes, in repunit format, are:"; for ( my $i = 7, my $count = 0; $count < 4; $i++ ) {
++$count and say "R($i)" if is_prime 1 x $i
}
say "\nRepunit testing:";
for (5003, 9887, 15073, 25031, 35317, 49081) {
say "R($_): Prime? " . (is_prime 1 x $_ ? 'True' : 'False');
}</lang>
- Output:
The first 19 circular primes are: 2 3 5 7 11 13 17 37 79 113 197 199 337 1193 3779 11939 19937 193939 199933 The next 4 circular primes, in repunit format, are: R(19) R(23) R(317) R(1031) Repunit testing: R(5003): Prime? False R(9887): Prime? False R(15073): Prime? False R(25031): Prime? False R(35317): Prime? False R(49081): Prime? True
Phix
<lang Phix>function circular(integer p)
integer len = length(sprintf("%d",p)), pow = power(10,len-1), p0 = p for i=1 to len-1 do p = pow*remainder(p,10)+floor(p/10) if p<p0 or not is_prime(p) then return false end if end for return true
end function
sequence c = {} integer n = 1 while length(c)<19 do
integer p = get_prime(n) if circular(p) then c &= p end if n += 1
end while printf(1,"The first 19 circular primes are:\n%v\n\n",{c})
include mpfr.e procedure repunit(mpz z, integer n)
mpz_set_str(z,repeat('1',n))
end procedure
c = {} n = 7 mpz z = mpz_init() randstate state = gmp_randinit_mt() while length(c)<4 do
repunit(z,n) if mpz_probable_prime_p(z,state) then c = append(c,sprintf("R(%d)",n)) end if n += 1
end while printf(1,"The next 4 circular primes, in repunit format, are:\n%s\n\n",{join(c)})</lang>
- Output:
The first 19 circular primes are: {2,3,5,7,11,13,17,37,79,113,197,199,337,1193,3779,11939,19937,193939,199933} The next 4 circular primes, in repunit format, are: R(19) R(23) R(317) R(1031)
stretch
<lang Phix>constant t = {5003, 9887, 15073, 25031, 35317, 49081} printf(1,"The following repunits are probably circular primes:\n") for i=1 to length(t) do
integer ti = t[i] atom t0 = time() repunit(z,ti) bool bPrime = mpz_probable_prime_p(z,state,1) printf(1,"R(%d) : %t (%s)\n", {ti, bPrime, elapsed(time()-t0)})
end for</lang>
- Output:
64-bit can only cope with the first five (it terminates abruptly on the sixth)
For comparison, the above Julia code (8/4/20, 64 bit) manages 1s, 5.6s, 15s, 50s, 1 min 50s (and 1 hour 45 min 40s) on the same box.
And Perl (somehow) manages 0/0/0/55s/0/21 mins 35 secs...
The following repunits are probably circular primes: R(5003) : false (2.0s) R(9887) : false (13.5s) R(15073) : false (45.9s) R(25031) : false (1 minute and 19s) R(35317) : false (3 minutes and 04s)
32-bit is much slower and can only cope with the first four
The following repunits are probably circular primes: R(5003) : false (10.2s) R(9887) : false (54.9s) R(15073) : false (2 minutes and 22s) R(25031) : false (7 minutes and 45s) diag looping, error code is 1, era is #00644651
Raku
Most of the repunit testing is relatively speedy using the ntheory library. The really slow ones are R(25031), at ~42 seconds and R(49081) at 922(!!) seconds.
<lang perl6>#!/usr/bin/env raku
- 20200406 Raku programming solution
sub isCircular(\n) {
return False unless n.is-prime; my @circular = n.comb; return False if @circular.min < @circular[0]; for 1 ..^ @circular -> $i { return False unless .is-prime and $_ >= n given @circular.rotate($i).join; } True
}
say "The first 19 circular primes are:"; say ((2..*).hyper.grep: { isCircular $_ })[^19];
say "\nThe next 4 circular primes, in repunit format, are:"; loop ( my $i = 7, my $count = 0; $count < 4; $i++ ) {
++$count, say "R($i)" if (1 x $i).is-prime
}
use ntheory:from<Perl5> qw[is_prime];
say "\nRepunit testing:";
(5003, 9887, 15073, 25031, 35317, 49081).map: {
my $now = now; say "R($_): Prime? ", ?is_prime("{1 x $_}"), " {(now - $now).fmt: '%.2f'}"
}</lang>
- Output:
The first 19 circular primes are: (2 3 5 7 11 13 17 37 79 113 197 199 337 1193 3779 11939 19937 193939 199933) The next 4 circular primes, in repunit format, are: R(19) R(23) R(317) R(1031) Repunit testing: R(5003): Prime? False 0.00 R(9887): Prime? False 0.01 R(15073): Prime? False 0.02 R(25031): Prime? False 41.40 R(35317): Prime? False 0.32 R(49081): Prime? True 921.73
REXX
By far, the greatest cost of CPU time is in the generation of a suitable number of primes. <lang rexx>/*REXX program finds & displays circular primes (with a title & in a horizontal format).*/ parse arg N hp . /*obtain optional arguments from the CL*/ if N== | N=="," then N= 19 /* " " " " " " */ if hp== | hp=="," then hp= 1000000 /* " " " " " " */ call genP /*gen primes up to hp (200,000). */ q= 024568 /*digs that most circular P can't have.*/
- = 0; $= /*#: circular P count; $: is a list.*/
do j=2 until #==N /* [↓] traipse through the number(s).*/ if \!.j then iterate /*Is J not prime? Then skip number.*/ if j>9 & verify(j, q, 'M')>0 then iterate /*Does J contain forbidden digs? Skip.*/ if \circP(j) then iterate /*Not circular? Then skip this number.*/ #= # + 1 /*bump the count of circular primes. */ $= $ j /*add this prime number ──► $ list. */ end /*j*/ /*at this point, $ has a leading blank.*/
say center(' first' # "circular primes ", 79, '─') say strip($) exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ circP: procedure expose @. !.; parse arg x 1 ox /*obtain a prime number to be examined.*/
do length(x)-1 /*parse X number, rotating the digits*/ parse var x f 2 y /*get 1st digit & the rightmost digits.*/ x= y || f /*construct a new possible circular P. */ if x<ox then return 0 /*is number < the original? ¬ circular*/ if \!.x then return 0 /* " " not prime? ¬ circular*/ end /*len···*/ return 1 /*passed all tests, X is a circular P.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/ genP: @.1=2; @.2=3; @.3=5; @.4=7; @.5=11; @.6= 13; nP=6 /*assign low primes; # primes. */
!.= 0; !.2=1; !.3=1; !.5=1; !.7=1; !.11=1; !.13=1 /*assign primality to numbers. */ do j=@.nP+4 by 2 to hp /*only find odd primes from here on. */ if j// 3==0 then iterate /*is J divisible by #3 Then not prime*/ parse var j -1 _;if _==5 then iterate /*Is last digit a "5"? " " " */ if j// 7==0 then iterate /*is J divisible by 7? " " " */ if j//11==0 then iterate /* " " " " 11? " " " */ if j//13==0 then iterate /*is " " " 13? " " " */ do k=7 while k*k<=j /*divide by some generated odd primes. */ if j // @.k==0 then iterate j /*Is J divisible by P? Then not prime*/ end /*k*/ /* [↓] a prime (J) has been found. */ nP= nP+1; !.j=1; @.nP= j /*bump P cnt; assign P to @. and !. */ end /*j*/; return</lang>
- output when using the default inputs:
────────────────────────── first 19 circular primes ─────────────────────────── 2 3 5 7 11 13 17 37 79 113 197 199 337 1193 3779 11939 19937 193939 199933
Rust
<lang rust>// [dependencies] // rug = "1.8"
fn is_prime(n : u32) -> bool {
if n < 2 { return false; } if n % 2 == 0 { return n == 2; } if n % 3 == 0 { return n == 3; } let mut p = 5; while p * p <= n { if n % p == 0 { return false; } p += 2; if n % p == 0 { return false; } p += 4; } true
}
fn cycle(n : u32) -> u32 {
let mut m : u32 = n; let mut p : u32 = 1; while m >= 10 { p *= 10; m /= 10; } m + 10 * (n % p)
}
fn is_circular_prime(p : u32) -> bool {
if !is_prime(p) { return false; } let mut p2 : u32 = cycle(p); while p2 != p { if p2 < p || !is_prime(p2) { return false; } p2 = cycle(p2); } true
}
fn test_repunit(digits : usize) {
use rug::{Integer, integer::IsPrime}; let repunit = "1".repeat(digits); let bignum = Integer::from_str_radix(&repunit, 10).unwrap(); if bignum.is_probably_prime(10) != IsPrime::No { println!("R({}) is probably prime.", digits); } else { println!("R({}) is not prime.", digits); }
}
fn main() {
use rug::{Integer, integer::IsPrime}; println!("First 19 circular primes:"); let mut count = 0; let mut p : u32 = 2; while count < 19 { if is_circular_prime(p) { if count > 0 { print!(", "); } print!("{}", p); count += 1; } p += 1; } println!(); println!("Next 4 circular primes:"); let mut repunit : u32 = 1; let mut digits : usize = 1; while repunit < p { repunit = 10 * repunit + 1; digits += 1; } let mut bignum = Integer::from(repunit); count = 0; while count < 4 { if bignum.is_probably_prime(15) != IsPrime::No { if count > 0 { print!(", "); } print!("R({})", digits); count += 1; } digits += 1; bignum = bignum * 10 + 1; } println!(); test_repunit(5003); test_repunit(9887); test_repunit(15073); test_repunit(25031); test_repunit(35317); test_repunit(49081);
}</lang>
- Output:
Execution time is about 805 seconds on my system (3.2 GHz Quad-Core Intel Core i5, macOS 10.15.4).
First 19 circular primes: 2, 3, 5, 7, 11, 13, 17, 37, 79, 113, 197, 199, 337, 1193, 3779, 11939, 19937, 193939, 199933 Next 4 circular primes: R(19), R(23), R(317), R(1031) R(5003) is not prime. R(9887) is not prime. R(15073) is not prime. R(25031) is not prime. R(35317) is not prime. R(49081) is probably prime.
Wren
Just the first part as no 'big integer' support. <lang ecmascript>var isPrime = Fn.new { |n|
if (n < 2 || !n.isInteger) return false if (n%2 == 0) return n == 2 if (n%3 == 0) return n == 3 var d = 5 while (d*d <= n) { if (n%d == 0) return false d = d + 2 if (n%d == 0) return false d = d + 4 } return true
}
var circs = []
var isCircular = Fn.new { |n|
var nn = n var pow = 1 // will eventually contain 10 ^ d where d is number of digits in n while (nn > 0) { pow = pow * 10 nn = (nn/10).floor } nn = n while (true) { nn = nn * 10 var f = (nn/pow).floor // first digit nn = nn + f * (1 - pow) if (circs.contains(nn)) return false if (nn == n) break if (!isPrime.call(nn)) return false } return true
}
System.print("The first 19 circular primes are:") var digits = [1, 3, 7, 9] var q = [1, 2, 3, 5, 7, 9] // queue the numbers to be examined var fq = [1, 2, 3, 5, 7, 9] // also queue the corresponding first digits var count = 0 while (true) {
var f = q[0] // peek first element var fd = fq[0] // peek first digit if (isPrime.call(f) && isCircular.call(f)) { circs.add(f) count = count + 1 if (count == 19) break } q.removeAt(0) // pop first element fq.removeAt(0) // ditto for first digit queue if (f != 2 && f != 5) { // if digits > 1 can't contain a 2 or 5 // add numbers with one more digit to queue // only numbers whose last digit >= first digit need be added for (d in digits) { if (d >= fd) { q.add(f*10+d) fq.add(fd) } } }
} System.print(circs)</lang>
- Output:
The first 19 circular primes are: [2, 3, 5, 7, 11, 13, 17, 37, 79, 113, 197, 199, 337, 1193, 3779, 11939, 19937, 193939, 199933]