Long primes
You are encouraged to solve this task according to the task description, using any language you may know.
A long prime (as defined here) is a prime number whose reciprocal (in decimal) has
a period length of one less than the prime number.
Long primes are also known as:
- base ten cyclic numbers
- full reptend primes
- golden primes
- long period primes
- maximal period primes
- proper primes
Another definition: primes p such that the decimal expansion of 1/p has period p-1, which is the greatest period possible for any integer.
- Example
7 is the first long prime, the reciprocal of seven is 1/7, which is equal to the repeating decimal fraction 0.142857142857···
The length of the repeating part of the decimal fraction
is six, (the underlined part) which is one less
than the (decimal) prime number 7.
Thus 7 is a long prime.
There are other (more) general definitions of a long prime which
include wording/verbiage for bases other than ten.
- Task
-
- Show all long primes up to 500 (preferably on one line).
- Show the number of long primes up to 500
- Show the number of long primes up to 1,000
- Show the number of long primes up to 2,000
- Show the number of long primes up to 4,000
- Show the number of long primes up to 8,000
- Show the number of long primes up to 16,000
- Show the number of long primes up to 32,000
- Show the number of long primes up to 64,000 (optional)
- Show all output here.
- Also see
11l
<lang 11l>F sieve(limit)
[Int] primes V c = [0B] * (limit + 1) V p = 3 L V p2 = p * p I p2 > limit L.break L(i) (p2 .< limit).step(2 * p) c[i] = 1B L p += 2 I !c[p] L.break
L(i) (3 .< limit).step(2) I !(c[i]) primes.append(i) R primes
F findPeriod(n)
V r = 1 L(i) 1 .< n r = (10 * r) % n V rr = r V period = 0 L r = (10 * r) % n period++ I r == rr L.break R period
V primes = sieve(64000) [Int] longPrimes L(prime) primes
I findPeriod(prime) == prime - 1 longPrimes.append(prime)
V numbers = [500, 1000, 2000, 4000, 8000, 16000, 32000, 64000] V count = 0 V index = 0 V totals = [0] * numbers.len L(longPrime) longPrimes
I longPrime > numbers[index] totals[index] = count index++ count++
totals.last = count print(‘The long primes up to 500 are:’) print(String(longPrimes[0 .< totals[0]]).replace(‘,’, ‘’)) print("\nThe number of long primes up to:") L(total) totals
print(‘ #5 is #.’.format(numbers[L.index], total))</lang>
- Output:
The long primes up to 500 are: [7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499] The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
AppleScript
The isLongPrime(n) handler here is a translation of the faster Phix one.
<lang applescript>on sieveOfEratosthenes(limit)
script o property numberList : {missing value} end script repeat with n from 2 to limit set end of o's numberList to n end repeat repeat with n from 2 to (limit ^ 0.5 div 1) if (item n of o's numberList is n) then repeat with multiple from (n * n) to limit by n set item multiple of o's numberList to missing value end repeat end if end repeat return o's numberList's numbers
end sieveOfEratosthenes
on factors(n)
set output to {} if (n < 0) then set n to -n set sqrt to n ^ 0.5 set limit to sqrt div 1 if (limit = sqrt) then set end of output to limit set limit to limit - 1 end if repeat with i from limit to 1 by -1 if (n mod i is 0) then set beginning of output to i set end of output to n div i end if end repeat return output
end factors
on isLongPrime(n)
if (n < 3) then return false script o property f : factors(n - 1) end script set counter to 0 repeat with fi in o's f set fi to fi's contents set e to 1 set base to 10 repeat until (fi = 0) if (fi mod 2 = 1) then set e to e * base mod n set base to base * base mod n set fi to fi div 2 end repeat if (e = 1) then set counter to counter + 1 if (counter > 1) then exit repeat end if end repeat return (counter = 1)
end isLongPrime
-- Task code: on longPrimesTask()
script o -- The isLongPrime() handler above returns the correct result for any number -- passed to it, but feeeding it only primes in the first place speeds things up. property primes : sieveOfEratosthenes(64000) property longs : {} end script set output to {} set counter to 0 set mileposts to {500, 1000, 2000, 4000, 8000, 16000, 32000, 64000} set m to 1 set nextMilepost to beginning of mileposts set astid to AppleScript's text item delimiters repeat with p in o's primes set p to p's contents if (isLongPrime(p)) then -- p being odd, it's never exactly one of the even mileposts. if (p < 500) then set end of o's longs to p else if (p > nextMilepost) then if (nextMilepost = 500) then set AppleScript's text item delimiters to space set end of output to "Long primes up to 500:" set end of output to o's longs as text end if set end of output to "Number of long primes up to " & nextMilepost & ": " & counter set m to m + 1 set nextMilepost to item m of mileposts end if set counter to counter + 1 end if end repeat set end of output to "Number of long primes up to " & nextMilepost & ": " & counter set AppleScript's text item delimiters to linefeed set output to output as text set AppleScript's text item delimiters to astid return output
end longPrimesTask
longPrimesTask()</lang>
- Output:
"Long primes up to 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes up to 500: 35 Number of long primes up to 1000: 60 Number of long primes up to 2000: 116 Number of long primes up to 4000: 218 Number of long primes up to 8000: 390 Number of long primes up to 16000: 716 Number of long primes up to 32000: 1300 Number of long primes up to 64000: 2430"
C
<lang c>#include <stdio.h>
- include <stdlib.h>
- define TRUE 1
- define FALSE 0
typedef int bool;
void sieve(int limit, int primes[], int *count) {
bool *c = calloc(limit + 1, sizeof(bool)); /* composite = TRUE */ /* no need to process even numbers */ int i, p = 3, p2, n = 0; p2 = p * p; while (p2 <= limit) { for (i = p2; i <= limit; i += 2 * p) c[i] = TRUE; do { p += 2; } while (c[p]); p2 = p * p; } for (i = 3; i <= limit; i += 2) { if (!c[i]) primes[n++] = i; } *count = n; free(c);
}
/* finds the period of the reciprocal of n */ int findPeriod(int n) {
int i, r = 1, rr, period = 0; for (i = 1; i <= n + 1; ++i) { r = (10 * r) % n; } rr = r; do { r = (10 * r) % n; period++; } while (r != rr); return period;
}
int main() {
int i, prime, count = 0, index = 0, primeCount, longCount = 0, numberCount; int *primes, *longPrimes, *totals; int numbers[] = {500, 1000, 2000, 4000, 8000, 16000, 32000, 64000}; primes = calloc(6500, sizeof(int)); numberCount = sizeof(numbers) / sizeof(int); totals = calloc(numberCount, sizeof(int)); sieve(64000, primes, &primeCount); longPrimes = calloc(primeCount, sizeof(int)); /* Surely longCount < primeCount */ for (i = 0; i < primeCount; ++i) { prime = primes[i]; if (findPeriod(prime) == prime - 1) { longPrimes[longCount++] = prime; } } for (i = 0; i < longCount; ++i, ++count) { if (longPrimes[i] > numbers[index]) { totals[index++] = count; } } totals[numberCount - 1] = count; printf("The long primes up to %d are:\n", numbers[0]); printf("["); for (i = 0; i < totals[0]; ++i) { printf("%d ", longPrimes[i]); } printf("\b]\n"); printf("\nThe number of long primes up to:\n"); for (i = 0; i < 8; ++i) { printf(" %5d is %d\n", numbers[i], totals[i]); } free(totals); free(longPrimes); free(primes); return 0;
}</lang>
- Output:
The long primes up to 500 are: [7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499] The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
C#
<lang csharp>using System; using System.Collections.Generic; using System.Linq;
public static class LongPrimes {
public static void Main() { var primes = SomePrimeGenerator.Primes(64000).Skip(1).Where(p => Period(p) == p - 1).Append(99999); Console.WriteLine(string.Join(" ", primes.TakeWhile(p => p <= 500))); int count = 0, limit = 500; foreach (int prime in primes) { if (prime > limit) { Console.WriteLine($"There are {count} long primes below {limit}"); limit *= 2; } count++; }
int Period(int n) { int r = 1, rr; for (int i = 0; i <= n; i++) r = 10 * r % n; rr = r; for (int period = 1;; period++) { r = (10 * r) % n; if (r == rr) return period; } } }
}
static class SomePrimeGenerator {
public static IEnumerable<int> Primes(int lim) { bool [] flags = new bool[lim + 1]; int j = 2; for (int d = 3, sq = 4; sq <= lim; j++, sq += d += 2) if (!flags[j]) { yield return j; for (int k = sq; k <= lim; k += j) flags[k] = true; } for (; j<= lim; j++) if (!flags[j]) yield return j; }
}</lang>
- Output:
7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 There are 35 long primes below 500 There are 60 long primes below 1000 There are 116 long primes below 2000 There are 218 long primes below 4000 There are 390 long primes below 8000 There are 716 long primes below 16000 There are 1300 long primes below 32000 There are 2430 long primes below 64000
C++
<lang cpp>
- include <iomanip>
- include <iostream>
- include <list>
using namespace std;
void sieve(int limit, list<int> &primes) {
bool *c = new bool[limit + 1]; for (int i = 0; i <= limit; i++) c[i] = false; // No need to process even numbers int p = 3, n = 0; int p2 = p * p; while (p2 <= limit) { for (int i = p2; i <= limit; i += 2 * p) c[i] = true; do p += 2; while (c[p]); p2 = p * p; } for (int i = 3; i <= limit; i += 2) if (!c[i]) primes.push_back(i); delete [] c;
}
// Finds the period of the reciprocal of n int findPeriod(int n) {
int r = 1, rr, period = 0; for (int i = 1; i <= n + 1; ++i) r = (10 * r) % n; rr = r; do { r = (10 * r) % n; period++; } while (r != rr); return period;
}
int main() {
int count = 0, index = 0; int numbers[] = {500, 1000, 2000, 4000, 8000, 16000, 32000, 64000}; list<int> primes; list<int> longPrimes; int numberCount = sizeof(numbers) / sizeof(int); int *totals = new int[numberCount]; cout << "Please wait." << endl << endl; sieve(64000, primes); for (list<int>::iterator iterPrime = primes.begin(); iterPrime != primes.end(); iterPrime++) { if (findPeriod(*iterPrime) == *iterPrime - 1) longPrimes.push_back(*iterPrime); } for (list<int>::iterator iterLongPrime = longPrimes.begin(); iterLongPrime != longPrimes.end(); iterLongPrime++) { if (*iterLongPrime > numbers[index]) totals[index++] = count; ++count; } totals[numberCount - 1] = count; cout << "The long primes up to " << totals[0] << " are:" << endl; cout << "["; int i = 0; for (list<int>::iterator iterLongPrime = longPrimes.begin(); iterLongPrime != longPrimes.end() && i < totals[0]; iterLongPrime++, i++) { cout << *iterLongPrime << " "; } cout << "\b]" << endl; cout << endl << "The number of long primes up to:" << endl; for (int i = 0; i < 8; ++i) cout << " " << setw(5) << numbers[i] << " is " << totals[i] << endl; delete [] totals;
} </lang>
- Output:
Please wait. The long primes up to 35 are: [7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499] The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
Common Lisp
<lang lisp>; Primality test using the Sieve of Eratosthenes with a couple minor optimizations (defun primep (n)
(cond ((and (<= n 3) (> n 1)) t) ((some #'zerop (mapcar (lambda (d) (mod n d)) '(2 3))) nil) (t (loop for i = 5 then (+ i 6) while (<= (* i i) n) when (some #'zerop (mapcar (lambda (d) (mod n (+ i d))) '(0 2))) return nil finally (return t)))))
- Translation of the long-prime algorithm from the Raku solution
(defun long-prime-p (n)
(cond ((< n 3) nil) ((not (primep n)) nil) (t (let* ((rr (loop repeat (1+ n) for r = 1 then (mod (* 10 r) n) finally (return r)))
(period (loop for p = 0 then (1+ p) for r = (mod (* 10 rr) n) then (mod (* 10 r) n) while (and (< p n) (/= r rr)) finally (return (1+ p)))))
(= period (1- n))))))
(format t "~{~a~^, ~}" (loop for n from 1 to 500 if (long-prime-p n) collect n)) </lang>
- Output:
7, 17, 19, 23, 29, 47, 59, 61, 97, 109, 113, 131, 149, 167, 179, 181, 193, 223, 229, 233, 257, 263, 269, 313, 337, 367, 379, 383, 389, 419, 433, 461, 487, 491, 499
Crystal
Simpler but slower. <lang ruby>require "big"
def prime?(n) # P3 Prime Generator primality test
return n | 1 == 3 if n < 5 # n: 2,3|true; 0,1,4|false return false if n.gcd(6) != 1 # this filters out 2/3 of all integers pc = typeof(n).new(5) # first P3 prime candidates sequence value until pc*pc > n return false if n % pc == 0 || n % (pc + 2) == 0 # if n is composite pc += 6 # 1st prime candidate for next residues group end true
end
- The smallest divisor d of p-1 such that 10^d = 1 (mod p),
- is the length of the period of the decimal expansion of 1/p.
def long_prime?(p)
return false unless prime? p (2...p).each do |d| return d == (p - 1) if (p - 1) % d == 0 && (10.to_big_i ** d) % p == 1 end false
end
start = Time.monotonic # time of starting puts "Long primes ≤ 500:" (2..500).each { |pc| print "#{pc} " if long_prime? pc } puts [500, 1000, 2000, 4000, 8000, 16000, 32000, 64000].each do |n|
puts "Number of long primes ≤ #{n}: #{(7..n).count { |pc| long_prime? pc }}"
end puts "\nTime: #{(Time.monotonic - start).total_seconds} secs"</lang>
- Output:
Long primes ≤ 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes ≤ 500: 35 Number of long primes ≤ 1000: 60 Number of long primes ≤ 2000: 116 Number of long primes ≤ 4000: 218 Number of long primes ≤ 8000: 390 Number of long primes ≤ 16000: 716 Number of long primes ≤ 32000: 1300 Number of long primes ≤ 64000: 2430
System: I7-6700HQ, 3.5 GHz, Linux Kernel 5.6.17, Crystal 0.35 Run as: $ crystal run longprimes.cr --release Time: 1.090748711 secs
Faster: using divisors of (p - 1) and powmod(). <lang ruby>require "big"
def prime?(n) # P3 Prime Generator primality test
n = n.to_big_i return n | 1 == 3 if n < 5 # n: 0,1,4|false, 2,3|true return false if n.gcd(6) != 1 # 1/3 (2/6) of integers are P3 pc p = typeof(n).new(5) # first P3 sequence value until p*p > n return false if n % p == 0 || n % (p + 2) == 0 # if n is composite p += 6 # first prime candidate for next kth residues group end true
end
def powmod(b, e, m) # Compute b**e mod m
r, b = 1, b.to_big_i while e > 0 r = (b * r) % m if e.odd? b = (b * b) % m e >>= 1 end r
end
def divisors(n) # divisors of n -> [1,..,n]
f = [] of Int32 (1..Math.sqrt(n)).each { |i| (n % i).zero? && (f << i; f << n // i if n // i != i) } f.sort
end
- The smallest divisor d of p-1 such that 10^d = 1 (mod p),
- is the length of the period of the decimal expansion of 1/p.
def long_prime?(p)
return false unless prime? p divisors(p - 1).each { |d| return d == (p - 1) if powmod(10, d, p) == 1 } false
end
start = Time.monotonic # time of starting puts "Long primes ≤ 500:" (7..500).each { |pc| print "#{pc} " if long_prime? pc } puts [500, 1000, 2000, 4000, 8000, 16000, 32000, 64000].each do |n|
puts "Number of long primes ≤ #{n}: #{(7..n).count { |pc| long_prime? pc }}"
end puts "\nTime: #{(Time.monotonic - start).total_seconds} secs"</lang>
- Output:
Long primes ≤ 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes ≤ 500: 35 Number of long primes ≤ 1000: 60 Number of long primes ≤ 2000: 116 Number of long primes ≤ 4000: 218 Number of long primes ≤ 8000: 390 Number of long primes ≤ 16000: 716 Number of long primes ≤ 32000: 1300 Number of long primes ≤ 64000: 2430
System: I7-6700HQ, 3.5 GHz, Linux Kernel 5.6.17, Crystal 0.35 Run as: $ crystal run longprimes.cr --release Time: 0.28927228 secs
Delphi
See Pascal.
F#
The Functions
This task uses Extensible Prime Generator (F#)
This task uses Factors_of_an_integer#F.23
<lang fsharp>
// Return true if prime n is a long prime. Nigel Galloway: September 25th., 2018
let fN n g = let rec fN i g e l = match e with | 0UL -> i
| _ when e%2UL = 1UL -> fN ((i*g)%l) ((g*g)%l) (e/2UL) l | _ -> fN i ((g*g)%l) (e/2UL) l fN 1UL 10UL (uint64 g) (uint64 n)
let isLongPrime n=Seq.length (factors (n-1) |> Seq.filter(fun g->(fN n g)=1UL))=1 </lang>
The Task
<lang fsharp> primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<500) |> Seq.filter isLongPrime |> Seq.iter(printf "%d ") </lang>
- Output:
7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499
<lang fsharp> printfn "There are %d long primes less than 500" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<500) |> Seq.filter isLongPrime |> Seq.length) </lang>
- Output:
There are 35 long primes less than 500
<lang fsharp> printfn "There are %d long primes less than 1000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<1000) |> Seq.filter isLongPrime |> Seq.length) </lang>
- Output:
There are 60 long primes less than 1000
<lang fsharp> printfn "There are %d long primes less than 2000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<2000) |> Seq.filter isLongPrime |> Seq.length) </lang>
- Output:
There are 116 long primes less than 2000
<lang fsharp> printfn "There are %d long primes less than 4000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<4000) |> Seq.filter isLongPrime|> Seq.length) </lang>
- Output:
There are 218 long primes less than 4000
<lang fsharp> printfn "There are %d long primes less than 8000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<8000) |> Seq.filter isLongPrime |> Seq.length) </lang>
- Output:
There are 390 long primes less than 8000
<lang fsharp> printfn "There are %d long primes less than 16000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<16000) |> Seq.filter isLongPrime |> Seq.length) </lang>
- Output:
There are 716 long primes less than 16000
<lang fsharp> printfn "There are %d long primes less than 32000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<32000) |> Seq.filter isLongPrime |> Seq.length) </lang>
- Output:
There are 1300 long primes less than 32000
<lang fsharp> printfn "There are %d long primes less than 64000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<64000) |> Seq.filter isLongPrime|> Seq.length) </lang>
- Output:
There are 2430 long primes less than 64000
<lang fsharp> printfn "There are %d long primes less than 128000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<128000) |> Seq.filter isLongPrime|> Seq.length) </lang>
- Output:
There are 4498 long primes less than 128000 Real: 00:00:01.294, CPU: 00:00:01.300, GC gen0: 27, gen1: 0
<lang fsharp> printfn "There are %d long primes less than 256000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<256000) |> Seq.filter isLongPrime|> Seq.length) </lang>
- Output:
There are 8434 long primes less than 256000 Real: 00:00:03.434, CPU: 00:00:03.440, GC gen0: 58, gen1: 0
<lang fsharp> printfn "There are %d long primes less than 512000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<512000) |> Seq.filter isLongPrime|> Seq.length) </lang>
- Output:
There are 15920 long primes less than 512000 Real: 00:00:09.248, CPU: 00:00:09.260, GC gen0: 128, gen1: 0
<lang fsharp> printfn "There are %d long primes less than 1024000" (primes |> Seq.skip 3 |> Seq.takeWhile(fun n->n<1024000) |> Seq.filter isLongPrime|> Seq.length) </lang>
- Output:
There are 30171 long primes less than 1024000 Real: 00:00:24.959, CPU: 00:00:25.020, GC gen0: 278, gen1: 1
Factor
<lang factor>USING: formatting fry io kernel math math.functions math.primes math.primes.factors memoize prettyprint sequences ; IN: rosetta-code.long-primes
- period-length ( p -- len )
[ 1 - divisors ] [ '[ 10 swap _ ^mod 1 = ] ] bi find nip ;
MEMO: long-prime? ( p -- ? ) [ period-length ] [ 1 - ] bi = ;
- .lp<=500 ( -- )
500 primes-upto [ long-prime? ] filter "Long primes <= 500:" print [ pprint bl ] each nl ;
- .#lp<=n ( n -- )
dup primes-upto [ long-prime? t = ] count swap "%-4d long primes <= %d\n" printf ;
- long-primes-demo ( -- )
.lp<=500 nl { 500 1,000 2,000 4,000 8,000 16,000 32,000 64,000 } [ .#lp<=n ] each ;
MAIN: long-primes-demo</lang>
- Output:
Long primes <= 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 35 long primes <= 500 60 long primes <= 1000 116 long primes <= 2000 218 long primes <= 4000 390 long primes <= 8000 716 long primes <= 16000 1300 long primes <= 32000 2430 long primes <= 64000
Forth
The prime sieve code was borrowed from Sieve of Eratosthenes#Forth. <lang forth>: prime? ( n -- ? ) here + c@ 0= ;
- notprime! ( n -- ) here + 1 swap c! ;
- sieve ( n -- )
here over erase 0 notprime! 1 notprime! 2 begin 2dup dup * > while dup prime? if 2dup dup * do i notprime! dup +loop then 1+ repeat 2drop ;
- modpow { c b a -- a^b mod c }
c 1 = if 0 exit then 1 a c mod to a begin b 0> while b 1 and 1 = if a * c mod then a a * c mod to a b 2/ to b repeat ;
- divide_out ( n1 n2 -- n )
begin 2dup mod 0= while tuck / swap repeat drop ;
- long_prime? ( n -- ? )
dup prime? invert if drop false exit then 10 over mod 0= if drop false exit then dup 1- 2 >r begin over r@ dup * > while r@ prime? if dup r@ mod 0= if over dup 1- r@ / 10 modpow 1 = if 2drop rdrop false exit then r@ divide_out then then r> 1+ >r repeat rdrop dup 1 = if 2drop true exit then over 1- swap / 10 modpow 1 <> ;
- next_long_prime ( n -- n )
begin 2 + dup long_prime? until ;
500 constant limit1 512000 constant limit2
- main
limit2 1+ sieve limit2 limit1 3 0 >r ." Long primes up to " over 1 .r ." :" cr begin 2 pick over > while next_long_prime dup limit1 < if dup . then dup 2 pick > if over limit1 = if cr then ." Number of long primes up to " over 6 .r ." : " r@ 5 .r cr swap 2* swap then r> 1+ >r repeat 2drop drop rdrop ;
main bye</lang>
- Output:
Execution time is about 1.1 seconds on my machine (3.2GHz Quad-Core Intel Core i5).
Long primes up to 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes up to 500: 35 Number of long primes up to 1000: 60 Number of long primes up to 2000: 116 Number of long primes up to 4000: 218 Number of long primes up to 8000: 390 Number of long primes up to 16000: 716 Number of long primes up to 32000: 1300 Number of long primes up to 64000: 2430 Number of long primes up to 128000: 4498 Number of long primes up to 256000: 8434 Number of long primes up to 512000: 15920
FreeBASIC
<lang freebasic>' version 01-02-2019 ' compile with: fbc -s console
Dim Shared As UByte prime()
Sub find_primes(n As UInteger)
ReDim prime(n) Dim As UInteger i, k
' need only to consider odd primes, 2 has no repetion For i = 3 To n Step 2 If prime(i) = 0 Then For k = i * i To n Step i + i prime(k) = 1 Next End If Next
End Sub
Function find_period(p As UInteger) As UInteger
' finds period for every positive number Dim As UInteger period, r = 1
Do r = (r * 10) Mod p period += 1 If r <= 1 Then Return period Loop
End Function
' ------=< MAIN >=------
- Define max 64000
Dim As UInteger p = 3, n1 = 3, n2 = 500, i, n50, count
find_primes(max) Print "Long primes upto 500 are ";
For i = n1 To n2 Step 2
If prime(i) = 0 Then If i -1 = find_period(i) Then If n50 <= 50 Then Print Str(i); " "; End If count += 1 End If End If
Next
Print : Print
Do
Print "There are "; Str(count); " long primes upto "; Str(n2)
n1 = n2 +1 n2 += n2 If n1 > max Then Exit Do
For i = n1 To n2 Step 2 If prime(i) = 0 Then If i -1 = find_period(i) Then count += 1 End If End If Next
Loop
' empty keyboard buffer While Inkey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>
- Output:
Long primes upto 500 are 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 There are 35 long primes upto 500 There are 60 long primes upto 1000 There are 116 long primes upto 2000 There are 218 long primes upto 4000 There are 390 long primes upto 8000 There are 716 long primes upto 16000 There are 1300 long primes upto 32000 There are 2430 long primes upto 64000
Go
<lang go>package main
import "fmt"
func sieve(limit int) []int {
var primes []int c := make([]bool, limit + 1) // composite = true // no need to process even numbers p := 3 p2 := p * p for p2 <= limit { for i := p2; i <= limit; i += 2 * p { c[i] = true } for ok := true; ok; ok = c[p] { p += 2 } p2 = p * p } for i := 3; i <= limit; i += 2 { if !c[i] { primes = append(primes, i) } } return primes
}
// finds the period of the reciprocal of n func findPeriod(n int) int {
r := 1 for i := 1; i <= n + 1; i++ { r = (10 * r) % n } rr := r period := 0 for ok := true; ok; ok = r != rr { r = (10 * r) % n period++ } return period
}
func main() {
primes := sieve(64000) var longPrimes []int for _, prime := range primes { if findPeriod(prime) == prime - 1 { longPrimes = append(longPrimes, prime) } } numbers := []int{500, 1000, 2000, 4000, 8000, 16000, 32000, 64000} index := 0 count := 0 totals := make([]int, len(numbers)) for _, longPrime := range longPrimes { if longPrime > numbers[index] { totals[index] = count index++ } count++ } totals[len(numbers)-1] = count fmt.Println("The long primes up to", numbers[0], "are: ") fmt.Println(longPrimes[:totals[0]]) fmt.Println("\nThe number of long primes up to: ") for i, total := range totals { fmt.Printf(" %5d is %d\n", numbers[i], total) }
}</lang>
- Output:
The long primes up to 500 are: [7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499] The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
Haskell
<lang Haskell>import Data.List (elemIndex)
longPrimesUpTo :: Int -> [Int] longPrimesUpTo n =
filter isLongPrime $ takeWhile (< n) primes where sieve (p : xs) = p : sieve [x | x <- xs, x `mod` p /= 0] primes = sieve [2 ..] isLongPrime n = found where cycles = take n (iterate ((`mod` n) . (10 *)) 1) index = elemIndex (head cycles) $ tail cycles found = case index of (Just i) -> n - i == 2 _ -> False
display :: Int -> IO () display n =
if n <= 64000 then do putStrLn ( show n <> " is " <> show (length $ longPrimesUpTo n) ) display (n * 2) else pure ()
main :: IO () main = do
let fiveHundred = longPrimesUpTo 500 putStrLn ( "The long primes up to 35 are:\n" <> show fiveHundred <> "\n" ) putStrLn ("500 is " <> show (length fiveHundred)) display 1000</lang>
- Output:
The long primes up to 35 are: [7,17,19,23,29,47,59,61,97,109,113,131,149,167,179,181,193,223,229,233,257,263,269,313,337,367,379,383,389,419,433,461,487,491,499] 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
J
NB. define the verb long NB. long is true iff the prime input greater than 2 NB. is a rosettacode long prime. NB. 0 is false, 1 is true. long =: ( <:@:[ = #@~.@( [: }. ( | 10&* )^:( <@[ ) ) )&1&> NB. demonstration of the long verb NB. long applied to integers 3 through 9 inclusively (,: long) 3 4 5 6 7 8 9 3 4 5 6 7 8 9 0 0 0 0 1 0 0 NB. find the number of primes through 64000 [ N =: p:^:_1 ] 64000 6413 NB. copy the long primes, excluding 2, the first. LONG_PRIMES =: (#~ long) p: >: i. N NB. those less than 500 ( #~ <&500) LONG_PRIMES 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 NB. counts [ MEASURE =: 500 * 2 ^ i. 8 500 1000 2000 4000 8000 16000 32000 64000 LONG_PRIMES ( ] ,: [: +/ </ ) MEASURE 500 1000 2000 4000 8000 16000 32000 64000 35 60 116 218 390 716 1300 2430
Java
<lang java> import java.util.LinkedList; import java.util.List;
public class LongPrimes {
private static void sieve(int limit, List<Integer> primes) { boolean[] c = new boolean[limit]; for (int i = 0; i < limit; i++) c[i] = false; // No need to process even numbers int p = 3, n = 0; int p2 = p * p; while (p2 <= limit) { for (int i = p2; i <= limit; i += 2 * p) c[i] = true; do p += 2; while (c[p]); p2 = p * p; } for (int i = 3; i <= limit; i += 2) if (!c[i]) primes.add(i); }
// Finds the period of the reciprocal of n private static int findPeriod(int n) { int r = 1, period = 0; for (int i = 1; i < n; i++) r = (10 * r) % n; int rr = r; do { r = (10 * r) % n; ++period; } while (r != rr); return period; } public static void main(String[] args) { int[] numbers = new int[]{500, 1000, 2000, 4000, 8000, 16000, 32000, 64000}; int[] totals = new int[numbers.length]; List<Integer> primes = new LinkedList<Integer>(); List<Integer> longPrimes = new LinkedList<Integer>(); sieve(64000, primes); for (int prime : primes) if (findPeriod(prime) == prime - 1) longPrimes.add(prime); int count = 0, index = 0; for (int longPrime : longPrimes) { if (longPrime > numbers[index]) totals[index++] = count; ++count; } totals[numbers.length - 1] = count; System.out.println("The long primes up to " + numbers[0] + " are:"); System.out.println(longPrimes.subList(0, totals[0])); System.out.println(); System.out.println("The number of long primes up to:"); for (int i = 0; i <= 7; i++) System.out.printf(" %5d is %d\n", numbers[i], totals[i]); }
} </lang>
- Output:
The long primes up to 500 are: [7, 17, 19, 23, 29, 47, 59, 61, 97, 109, 113, 131, 149, 167, 179, 181, 193, 223, 229, 233, 257, 263, 269, 313, 337, 367, 379, 383, 389, 419, 433, 461, 487, 491, 499] The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
jq
Adapted from Wren
Works with gojq, the Go implementation of jq (*)
This entry does not attempt to avoid the redundancy involved in the subtasks, but instead includes a prime number generator intended for efficiently generating large numbers of primes.
(*) For the computationally intensive subtasks, gojq will require too much memory.
Preliminaries <lang jq> def count(s): reduce s as $x (0; .+1);
- Is the input integer a prime?
- "previous" should be a sorted array of consecutive primes
- from 2 on that includes the greatest prime less than (.|sqrt)
def is_prime(previous):
. as $in | (($in + 1) | sqrt) as $sqrt | first(previous[] | if . > $sqrt then 1 elif 0 == ($in % .) then 0 else empty end) // 1 | . == 1;
- This assumes . is an array of consecutive primes beginning with [2,3]
def next_prime:
. as $previous | (2 + .[-1] ) | until(is_prime($previous); . + 2) ;
- Emit primes from 2 up
def primes:
# The helper function has arity 0 for TCO # It expects its input to be an array of previously found primes, in order: def next: . as $previous | ($previous|next_prime) as $next | $next, (($previous + [$next]) | next) ; 2, 3, ([2,3] | next);
</lang> Long Primes <lang jq>
- finds the period of the reciprocal of .
- (The following definition does not make a special case of 2
- but yields a justifiable result for 2, namely 1.)
def findPeriod:
. as $n | (reduce range(1; $n+2) as $i (1; (. * 10) % $n)) as $rr | {r: $rr, period:0, ok:true} | until( .ok|not; .r = (10 * .r) % $n | .period += 1 | .ok = (.r != $rr) ) | .period ;
- This definition takes into account the
- claim in the preamble that the first long prime is 7:
def long_primes_less_than($n):
label $out | primes | if . >= $n then break $out else . end | select(. > 2 and (findPeriod == . - 1));
def count_long_primes:
count(long_primes_less_than(.));
- Since 2 is not a "long prime" for the purposes of this
- article, we can begin searching at 3:
"Long primes ≤ 500: ", long_primes_less_than(500),
"\n",
(500,1000, 2000, 4000, 8000, 16000, 32000, 64000 | "Number of long primes ≤ \(.): \(count_long_primes)" )
</lang>
- Output:
Long primes ≤ 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes ≤ 500: 35 Number of long primes ≤ 1000: 60 Number of long primes ≤ 2000: 116 Number of long primes ≤ 4000: 218 Number of long primes ≤ 8000: 390 Number of long primes ≤ 16000: 716 Number of long primes ≤ 32000: 1300 Number of long primes ≤ 64000: 2430
Julia
<lang julia> using Primes
function divisors(n)
f = [one(n)] for (p,e) in factor(n) f = reduce(vcat, [f*p^j for j in 1:e], init=f) end return length(f) == 1 ? [one(n), n] : sort!(f)
end
function islongprime(p)
for i in divisors(p-1) if powermod(10, i, p) == 1 return i + 1 == p end end false
end
println("Long primes ≤ 500: ") for i in 2:500
if islongprime(i) i == 229 ? println(i) : print(i, " ") end
end print("\n\n")
for i in [500, 1000, 2000, 4000, 8000, 16000, 32000, 64000]
println("Number of long primes ≤ $i: $(sum(map(x->islongprime(x), 1:i)))")
end </lang>
- Output:
Long primes ≤ 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499
Number of long primes ≤ 500: 35 Number of long primes ≤ 1000: 60 Number of long primes ≤ 2000: 116 Number of long primes ≤ 4000: 218 Number of long primes ≤ 8000: 390 Number of long primes ≤ 16000: 716 Number of long primes ≤ 32000: 1300 Number of long primes ≤ 64000: 2430
Kotlin
<lang scala>// Version 1.2.60
fun sieve(limit: Int): List<Int> {
val primes = mutableListOf<Int>() val c = BooleanArray(limit + 1) // composite = true // no need to process even numbers var p = 3 var p2 = p * p while (p2 <= limit) { for (i in p2..limit step 2 * p) c[i] = true do { p += 2 } while (c[p]) p2 = p * p } for (i in 3..limit step 2) { if (!c[i]) primes.add(i) } return primes
}
// finds the period of the reciprocal of n fun findPeriod(n: Int): Int {
var r = 1 for (i in 1..n + 1) r = (10 * r) % n val rr = r var period = 0 do { r = (10 * r) % n period++ } while (r != rr) return period
}
fun main(args: Array<String>) {
val primes = sieve(64000) val longPrimes = mutableListOf<Int>() for (prime in primes) { if (findPeriod(prime) == prime - 1) { longPrimes.add(prime) } } val numbers = listOf(500, 1000, 2000, 4000, 8000, 16000, 32000, 64000) var index = 0 var count = 0 val totals = IntArray(numbers.size) for (longPrime in longPrimes) { if (longPrime > numbers[index]) { totals[index++] = count } count++ } totals[numbers.lastIndex] = count println("The long primes up to " + numbers[0] + " are:") println(longPrimes.take(totals[0])) println("\nThe number of long primes up to:") for ((i, total) in totals.withIndex()) { System.out.printf(" %5d is %d\n", numbers[i], total) }
}</lang>
- Output:
The long primes up to 500 are: [7, 17, 19, 23, 29, 47, 59, 61, 97, 109, 113, 131, 149, 167, 179, 181, 193, 223, 229, 233, 257, 263, 269, 313, 337, 367, 379, 383, 389, 419, 433, 461, 487, 491, 499] The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
M2000 Interpreter
Sieve leave to stack of values primes. This happen because we call the function as a module, so we pass the current stack (modules call modules passing own stack of values). We can place value to stack using Push to the top (as LIFO) or using Data to bottom (as FIFO). Variable Number read a number from stack and drop it.
<lang M2000 Interpreter> Module LongPrimes {
Sieve=lambda (limit)->{ Flush Buffer clear c as byte*limit+1 \\ no need to process even numbers p=3 do p2=p^2 if p2>limit then exit i=p2 while i<=limit Return c, i:=1 i+=2*p end While do p+=2 Until not eval(c,p) always for i = 3 to limit step 2 if eval(c,i) else data i next i } findPeriod=lambda (n) -> { r = 1 for i = 1 to n+1 {r = (10 * r) mod n} rr = r : period = 0 do r = (10 * r) mod n period++ if r == rr then exit always =period } Call sieve(64000) ' leave stack with primes stops=(500,1000,2000,4000,8000,16000,32000,64000) acc=0 stp=0 limit=array(stops, stp) p=number ' pop one Print "Long primes up to 500:" document lp500$ for i=1 to 500 if i=p then if findPeriod(i)=i-1 then acc++ :lp500$=str$(i) p=number end if if empty then exit for next i lp500$="]" insert 1,1 lp500$="[" Print lp500$ Print i=500 Print "The number of long primes up to:" print i," is ";acc stp++ m=each(stops,1,-2) while m for i=array(m)+1 to array(m,m^+1) if i=p then if findPeriod(i)=i-1 then acc++ p=number end if if empty then exit for next i print array(m,m^+1)," is ";acc end While
} LongPrimes </lang>
- Output:
The long primes up to 500 are: [7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499] The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
Maple
<lang Maple> with(NumberTheory): with(ArrayTools):
isLong := proc(x::integer)
if irem(10^(x - 1) - 1, x) = 0 then for local count from 1 to x - 2 do if irem(10^(count) - 1, x) = 0 then return false; end if; end do; else return false; end if; return true;
end proc:
longPrimes := Array([]):
for number from 1 to PrimeCounting(500) do
if isLong(ithprime(number)) then Append(longPrimes, ithprime(number)): end if:
end:
longPrimes; lpcount := ArrayNumElems(longPrimes): numOfLongPrimes := Array([lpcount]): for expon from 1 to 7 do
for number from PrimeCounting(500 * 2^(expon - 1)) + 1 to PrimeCounting(500 * 2^expon) do if isLong(ithprime(number)) then lpcount += 1: end if: end: Append(numOfLongPrimes, lpcount):
end:
numOfLongPrimes;
</lang>
- Output:
[7, 17, 19, 23, 29, 47, 59, 61, 97, 109, 113, 131, 149, 167, 179, 181, 193,
223, 229, 233, 257, 263, 269, 313, 337, 367, 379, 383, 389, 419, 433, 461,487, 491, 499][35, 60, 116, 218, 390, 716, 1300, 2430]
Mathematica/Wolfram Language
<lang Mathematica>lPrimes[n_] := Select[Range[2, n], Length[RealDigits[1/#]1, 1] == # - 1 &]; lPrimes[500] Length /@ lPrimes /@ ( 250*2^Range[8])</lang>
- Output:
{7, 17, 19, 23, 29, 47, 59, 61, 97, 109, 113, 131, 149, 167, 179, 181, 193, 223, 229, 233, 257, 263, 269, 313, 337, 367, 379, 383, 389, 419, 433, 461, 487, 491, 499} {35, 60, 116, 218, 390, 716, 1300, 2430}
NewLISP
<lang NewLISP>
- Using the fact that 10 has to be a primitive root mod p
- for p to be a reptend/long prime.
- p supposed prime and >= 7
(define (cycle-mod p) (let (n 10 tally 1) (while (!= n 1) (++ tally) (setq n (% (* n 10) p)) tally)))
- Primality test
(define (prime? n) (= (length (factor n)) 1))
- Reptend test (p >= 7)
(define (reptend? p) (if (prime? p) (= (- p (cycle-mod p)) 1) false))
- Find reptends in interval 7 .. n
(define (find-reptends n) (filter reptend? (sequence 7 n)))
- Task
(println (find-reptends 500)) (println (map (fn(n) (println n " --> " (length (find-reptends n)))) '(500 1000 2000 4000 8000 16000 32000 64000))) </lang>
- Output:
(7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499) 500 --> 35 1000 --> 60 2000 --> 116 4000 --> 218 8000 --> 390 16000 --> 716 32000 --> 1300 64000 --> 2430
Nim
<lang Nim>import strformat
func sieve(limit: int): seq[int] =
var composite = newSeq[bool](limit + 1) var p = 3 var p2 = p * p while p2 < limit: if not composite[p]: for n in countup(p2, limit, 2 * p): composite[n] = true inc p, 2 p2 = p * p
for n in countup(3, limit, 2): if not composite[n]: result.add n
func period(n: int): int =
## Find the period of the reciprocal of "n". var r = 1 for i in 1..(n + 1): r = 10 * r mod n let r1 = r while true: r = 10 * r mod n inc result if r == r1: break
let primes = sieve(64000)
var longPrimes: seq[int]
for prime in primes:
if prime.period() == prime - 1: longPrimes.add prime
const Numbers = [500, 1000, 2000, 4000, 8000, 16000, 32000, 64000] var index, count = 0 var totals = newSeq[int](Numbers.len) for longPrime in longPrimes:
if longPrime > Numbers[index]: totals[index] = count inc index inc count
totals[^1] = count
echo &"The long primes up to {Numbers[0]} are:" for i in 0..<totals[0]:
stdout.write ' ', longPrimes[i]
stdout.write '\n'
echo "\nThe number of long primes up to:" for i, total in totals:
echo &" {Numbers[i]:>5} is {total}"</lang>
- Output:
The long primes up to 500 are: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
Pascal
first post.old program modified. Using Euler Phi
www . arndt-bruenner.de/mathe/scripts/periodenlaenge.htm
<lang pascal> program Periode;
{$IFDEF FPC}
{$MODE Delphi} {$OPTIMIZATION ON} {$OPTIMIZATION Regvar} {$OPTIMIZATION Peephole} {$OPTIMIZATION cse} {$OPTIMIZATION asmcse}
{$else}
{$Apptype Console}
{$ENDIF}
uses
sysutils;
const
cBASIS = 10; PRIMFELDOBERGRENZE = 6542; {Das sind alle Primzahlen bis 2^16} {Das reicht fuer al8le Primzahlen bis 2^32} TESTZAHL = 500; //429496709;//High(Cardinal) DIV cBasis;
type
tPrimFeld = array[1..PRIMFELDOBERGRENZE] of Word;
tFaktorPotenz = record Faktor, Potenz: Cardinal; end; //2*3*5*7*11*13*17*19*23 *29 > Cardinal also maximal 9 Faktoren
tFaktorFeld = array[1..9] of TFaktorPotenz; //Cardinal
// tFaktorFeld = array [1..15] of TFaktorPotenz;//QWord
tFaktorisieren = class(TObject) private fFakZahl: Cardinal; fFakBasis: Cardinal; fFakAnzahl: Cardinal; fAnzahlMoeglicherTeiler: Cardinal; fEulerPhi: Cardinal; fStartPeriode: Cardinal; fPeriodenLaenge: Cardinal; fTeiler: array of Cardinal; fFaktoren: tFaktorFeld; fBasFakt: tFaktorFeld; fPrimfeld: tPrimFeld; procedure PrimFeldAufbauen; procedure Fakteinfuegen(var Zahl: Cardinal; inFak: Cardinal); function BasisPeriodeExtrahieren(var inZahl: Cardinal): Cardinal; procedure NachkommaPeriode(var OutText: string); public constructor create; overload; function Prim(inZahl: Cardinal): Boolean; procedure AusgabeFaktorfeld(n: Cardinal); procedure Faktorisierung(inZahl: Cardinal); procedure TeilerErmitteln; procedure PeriodeErmitteln(inZahl: Cardinal); function BasExpMod(b, e, m: Cardinal): Cardinal; property EulerPhi: Cardinal read fEulerPhi; property PeriodenLaenge: Cardinal read fPeriodenLaenge; property StartPeriode: Cardinal read fStartPeriode; end;
constructor tFaktorisieren.create; begin
inherited; PrimFeldAufbauen;
fFakZahl := 0; fFakBasis := cBASIS; Faktorisierung(fFakBasis); fBasFakt := fFaktoren;
fFakZahl := 0; fEulerPhi := 1; fPeriodenLaenge := 0; fFakZahl := 0; fFakAnzahl := 0; fAnzahlMoeglicherTeiler := 0;
end;
function tFaktorisieren.Prim(inZahl: Cardinal): Boolean; {Testet auf PrimZahl} var
Wurzel, pos: Cardinal;
begin
if fFakZahl = inZahl then begin result := (fAnzahlMoeglicherTeiler = 2); exit; end; result := false; if inZahl > 1 then begin result := true; pos := 1; Wurzel := trunc(sqrt(inZahl)); while fPrimFeld[pos] <= Wurzel do begin if (inZahl mod fPrimFeld[pos]) = 0 then begin result := false; break; end; inc(pos); if pos > High(fPrimFeld) then break; end; end;
end;
procedure tFaktorisieren.PrimFeldAufbauen; {Baut die Liste der Primzahlen bis Obergrenze auf} const
MAX = 65536;
var
TestaufPrim, Zaehler, delta: Cardinal;
begin
Zaehler := 1; fPrimFeld[Zaehler] := 2; inc(Zaehler); fPrimFeld[Zaehler] := 3;
delta := 2; TestaufPrim := 5; repeat if prim(TestaufPrim) then begin inc(Zaehler); fPrimFeld[Zaehler] := TestaufPrim; end; inc(TestaufPrim, delta); delta := 6 - delta; // 2,4,2,4,2,4,2, until (TestaufPrim >= MAX);
end; {PrimfeldAufbauen}
procedure tFaktorisieren.Fakteinfuegen(var Zahl: Cardinal; inFak: Cardinal); var
i: Cardinal;
begin
inc(fFakAnzahl); with fFaktoren[fFakAnzahl] do begin fEulerPhi := fEulerPhi * (inFak - 1); Faktor := inFak; Potenz := 0; while (Zahl mod inFak) = 0 do begin Zahl := Zahl div inFak; inc(Potenz); end; for i := 2 to Potenz do fEulerPhi := fEulerPhi * inFak; end; fAnzahlMoeglicherTeiler := fAnzahlMoeglicherTeiler * (1 + fFaktoren[fFakAnzahl].Potenz);
end;
procedure tFaktorisieren.Faktorisierung(inZahl: Cardinal); var
j, og: longint;
begin
if fFakZahl = inZahl then exit;
fPeriodenLaenge := 0; fFakZahl := inZahl; fEulerPhi := 1; fFakAnzahl := 0; fAnzahlMoeglicherTeiler := 1; setlength(fTeiler, 0);
if inZahl < 2 then exit; og := round(sqrt(inZahl) + 1.0);
{Suche Teiler von inZahl}
for j := 1 to High(fPrimfeld) do begin if fPrimfeld[j] > og then Break; if (inZahl mod fPrimfeld[j]) = 0 then Fakteinfuegen(inZahl, fPrimfeld[j]); end; if inZahl > 1 then Fakteinfuegen(inZahl, inZahl); TeilerErmitteln;
end; {Faktorisierung}
procedure tFaktorisieren.AusgabeFaktorfeld(n: Cardinal); var
i: integer;
begin
if fFakZahl <> n then Faktorisierung(n); write(fAnzahlMoeglicherTeiler: 5, ' Faktoren ');
for i := 1 to fFakAnzahl - 1 do with fFaktoren[i] do if potenz > 1 then write(Faktor, '^', Potenz, '*') else write(Faktor, '*'); with fFaktoren[fFakAnzahl] do if potenz > 1 then write(Faktor, '^', Potenz) else write(Faktor);
writeln(' Euler Phi: ', fEulerPhi: 12, PeriodenLaenge: 12);
end;
procedure tFaktorisieren.TeilerErmitteln; var
Position: Cardinal; i, j: Cardinal;
procedure FaktorAufbauen(Faktor: Cardinal; n: Cardinal); var i, Pot: Cardinal; begin Pot := 1; i := 0; repeat if n > Low(fFaktoren) then FaktorAufbauen(Pot * Faktor, n - 1) else begin FTeiler[Position] := Pot * Faktor; inc(Position); end; Pot := Pot * fFaktoren[n].Faktor; inc(i); until i > fFaktoren[n].Potenz; end;
begin
Position := 0; setlength(FTeiler, fAnzahlMoeglicherTeiler); FaktorAufbauen(1, fFakAnzahl); //Sortieren for i := Low(fTeiler) to fAnzahlMoeglicherTeiler - 2 do begin j := i; while (j >= Low(fTeiler)) and (fTeiler[j] > fTeiler[j + 1]) do begin Position := fTeiler[j]; fTeiler[j] := fTeiler[j + 1]; fTeiler[j + 1] := Position; dec(j); end; end;
end;
function tFaktorisieren.BasisPeriodeExtrahieren(var inZahl: Cardinal): Cardinal; var
i, cnt, Teiler: Cardinal;
begin
cnt := 0; result := 0; for i := Low(fBasFakt) to High(fBasFakt) do begin with fBasFakt[i] do begin if Faktor = 0 then BREAK; Teiler := Faktor; for cnt := 2 to Potenz do Teiler := Teiler * Faktor; end; cnt := 0; while (inZahl <> 0) and (inZahl mod Teiler = 0) do begin inZahl := inZahl div Teiler; inc(cnt); end; if cnt > result then result := cnt; end;
end;
procedure tFaktorisieren.PeriodeErmitteln(inZahl: Cardinal); var
i, TempZahl, TempPhi, TempPer, TempBasPer: Cardinal;
begin
Faktorisierung(inZahl); TempZahl := inZahl; //Die Basis_Nicht_Periode ermitteln TempBasPer := BasisPeriodeExtrahieren(TempZahl); TempPer := 0; if TempZahl > 1 then begin Faktorisierung(TempZahl); TempPhi := fEulerPhi; if (TempPhi > 1) then begin Faktorisierung(TempPhi); i := 0; repeat TempPer := fTeiler[i]; if BasExpMod(fFakBasis, TempPer, TempZahl) = 1 then Break; inc(i); until i >= Length(fTeiler); if i >= Length(fTeiler) then TempPer := inZahl - 1; end; end;
Faktorisierung(inZahl); fPeriodenlaenge := TempPer; fStartPeriode := TempBasPer;
end;
procedure tFaktorisieren.NachkommaPeriode(var OutText: string); var
i, limit: integer; Rest, Rest1, Divi, basis: Cardinal; pText: pChar;
procedure Ziffernfolge(Ende: longint); var j: longint; begin j := i - Ende;
while j < 0 do begin Rest := Rest * basis; Rest1 := Rest div Divi; Rest := Rest - Rest1 * Divi; //== Rest1 Mod Divi
pText^ := chr(Rest1 + Ord('0')); inc(pText);
inc(j); end;
i := Ende; end;
begin
limit := fStartPeriode + fPeriodenlaenge;
setlength(OutText, limit + 2 + 2 + 5); OutText[1] := '0'; OutText[2] := '.'; pText := @OutText[3];
Rest := 1; Divi := fFakZahl; basis := fFakBasis;
i := 0; Ziffernfolge(fStartPeriode); if fPeriodenlaenge = 0 then begin setlength(OutText, fStartPeriode + 2); EXIT; end;
pText^ := '_'; inc(pText); Ziffernfolge(limit); pText^ := '_'; inc(pText);
Ziffernfolge(limit + 5);
end;
type
tZahl = integer;
tRestFeld = array[0..31] of integer;
var
F: tFaktorisieren;
function tFaktorisieren.BasExpMod(b, e, m: Cardinal): Cardinal; begin
Result := 1; if m = 0 then exit; Result := 1; while (e > 0) do begin if (e and 1) <> 0 then Result := (Result * int64(b)) mod m; b := (int64(b) * b) mod m; e := e shr 1; end;
end;
procedure start; var
Limit, Testzahl: Cardinal; longPrimCount: int64; t1, t0: TDateTime;
begin
Limit := 500; Testzahl := 2; longPrimCount := 0; t0 := time;
repeat write(Limit: 8, ': '); repeat if F.Prim(Testzahl) then begin F.PeriodeErmitteln(Testzahl); if F.PeriodenLaenge = Testzahl - 1 then begin inc(longPrimCount); if Limit = 500 then write(Testzahl, ','); end end; inc(Testzahl); until Testzahl = Limit; inc(Limit, Limit); write(' .. count ', longPrimCount: 8, ' '); t1 := time; if (t1 - t0) > 1 / 864000 then write(FormatDateTime('HH:NN:SS.ZZZ', t1 - t0)); writeln; until Limit > 10 * 1000 * 1000;
t1 := time; writeln; writeln('count of long primes ', longPrimCount); writeln('Benoetigte Zeit ', FormatDateTime('HH:NN:SS.ZZZ', t1 - t0));
end;
begin
F := tFaktorisieren.create; writeln('Start'); start; writeln('Fertig.'); F.free; readln;
end.</lang>
- Output:
sh-4.4# ./Periode Start 500: 7,17,19,23,29,47,59,61,97,109,113,131,149,167,179,181,193,223,229,233,257,263,269,313,337,367,379,383,389,419,433,461,487,491,499, .. count 35 1000: .. count 60 2000: .. count 116 4000: .. count 218 8000: .. count 390 16000: .. count 716 32000: .. count 1300 64000: .. count 2430 128000: .. count 4498 256000: .. count 8434 00:00:00.100 512000: .. count 15920 00:00:00.220 1024000: .. count 30171 00:00:00.494 2048000: .. count 57115 00:00:01.140 4096000: .. count 108381 00:00:02.578 8192000: .. count 206594 00:00:06.073 count of long primes 206594 Benoetigte Zeit 00:00:06.073 Fertig.
Perl
<lang perl>use ntheory qw/divisors powmod is_prime/;
sub is_long_prime {
my($p) = @_; return 0 unless is_prime($p); for my $d (divisors($p-1)) { return $d+1 == $p if powmod(10, $d, $p) == 1; } 0;
}
print "Long primes ≤ 500:\n"; print join(' ', grep {is_long_prime($_) } 1 .. 500), "\n\n";
for my $n (500, 1000, 2000, 4000, 8000, 16000, 32000, 64000) {
printf "Number of long primes ≤ $n: %d\n", scalar grep { is_long_prime($_) } 1 .. $n;
}</lang>
- Output:
Long primes ≤ 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes ≤ 500: 35 Number of long primes ≤ 1000: 60 Number of long primes ≤ 2000: 116 Number of long primes ≤ 4000: 218 Number of long primes ≤ 8000: 390 Number of long primes ≤ 16000: 716 Number of long primes ≤ 32000: 1300 Number of long primes ≤ 64000: 2430
Using znorder
Faster due to going directly over primes and using znorder. Takes one second to count to 8,192,000. <lang perl>use ntheory qw/forprimes znorder/; my($t,$z)=(0,0); forprimes {
$z = znorder(10, $_); $t++ if defined $z && $z+1 == $_;
} 8192000; print "$t\n";</lang>
- Output:
206594
Phix
Slow version:
function is_long_prime(integer n) integer r = 1, rr, period = 0 for i=1 to n+1 do r = mod(10*r,n) end for rr = r while true do r = mod(10*r,n) period += 1 if period>=n then return false end if if r=rr then exit end if end while return period=n-1 end function
(use the same main() as below but limit maxN to 8 iterations)
Much faster version:
function is_long_prime(integer n) sequence f = factors(n-1,1) integer count = 0 for i=1 to length(f) do integer fi = f[i], e=1, base=10 while fi!=0 do if mod(fi,2)=1 then e = mod(e*base,n) end if base = mod(base*base,n) fi = floor(fi/2) end while if e=1 then count += 1 if count>1 then exit end if end if end for return count=1 end function procedure main() atom t0 = time() integer maxN = 500*power(2,14) --integer maxN = 500*power(2,7) -- (slow version) sequence long_primes = {} integer count = 0, n = 500, i = 2 while true do integer prime = get_prime(i) if is_long_prime(prime) then if prime<500 then long_primes &= prime end if if prime>n then if n=500 then printf(1,"The long primes up to 500 are:\n %V\n",{long_primes}) printf(1,"\nThe number of long primes up to:\n") end if printf(1," %7d is %d (%s)\n", {n, count, elapsed(time()-t0)}) if n=maxN then exit end if n *= 2 end if count += 1 end if i += 1 end while end procedure main()
- Output:
slow version:
The long primes up to 500 are: {7,17,19,23,29,47,59,61,97,109,113,131,149,167,179,181,193,223,229,233,257,263,269,313,337,367,379,383,389,419,433,461,487,491,499} The number of long primes up to: 500 is 35 (0.2s) 1000 is 60 (0.2s) 2000 is 116 (0.3s) 4000 is 218 (0.5s) 8000 is 390 (1.4s) 16000 is 716 (4.5s) 32000 is 1300 (16.0s) 64000 is 2430 (59.5s)
fast version:
The long primes up to 500 are: {7,17,19,23,29,47,59,61,97,109,113,131,149,167,179,181,193,223,229,233,257,263,269,313,337,367,379,383,389,419,433,461,487,491,499} The number of long primes up to: 500 is 35 (0.2s) 1000 is 60 (0.3s) 2000 is 116 (0.3s) 4000 is 218 (0.3s) 8000 is 390 (0.3s) 16000 is 716 (0.4s) 32000 is 1300 (0.5s) 64000 is 2430 (0.7s) 128000 is 4498 (1.4s) 256000 is 8434 (2.7s) 512000 is 15920 (5.8s) 1024000 is 30171 (12.3s) 2048000 is 57115 (26.3s) 4096000 is 108381 (56.5s) 8192000 is 206594 (1 minute and 60s)
Prolog
<lang prolog>% See https://en.wikipedia.org/wiki/Full_reptend_prime long_prime(Prime):-
is_prime(Prime), M is 10 mod Prime, M > 1, primitive_root(10, Prime).
% See https://en.wikipedia.org/wiki/Primitive_root_modulo_n#Finding_primitive_roots primitive_root(Base, Prime):-
Phi is Prime - 1, primitive_root(Phi, 2, Base, Prime).
primitive_root(1, _, _, _):-!. primitive_root(N, P, Base, Prime):-
is_prime(P), 0 is N mod P, !, X is (Prime - 1) // P, R is powm(Base, X, Prime), R \= 1, divide_out(N, P, M), Q is P + 1, primitive_root(M, Q, Base, Prime).
primitive_root(N, P, Base, Prime):-
Q is P + 1, Q * Q < Prime, !, primitive_root(N, Q, Base, Prime).
primitive_root(N, _, Base, Prime):-
X is (Prime - 1) // N, R is powm(Base, X, Prime), R \= 1.
divide_out(N, P, M):-
divmod(N, P, Q, 0), !, divide_out(Q, P, M).
divide_out(N, _, N).
print_long_primes([], _):-
!, nl.
print_long_primes([Prime|_], Limit):-
Prime > Limit, !, nl.
print_long_primes([Prime|Primes], Limit):-
writef('%w ', [Prime]), print_long_primes(Primes, Limit).
count_long_primes(_, L, Limit, _):-
L > Limit, !.
count_long_primes([], Limit, _, Count):-
writef('Number of long primes up to %w: %w\n', [Limit, Count]), !.
count_long_primes([Prime|Primes], L, Limit, Count):-
Prime > L, !, writef('Number of long primes up to %w: %w\n', [L, Count]), Count1 is Count + 1, L1 is L * 2, count_long_primes(Primes, L1, Limit, Count1).
count_long_primes([_|Primes], L, Limit, Count):-
Count1 is Count + 1, count_long_primes(Primes, L, Limit, Count1).
main(Limit):-
find_prime_numbers(Limit), findall(Prime, long_prime(Prime), Primes), writef('Long primes up to 500:\n'), print_long_primes(Primes, 500), count_long_primes(Primes, 500, Limit, 0).
main:-
main(256000).</lang>
Module for finding prime numbers up to some limit: <lang prolog>:- module(prime_numbers, [find_prime_numbers/1, is_prime/1]).
- - dynamic is_prime/1.
find_prime_numbers(N):-
retractall(is_prime(_)), assertz(is_prime(2)), init_sieve(N, 3), sieve(N, 3).
init_sieve(N, P):-
P > N, !.
init_sieve(N, P):-
assertz(is_prime(P)), Q is P + 2, init_sieve(N, Q).
sieve(N, P):-
P * P > N, !.
sieve(N, P):-
is_prime(P), !, S is P * P, cross_out(S, N, P), Q is P + 2, sieve(N, Q).
sieve(N, P):-
Q is P + 2, sieve(N, Q).
cross_out(S, N, _):-
S > N, !.
cross_out(S, N, P):-
retract(is_prime(S)), !, Q is S + 2 * P, cross_out(Q, N, P).
cross_out(S, N, P):-
Q is S + 2 * P, cross_out(Q, N, P).</lang>
- Output:
?- time(main(256000)). Long primes up to 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes up to 500: 35 Number of long primes up to 1000: 60 Number of long primes up to 2000: 116 Number of long primes up to 4000: 218 Number of long primes up to 8000: 390 Number of long primes up to 16000: 716 Number of long primes up to 32000: 1300 Number of long primes up to 64000: 2430 Number of long primes up to 128000: 4498 Number of long primes up to 256000: 8434 % 8,564,024 inferences, 0.991 CPU in 1.040 seconds (95% CPU, 8641390 Lips) true.
PureBasic
<lang PureBasic>#MAX=64000 If OpenConsole()=0 : End 1 : EndIf
Dim p.b(#MAX) : FillMemory(@p(),#MAX,#True,#PB_Byte) For n=2 To Int(Sqr(#MAX))+1 : If p(n) : m=n*n : While m<=#MAX : p(m)=#False : m+n : Wend : EndIf : Next
Procedure.i periodic(v.i)
r=1 : Repeat : r=(r*10)%v : c+1 : If r<=1 : ProcedureReturn c : EndIf : ForEver
EndProcedure
n=500 PrintN(LSet("_",15,"_")+"Long primes upto "+Str(n)+LSet("_",15,"_")) For i=3 To 500 Step 2
If p(i) And (i-1)=periodic(i) Print(RSet(Str(i),5)) : c+1 : If c%10=0 : PrintN("") : EndIf EndIf
Next
PrintN(~"\n") PrintN("The number of long primes up to:") PrintN(RSet(Str(n),8)+" is "+Str(c)) : n+n For i=501 To #MAX+1 Step 2
If p(i) And (i-1)=periodic(i) : c+1 : EndIf If i>n : PrintN(RSet(Str(n),8)+" is "+Str(c)) : n+n : EndIf
Next Input()</lang>
- Output:
_______________Long primes upto 500_______________ 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
Python
<lang python>def sieve(limit):
primes = [] c = [False] * (limit + 1) # composite = true # no need to process even numbers p = 3 while True: p2 = p * p if p2 > limit: break for i in range(p2, limit, 2 * p): c[i] = True while True: p += 2 if not c[p]: break
for i in range(3, limit, 2): if not c[i]: primes.append(i) return primes
- finds the period of the reciprocal of n
def findPeriod(n):
r = 1 for i in range(1, n): r = (10 * r) % n rr = r period = 0 while True: r = (10 * r) % n period += 1 if r == rr: break return period
primes = sieve(64000) longPrimes = [] for prime in primes:
if findPeriod(prime) == prime - 1: longPrimes.append(prime)
numbers = [500, 1000, 2000, 4000, 8000, 16000, 32000, 64000] count = 0 index = 0 totals = [0] * len(numbers) for longPrime in longPrimes:
if longPrime > numbers[index]: totals[index] = count index += 1 count += 1
totals[-1] = count print('The long primes up to 500 are:') print(str(longPrimes[:totals[0]]).replace(',', )) print('\nThe number of long primes up to:') for (i, total) in enumerate(totals):
print(' %5d is %d' % (numbers[i], total))</lang>
- Output:
The long primes up to 500 are: [7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499] The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
Quackery
eratosthenes
and isprime
are defined at Sieve of Eratosthenes#Quackery.
bsearchwith
is defined at Binary search#Quackery.
<lang Quackery> [ over size 0 swap 2swap
bsearchwith < drop ] is search ( [ --> n )
[ 1 over 1 - times [ 10 * over mod ] tuck 0 temp put [ 10 * over mod 1 temp tally rot 2dup != while unrot again ] 2drop drop temp take ] is period ( n --> n )
[ dup isprime not iff [ drop false ] done dup period 1+ = ] is islongprime ( n --> b ) 64000 eratosthenes
[] 64000 times [ i^ islongprime if [ i^ join ] ] behead drop dup dup 500 search split drop echo cr cr ' [ 500 1000 2000 4000 8000 16000 32000 64000 ] witheach [ dup echo say " --> " dip dup search echo cr ] drop</lang>
- Output:
[ 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 ] 500 --> 35 1000 --> 60 2000 --> 116 4000 --> 218 8000 --> 390 16000 --> 716 32000 --> 1300 64000 --> 2430
Racket
(at least find-period)
<lang racket>#lang racket (require math/number-theory)
(define (find-period n)
(let ((rr (for/fold ((r 1)) ((i (in-range 1 (+ n 2)))) (modulo (* 10 r) n)))) (let period-loop ((r rr) (p 1)) (let ((r′ (modulo (* 10 r) n))) (if (= r′ rr) p (period-loop r′ (add1 p)))))))
(define (long-prime? n)
(and (prime? n) (= (find-period n) (sub1 n))))
(define memoised-long-prime? (let ((h# (make-hash))) (λ (n) (hash-ref! h# n (λ () (long-prime? n))))))
(module+ main
;; strictly, won't test 500 itself... but does it look prime to you? (filter memoised-long-prime? (range 7 500 2)) (for-each (λ (n) (displayln (cons n (for/sum ((i (in-range 7 n 2))) (if (memoised-long-prime? i) 1 0))))) '(500 1000 2000 4000 8000 16000 32000 64000)))
(module+ test
(require rackunit) (check-equal? (map find-period '(7 11 977)) '(6 2 976)))</lang>
- Output:
'(7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499) (500 . 35) (1000 . 60) (2000 . 116) (4000 . 218) (8000 . 390) (16000 . 716) (32000 . 1300) (64000 . 2430)
Raku
(formerly Perl 6)
Not very fast as the numbers get larger. <lang perl6>use Math::Primesieve; my $sieve = Math::Primesieve.new;
sub is-long (Int $p) {
my $r = 1; my $rr = $r = (10 * $r) % $p for ^$p; my $period; loop { $r = (10 * $r) % $p; ++$period; last if $period >= $p or $r == $rr; } $period == $p - 1 and $p > 2;
}
my @primes = $sieve.primes(500); my @long-primes = @primes.grep: {.&is-long};
put "Long primes ≤ 500:\n", @long-primes;
@long-primes = ();
for 500, 1000, 2000, 4000, 8000, 16000, 32000, 64000 -> $upto {
state $from = 0; my @extend = $sieve.primes($from, $upto); @long-primes.append: @extend.hyper(:8degree).grep: {.&is-long}; say "\nNumber of long primes ≤ $upto: ", +@long-primes; $from = $upto;
}</lang>
- Output:
Long primes ≤ 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes ≤ 500: 35 Number of long primes ≤ 1000: 60 Number of long primes ≤ 2000: 116 Number of long primes ≤ 4000: 218 Number of long primes ≤ 8000: 390 Number of long primes ≤ 16000: 716 Number of long primes ≤ 32000: 1300 Number of long primes ≤ 64000: 2430
REXX
For every doubling of the limit, it takes about roughly 5 times longer to compute the long primes.
uses odd numbers
<lang rexx>/*REXX pgm calculates/displays base ten long primes (AKA golden primes, proper primes,*/ /*───────────────────── maximal period primes, long period primes, full reptend primes).*/ parse arg a /*obtain optional argument from the CL.*/ if a= | a="," then a= '500 -500 -1000 -2000 -4000 -8000 -16000' , /*Not specified? */
'-32000 -64000 -128000 -512000 -1024000' /*Then use default*/ do k=1 for words(a); H=word(a, k) /*step through the list of high limits.*/ neg= H<1 /*used as an indicator to display count*/ H= abs(H) /*obtain the absolute value of H. */ $= /*the list of long primes (so far). */ do j=7 to H by 2 /*start with 7, just use odd integers.*/ if .len(j) + 1 \== j then iterate /*Period length wrong? Then skip it. */ $=$ j /*add the long prime to the $ list.*/ end /*j*/ say if neg then do; say 'number of long primes ≤ ' H " is: " words($); end else do; say 'list of long primes ≤ ' H":"; say strip($); end end /*k*/
exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ .len: procedure; parse arg x; r=1; do x; r= 10*r // x; end /*x*/
rr=r; do p=1 until r==rr; r= 10*r // x; end /*p*/ return p</lang>
- output when using the internal default inputs:
list of long primes ≤ 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 number of long primes ≤ 500 is: 35 number of long primes ≤ 1000 is: 60 number of long primes ≤ 2000 is: 116 number of long primes ≤ 4000 is: 218 number of long primes ≤ 8000 is: 390 number of long primes ≤ 16000 is: 716 number of long primes ≤ 32000 is: 1300 number of long primes ≤ 64000 is: 2430 number of long primes ≤ 128000 is: 4498 number of long primes ≤ 512000 is: 15920 number of long primes ≤ 1024000 is: 30171
uses odd numbers, some prime tests
This REXX version is about 2 times faster than the 1st REXX version (because it does some primality testing). <lang rexx>/*REXX pgm calculates/displays base ten long primes (AKA golden primes, proper primes,*/ /*───────────────────── maximal period primes, long period primes, full reptend primes).*/ parse arg a /*obtain optional argument from the CL.*/ if a= | a="," then a= '500 -500 -1000 -2000 -4000 -8000 -16000' , /*Not specified? */
'-32000 -64000 -128000 -512000 -1024000' /*Then use default*/ do k=1 for words(a); H=word(a, k) /*step through the list of high limits.*/ neg= H<1 /*used as an indicator to display count*/ H= abs(H) /*obtain the absolute value of H. */ $= /*the list of long primes (so far). */ do j=7 to H by 2; parse var j -1 _ /*start with 7, just use odd integers.*/ if _==5 then iterate /*last digit a five? Then not a prime.*/ if j// 3==0 then iterate /*Is divisible by 3? " " " " */ if j\==11 then if j//11==0 then iterate /* " " " 11? " " " " */ if j\==13 then if j//13==0 then iterate /* " " " 13? " " " " */ if j\==17 then if j//17==0 then iterate /* " " " 17? " " " " */ if j\==19 then if j//19==0 then iterate /* " " " 19? " " " " */ if .len(j) + 1 \== j then iterate /*Period length wrong? Then skip it. */ $=$ j /*add the long prime to the $ list.*/ end /*j*/ say if neg then do; say 'number of long primes ≤ ' H " is: " words($); end else do; say 'list of long primes ≤ ' H":"; say strip($); end end /*k*/
exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ .len: procedure; parse arg x; r=1; do x; r= 10*r // x; end /*x*/
rr=r; do p=1 until r==rr; r= 10*r // x; end /*p*/ return p</lang>
- output is identical to the 1st REXX version.
uses primes
This REXX version is about 5 times faster than the 1st REXX version (because it only tests primes). <lang rexx>/*REXX pgm calculates/displays base ten long primes (AKA golden primes, proper primes,*/ /*───────────────────── maximal period primes, long period primes, full reptend primes).*/ parse arg a /*obtain optional argument from the CL.*/ if a= | a="," then a= '500 -500 -1000 -2000 -4000 -8000 -16000' , /*Not specified? */
'-32000 -64000 -128000 -512000 -1024000' /*Then use default*/
m=0; aa= words(a) /* [↑] two list types of low primes. */
do j=1 for aa; m= max(m, abs(word(a, j))) /*find the maximum argument in the list*/ end /*j*/
call genP /*go and generate some primes. */
do k=1 for aa; H= word(a, k) /*step through the list of high limits.*/ neg= H<1 /*used as an indicator to display count*/ H= abs(H) /*obtain the absolute value of H. */ $= /*the list of long primes (so far). */ do j=7 to H by 2 if \@.j then iterate /*Is J not a prime? Then skip it. */ if .len(j) + 1 \== j then iterate /*Period length wrong? " " " */ $= $ j /*add the long prime to the $ list.*/ end /*j*/ /* [↑] some pretty weak prime testing.*/ say if neg then say 'number of long primes ≤ ' H " is: " words($) else do; say 'list of long primes ≤ ' H":"; say strip($); end end /*k*/
exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ genP: @.=0; @.2=1; @.3=1; @.5=1; @.7=1; @.11=1; !.=0; !.1=2; !.2=3; !.3=5; !.4=7; !.5=11
#= 5 /*the number of primes (so far). */ do g=!.#+2 by 2 until g>=m /*gen enough primes to satisfy max A. */ if @.g\==0 then iterate /*Is it not a prime? Then skip it. */ do d=2 until !.d**2>g /*only divide up to square root of X. */ if g//!.d==0 then iterate g /*Divisible? Then skip this integer. */ end /*d*/ /* [↓] a spanking new prime was found.*/ #= #+1 @.g= 1; !.#= g /*bump P counter; assign P, add to P's.*/ end /*g*/ return
/*──────────────────────────────────────────────────────────────────────────────────────*/ .len: procedure; parse arg x; r=1; do x; r= 10*r // x; end /*x*/
rr=r; do p=1 until r==rr; r= 10*r // x; end /*p*/ return p</lang>
- output is identical to the 1st REXX version.
Ruby
System: I7-6700HQ, 3.5 GHz, Linux Kernel 5.6.17 Run as: $ ruby longprimes.rb
Finding long prime numbers using finding period location (translation of Python's module def findPeriod(n)) <lang Ruby>require 'prime'
batas = 64_000 # limit number start = Time.now # time of starting lp_array = [] # array of long-prime numbers
def find_period(n)
r, period = 1, 0 (1...n).each {r = (10 * r) % n} rr = r loop do r = (10 * r) % n period += 1 break if r == rr end return period
end
Prime.each(batas).each do |prime|
lp_array.push(prime) if find_period(prime) == prime-1 && prime != 2
end
[500, 1000, 2000, 4000, 8000, 16000, 32000, 64000].each do |s|
if s == 500 puts "\nAll long primes up to #{s} are: #{lp_array.count {|x| x < s}}. They are:" lp_array.each {|x| print x, " " if x < s} else print "\nAll long primes up to #{s} are: #{lp_array.count {|x| x < s}}" end
end
puts "\n\nTime: #{Time.now - start}"</lang>
- Output:
All long primes up to 500 are: 35. They are: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 All long primes up to 1000 are: 60 All long primes up to 2000 are: 116 All long primes up to 4000 are: 218 All long primes up to 8000 are: 390 All long primes up to 16000 are: 716 All long primes up to 32000 are: 1300 All long primes up to 64000 are: 2430
Time: 16.212039 # Ruby 2.7.1 Time: 18.664795 # JRuby 9.2.11.1 Time: 3.198 # Truffleruby 20.1.0
Alternatively, by using primitive way: converting value into string and make assessment for proper repetend position. Indeed produce same information, but with longer time. <lang Ruby>require 'prime' require 'bigdecimal' require 'strscan'
batas = 64_000 # limit number start = Time.now # time of starting lp_array = [] # array of long-prime numbers a = BigDecimal.("1") # number being divided, that is 1.
Prime.each(batas).each do |prime|
cek = a.div(prime, (prime-1)*2).truncate((prime-1)*2).to_s('F')[2..-1] # Dividing 1 with prime and take its value as string. if (cek[0, prime-1] == cek[prime-1, prime-1]) i = prime-2 until i < 5 break if cek[0, i] == cek[i, i] i-=1 cek.slice!(-2, 2) # Shortening checked string to reduce checking process load end until i == 0 break if cek[0, (cek.size/i)*i].scan(/.{#{i}}/).uniq.length == 1 i-=1 end
lp_array.push(prime) if i == 0 end
end
[500, 1000, 2000, 4000, 8000, 16000, 32000, 64000].each do |s|
if s == 500 puts "\nAll long primes up to #{s} are: #{lp_array.count {|x| x < s}}. They are:" lp_array.each {|x| print x, " " if x < s} else print "\nAll long primes up to #{s} are: #{lp_array.count {|x| x < s}}" end
end
puts "\n\nTime: #{Time.now - start}"</lang>
- Output:
(same output with previous version, but longer time elapse) Time: 629.360075011 secs # Ruby 2.7.1
Fastest version. <lang ruby>def prime?(n) # P3 Prime Generator primality test
return n | 1 == 3 if n < 5 # n: 2,3|true; 0,1,4|false return false if n.gcd(6) != 1 # this filters out 2/3 of all integers pc, sqrtn = 5, Integer.sqrt(n) # first P3 prime candidates sequence value until pc > sqrtn return false if n % pc == 0 || n % (pc + 2) == 0 # if n is composite pc += 6 # 1st prime candidate for next residues group end true
end
def divisors(n) # divisors of n -> [1,..,n]
f = [] (1..Integer.sqrt(n)).each { |i| (n % i).zero? && (f << i; f << n / i if n / i != i) } f.sort
end
- The smallest divisor d of p-1 such that 10^d = 1 (mod p),
- is the length of the period of the decimal expansion of 1/p.
def long_prime?(p)
return false unless prime? p divisors(p - 1).each { |d| return d == (p - 1) if 10.pow(d, p) == 1 } false
end
start = Time.now puts "Long primes ≤ 500:" (7..500).each { |pc| print "#{pc} " if long_prime? pc } puts [500, 1000, 2000, 4000, 8000, 16000, 32000, 64000].each do |n|
puts "Number of long primes ≤ #{n}: #{(7..n).count { |pc| long_prime? pc }}"
end puts "\nTime: #{(Time.now - start)} secs"</lang>
- Output:
Long primes ≤ 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes ≤ 500: 35 Number of long primes ≤ 1000: 60 Number of long primes ≤ 2000: 116 Number of long primes ≤ 4000: 218 Number of long primes ≤ 8000: 390 Number of long primes ≤ 16000: 716 Number of long primes ≤ 32000: 1300 Number of long primes ≤ 64000: 2430
Time: 0.228912772 secs # Ruby 2.7.1 Time: 0.544648 secs # JRuby 9.2.11.1 Time: 11.985 secs # Truffleruby 20.1.0
Rust
<lang rust>// main.rs // References: // https://en.wikipedia.org/wiki/Full_reptend_prime // https://en.wikipedia.org/wiki/Primitive_root_modulo_n#Finding_primitive_roots
mod bit_array; mod prime_sieve;
use prime_sieve::PrimeSieve;
fn modpow(mut base: usize, mut exp: usize, n: usize) -> usize {
if n == 1 { return 0; } let mut result = 1; base %= n; while exp > 0 { if (exp & 1) == 1 { result = (result * base) % n; } base = (base * base) % n; exp >>= 1; } result
}
fn is_long_prime(sieve: &PrimeSieve, prime: usize) -> bool {
if !sieve.is_prime(prime) { return false; } if 10 % prime == 0 { return false; } let n = prime - 1; let mut m = n; let mut p = 2; while p * p <= n { if sieve.is_prime(p) && m % p == 0 { if modpow(10, n / p, prime) == 1 { return false; } while m % p == 0 { m /= p; } } p += 1; } if m == 1 { return true; } modpow(10, n / m, prime) != 1
}
fn long_primes(limit1: usize, limit2: usize) {
let sieve = PrimeSieve::new(limit2); let mut count = 0; let mut limit = limit1; let mut prime = 3; while prime < limit2 { if is_long_prime(&sieve, prime) { if prime < limit1 { print!("{} ", prime); } if prime > limit { print!("\nNumber of long primes up to {}: {}", limit, count); limit *= 2; } count += 1; } prime += 2; } println!("\nNumber of long primes up to {}: {}", limit, count);
}
fn main() {
long_primes(500, 8192000);
}</lang>
<lang rust>// prime_sieve.rs use crate::bit_array;
pub struct PrimeSieve {
composite: bit_array::BitArray,
}
impl PrimeSieve {
pub fn new(limit: usize) -> PrimeSieve { let mut sieve = PrimeSieve { composite: bit_array::BitArray::new(limit / 2), }; let mut p = 3; while p * p <= limit { if !sieve.composite.get(p / 2 - 1) { let inc = p * 2; let mut q = p * p; while q <= limit { sieve.composite.set(q / 2 - 1, true); q += inc; } } p += 2; } sieve } pub fn is_prime(&self, n: usize) -> bool { if n < 2 { return false; } if n % 2 == 0 { return n == 2; } !self.composite.get(n / 2 - 1) }
}</lang>
<lang rust>// bit_array.rs pub struct BitArray {
array: Vec<u32>,
}
impl BitArray {
pub fn new(size: usize) -> BitArray { BitArray { array: vec![0; (size + 31) / 32], } } pub fn get(&self, index: usize) -> bool { let bit = 1 << (index & 31); (self.array[index >> 5] & bit) != 0 } pub fn set(&mut self, index: usize, new_val: bool) { let bit = 1 << (index & 31); if new_val { self.array[index >> 5] |= bit; } else { self.array[index >> 5] &= !bit; } }
}</lang>
- Output:
Execution time is just over 1.5 seconds on my system (macOS 10.15, 3.2 GHz Quad-Core Intel Core i5)
7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes up to 500: 35 Number of long primes up to 1000: 60 Number of long primes up to 2000: 116 Number of long primes up to 4000: 218 Number of long primes up to 8000: 390 Number of long primes up to 16000: 716 Number of long primes up to 32000: 1300 Number of long primes up to 64000: 2430 Number of long primes up to 128000: 4498 Number of long primes up to 256000: 8434 Number of long primes up to 512000: 15920 Number of long primes up to 1024000: 30171 Number of long primes up to 2048000: 57115 Number of long primes up to 4096000: 108381 Number of long primes up to 8192000: 206594
Rust FP
<lang rust>fn is_oddprime(n: u64) -> bool {
let limit = (n as f64).sqrt().ceil() as u64; (3..=limit).step_by(2).all(|a| n % a > 0)
}
fn divisors(n: u64) -> Vec<u64> {
let list1: Vec<u64> = (1..=(n as f64).sqrt().floor() as u64) .filter(|d| n % d == 0).collect(); let list2: Vec<u64> = list1.iter().rev() .skip_while(|&d| d * d == n).map(|d| n / d).collect(); [list1, list2].concat()
}
fn power_mod(base: u64, exp: u64, modulo: u64) -> u64 {
fn iter(base: u64, modu: &u64, exp: u64, res: u64) -> u64 { if exp > 0 { let base1 = (base * base) % modu; let res1 = if exp & 1 > 0 {(base * res) % modu} else {res}; iter(base1, modu, exp >> 1, res1) } else {res} } iter(base, &modulo,exp - 1,base % modulo)
}
// the smallest divisor d of p-1 such that 10^d = 1 (mod p) // is the length of the period of the decimal expansion of 1/p fn is_longprime(p: u64) -> bool {
match divisors(p - 1).into_iter() .skip_while(|&d| power_mod(10, d, p) != 1) .next() { Some(d) => d + 1 == p, None => false }
}
fn long_primes() -> impl Iterator<Item = u64> {
(7..).step_by(2).filter(|&p|is_oddprime(p)) .filter(|&p| is_longprime(p))
}
fn main() {
println!("long primes up to 500:"); let list500: Vec<u64> = long_primes() .take_while(|&p| p <= 500) .collect(); println!("{:?}\n", list500); let limits: Vec<u64> = (0..8).map(|n| 2u64.pow(n) * 500).collect(); for limit in limits { let start = std::time::Instant::now(); let count = long_primes().take_while(|&p| p <= limit).count(); let duration = start.elapsed().as_millis(); println!("there are {:4} long primes up to {:5} [time(ms) {:3}]", count, limit, duration); }
}</lang>
Scala
<lang scala>object LongPrimes extends App {
def primeStream = LazyList.from(3, 2) .filter(p => (3 to math.sqrt(p).ceil.toInt by 2).forall(p % _ > 0)) def longPeriod(p: Int): Boolean = { val mstart = 10 % p @annotation.tailrec def iter(mod: Int, period: Int): Int = { val mod1 = (10 * mod) % p if (mod1 == mstart) period else iter(mod1, period + 1) } iter(mstart, 1) == p - 1 }
val longPrimes = primeStream.filter(longPeriod(_)) println("long primes up to 500:") println(longPrimes.takeWhile(_ <= 500).mkString(" ")) println val limitList = Seq.tabulate(8)(math.pow(2, _).toInt * 500) for (limit <- limitList) { val count = longPrimes.takeWhile(_ <= limit).length println(f"there are $count%4d long primes up to $limit%5d") }
}</lang>
- Output:
long primes up to 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 there are 35 long primes up to 500 there are 60 long primes up to 1000 there are 116 long primes up to 2000 there are 218 long primes up to 4000 there are 390 long primes up to 8000 there are 716 long primes up to 16000 there are 1300 long primes up to 32000 there are 2430 long primes up to 64000
Sidef
The smallest divisor d of p-1 such that 10^d = 1 (mod p), is the length of the period of the decimal expansion of 1/p. <lang ruby>func is_long_prime(p) {
for d in (divisors(p-1)) { if (powmod(10, d, p) == 1) { return (d+1 == p) } }
return false
}
say "Long primes ≤ 500:" say primes(500).grep(is_long_prime).join(' ')
for n in ([500, 1000, 2000, 4000, 8000, 16000, 32000, 64000]) {
say ("Number of long primes ≤ #{n}: ", primes(n).count_by(is_long_prime))
}</lang>
- Output:
Long primes ≤ 500: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 Number of long primes ≤ 500: 35 Number of long primes ≤ 1000: 60 Number of long primes ≤ 2000: 116 Number of long primes ≤ 4000: 218 Number of long primes ≤ 8000: 390 Number of long primes ≤ 16000: 716 Number of long primes ≤ 32000: 1300 Number of long primes ≤ 64000: 2430
Alternatively, we can implement the is_long_prime(p) function using the `znorder(a, p)` built-in method, which is considerably faster:
<lang ruby>func is_long_prime(p) {
znorder(10, p) == p-1
}</lang>
Swift
<lang swift>public struct Eratosthenes: Sequence, IteratorProtocol {
private let n: Int private let limit: Int
private var i = 2 private var sieve: [Int]
public init(upTo: Int) { if upTo <= 1 { self.n = 0 self.limit = -1 self.sieve = [] } else { self.n = upTo self.limit = Int(Double(n).squareRoot()) self.sieve = Array(0...n) } }
public mutating func next() -> Int? { while i < n { defer { i += 1 }
if sieve[i] != 0 { if i <= limit { for notPrime in stride(from: i * i, through: n, by: i) { sieve[notPrime] = 0 } }
return i } }
return nil }
}
func findPeriod(n: Int) -> Int {
let r = (1...n+1).reduce(1, {res, _ in (10 * res) % n }) var rr = r var period = 0
repeat { rr = (10 * rr) % n period += 1 } while r != rr
return period
}
let longPrimes = Eratosthenes(upTo: 64000).dropFirst().lazy.filter({ findPeriod(n: $0) == $0 - 1 })
print("Long primes less than 500: \(Array(longPrimes.prefix(while: { $0 <= 500 })))")
let counts =
longPrimes.reduce(into: [500: 0, 1000: 0, 2000: 0, 4000: 0, 8000: 0, 16000: 0, 32000: 0, 64000: 0], {counts, n in for key in counts.keys where n < key { counts[key]! += 1 } })
for key in counts.keys.sorted() {
print("There are \(counts[key]!) long primes less than \(key)")
}</lang>
- Output:
Long primes less than 500: [7, 17, 19, 23, 29, 47, 59, 61, 97, 109, 113, 131, 149, 167, 179, 181, 193, 223, 229, 233, 257, 263, 269, 313, 337, 367, 379, 383, 389, 419, 433, 461, 487, 491, 499] There are 35 long primes less than 500 There are 60 long primes less than 1000 There are 116 long primes less than 2000 There are 218 long primes less than 4000 There are 390 long primes less than 8000 There are 716 long primes less than 16000 There are 1300 long primes less than 32000 There are 2430 long primes less than 64000
Visual Basic .NET
<lang vbnet>Imports System, System.Collections.Generic, System.Linq, System.Console
Module LongPrimes
Function Period(ByVal n As Integer) As Integer Dim m As Integer, r As Integer = 1 For i As Integer = 0 To n : r = 10 * r Mod n : Next m = r : Period = 1 : While True r = (10 * r) Mod n : If r = m Then Return Period Period += 1 : End While End Function
Sub Main() Dim primes As IEnumerable(Of Integer) = SomePrimeGenerator.Primes(64000).Skip(1).Where(Function(p) Period(p) = p - 1).Append(99999) Dim count As Integer = 0, limit As Integer = 500 WriteLine(String.Join(" ", primes.TakeWhile(Function(p) p <= limit))) For Each prime As Integer In primes If prime > limit Then WriteLine($"There are {count} long primes below {limit}") limit <<= 1 : End If : count += 1 : Next End Sub
End Module
Module SomePrimeGenerator
Iterator Function Primes(lim As Integer) As IEnumerable(Of Integer) Dim flags As Boolean() = New Boolean(lim) {}, j As Integer = 2, d As Integer = 3, sq As Integer = 4 While sq <= lim If Not flags(j) Then Yield j : For k As Integer = sq To lim step j flags(k) = True : Next End If : j += 1 : d += 2 : sq += d End While : While j <= lim If Not flags(j) Then Yield j j += 1 : End While End Function
End Module</lang>
- Output:
Same output as C#.
Wren
<lang ecmascript>import "/fmt" for Fmt import "/math" for Int
// finds the period of the reciprocal of n var findPeriod = Fn.new { |n|
var r = 1 for (i in 1..n+1) r = (10*r) % n var rr = r var period = 0 var ok = true while (ok) { r = (10*r) % n period = period + 1 ok = (r != rr) } return period
}
var primes = Int.primeSieve(64000).skip(1) var longPrimes = [] for (prime in primes) {
if (findPeriod.call(prime) == prime - 1) longPrimes.add(prime)
} var numbers = [500, 1000, 2000, 4000, 8000, 16000, 32000, 64000] var index = 0 var count = 0 var totals = List.filled(numbers.count, 0) for (longPrime in longPrimes) {
if (longPrime > numbers[index]) { totals[index] = count index = index + 1 } count = count + 1
} totals[-1] = count System.print("The long primes up to %(numbers[0]) are: ") System.print(longPrimes[0...totals[0]].join(" "))
System.print("\nThe number of long primes up to: ") var i = 0 for (total in totals) {
System.print(" %(Fmt.d(5, numbers[i])) is %(total)") i = i + 1
}</lang>
- Output:
The long primes up to 500 are: 7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499 The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
XBasic
<lang xbasic> PROGRAM "longprimes" VERSION "0.0002"
DECLARE FUNCTION Entry() INTERNAL FUNCTION Sieve(limit&, primes&[], count%) INTERNAL FUNCTION FindPeriod(n&)
FUNCTION Entry()
DIM numbers&[7] numbers&[0] = 500 numbers&[1] = 1000 numbers&[2] = 2000 numbers&[3] = 4000 numbers&[4] = 8000 numbers&[5] = 16000 numbers&[6] = 32000 numbers&[7] = 64000 numberUpperBound% = UBOUND(numbers&[]) DIM totals%[numberUpperBound%] DIM primes&[6499] PRINT "Please wait." PRINT Sieve(64000, @primes&[], @primeCount%) DIM longPrimes&[primeCount% - 1] ' Surely longCount% < primeCount% longCount% = 0 FOR i% = 0 TO primeCount% - 1 prime& = primes&[i%] IF FindPeriod(prime&) = prime& - 1 THEN longPrimes&[longCount%] = prime& INC longCount% END IF NEXT i% count% = 0 index% = 0 FOR i% = 0 TO longCount% - 1 IF longPrimes&[i%] > numbers&[index%] THEN totals%[index%] = count% INC index% END IF INC count% NEXT i% totals%[numberUpperBound%] = count% PRINT "The long primes up to"; numbers&[0]; " are:" PRINT "["; FOR i% = 0 TO totals%[0] - 2 PRINT STRING$(longPrimes&[i%]); " "; NEXT i% IF totals%[0] > 0 THEN PRINT STRING$(longPrimes&[totals%[0] - 1]); END IF PRINT "]" PRINT PRINT "The number of long primes up to:" FOR i% = 0 TO numberUpperBound% PRINT FORMAT$(" #####", numbers&[i%]); " is"; totals%[i%] NEXT i%
END FUNCTION
FUNCTION Sieve(limit&, primes&[], count%)
DIM c@[limit&] FOR i& = 0 TO limit& c@[i&] = 0 NEXT i& ' No need to process even numbers p% = 3 n% = 0 p2& = p% * p% DO WHILE p2& <= limit& FOR i& = p2& TO limit& STEP 2 * p% c@[i&] = 1 NEXT i& DO p% = p% + 2 LOOP UNTIL !c@[p%] p2& = p% * p% LOOP FOR i& = 3 TO limit& STEP 2 IFZ c@[i&] THEN primes&[n%] = i& INC n% END IF NEXT i& count% = n%
END FUNCTION
' Finds the period of the reciprocal of n& FUNCTION FindPeriod(n&)
r& = 1 period& = 0 FOR i& = 1 TO n& + 1 r& = (10 * r&) MOD n& NEXT i& rr& = r& DO r& = (10 * r&) MOD n& INC period& LOOP UNTIL r& = rr&
END FUNCTION period&
END PROGRAM </lang>
- Output:
Please wait. The long primes up to 500 are: [7 17 19 23 29 47 59 61 97 109 113 131 149 167 179 181 193 223 229 233 257 263 269 313 337 367 379 383 389 419 433 461 487 491 499] The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
zkl
Using GMP (GNU Multiple Precision Arithmetic Library, probabilistic primes), because it is easy and fast to generate primes. <lang zkl>var [const] BN=Import("zklBigNum"); // libGMP primes,p := List.createLong(7_000), BN(3); // one big alloc vs lots of allocs while(p.nextPrime()<=64_000){ primes.append(p.toInt()) } // 6412 of them, skipped 2 primes.append(p.toInt()); // and one more so tail prime is >64_000
longPrimes:=primes.filter(fcn(p){ findPeriod(p)==p-1 }); // yawn fcn findPeriod(n){
r,period := 1,0; do(n){ r=(10*r)%n } rr:=r; while(True){ // reduce is more concise but 2.5 times slower r=(10*r)%n; period+=1; if(r==rr) break; } period
}</lang>
<lang zkl>fiveHundred:=longPrimes.filter('<(500)); println("The long primes up to 500 are:\n",longPrimes.filter('<(500)).concat(","));
println("\nThe number of long primes up to:"); foreach n in (T(500, 1000, 2000, 4000, 8000, 16000, 32000, 64000)){
println(" %5d is %d".fmt( n, longPrimes.filter1n('>(n)) ));
}</lang>
- Output:
The long primes up to 500 are: 7,17,19,23,29,47,59,61,97,109,113,131,149,167,179,181,193,223,229,233,257,263,269,313,337,367,379,383,389,419,433,461,487,491,499 The number of long primes up to: 500 is 35 1000 is 60 2000 is 116 4000 is 218 8000 is 390 16000 is 716 32000 is 1300 64000 is 2430
- Programming Tasks
- Prime Numbers
- 11l
- AppleScript
- C
- C sharp
- C++
- Common Lisp
- Crystal
- Delphi
- F Sharp
- Factor
- Forth
- FreeBASIC
- Go
- Haskell
- J
- Java
- Jq
- Julia
- Kotlin
- M2000 Interpreter
- Maple
- Mathematica
- Wolfram Language
- NewLISP
- Nim
- Pascal
- Perl
- Ntheory
- Phix
- Prolog
- PureBasic
- Python
- Quackery
- Racket
- Raku
- REXX
- Ruby
- Rust
- Scala
- Sidef
- Swift
- Visual Basic .NET
- Wren
- Wren-fmt
- Wren-math
- XBasic
- Zkl