[ANN] nan_numerics_prime (1.2-beta)

[Julio Di Egidio]

Hello all,

I am pleased to announce release of the library:
nan_numerics_prime - A simple prime number library

Source: <https://github.com/jp-diegidio/Nan.Numerics.Prime-Prolog>
Download: <https://github.com/jp-diegidio/Nan.Numerics.Prime-Prolog/tree/1.2-beta>

This library was developed and tested with:
SWI-Prolog 7.3.24: http://www.swi-prolog.org/

NOTE: Please do not install from the SWI repository as that is broken at the moment.
Instead, download the .zip file from the Download page to your machine, then run:

?- pack_install('/<your_path>/nan_numerics_prime-1.2.zip').

Here follow contents of the README file:

============================================================
Nan.Numerics.Prime 1.2-beta
A simple prime number library
Copyright 2016 Julio P. Di Egidio
Licensed under GNU GPLv3.
http://julio.diegidio.name/Projects/Nan.Numerics.Prime/
https://github.com/jp-diegidio/Nan.Numerics.Prime-Prolog/

The module =prime= provides predicates to test (positive integer) numbers
for primality, find divisors and factor numbers, generate prime numbers in
some interval, find consecutive prime numbers, and save/load all prime
numbers up to some value to/from a file or stream.

Implements a variant of the *Miller-Rabin* primality test that is
_deterministic_ for numbers up to =3317044064679887385961980=, otherwise
it is _probabilistic_ with the number of iterations fixed at =20=. For
better performance, leverages a prime wheel of level =4= and memoization.

All predicates in module =prime= are _safe_, i.e. validate input arguments
and ensure steadfastness. For maximum performance, user code can directly
call the _unsafe_ =public= (not exported) predicates in sub-module
=prime_lgc=.

*NOTE*: Since the primality test in use is _probabilistic_ in general, this
library is not suitable for cryptographic applications.

This library was developed and tested with:
SWI-Prolog 7.3.24: http://www.swi-prolog.org/

Example:

?- pack_install('nan_numerics_prime-1.2.zip').
true.

?- use_module(library(nan_numerics_prime)).
true.

?- time(prime_right(1234567891012345678901234567890123456789011111,P)).
% 1,205 inferences, 0.000 CPU in 0.000 seconds (?% CPU, Infinite Lips)
P = 1234567891012345678901234567890123456789011139.

To be done: Implement prime counting/n-th prime functions.
To be done: Implement probabilitic test error estimates?
To be done: Implement deterministic tests (elliptic curves)?
============================================================

Julio

[Markus Triska]

Hi Julio,

Julio Di Egidio <ju...@diegidio.name> writes:

> nan_numerics_prime - A simple prime number library

Nice, thank you for sharing!

I have already spotted a use case of your library "in the wild":

http://stackoverflow.com/a/39567840

This seems to blend in nicely with CLP(FD) constraints at least in this
specific example.

All the best,
Markus

--
comp.lang.prolog FAQ: http://www.logic.at/prolog/faq/

[Julio Di Egidio]

On Monday, September 19, 2016 at 7:01:48 PM UTC+2, Markus Triska wrote:
> Julio Di Egidio <ju...@diegidio.name> writes:
>
> > nan_numerics_prime - A simple prime number library
>
> Nice, thank you for sharing!
>
> I have already spotted a use case of your library "in the wild":
>
> http://stackoverflow.com/a/39567840
>
> This seems to blend in nicely with CLP(FD) constraints at least in this
> specific example.

Hi Markus, as a numeric library it's still in its infancy, but I'm very glad
it was useful: indeed at least for didactic purposes, about Prolog and some
of the common problems/patterns one encounters in Prolog.

Cheers,

Julio

[burs...@gmail.com]

If I had time, I would really like to parallelize the
algorithm, like is done here, for different Prolog
systems (*):

http://primesieve.org/
(See CPU scaling)

But currently I am too stupid or to old, I even don't
know how wheel factorization works and how it would be
paralellized at all.

Bye

(*) Candidates currently SWI-Prolog and Jekejeke Prolog,
where I have recently conducted a lot of embarassingly
parallel experiment.

Are wheels also embarassingly parallelizable? Is the
web site humbug?

[j4n bur53]

BTW: If you only want to count prime numbers,
you can also do without constructing a sieve,

there are some recursive algorithms based on
number theory and that can do that, notes lost.

Something along this I guess:
https://en.wikipedia.org/wiki/Prime-counting_function#Algorithms_for_evaluating_.CF.80.28x.29

With very few computer resources you can go
very high up, if its only about counting...

burs...@gmail.com schrieb:

[Julio Di Egidio]

There is no point in belaborating prime wheels, a *fixed* wheel of order 4 to maximum 6 does the trick, which is just that of skipping numbers that are divisible by the first few prime numbers, to then be used by sieving and other prime testing or finding methods: and nothing more than order 4 to 6 because after that the memory footprint becomes anyway prohibitive while the gain in efficiency is practically negligible, as in something like 99.99999+% of the composite numbers have been captured already (sorry, cannot find a proper reference and the exact figures right now).

OTOH, keep in mind that with Miller-Rabin we have deterministic tests for numbers up to ~3*10^24, with performances that far exceed any sieve...
(that's implemented in my nan_numerics_prime_prb.pl).

There is in fact room for parallelisation, but that is about *finding divisors and factoring*, as nothing (I presently know of) exists much better than testing the prime numbers one by one in the range [2;sqrt(n)] for n the number to factor, i.e. trial division... (Of course one first checks that the number in question is not a prime number, as testing only is as cheap as I was saying above, then the wheel makes it so that we do not need to check divisibility by the first few prime numbers.)

Anyway it is certainly a broad and very technical topic, and I do not claim I master it all, on the contrary...

Julio

[j4n bur53]

Jokingly in the past, I was thinking about checking
for a prime number as follows:

is_prime(x) = count_to(x)-count_to(x-1).

So is_prime(x) would yield 1, if there is one more
prime between [2..x-1] and [2..x], and otherwise it
would yield 0.

Julio Di Egidio schrieb:

> There is in fact room for parallelisation, but that is

> about*finding divisors and factoring*, as nothing

> (I presently know of) exists much better than testing
> the prime numbers one by one in the range [2;sqrt(n)]
> for n the number to factor, i.e. trial division...
> (Of course one first checks that the number in question
> is not a prime number, as testing only is as cheap as
> I was saying above, then the wheel makes it so that
> we do not need to check divisibility by the first
> few prime numbers.)

Now I am perplexed, and this relates deeply into
logic, i.e. realization problem, witness problem,
Herbrand's theorem,

if you do your is prime test before hand, wouldn't
this is prime test, if x is not a prime, already
throw out a factor of x?

(For example in my is_prime/1 defined via
sieve-less count_to/1 it wouldn't throw out
a factor when is_prime(x) yields 0)

[j4n bur53]

j4n bur53 schrieb:

> (For example in my is_prime/1 defined via
> sieve-less count_to/1 it wouldn't throw out
> a factor when is_prime(x) yields 0)

Addendum: Not automatically...

[Julio Di Egidio]

On Thursday, September 22, 2016 at 7:06:57 AM UTC+2, j4n bur53 wrote:
> Jokingly in the past, I was thinking about checking
> for a prime number as follows:
>
> is_prime(x) = count_to(x)-count_to(x-1).
>
> So is_prime(x) would yield 1, if there is one more
> prime between [2..x-1] and [2..x], and otherwise it
> would yield 0.
>
> Julio Di Egidio schrieb:
> > There is in fact room for parallelisation, but that is
> > about*finding divisors and factoring*, as nothing
> > (I presently know of) exists much better than testing
> > the prime numbers one by one in the range [2;sqrt(n)]
> > for n the number to factor, i.e. trial division...
> > (Of course one first checks that the number in question
> > is not a prime number, as testing only is as cheap as
> > I was saying above, then the wheel makes it so that
> > we do not need to check divisibility by the first
> > few prime numbers.)
>
> Now I am perplexed, and this relates deeply into
> logic, i.e. realization problem, witness problem,
> Herbrand's theorem,

Number theory provides methods to test numbers for primality without finding
any factors. It's all about efficiency.

> if you do your is prime test before hand, wouldn't
> this is prime test, if x is not a prime, already
> throw out a factor of x?

There are prime testing algorithms that also spit out factors, but if there
was any *efficient* way to do that, we'd have solved the factoring problem...

Julio

[j4n bur53]

burs...@gmail.com schrieb:

> But currently I am too stupid or to old, I even don't

Anyway, you have to start when your 15 or 16 years old:
http://www.maths.ox.ac.uk/node/16561

[j4n bur53]

Lets win a prize:
http://www.math.unt.edu/~mauldin/beal.html

Yesterday, on another but similar line, just found
n and m such that (x^n)^2+(y^m)^2 = z^2, and I
used CLP(FD), with the new infinite set feature:

Jekejeke Prolog 2, Runtime Library 1.1.6
(c) 1985-2016, XLOG Technologies GmbH, Switzerland

?- use_module(library(finite/clpfd)).
% 20 konsultiert und 0 entladen in 1859 ms.
Ja

?- [X,Y,Z] ins 1..sup, X*X*X*X+Y*Y*Y*Y*Y*Y #= Z*Z, label([X,Y,Z]).
X = 6,
Y = 3,
Z = 45

?- 6^4+3^6 =:= 45^2.
Ja

But they are not co-prime. :-( And its not exactly
the Beal Conjecture, but it was also fun!

Can be also done with SWI-Prolog, but there the finite
set feature is currently missing, but one can directly
write X^4 + Y^6 #= Z^2. I do not yet support (^)/2.

I don't know how much CLP(FD) could profit from larger
or faster prime stuff. I am currently using a native
implementd gcd/2 evaluable function thats all.(*)

I don't know how much confidence can be gained by
checking a conjecture respectively try to find
counter examples up to 10^10 or so?

Maybe if we don't find a counter example in a very large
range, its time to give a try in proving the conjecture.
But will be able to prove it, even if its a true
conjecture?

Bye

(*)
https://en.wikipedia.org/wiki/Binary_GCD_algorithm


Julio Di Egidio schrieb:

[j4n bur53]

j4n bur53 schrieb:

> I don't know how much CLP(FD) could profit from larger
> or faster prime stuff. I am currently using a native
> implementd gcd/2 evaluable function thats all.(*)

Its is used for fast failure of linear
forms A1*X1+..+An*Xn #= B.

fast failure: B cannot be divided by gcd(A1,..,An)

If there is no fast failure, we can normalize
it, or its already in normal form.

Also applicable to inequalities to some extend,
for example A1*X1+..+An*Xn in G+1..2*G-1 is also
subject to fast failure.

Not an extremly complex heuristic, but a cheap one.
And an useful one, I guess this failure happens quite
often, but I don't have exact measures/figures here.

Also all this CLP(FD) stuff is trial and error,
the labeling method, probably a horror for a real
mathematicians.

[j4n bur53]

j4n bur53 schrieb:

>
> I don't know how much confidence can be gained by
> checking a conjecture respectively try to find
> counter examples up to 10^10 or so?

Oops, according to Peter Norvig:

"But Witold Jarnicki and David Konerding already did that: they wrote a
C++ program that built a table of Cz(modp) up to 50005000, and, in
parallel across thousands of machines, searched for A,B up to 200,000
and x,y up to 5,000, but found no counterexamples. On a smaller scale,
Edwin P. Berlin Jr. searched all Cz up to 1017 and also found nothing.
So I don't think it is worthwhile to continue on that path."
http://www.norvig.com/beal.html

[j4n bur53]

j4n bur53 schrieb:

> Can be also done with SWI-Prolog, but there the finite
> set feature is currently missing, but one can directly

Corr.: infinite set feature is miissing

But the follow works fine in SWI-Prolog:

?- time(([X,Y,Z] ins 1..10000, X^4+Y^6 #= Z^2, label([X,Y,Z]))).
% 124,509 inferences, 0.027 CPU in 0.027 seconds (99% CPU, 4563611 Lips)

X = 6,
Y = 3,
Z = 45

% 165,052 inferences, 0.035 CPU in 0.036 seconds (98% CPU, 4750928 Lips)
X = 15,
Y = 10,
Z = 1025
% 365,902 inferences, 0.082 CPU in 0.083 seconds (99% CPU, 4487943 Lips)
X = 36,
Y = 12,
Z = 2160
% 205,035 inferences, 0.049 CPU in 0.050 seconds (98% CPU, 4142372 Lips)
X = 48,
Y = 12,
Z = 2880
% 647,656 inferences, 0.138 CPU in 0.140 seconds (99% CPU, 4687114 Lips)
X = 90,
Y = 15,
Z = 8775
% 96,346 inferences, 0.023 CPU in 0.024 seconds (96% CPU, 4163433 Lips)
false.

If the intervsl of ins is chosen bigger, the thing gets slower:

time(([X,Y,Z] ins 1..100000000, X^4+Y^6 #= Z^2, label([X,Y,Z]))).
% 2,715,197 inferences, 0.893 CPU in 0.903 seconds (99% CPU, 3040102 Lips)

X = 6,
Y = 3,
Z = 45

% 3,941,142 inferences, 1.252 CPU in 1.262 seconds (99% CPU, 3147274 Lips)
X = 15,
Y = 10,
Z = 1025
% 8,883,966 inferences, 2.771 CPU in 2.785 seconds (99% CPU, 3206172 Lips)
X = 36,
Y = 12,
Z = 2160
% 5,166,670 inferences, 1.817 CPU in 1.824 seconds (100% CPU, 2843287 Lips)
X = 48,
Y = 12,
Z = 2880
% 17,499,301 inferences, 6.189 CPU in 6.212 seconds (100% CPU, 2827376 Lips)
X = 90,
Y = 15,
Z = 8775 .

But its the same solutions!!!!???? In my infinite search I don't have
an upper bound, so I cannot make this comparison. I didn't make much
measures yet, only made this strange observatiin of timing now.

Bye

[j4n bur53]

j4n bur53 schrieb:

> But its the same solutions!!!!???? In my infinite search I don't have
> an upper bound, so I cannot make this comparison. I didn't make much
> measures yet, only made this strange observatiin of timing now.

Probably not a defect, maybe an artefact of the solution structure,
Y is smaller than X, and there is only one Z. Which results for a
larger range in bigger unnessesary exhaust. Just a guess.

[Julio Di Egidio]

What are you, stupid compulsive that you have to spam all threads??

*Plonk*

Julio

[burs...@gmail.com]

> What are you, stupid compulsive that you have to spam all threads?

My spologies for also hitchhiking your threads
here and not only on SWI-Prolog forum.

Maybe I take usenet for twitter sometimes, although I feel
harassed myself when my threads are obviously abused.

Fortunately I wasnt victim for long, maybe I am too fast
with spamming right now, so harasser are stressed.

[j4n bur53]

Just research J.H.Lambert, and then this:
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.233.3230

Please enjoy this history bit.

Julio Di Egidio schrieb:

[burs...@gmail.com]

Are there some new additions to the library?

I would be interested in Euler’s totient function.
Interestingly this Python library has also prime wheels.

http://pythonhosted.org/eulerlib/eulerlib.html

Am Samstag, 27. August 2016 07:53:45 UTC+2 schrieb Julio Di Egidio:

[burs...@gmail.com]

They use Euler's product formula:

if(num < 1):
return 0
if(num == 1):
return 1
if(num in self.primes_table):
return num-1
pfs = self.prime_factors_only(num)
prod = num
for pfi in pfs:
prod = prod*(pfi-1)/pfi
return prod

http://pythonhosted.org/eulerlib/_modules/eulerlib/numtheory.html#Divisors.phi

[Julio Di Egidio]

On Friday, February 3, 2017 at 2:29:16 PM UTC+1, burs...@gmail.com wrote:

> Are there some new additions to the library?
>
> I would be interested in Euler’s totient function.

I am going to published version 1.3.0 shortly, but it is only going to be a
major clean up and consolidation. In the near future, I might just play with
parallelising the factoring function (using my answer sources), anyway this
remains a "simple" library for now.

For the future I don't know yet: for one thing, to progress further I'd really
need support from a professional mathematician (or I should finally subscribe
to math overflow, I guess...). But there is also the issue of the underlying
technology: to get a "professional" product I'd want to rethink the "Prolog
engine" from scratch, which might in fact be my next "big" project: or maybe
closed numerics in C on top of GMP... I just don't know yet, OTOH, working on
this library does not look like a priority: there are several professional
packages available already. As a "simple" library I think it makes sense
because, for simple usage, on one hand, it provides pretty decent functionality,
on the other hand, one does not want the burden of the complexity involved in
deploying and using more professional solutions.

> Interestingly this Python library has also prime wheels.

Using prime wheels to at least prune the search space is pretty common.

Julio