OK, I think I get it.
The problem lies in:

                "Correct overestimate of q.               
                Max of 2 iterations through loop -- see Knuth vol. 2"
                j < 3
                    ifTrue: [r3 := 0]
                    ifFalse: [r3 := pRem at: j - 3].
               
                [(t < hi
                    or: [t = hi and: [r3 < lo]])
                    ifTrue:
                        ["i.e. (t,r3) < (hi,lo)"
                        q := q - 1.
                        lo < dnh
                            ifTrue:
                                [hi := hi - 1.
                                lo := lo + 16r100 - dnh]
                            ifFalse:
                                [lo := lo - dnh].
                        cond := hi >= dh]
                    ifFalse: [cond := false].
                cond]
                    whileTrue: [hi := hi - dh]].

What happens here is that hi may become negative.
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...

Here we divide (r1,r2,r3,...) by (dh,dnh,...)
we are trying to guess q based on 2 limbs
(r1,r2) = q*dh + t    that is   (r1,r2-t) = q*dh
(hi,lo) =q*dnh

Thus q*(dh,dnh)=(r1,r2-t+hi,lo)

If t<hi then q*(dh,dnh) is greater than (r1,r2) and we have an overestimate of q
Same if t=hi and lo<r3.

Here we have the tricky case when t=0, r3=0 and lo>0, and hi=dnh+1
In the first loop, h gets decremented, thus hi=dnh.
Then it goes to 0, and gets decremented a second time...

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).

I will fix it tonight or tomorrow.


2016-03-27 19:03 GMT+02:00 Nicolas Cellier <nicolas.cellier.aka.nice@gmail.com>:
Thanks David.

By running the snippet, instrumenting erf, printing each and every intermediate value to a text file, and diffing:
{
(10 asArbitraryPrecisionFloatNumBits: 33) erf.
(10 asArbitraryPrecisionFloatNumBits: 34) erf.
}

I reduced the failing test to:

 a := 1598335257761788022467377781654101148543282249044465229239888363328190330275719997501596724768507889233831388734160190922469363547795602076820594918.
b := 49612.
self assert: a - ((a quo: b)*b) < b

So it must be some nasty sign promotion (or lack of?) in division.
Normally, this could be simulated, but with signedness, I'm never sure, maybe my next step will be a breakpoint in xcode/gdb...


Curiously, (a rem: b) is correct !

2016-03-27 17:36 GMT+02:00 David T. Lewis <lewis@mail.msen.com>:

I confirm the same results for your ArbitraryPrecisionFloat example below,
testing on Linux with a 64 bit interpreter VM and 32 bit image.

  Image
  -----
  /home/lewis/squeak/Squeak4.6/squeak.37.image
  Squeak4.6
  latest update: #15723
  Current Change Set: Unnamed2
  Image format 6504 (32 bit)

  Virtual Machine
  ---------------
  /usr/local/lib/squeak/4.15.4-3634/squeakvm
  Squeak4.5 of 10 December 2015 [latest update: #1195]
  Unix built on Mar 25 2016 13:27:51 Compiler: 4.9.2
  platform sources revision 3634
  VMMaker versionString 4.15.4

Dave


On Sun, Mar 27, 2016 at 03:25:58PM +0200, Nicolas Cellier wrote:
>
> Hi Levente, thank you.
> Here is a follow up:
>
> I was able to reproduce the same symptom on the classical interpreter VM.
> <rant> It's just that for MacOSX, there is no recent VM available, and that
> the MacOS build is somehow abandonware, so it took much more time than
> necessary :( </rant>
>
> So the changes that have been integrated on Interpreter Vm in july 2014 and
> that I thought correct were not.
> It's certainly related to a LargeInteger plugin primitive.
> Those are triggered when nbits of magnitude > 64.
> (otherwise, <= 64 bits, the LargeInteger numbered primitives do their job).
>
> Now the job is to produce/isolate an elementary failing LargeInteger test.
> But you need to generate the VM from VMMaker.oscog head for COG/STACK
> V3/SPUR 32/64 bits
> or use a sufficiently recent interpreter VM (> july 2014) to trigger such
> failure.
>
>
> --------------------------------------------
>
> Here is my starting point:
>
> If you load ArbitraryPrecisionFloat and Tests from
> http://www.squeaksource.com/ArbitraryPrecisionFl
> then you'll see this strange pattern:
>
> (53 to: 70) collect: [:i |
>     i -> (10 asArbitraryPrecisionFloatNumBits: i) erf ]
>
>  {53->(ArbitraryPrecisionFloat readFrom: '0.9999998981732067' readStream
> numBits: 53) .
>  54->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 54) .
>  55->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 55) .
>  56->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 56) .
>  57->(ArbitraryPrecisionFloat readFrom: '0.99999989817320666' readStream
> numBits: 57) .
>  58->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 58) .
>  59->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 59) .
>  60->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 60) .
>  61->(ArbitraryPrecisionFloat readFrom: '0.9999998981732066655' readStream
> numBits: 61) .
>  62->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 62) .
>  63->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 63) .
>  64->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 64) .
>  65->(ArbitraryPrecisionFloat readFrom: '0.99999989817320666564' readStream
> numBits: 65) .
>  66->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 66) .
>  67->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 67) .
>  68->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 68) .
>  69->(ArbitraryPrecisionFloat readFrom: '0.999999898173206665646'
> readStream numBits: 69) .
>  70->(ArbitraryPrecisionFloat readFrom: '1.0' readStream numBits: 70)}
>
> The value seems incorrect every four numbits
> The exact value is near 1- 2.0e-45, and given the insufficient precision
> requested, the result should be rounded to 1.
> We can confirm with online high precision erf calculator like found for
> example with http://keisan.casio.com/exec/system/1180573449
>
> Or we can check with
> (10 asArbitraryPrecisionFloatNumBits: 200) erf (ArbitraryPrecisionFloat
> readFrom: '0.999999999999999999999999999999999999999999997911512416237455'
> readStream numBits: 200)
>
> -----------------------------
>
> Note that when I change the LargeInteger digitLength to 32bits instead of 8
> (my own VM branch) then the problem magically disappear.
>
> I plan to introduce these change in the VMMaker.oscog, so I can as well
> ignore the problem and proceed, however I much prefer to isolate it, and
> understand it for gain of knowledge :)
>
>
> 2016-03-27 4:27 GMT+02:00 Levente Uzonyi <leves@caesar.elte.hu>:
>
> >
> > Hi Nicolas,
> >
> > If you give us a bit more detailed instructions about what to test, then
> > we might be able to help you get this sorted out.
> >
> > Levente
> >
> >