1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    2/*	Nan.Numerics.Prime
    3	A simple prime number library
    4	Copyright 2016 Julio P. Di Egidio
    5	<mailto:julio@diegidio.name>
    6	<http://julio.diegidio.name/Projects/Nan.Numerics.Prime/>
    7	
    8	This file is part of Nan.Numerics.Prime.
    9	
   10	Nan.Numerics.Prime is free software: you can redistribute it and/or modify
   11	it under the terms of the GNU General Public License as published by
   12	the Free Software Foundation, either version 3 of the License, or
   13	(at your option) any later version.
   14	
   15	Nan.Numerics.Prime is distributed in the hope that it will be useful,
   16	but WITHOUT ANY WARRANTY; without even the implied warranty of
   17	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   18	GNU General Public License for more details.
   19	
   20	You should have received a copy of the GNU General Public License
   21	along with Nan.Numerics.Prime.  If not, see <http://www.gnu.org/licenses/>.
   22*/

   23%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   24
   25% (SWI-Prolog 7.3.25)

   26
   27:- module(prime_lgc, []).
   28
   29:- public
   30	test_/1,		% +N:posint

   31	div_/2,			% +N:posint, -P:prime

   32	div_rev_/2,		% +N:posint, -P:prime

   33	fact_/2,		% +N:posint, -PFs:list(pfact)

   34	gen_/2,			% +Inf:posint, -P:prime

   35	gen_/3,			% +Inf:posint, +Sup:posint, -P:prime

   36	gen_p_/2,		% +L:prime, -P:prime

   37	gen_p_/3,		% +L:prime, +H:prime, -P:prime

   38	gen_rev_/2,		% +Sup:posint, -P:prime

   39	gen_rev_/3,		% +Inf:posint, +Sup:posint, -P:prime

   40	gen_rev_p_/2,	% +H:prime, -P:prime

   41	gen_rev_p_/3,	% +L:prime, +H:prime, -P:prime

   42	next_/2,		% +N:posint, -P:prime

   43	next_p_/2,		% +P0:prime, -P:prime

   44	prev_/2,		% +N:posint, -P:prime

   45	prev_p_/2,		% +P0:prime, -P:prime

   46	right_/2,		% +N:posint, -P:prime

   47	left_/2.		% +N:posint, -P:prime

   68:- use_module(nan_numerics_prime_mem).
   69:- use_module(nan_numerics_prime_whl).
   70:- use_module(nan_numerics_prime_prb).
   71:- use_module(library(lists)).
   72
   73%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 test_(+N:posint) is semidet

   76%

   77%	True if N is a prime number.

   78
   79test_(N) :-
   80	prime_mem:get_(N), !.
   81test_(N) :-
   82	prime_whl:test_(N, Cert),
   83	test__do(Cert, N).

 div_(+N:posint, -P:prime) is semidet

   86%

   87%	True if N is a composite number with P its smallest prime divisor.

   88
   89div_(N, P) :-
   90	N > 3,
   91	\+ test_(N),
   92	div__do(N, P).

 div_rev_(+N:posint, -P:prime) is semidet

   95%

   96%	True if N is a composite number with P its greatest prime divisor.

   97
   98div_rev_(N, P) :-
   99	N > 3,
  100	\+ test_(N),
  101	div_rev__do(N, P).

 fact_(+N:posint, -PFs:list(pfact)) is det

  104%

  105%	PFs is the list of all prime factors of N in _ascending_ order of the

  106%	prime divisors.

  107%

  108%	Elements of PFs are of the form =|P^F|= with _P_ the prime divisor and

  109%	_F_ the corresponding power.

  110%

  111%	If N is equal to =1= or if N is a prime number, PFs is =|[N^1]|=.

  112
  113fact_(N, PFs) :-
  114	fact__sel_rev(N, [], PFsRev),
  115	reverse(PFsRev, PFs).
  116
  117fact__sel_rev(N, PFs0, PFs) :-
  118	div_(N, P), !,
  119	N1 is N // P,
  120	fact__pfs(P, PFs0, PFs1),
  121	fact__sel_rev(N1, PFs1, PFs).
  122fact__sel_rev(P, PFs0, PFs) :-
  123	fact__pfs(P, PFs0, PFs).
  124
  125fact__pfs(P, [P^F| PFs], [P^F1| PFs]) :- !, F1 is F + 1.
  126fact__pfs(P, PFs, [P^1| PFs]).

 gen_(+Inf:posint, -P:prime) is multi

 gen_(+Inf:posint, +Sup:posint, -P:prime) is nondet

  130%

  131%	Generates in _ascending_ order all prime numbers P greater than or

  132%	equal to Inf, and less than or equal to Sup in the variant with arity

  133%	=3=.  Fails if the prime to the left of Sup is less than the prime to

  134%	the right of Inf.

  135
  136gen_(Inf, P) :-
  137	right_(Inf, L),
  138	gen_p_(L, P).
  139
  140gen_(Inf, Sup, P) :-
  141	gen__lh(Inf, Sup, L, H),
  142	gen_p_(L, H, P).
  143
  144gen__lh(Inf, Sup, L, H) :-
  145	Sup >= Inf,
  146	right_(Inf, L),
  147	left_(Sup, H),
  148	H >= L.

 gen_p_(+L:prime, -P:prime) is multi

 gen_p_(+L:prime, +H:prime, -P:prime) is nondet

  152%

  153%	Generates in _ascending_ order all prime numbers P starting from L, and

  154%	up to H in the variant with arity =3=.  Fails if H is less than L.

  155
  156gen_p_(P, P).
  157gen_p_(L, P) :-
  158	next_p_(L, L1),
  159	gen_p_(L1, P).
  160
  161gen_p_(P, H, P) :- P >= H, !.
  162gen_p_(P, _, P).
  163gen_p_(L, H, P) :-
  164	next_p_(L, L1),
  165	gen_p_(L1, H, P).

 gen_rev_(+Sup:prime, -P:prime) is nondet

 gen_rev_(+Inf:posint, +Sup:posint, -P:prime) is nondet

  169%

  170%	Generates in _descending_ order all prime numbers P less than or equal

  171%	to Sup, and greater than or equal to Inf in the variant with arity =3=.

  172%	Fails if Sup is equal to =1= or if the prime to the left of Sup is less

  173%	than the prime to the right of Inf.

  174
  175gen_rev_(Sup, P) :-
  176	left_(Sup, H),
  177	gen_rev_p_(H, P).
  178
  179gen_rev_(Inf, Sup, P) :-
  180	gen__lh(Inf, Sup, L, H),
  181	gen_rev_p_(L, H, P).

 gen_rev_p_(+H:prime, -P:prime) is multi

 gen_rev_p_(+L:prime, +H:prime, -P:prime) is nondet

  185%

  186%	Generates in _descending_ order all prime numbers P starting from H,

  187%	and down to L in the variant with arity =3=.  Fails if H is less than

  188%	L.

  189
  190gen_rev_p_(2, 2) :- !.
  191gen_rev_p_(P, P).
  192gen_rev_p_(H, P) :-
  193	prev_p_(H, H1),
  194	gen_rev_p_(H1, P).
  195
  196gen_rev_p_(L, P, P) :- P =< L, !.
  197gen_rev_p_(_, P, P).
  198gen_rev_p_(L, H, P) :-
  199	prev_p_(H, H1),
  200	gen_rev_p_(L, H1, P).

 next_(+N:posint, -P:prime) is det

  203%

  204%	P is the smallest prime number greater than N.

  205
  206next_(N, P) :-
  207	right_(N, P0),
  208	next__sel(N, P0, P).
  209
  210next__sel(P0, P0, P) :- !,
  211	next_p_(P0, P).
  212next__sel(_, P, P).

 next_p_(+P0:prime, -P:prime) is det

  215%

  216%	P is the smallest prime number greater than P0.

  217
  218next_p_(P0, P) :-
  219	prime_mem:get_(P0, P), !.
  220next_p_(P0, P) :-
  221	N is P0 + 1,
  222	right_(N, P),
  223	prime_mem:add_(P0, P).

 prev_(+N:posint, -P:prime) is semidet

  226%

  227%	P is the greatest prime number less than N.  Fails if N is less than or

  228%	equal to =2=.

  229
  230prev_(N, P) :-
  231	left_(N, P0),
  232	prev__sel(N, P0, P).
  233
  234prev__sel(P0, P0, P) :- !,
  235	prev_p_(P0, P).
  236prev__sel(_, P, P).

 prev_p_(+P0:prime, -P:prime) is semidet

  239%

  240%	P is the greatest prime number less than P0.  Fails if P is equal to

  241%	=2=.

  242
  243prev_p_(P0, P) :-
  244	prime_mem:get_(P, P0), !.
  245prev_p_(P0, P) :-
  246	N is P0 - 1,
  247	left_(N, P),
  248	prime_mem:add_(P, P0).

 right_(+N:posint, -P:prime) is det

  251%

  252%	P is the smallest prime number greater than or equal to N.

  253
  254right_(N, P) :-
  255	prime_whl:right_(N, W, Cert),
  256	right__sel(Cert, W, P).
  257
  258right__sel(Cert, P, P) :-
  259	test__do(Cert, P), !.
  260right__sel(_, W, P) :-
  261	W1 is W + 2,
  262	right_(W1, P).

 left_(+N:posint, -P:prime) is semidet

  265%

  266%	P is the greatest prime number less than or equal to N.  Fails if N is

  267%	equal to =1=.

  268
  269left_(N, P) :-
  270	prime_whl:left_(N, W, Cert),
  271	left__sel(Cert, W, P).
  272
  273left__sel(Cert, P, P) :-
  274	test__do(Cert, P), !.
  275left__sel(_, W, P) :-
  276	W1 is W - 2,
  277	left_(W1, P).
  278
  279%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  280
  281%	test__do(+Cert:boolean, +N:posint) is semidet.

  282
  283test__do(true, _) :- !.
  284test__do(_, N) :-
  285	prime_prb:test_(N, _).
  286
  287%	div__do(+N:posint, -P:prime) is semidet.

  288
  289div__do(N, P) :-
  290	% Sup is floor(sqrt(N)),		% TODO: Integrate =isqrt= function from GMP?

  291	gen_p_(2, P),
  292	% (	P > Sup, !, fail			% 

  293	(	P * P > N, !, fail			% 

  294	;	N mod P =:= 0
  295	), !.
  296
  297%	div_rev__do(+N:posint, -P:prime) is semidet.

  298
  299div_rev__do(N, P) :-
  300	Sup is N // 2,
  301	left_(Sup, H),
  302	gen_rev_p_(H, P),
  303	N mod P =:= 0, !.
  304
  305%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%