<div dir="ltr"><div><div><div><div><div><div><div>OK, I think I get it.<br></div>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></div>What happens here is that hi may become negative.<br></div>But since I declared it unsigned, t<hi and cond := hi >= dh remains false quite a long time and the quotient q is decremented way too many times...<br><br></div></div>Here we divide (r1,r2,r3,...) by (dh,dnh,...)<br>we are trying to guess q based on 2 limbs<br></div>(r1,r2) = q*dh + t that is (r1,r2-t) = q*dh</div><div>(hi,lo) =q*dnh<br><br></div><div>Thus q*(dh,dnh)=(r1,r2-t+hi,lo)<br><br></div><div>If t<hi then q*(dh,dnh) is greater than (r1,r2) and we have an overestimate of q<br></div><div>Same if t=hi and lo<r3.<br><br></div><div>Here we have the tricky case when t=0, r3=0 and lo>0, and hi=dnh+1<br></div><div>In the first loop, h gets decremented, thus hi=dnh.<br></div><div>Then it goes to 0, and gets decremented a second time...<br></div><div><br>The same algorithm is used with 32bits limbs, so it
may also happen in my own VM branch. I believe it must be even rarest an
event though (2^32 possible values for hi makes hi=dnh+1 less probable than with 2^8 combinations).<br><br></div><div>I will fix it tonight or tomorrow.<br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2016-03-27 19:03 GMT+02:00 Nicolas Cellier <span dir="ltr"><<a href="mailto:nicolas.cellier.aka.nice@gmail.com" target="_blank">nicolas.cellier.aka.nice@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div>Thanks David.<br><br></div>By running the snippet, instrumenting erf, printing each and every intermediate value to a text file, and diffing:<br>{<br>(10 asArbitraryPrecisionFloatNumBits: 33) erf.<br>(10 asArbitraryPrecisionFloatNumBits: 34) erf.<br>}<br><br></div>I reduced the failing test to:<br><br> a := 1598335257761788022467377781654101148543282249044465229239888363328190330275719997501596724768507889233831388734160190922469363547795602076820594918.<br>b := 49612.<br>self assert: a - ((a quo: b)*b) < b<br><br></div><div>So it must be some nasty sign promotion (or lack of?) in division.<br></div><div>Normally, this could be simulated, but with signedness, I'm never sure, maybe my next step will be a breakpoint in xcode/gdb...<br><br><br>Curiously, (a rem: b) is correct !<br></div></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">2016-03-27 17:36 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>
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>
<div><div><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 VM.<br>
> <rant> It's just that for MacOSX, there is no recent VM available, and 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 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 job).<br>
><br>
> Now the job is to produce/isolate an elementary failing LargeInteger 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 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' 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' 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: '0.999999999999999999999999999999999999999999997911512416237455'<br>
> readStream numBits: 200)<br>
><br>
> -----------------------------<br>
><br>
> Note that when I change the LargeInteger digitLength to 32bits instead 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" target="_blank">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, then<br>
> > we might be able to help you get this sorted out.<br>
> ><br>
> > Levente<br>
> ><br>
> ><br>
<br>
</div></div></blockquote></div><br></div>
</div></div></blockquote></div><br></div>