Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Replace isprime #50

Closed
wants to merge 2 commits into from
Closed

Replace isprime #50

wants to merge 2 commits into from

Conversation

enedil
Copy link

@enedil enedil commented Jun 21, 2017

Replace isprime with weaker, but orders of magnitude faster test. Better to finish testing primes up to 2^16, then performing Miller-Rabin as many times as different prime factors.

Related issue: #49

Replace isprime with weaker, but orders of magnitude faster test. Better to finish testing primes up to 2^16, then performing Miller-Rabin as many times as different prime factors.
@enedil enedil mentioned this pull request Jun 21, 2017
that would do
@@ -264,9 +264,10 @@ function factor!{T<:Integer,K<:Integer}(n::T, h::Associative{K,Int})
n = div(n, p)
end
n == 1 && return h
isprime(n) && (h[n] = 1; return h)
p^2 >= n && (h[n] = 1; return h)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If p^2 >= n, it means that n is the factor of primes less than or equal to p: hence n should have already been factored out by that time, and the condition n == 1 already met!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A more sensible change would be to test PRIMES[end]^2 >= n as a condition to break out of the loop.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rfourquet no, that's exactly what I meant. Either n is 1, or n is prime by this moment.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But can you give me an exampleof initial value of n such that p^2 >= n > 1 at this point?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One example is 63. After factoring out 3*3, we are left with 7, which satisfies conditions : 1< 7 < 3^2. As I said, it is a prime. Also, equal sign is here redundant, but the point still holds.

Let's put it this way:
You might want to check PRIMES[end]^2 >= n. Sure, if it is met, then n is at that point either prime or 1. But we could test something similar - p^2 >= n (p being the rithg now excluded factor). If this condition is met, then suppose n=ab with a, b proper divisors of n. Bot a and b are at this point greater then p. So n=ab>p*p>=n, a contradiction. So n has only trivial factorisation, hence it's a prime.

But it is a considerably stronger test! For if n falls into the latter, it falls also into the former, but nor in converse. That is why I propose such solution.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes of course, thank you and sorry, I should have analysed your change more carefully!

@rfourquet
Copy link
Member

So you are proposing a different strategy: can you show some timings, both what speed-up or slow-down can be expected?

@enedil
Copy link
Author

enedil commented Jun 22, 2017

@rfourquet i'll reply when I'll get to my notebook

@enedil
Copy link
Author

enedil commented Jun 22, 2017

@rfourquet
Ok - some tests:

test_numbers = [#Primes.PRIMES,
                [921236161, 931235651, 941236273, 951236179, 901236131, 961236149, 911236121],
                [2, 3, 5, 13, 17, 61, 623341, 252097800623, 252097800629, 252097800667, 2760727302517],
                [63721832374, 13642987325257, 13789015084387, 25092045343291, 18748357090643, 19019197924633, 1927861641436, 115215177288657],
                [2, 23]*1000000000]
for n in test_numbers
    num = reduce(*, big(1), n)
    tic(); factor(num); toc()
end

And results:

~/c/g/P/src ❯❯❯ julia old.jl 
elapsed time: 0.892513902 seconds
elapsed time: 5.628045027 seconds
elapsed time: 40.111286483 seconds
elapsed time: 5.1307e-5 seconds
~/c/g/P/src ❯❯❯ julia new.jl 
elapsed time: 0.215181272 seconds
elapsed time: 1.384578026 seconds
elapsed time: 43.26424359 seconds
elapsed time: 4.0111e-5 seconds

Every factor in the lists is prime. I commented out Primes.PRIMES array, so that I can ever finish factoring in a reasonable time. As seen, old implementation is rather generally slower (although in some cases another test run might change the output).

As for Primes.PRIMES array, time for new.jl is:

elapsed time: 0.553304032 seconds

This is considerably faster that at least 15 minutes.

@rfourquet
Copy link
Member

OK thanks for benchmarking. What about a compromise: adding an isprime test every 10 or 100 iteration?

@enedil
Copy link
Author

enedil commented Jun 23, 2017

The point is - numbers with large number of prime factors and numbers with very little prime factors, both are rather rare.

You may be right. In fact, Primes.PRIMES table is about 6300 long. Doing Rabin-Miller after every iteration is not necessary, but if it was only 63 times, maybe it would make sense? I'd have to do timings.

@rfourquet
Copy link
Member

When there are very few factors from PRIMES, maybe the test PRIMES[end]^2 > n can be used to break out of the loop early (meaning that either n is prime, or at least has factors greater than PRIMES[end]) ?

@enedil
Copy link
Author

enedil commented Jun 23, 2017

That is not true. What if n = 2039^3? 2039 is prime. We're left with 2039^3 = 8477185319 > 4293001441 = 65521^2 = PRIMES[end]^2. So in fact, n is composite and all its' factors are smaller than 2^16.

@enedil
Copy link
Author

enedil commented Jun 23, 2017

Maybe there should be an optional argument which is how often make tests and default it to about 100? And make docs that say what is the specifics of such argument? And furthermore, don't perform isprime, which tests for first 9 primes, maybe isprime should be splitted so every phase has its own function. How do you see such design?

Also, excuse me if my code is poor, this pr is literally first 3 lines I ever wrote in Julia.

@rfourquet
Copy link
Member

Maybe there should be an optional argument which is how often make tests

I would say it's a bit premature to offer an optional argument for that. And if someone has a wish that this option is available, meaning she has an insight how to optimize this value, maybe would be good then to include this heuristic in this package!

maybe isprime should be splitted so every phase has its own function. How do you see such design?

It seems like a good idea, I thought about that too. Whether you want to implement that in this PR or in another one is your call!

@enedil
Copy link
Author

enedil commented Jun 24, 2017

I'd better if we finished this PR (separation of concerns). Then I will start with clean state. And as far as optional arg is concerned, I found that in very most cases, the less prime tests, the better. So I'd stick with about 200 probably.

Btw, do you have any idea why "AppVeyor build failed"? Logs mention reducing over empty collection, I cannot locate where this function could yield empty Dict.

@rfourquet
Copy link
Member

So I'd stick with about 200 probably.

OK.

Btw, do you have any idea why "AppVeyor build failed"? Logs mention reducing over empty collection, I cannot locate where this function could yield empty Dict.

Try factor(Dict, 1). I will open a PR for that.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants