<div dir="ltr"><div><div><div><div><div><div>Thanks David, good!<br></div>I'm really happy about the decision to introduce changes by small steps.<br>And we are lucky that I recently added the erf function in ArbitraryPrecisionFloat.<br>Otherwise, with 32 bits limbs, the division bug could have been un-noticed for years.<br><br></div>When I contemplate the division code, I wonder how to test such beast...<br></div>Long method, high cyclomatic complexity...<br></div><div>Given the very low probability of failure, testing with random values is not appropriate.<br></div>I wonder if some kind of formal proof would be possible.<br></div>The first thing would be to describe precisely the invariants.<br></div>If I have an occasion, I'll try to submit the case to the frama-C guys (<a href="http://frama-c.com/about.html">http://frama-c.com/about.html</a>).<br><br></div><div class="gmail_extra"><br><div class="gmail_quote">2016-03-28 4:05 GMT+02:00 David T. Lewis <span dir="ltr"><<a href="mailto:lewis@mail.msen.com" target="_blank">lewis@mail.msen.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
Nicolas,<br>
<br>
I updated VMM trunk with all of your changes for LargeIntegersPlugin copied from<br>
the latest oscog branch. All number tests pass, and your ArbitraryPrecisionFloat<br>
example works properly now.<br>
<br>
Dave<br>
<div class="HOEnZb"><div class="h5"><br>
<br>
On Mon, Mar 28, 2016 at 12:09:37AM +0200, Nicolas Cellier wrote:<br>
><br>
> OK, I think I get it.<br>
> The problem lies in:<br>
><br>
> "Correct overestimate of q.<br>
> Max of 2 iterations through loop -- see Knuth vol. 2"<br>
> j < 3<br>
> ifTrue: [r3 := 0]<br>
> ifFalse: [r3 := pRem at: j - 3].<br>
><br>
> [(t < hi<br>
> or: [t = hi and: [r3 < lo]])<br>
> ifTrue:<br>
> ["i.e. (t,r3) < (hi,lo)"<br>
> q := q - 1.<br>
> lo < dnh<br>
> ifTrue:<br>
> [hi := hi - 1.<br>
> lo := lo + 16r100 - dnh]<br>
> ifFalse:<br>
> [lo := lo - dnh].<br>
> cond := hi >= dh]<br>
> ifFalse: [cond := false].<br>
> cond]<br>
> whileTrue: [hi := hi - dh]].<br>
><br>
> What happens here is that hi may become negative.<br>
> But since I declared it unsigned, t<hi and cond := hi >= dh remains false<br>
> quite a long time and the quotient q is decremented way too many times...<br>
><br>
> Here we divide (r1,r2,r3,...) by (dh,dnh,...)<br>
> we are trying to guess q based on 2 limbs<br>
> (r1,r2) = q*dh + t that is (r1,r2-t) = q*dh<br>
> (hi,lo) =q*dnh<br>
><br>
> Thus q*(dh,dnh)=(r1,r2-t+hi,lo)<br>
><br>
> If t<hi then q*(dh,dnh) is greater than (r1,r2) and we have an overestimate<br>
> of q<br>
> Same if t=hi and lo<r3.<br>
><br>
> Here we have the tricky case when t=0, r3=0 and lo>0, and hi=dnh+1<br>
> In the first loop, h gets decremented, thus hi=dnh.<br>
> Then it goes to 0, and gets decremented a second time...<br>
><br>
> The same algorithm is used with 32bits limbs, so it may also happen in my<br>
> own VM branch. I believe it must be even rarest an event though (2^32<br>
> possible values for hi makes hi=dnh+1 less probable than with 2^8<br>
> combinations).<br>
><br>
> I will fix it tonight or tomorrow.<br>
><br>
><br>
> 2016-03-27 19:03 GMT+02:00 Nicolas Cellier <<br>
> <a href="mailto:nicolas.cellier.aka.nice@gmail.com">nicolas.cellier.aka.nice@gmail.com</a>>:<br>
><br>
> > Thanks David.<br>
> ><br>
> > By running the snippet, instrumenting erf, printing each and every<br>
> > intermediate value to a text file, and diffing:<br>
> > {<br>
> > (10 asArbitraryPrecisionFloatNumBits: 33) erf.<br>
> > (10 asArbitraryPrecisionFloatNumBits: 34) erf.<br>
> > }<br>
> ><br>
> > I reduced the failing test to:<br>
> ><br>
> > a :=<br>
> > 1598335257761788022467377781654101148543282249044465229239888363328190330275719997501596724768507889233831388734160190922469363547795602076820594918.<br>
> > b := 49612.<br>
> > self assert: a - ((a quo: b)*b) < b<br>
> ><br>
> > So it must be some nasty sign promotion (or lack of?) in division.<br>
> > Normally, this could be simulated, but with signedness, I'm never sure,<br>
> > maybe my next step will be a breakpoint in xcode/gdb...<br>
> ><br>
> ><br>
> > Curiously, (a rem: b) is correct !<br>
> ><br>
> > 2016-03-27 17:36 GMT+02:00 David T. Lewis <<a href="mailto:lewis@mail.msen.com">lewis@mail.msen.com</a>>:<br>
> ><br>
> >><br>
> >> I confirm the same results for your ArbitraryPrecisionFloat example below,<br>
> >> testing on Linux with a 64 bit interpreter VM and 32 bit image.<br>
> >><br>
> >> Image<br>
> >> -----<br>
> >> /home/lewis/squeak/Squeak4.6/squeak.37.image<br>
> >> Squeak4.6<br>
> >> latest update: #15723<br>
> >> Current Change Set: Unnamed2<br>
> >> Image format 6504 (32 bit)<br>
> >><br>
> >> Virtual Machine<br>
> >> ---------------<br>
> >> /usr/local/lib/squeak/4.15.4-3634/squeakvm<br>
> >> Squeak4.5 of 10 December 2015 [latest update: #1195]<br>
> >> Unix built on Mar 25 2016 13:27:51 Compiler: 4.9.2<br>
> >> platform sources revision 3634<br>
> >> VMMaker versionString 4.15.4<br>
> >><br>
> >> Dave<br>
> >><br>
> >><br>
> >> On Sun, Mar 27, 2016 at 03:25:58PM +0200, Nicolas Cellier wrote:<br>
> >> ><br>
> >> > Hi Levente, thank you.<br>
> >> > Here is a follow up:<br>
> >> ><br>
> >> > I was able to reproduce the same symptom on the classical interpreter<br>
> >> VM.<br>
> >> > <rant> It's just that for MacOSX, there is no recent VM available, and<br>
> >> that<br>
> >> > the MacOS build is somehow abandonware, so it took much more time than<br>
> >> > necessary :( </rant><br>
> >> ><br>
> >> > So the changes that have been integrated on Interpreter Vm in july 2014<br>
> >> and<br>
> >> > that I thought correct were not.<br>
> >> > It's certainly related to a LargeInteger plugin primitive.<br>
> >> > Those are triggered when nbits of magnitude > 64.<br>
> >> > (otherwise, <= 64 bits, the LargeInteger numbered primitives do their<br>
> >> job).<br>
> >> ><br>
> >> > Now the job is to produce/isolate an elementary failing LargeInteger<br>
> >> test.<br>
> >> > But you need to generate the VM from VMMaker.oscog head for COG/STACK<br>
> >> > V3/SPUR 32/64 bits<br>
> >> > or use a sufficiently recent interpreter VM (> july 2014) to trigger<br>
> >> such<br>
> >> > failure.<br>
> >> ><br>
> >> ><br>
> >> > --------------------------------------------<br>
> >> ><br>
> >> > Here is my starting point:<br>
> >> ><br>
> >> > If you load ArbitraryPrecisionFloat and Tests from<br>
> >> > <a href="http://www.squeaksource.com/ArbitraryPrecisionFl" rel="noreferrer" target="_blank">http://www.squeaksource.com/ArbitraryPrecisionFl</a><br>
> >> > then you'll see this strange pattern:<br>
> >> ><br>
> >> > (53 to: 70) collect: [:i |<br>
> >> > i -> (10 asArbitraryPrecisionFloatNumBits: i) erf ]<br>
> >> ><br>
> >> > {53->(ArbitraryPrecisionFloat readFrom: '0.9999998981732067' readStream<br>
> >> > numBits: 53) .<br>
> >> > 54->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 54) .<br>
> >> > 55->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 55) .<br>
> >> > 56->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 56) .<br>
> >> > 57->(ArbitraryPrecisionFloat readFrom: '0.99999989817320666' readStream<br>
> >> > numBits: 57) .<br>
> >> > 58->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 58) .<br>
> >> > 59->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 59) .<br>
> >> > 60->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 60) .<br>
> >> > 61->(ArbitraryPrecisionFloat readFrom: '0.9999998981732066655'<br>
> >> readStream<br>
> >> > numBits: 61) .<br>
> >> > 62->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 62) .<br>
> >> > 63->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 63) .<br>
> >> > 64->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 64) .<br>
> >> > 65->(ArbitraryPrecisionFloat readFrom: '0.99999989817320666564'<br>
> >> readStream<br>
> >> > numBits: 65) .<br>
> >> > 66->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 66) .<br>
> >> > 67->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 67) .<br>
> >> > 68->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 68) .<br>
> >> > 69->(ArbitraryPrecisionFloat readFrom: '0.999999898173206665646'<br>
> >> > readStream numBits: 69) .<br>
> >> > 70->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 70)}<br>
> >> ><br>
> >> > The value seems incorrect every four numbits<br>
> >> > The exact value is near 1- 2.0e-45, and given the insufficient precision<br>
> >> > requested, the result should be rounded to 1.<br>
> >> > We can confirm with online high precision erf calculator like found for<br>
> >> > example with <a href="http://keisan.casio.com/exec/system/1180573449" rel="noreferrer" target="_blank">http://keisan.casio.com/exec/system/1180573449</a><br>
> >> ><br>
> >> > Or we can check with<br>
> >> > (10 asArbitraryPrecisionFloatNumBits: 200) erf (ArbitraryPrecisionFloat<br>
> >> > readFrom:<br>
> >> '0.999999999999999999999999999999999999999999997911512416237455'<br>
> >> > readStream numBits: 200)<br>
> >> ><br>
> >> > -----------------------------<br>
> >> ><br>
> >> > Note that when I change the LargeInteger digitLength to 32bits instead<br>
> >> of 8<br>
> >> > (my own VM branch) then the problem magically disappear.<br>
> >> ><br>
> >> > I plan to introduce these change in the VMMaker.oscog, so I can as well<br>
> >> > ignore the problem and proceed, however I much prefer to isolate it, and<br>
> >> > understand it for gain of knowledge :)<br>
> >> ><br>
> >> ><br>
> >> > 2016-03-27 4:27 GMT+02:00 Levente Uzonyi <<a href="mailto:leves@caesar.elte.hu">leves@caesar.elte.hu</a>>:<br>
> >> ><br>
> >> > ><br>
> >> > > Hi Nicolas,<br>
> >> > ><br>
> >> > > If you give us a bit more detailed instructions about what to test,<br>
> >> then<br>
> >> > > we might be able to help you get this sorted out.<br>
> >> > ><br>
> >> > > Levente<br>
> >> > ><br>
> >> > ><br>
> >><br>
> >><br>
> ><br>
<br>
</div></div></blockquote></div><br></div>