<div dir="ltr">Sorry, I think this does not work for the numbers 3, 5, 7, 11, 13, 17, 19, 23, 29 and 31.</div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">El mié., 9 oct. 2019 a las 0:34, Andres Valloud via Cuis-dev (<<a href="mailto:cuis-dev@lists.cuis.st">cuis-dev@lists.cuis.st</a>>) escribió:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">I played a bit with the guard clauses and found the attached one is <br>
simpler yet just as fast.<br>
<br>
On 10/8/19 20:11, Andres Valloud via Cuis-dev wrote:<br>
> Regarding each+31, sure, 30*k+1 comes first, except when k = 0 because <br>
> why would anyone try dividing by 1. So this is why that case is shifted <br>
> by 30. However, when written this way, the actual divisor evaluation <br>
> order is 31, 7, 11, and so on. It's more likely that a random integer <br>
> is 0 mod 7 than 0 mod 31, and the sooner one detects exact division, the <br>
> sooner the computation can stop. Because of that, I think the each+31 <br>
> case should be the last one in the division loop.<br>
> <br>
> On 10/8/19 19:17, Agustín Sansone via Cuis-dev wrote:<br>
>> Hello!<br>
>><br>
>> I agree with you. I don't think isPrime should send isProbablyPrime <br>
>> because it could fail in the future.<br>
>> I leave you here the implementation with this taken care of.<br>
>> I wrote the (each+31) case first in the trial division loop, because <br>
>> it is testing the 30*k+1 case, wich I also wrote first in the comment.<br>
>><br>
>> Thanks,<br>
>> Agustín<br>
>><br>
>> El mar., 8 oct. 2019 a las 8:11, Juan Vuletich via Cuis-dev <br>
>> (<<a href="mailto:cuis-dev@lists.cuis.st" target="_blank">cuis-dev@lists.cuis.st</a> <mailto:<a href="mailto:cuis-dev@lists.cuis.st" target="_blank">cuis-dev@lists.cuis.st</a>>>) escribió:<br>
>><br>
>> Hi Folks,<br>
>><br>
>> I agree with Andrés comments, and will just focusing on the proposed<br>
>> changes.<br>
>> (snip)<br>
>><br>
>> On 10/8/2019 2:20 AM, Andres Valloud via Cuis-dev wrote:<br>
>> > Agustin, nice to see someone looking into these kinds of things <br>
>> :).<br>
>> > ...<br>
>> >> * The *raisedToInteger: exp modulo: m *method**in Integer has<br>
>> a very<br>
>> >> big problem. If we compute, for example, /"5 raisedTo: 0 <br>
>> modulo:<br>
>> >> 0"/, this returns 1. This means, that according to<br>
>> Smalltalk, the<br>
>> >> rest of the division by 0 of 1(=5^0) is equal to 1 (Yes,<br>
>> division by<br>
>> >> zero!!). I think you can see the problem. This is due the<br>
>> first line<br>
>> >> of the method, that says /"(exp = 0) ifTrue: [^ 1].", /does<br>
>> >> not check anything else. This problem can be easily fixed by<br>
>> >> checking if m=0 just before.<br>
>> ><br>
>> > I agree, the current code appears to be wrong. The initials on <br>
>> the<br>
>> > code belong to Juan Vuletich and Nicolas Cellier. Guys, is there<br>
>> > reason why e.g. 5 raisedTo: 0 modulo: 0 should answer 1 rather <br>
>> than<br>
>> > fail? I don't see any, but...<br>
>> ><br>
>> > Assuming the code is broken and needs to be fixed, <br>
>> alternatively one<br>
>> > could also write the initial guard clause like this:<br>
>> ><br>
>> > n = 0 ifTrue: [^1 \\ m].<br>
>> ><br>
>> > because the case m = 0 will fail.<br>
>> > ...<br>
>><br>
>> Just added this suggestion as an update to GitHub. Andrés, I did it<br>
>> with<br>
>> your author initials, it's your code!<br>
>><br>
>> > ...<br>
>> >> * The *isPrime *method in Integer makes some optimization in<br>
>> order to<br>
>> >> run the algorithm in O(sqrt(self)) instead of the naive <br>
>> way in<br>
>> >> O(self). This is very intelligent, but the constant factor<br>
>> of this<br>
>> >> method can be still improved significantly. I share with <br>
>> you my<br>
>> >> implementation of *isPrimeFast *with a small explanation. <br>
>> This<br>
>> >> implementation runs in general more than 3 times faster <br>
>> than the<br>
>> >> actual one. I leave you a test that checks the correctness<br>
>> of it as<br>
>> >> well, and some other tests that check this complexity I<br>
>> mentioned.<br>
>> ><br>
>> > I see what you did there, but I do not know how to reproduce the<br>
>> time<br>
>> > tests you mention. I built a sample of integers between 1 and<br>
>> 2^32 (I<br>
>> > didn't go up to 2^64 because that would require O(2^32) operations<br>
>> > each, and I want that to finish in reasonable time), and I get<br>
>> > something like a 2x performance improvement rather than 3x. This<br>
>> > seems to make sense because the approach you propose halves the \\<br>
>> > operations (8 remain out of the 16 the current code is doing, for<br>
>> > every batch of 30 potential divisors).<br>
>> ><br>
>> > slicer := 1024.<br>
>> > thickness := 255.<br>
>> > maxK := 1 bitShift: 32.<br>
>> > integers := 1 to: maxK by: maxK // slicer<br>
>> > :: inject: OrderedCollection new<br>
>> > into: [:t :x |<br>
>> > t add: x.<br>
>> > thickness timesRepeat: [t add: t last + 1].<br>
>> > t yourself]<br>
>> > :: asArray.<br>
>> > Time millisecondsToRun:<br>
>> > [1 to: integers size do:<br>
>> > [:x | (integers at: x) isPrime]]<br>
>> ><br>
>> > Using the above code (which I could not format more nicely in this<br>
>> > email), I get about 4.8s for isPrime, and about 2.4s for <br>
>> isPrimeFast.<br>
>> ><br>
>> > Generally, isPrime shouldn't send isProbablyPrime because <br>
>> isPrime is<br>
>> > meant to be deterministic, and one shouldn't assume that the<br>
>> > probabilistic algorithm today will happen to provide the correct<br>
>> > deterministic answer tomorrow.<br>
>> ><br>
>> > Why is the (each+31) case first in the trial division loop?<br>
>> ><br>
>> > Andres.<br>
>><br>
>> I'll wait for your consensus on what to do here. Making isPrime not<br>
>> send<br>
>> isProbablyPrime sounds reasonable to me, but folks, I prefer to wait<br>
>> for<br>
>> your suggestion.<br>
>><br>
>> Thanks,<br>
>><br>
>> -- Juan Vuletich<br>
>> <a href="http://www.cuis-smalltalk.org" rel="noreferrer" target="_blank">www.cuis-smalltalk.org</a> <<a href="http://www.cuis-smalltalk.org" rel="noreferrer" target="_blank">http://www.cuis-smalltalk.org</a>><br>
>> <a href="https://github.com/Cuis-Smalltalk/Cuis-Smalltalk-Dev" rel="noreferrer" target="_blank">https://github.com/Cuis-Smalltalk/Cuis-Smalltalk-Dev</a><br>
>> <a href="https://github.com/jvuletich" rel="noreferrer" target="_blank">https://github.com/jvuletich</a><br>
>> <a href="https://www.linkedin.com/in/juan-vuletich-75611b3" rel="noreferrer" target="_blank">https://www.linkedin.com/in/juan-vuletich-75611b3</a><br>
>> @JuanVuletich<br>
>><br>
>> -- Cuis-dev mailing list<br>
>> <a href="mailto:Cuis-dev@lists.cuis.st" target="_blank">Cuis-dev@lists.cuis.st</a> <mailto:<a href="mailto:Cuis-dev@lists.cuis.st" target="_blank">Cuis-dev@lists.cuis.st</a>><br>
>> <a href="https://lists.cuis.st/mailman/listinfo/cuis-dev" rel="noreferrer" target="_blank">https://lists.cuis.st/mailman/listinfo/cuis-dev</a><br>
>><br>
>><br>
-- <br>
Cuis-dev mailing list<br>
<a href="mailto:Cuis-dev@lists.cuis.st" target="_blank">Cuis-dev@lists.cuis.st</a><br>
<a href="https://lists.cuis.st/mailman/listinfo/cuis-dev" rel="noreferrer" target="_blank">https://lists.cuis.st/mailman/listinfo/cuis-dev</a><br>
</blockquote></div>