<div dir="ltr"><div dir="ltr"><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, Aug 6, 2020 at 8:38 PM Martin Budaj <<a href="mailto:m.budaj@gmail.com">m.budaj@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hi,<br>
<br>
there is a bug in arctime calculation when using numbersystems decimal<br>
and binary (double and scaled work fine). See the following example:<br>
<br>
def test(expr P) =<br>
pos := 0;<br>
len := arclength P;<br>
d := 10bp;<br>
forever:<br>
arc_t := arctime pos of P;<br>
draw ((point arc_t of P) --<br>
(point arc_t of P) + 1mm * unitvector(direction arc_t of P rotated 90))<br>
withcolor (1-pos/len);<br>
message decimal pos & " " & decimal arc_t;<br>
pos := pos + d;<br>
exitif pos > len;<br>
endfor;<br>
draw P;<br>
enddef;<br>
<br>
beginfig(1);<br>
test((0,0)--(100,0));<br>
endfig;<br>
<br>
end<br>
<br>
As the test is run on a simple line segment, arctime should increase<br>
linearly from 0 to 1. This is the case for scaled and double<br>
numbersystems:<br>
<br>
pos arc_t<br>
0 0<br>
10 0.09999<br>
20 0.2<br>
30 0.29999<br>
40 0.4<br>
50 0.5<br>
60 0.59999<br>
70 0.7<br>
80 0.79999<br>
90 0.9<br>
100 1<br>
<br>
But for decimal and binary, arctime is not only non-linear, it's even<br>
non-monotonous:<br>
<br>
pos arc_t<br>
0 0<br>
10 0.218357288800015125792849022379372<br>
20 0.338083509141544805383550503071265<br>
30 0.45780972290631604654943093962971<br>
40 0.038767960662661419366106762617287<br>
50 0.098631047814785205469201499670517<br>
60 0.718357199283058404647430603398735<br>
70 0.838083397702050143830649542333004<br>
80 0.957809611466813022214145595635334<br>
90 0.538767915904170514622883519402551<br>
100 0.598631024978811668703918040345015<br>
<br>
Could you please look into this issue?<br>
</blockquote></div><div><br></div><div>The bug is in mp_solve_rising_cubic:</div><div><br></div><div> @ Here is the |solve_rising_cubic| routine that finds the time~$t$ when<br> $$ B(0, a, a+b, a+b+c; t) = x. $$<br>This routine is based on |crossing_point| but is simplified by the<br>assumptions that $B(a,b,c;t)\ge0$ for $0\le t\le1$ and that |0<=x<=a+b+c|.<br>If rounding error causes this condition to be violated slightly, we just ignore<br>it and proceed with binary search. This finds a time when the function value<br>reaches |x| and the slope is positive.<br></div><div><br></div><div><br></div><div>Both binary and decimal uses and epsilon = 1E-52</div><div>buth double has epsilon = 2^-52 and this matches with the binary search.</div>Locally I have set epsilon = 2^-52 for both binary and decimal <div>and the results are now:<div><br><div># mpost --numbersystem=double <a href="http://test-b.mp">test-b.mp</a> <br>0 0<br>10 0.10000000335693371<br>20 0.20000000854492228<br>30 0.3000000082397456<br>40 0.39999999511718842<br>50 0.49999996185304463<br>60 0.59999992858888307<br>70 0.69999991546631324<br>80 0.79999991516112967<br>90 0.89999992034911669<br>100 0.99999992370605462 <br><br></div><div><br></div><div># mpost --numbersystem=binary <a href="http://test-b.mp">test-b.mp</a> <br>0 0<br>10 0.10000000335693370878686891956022009<br>20 0.20000000854492216983970820365357213<br>30 0.30000000823974559693851915653795004<br>40 0.39999999511718853373309912058175541<br>50 0.49999996185304473872434982695267536<br>60 0.5999999285888828470802991432719864<br>70 0.69999991546631312733239838053123094<br>80 0.79999991516112967104845665744505823<br>90 0.89999992034911668881136392883490771<br>100 0.99999992370605461644572642398998141 <br><br></div><div># mpost --numbersystem=decimal <a href="http://test-b.mp">test-b.mp</a> <br>0 0<br>10 0.10000000335693370878686891956022<br>20 0.200000008544922169839708203653571<br>30 0.30000000823974559693851915653795<br>40 0.399999995117188533733099120581755<br>50 0.499999961853044738724349826952676<br>60 0.599999928588882847080299143271985<br>70 0.699999915466313127332398380531229<br>80 0.799999915161129671048456657445057<br>90 0.899999920349116688811363928834907<br>100 0.99999992370605461644572642398998 <br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div>--</div><div>luigi<br></div></div></div></div>