For generating primes, the Julia interface to libprimesieve is 100 times
faster than the current Julia code, for every number of primes that I
tried. (not bad for your code considering Kim Walisch has poured his soul
into this for several years!) It also generates the nth prime, primes in an
interval, etc. Mathematica is very slow, at least using their example
code Table[Prime[i],{i,10^i}]; certainly limited by their List structure.
The upper limit for the argument to the prime counting function, PrimePi
(Mma 8), is lower by 5 orders of magnitude, than for libprimesieve (and the
tables in PrimeSieve.jl have the largest known value in them). Near the
top of PrimePi's domain and below 10^14, the worst case for countprimes()
is 30 times faster than Mma. (countprimes(10^14-1) is just before a table
entry). Above 10^14, a new, sparser table is used, and countprimes is only
faster, in the worst case, by a factor of about 4.
Of course, it's hard to understand what Mathematica is doing. They do
caching, and probably something similar to the Walisch code I mention
below. Also PrimePi can't compute PrimePi[n] - PrimePi[m] efficiently,
which libprimesieve does. This is why adding the table lookup gives an
effective worst-case speedup of 10^4. I tried Sage, for prime counting, but
it was much slower.
I should also note that Kim Walisch has a library just for the prime
counting function using algorithms that are orders of magnitude faster than
using libprimesieve... hmm probably easy to include that in the Julia
package.
PrimeSieve.jl's nthprime(x,1) gives the first prime larger than x. This
is very inefficient. Mathematica does it almost instantaneously with
NextPrime. The Maxima CAS primes code is very good for some things,
https://github.com/andrejv/maxima/blob/master/src/ifactor.lisp
Its nextprime() seems to be very fast. Its integer factoring routine easily
outperformed Mathematica the last time I checked ( a couple of years ago)
The current Julia integer factoring algorithm is of higher time complexity
so that eventually it just stops. The Maxima prime and factoring functions
are about 800 lines of self-contained common lisp (very little if any
reference to other Maxima code) Maxima's prime() function to fill an array
uses nextprime repeatedly, which is inefficient and should not be used if
other routines are available. But if someone wanted to translate the file
to Julia, it might be worth it for nextprime() and ifactors(). The authors
of the maxima code translated much of it from Gap. I compared the Gap
factoring algorithm to Maxima. It was orders of magnitude slower. It could
be that Maxima mixed in other algorithms; but the Gap code is written in
Gap, which may be interpreted, while the Maxima (that I used) was compiled
sbcl... or I didn't spend enough time with it.
Factoring speed depends on the number. Eg some algorithms are better for
larger numbers, others for larger factors. I tried the msieve program which
was much faster than Maxima (which is faster than Mma) for many numbers
(by faster, I mean sometimes I don't know if the slower one ever
finished). But Maxima was faster for some numbers. I suspect this is
because I did not compile msieve with a library giving "ecm" support, while
Maxima includes ecm algorithms. You can do even better than msieve, but
then the available code is messier and you are entering the realm of the
cutting edge factoring projects/theory. But, msieve seems to have a good
API and to be easy to compile. For factoring integers by end users in a
"general purpose" language like Julia, I think wrapping msieve+ecm would
give extremely good results. But, I'm not an expert and did not do an
exhaustive search. I could be missing something.
On Monday, December 22, 2014 4:08:39 PM UTC+1, Stefan Karpinski wrote:
>
> Very cool. I did the best I could with relatively simple code for prime
> stuff in Base, but this looks far more comprehensive. It would be
> interesting to see comparisons to the simple stuff in base and the
> sophisticated tools in other systems like Mathematica.
>
> On Sun, Dec 21, 2014 at 7:13 PM, <[email protected] <javascript:>>
> wrote:
>
>> Here are some packages for people to look at. In particular,
>> I don't want to duplicate effort.
>>
>> PrimeSieve: This is an interface to the (probably) fastest
>> opensource prime number sieve, and the (probably) most extensive
>> tables of the prime pi function (but only about 4 megabytes). It
>> also generates primes, prime twins, etc.
>>
>> https://github.com/jlapeyre/PrimeSieve.jl
>>
>> DeepConvert: This issue came up when playing with PrimeSieve. It
>> converts numbers in an expression, that will cause overflow in
>> intermediate expressions, to a different type, I experimented
>> with nonstandard string and macro interfaces.
>>
>> Eg, these do what you want:
>> a = bi"[2^63, 2^64]"
>> round(bf"((1+2e-50)-(1+1e-50))/(1e-50)")
>>
>> @bigint function g(x) # does what you want
>> return 2^64 * x
>> end
>>
>> https://github.com/jlapeyre/DeepConvert.jl
>>
>> ZChop: Like Mathematica's Chop.
>>
>> https://github.com/jlapeyre/ZChop.jl
>>
>> PermutationsA:
>>
>> This is bigger than Permutations.jl. I already contributed some
>> things to Permutations, as well.
>>
>> https://github.com/jlapeyre/PermutationsA.jl
>> https://github.com/jlapeyre/PermPlain.jl
>>
>> --John
>>
>>
>
