Re: bug in arctime calculation

On Thu, Aug 6, 2020 at 8:38 PM Martin Budaj <[email protected]> wrote:

> Hi,
>
> there is a bug in arctime calculation when using numbersystems decimal
> and binary (double and scaled work fine). See the following example:
>
> def test(expr P) =
>   pos := 0;
>   len := arclength P;
>   d := 10bp;
>   forever:
>     arc_t := arctime pos of P;
>     draw ((point arc_t of P) --
>          (point arc_t of P) + 1mm * unitvector(direction arc_t of P
> rotated 90))
>          withcolor (1-pos/len);
>     message decimal pos & " " & decimal arc_t;
>     pos := pos + d;
>     exitif pos > len;
>   endfor;
>   draw P;
> enddef;
>
> beginfig(1);
>   test((0,0)--(100,0));
> endfig;
>
> end
>
> As the test is run on a simple line segment, arctime should increase
> linearly from 0 to 1. This is the case for scaled and double
> numbersystems:
>
> pos arc_t
> 0 0
> 10 0.09999
> 20 0.2
> 30 0.29999
> 40 0.4
> 50 0.5
> 60 0.59999
> 70 0.7
> 80 0.79999
> 90 0.9
> 100 1
>
> But for decimal and binary, arctime is not only non-linear, it's even
> non-monotonous:
>
> pos arc_t
> 0 0
> 10 0.218357288800015125792849022379372
> 20 0.338083509141544805383550503071265
> 30 0.45780972290631604654943093962971
> 40 0.038767960662661419366106762617287
> 50 0.098631047814785205469201499670517
> 60 0.718357199283058404647430603398735
> 70 0.838083397702050143830649542333004
> 80 0.957809611466813022214145595635334
> 90 0.538767915904170514622883519402551
> 100 0.598631024978811668703918040345015
>
> Could you please look into this issue?
>

The bug is in mp_solve_rising_cubic:

 @ Here is the |solve_rising_cubic| routine that finds the time~$t$ when
 $$ B(0, a, a+b, a+b+c; t) = x. $$
This routine is based on |crossing_point| but is simplified by the
assumptions that $B(a,b,c;t)\ge0$ for $0\le t\le1$ and that |0<=x<=a+b+c|.
If rounding error causes this condition to be violated slightly, we just
ignore
it and proceed with binary search.  This finds a time when the function
value
reaches |x| and the slope is positive.


Both binary and decimal uses and epsilon = 1E-52
buth double has epsilon = 2^-52 and this matches with the binary search.
Locally I have set  epsilon = 2^-52 for both binary and decimal
and the results  are now:

# mpost --numbersystem=double test-b.mp
0 0
10 0.10000000335693371
20 0.20000000854492228
30 0.3000000082397456
40 0.39999999511718842
50 0.49999996185304463
60 0.59999992858888307
70 0.69999991546631324
80 0.79999991516112967
90 0.89999992034911669
100 0.99999992370605462


# mpost --numbersystem=binary test-b.mp
0 0
10 0.10000000335693370878686891956022009
20 0.20000000854492216983970820365357213
30 0.30000000823974559693851915653795004
40 0.39999999511718853373309912058175541
50 0.49999996185304473872434982695267536
60 0.5999999285888828470802991432719864
70 0.69999991546631312733239838053123094
80 0.79999991516112967104845665744505823
90 0.89999992034911668881136392883490771
100 0.99999992370605461644572642398998141

# mpost --numbersystem=decimal test-b.mp
0 0
10 0.10000000335693370878686891956022
20 0.200000008544922169839708203653571
30 0.30000000823974559693851915653795
40 0.399999995117188533733099120581755
50 0.499999961853044738724349826952676
60 0.599999928588882847080299143271985
70 0.699999915466313127332398380531229
80 0.799999915161129671048456657445057
90 0.899999920349116688811363928834907
100 0.99999992370605461644572642398998





--
luigi

--
http://tug.org/metapost/