Source code in the prolog programming language

% 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).


:- 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).


isPrime(A):-
  	A1 is ceil(sqrt(A)),
  	between(2, A1, N),
  	0 =:= A mod N,!,
  	false.
isPrime(_).

divisors(N, Dlist):-
	N1 is floor(sqrt(N)),
	numlist(1, N1, Ds0),
	include([D]>>(N mod D =:= 0), Ds0, Ds1),
	reverse(Ds1, [Dh|Dt]),
	( Dh * Dh < N
	-> Ds1a = [Dh|Dt]
	;  Ds1a = Dt
	),
	maplist([X,Y]>>(Y is N div X), Ds1a, Ds2),
	append(Ds1, Ds2, Dlist).

longPrime(P):-
	divisors(P - 1, Dlist),
	longPrime(P, Dlist).

longPrime(_,[]):- false.
longPrime(P, [D|Dtail]):-
	powm(10, D, P) =\= 1,!,
	longPrime(P, Dtail).
longPrime(P, [D|_]):-!,
	D =:= P - 1.

isLongPrime(N):-
	isPrime(N),
	longPrime(N).

longPrimes(N, LongPrimes):-
	numlist(7, N, List),
	include(isLongPrime, List, LongPrimes).

run([]):-!.
run([Limit|Tail]):-
	statistics(runtime,[Start|_]),
	longPrimes(Limit, LongPrimes),
	length(LongPrimes, Num),
	statistics(runtime,[Stop|_]),
	Runtime is Stop - Start,
	writef('there are%5r long primes up to%6r [time (ms)%5r]\n',[Num, Limit, Runtime]),
	run(Tail).

do:-	longPrimes(500, LongPrimes),
	writeln('long primes up to 500:'),
	writeln(LongPrimes),
	numlist(0, 7, List),
	maplist([X, Y]>>(Y is 500 * 2**X), List, LimitList),
	run(LimitList).