weak log calculations, possible to improve?

newbie nullzwei via gnumeric-list <[email protected]> Tue, 5 Jul 2022 19:33:56 +0200
Newsgroups gmane.comp.gnome.apps.gnumeric
Message-ID <trinity-965b6544-f8ff-49d8-8d84-b2bc582bab8d-1657042436945@3c-app-gmx-bap65> 
--===============7742339317074401511==
Content-Type: text/html; charset=UTF-8

<html><head></head><body><div style="font-family: Verdana;font-size: 12.0px;"><div style="font-family: Verdana;font-size: 12.0px;">
<div>&nbsp;<br/>
hello @all,&nbsp;<br/>
&nbsp;<br/>
IMHO it need not be that the result of &#39;=log( 8, 3 )&#39; is off by two ULP.&nbsp;<br/>
&nbsp;<br/>
( it is calculated to&nbsp;<br/>
1.892789260714371923910E00&nbsp; while &#39;3&#39; re-powered with that value becomes only&nbsp;<br/>
7.999999999999996447286E00&nbsp; , and even 3 powered to the next representable double&nbsp;<br/>
1.892789260714372145955E00&nbsp; does not! reach 8.0 but it needs the &#39;nextnext&#39; representable double&nbsp;<br/>
1.892789260714372368000E00&nbsp; as exponent to get 3 powered up to 8.0. )&nbsp;<br/>
&nbsp;<br/>
pls. check and correct if I&#39;m wrong ...&nbsp;<br/>
&nbsp;<br/>
( above may be one of few extreme cases? but there are more where gnumerics calculation is off by 1 ULP, e.g. log( 243, 3 ) results in &#39;4.999999999999999...&#39; despite the correct result &#39;5.0&#39; is! a representable double value. )&nbsp;<br/>
&nbsp;<br/>
( The imprecision is &#39;FP typical&#39; and - IMHO - results from calculating log( x, y ) as &#39;ln( x ) / ln( y )&#39;, ln( 8 ) and ln( 243 ) are undershot, ln( 3 ) is overshot, results are undershot. )&nbsp;<br/>
&nbsp;<br/>
We can&#39;t solve cases where we&#39;d need additional values between representable doubles, but we can try correction for cases which are uneccessary far off.&nbsp;<br/>
&nbsp;<br/>
As long as we don&#39;t find a better solution / calculation / libraries we can do the following in the gnumeric calculation:&nbsp;<br/>
&nbsp;<br/>
// edit b. - TESTING - 2022-07-03: correcting weak logs,&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; {&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; res = gnm_log (t) / gnm_log (base);&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if( gnm_abs( gnm_pow( base, nextafter( res, GNM_MAX ) ) - t ) &lt; gnm_abs( gnm_pow( base, res ) - t ) )&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; res = nextafter ( res, GNM_MAX );&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if( gnm_abs( gnm_pow( base, nextafter( res, GNM_MAX ) ) - t ) &lt; gnm_abs( gnm_pow( base, res ) - t ) )&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; res = nextafter ( res, GNM_MAX );&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if( gnm_abs( gnm_pow( base, nextafter( res, -GNM_MAX ) ) - t ) &lt; ( gnm_abs( gnm_pow( base, res ) ) - t ) )&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; res = nextafter ( res, -GNM_MAX );&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if( gnm_abs( gnm_pow( base, nextafter( res, -GNM_MAX ) ) - t ) &lt; ( gnm_abs( gnm_pow( base, res ) ) - t ) )&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; res = nextafter ( res, -GNM_MAX );&nbsp;<br/>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }&nbsp;<br/>
// edit b. - end edit,&nbsp;<br/>
&nbsp;<br/>
it is some overhead ... but IMHO tolerable, there is no footprint in timing the gnumeric test suite.&nbsp;<br/>
&nbsp;<br/>
With a simple test field of the integers from 2 .. 20 as base and argument, and gnumeric double &#39;standard&#39; we have a problem quota of about 61% where re-powering is imprecise, and can &#39;only&#39; reduce that to about 48% with above &#39;afterburner&#39;. That looks weak at a first glance, but it&#39;s eliminating all cases with &#39;big&#39; deviation where a better solution is available in doubles, thus it is an improvement ( from 15.7899&nbsp;average, 15.2028 min. &#39;quality&#39; to 15.8491&nbsp;average, 15.4337 min. quality ( &#39;quality&#39;: ~no. of correct digits ( in the re-powered value, I was too lazy to calculate 361 results with wolframalpha ) calculated similar to other tests by sthg. like &#39;-log10( abs( result - reference ) / abs( reference ) )&#39; ) )&nbsp;<
 br/>
&nbsp;<br/>
( &#39;long&#39; has ~45% problem quota in the above scope which can be reduced to ~33% with the patch, improving quality from 19.1354 avg. / 18.6639 min. to 19.1960 avg. / 18.7404 min., needs nextafter**l** instead of nextafter. )&nbsp;<br/>
&nbsp;<br/>
I&#39;m aware that such can&#39;t clean up FP-weaknesses in general, it&#39;s just one small step to avoid uneccessary imprecision.&nbsp;<br/>
I&#39;m aware that &#39;power&#39; is FP-math too and may have weaknesses ... WIP ... perhaps one will find improvements too ...&nbsp;<br/>
The patch above is most likely not well programmed, I do not! write here because I&#39;m good at programming, but because I&#39;m good at spotting weak points, even if I sometimes slip into mistakes ... thus: whoever can do, pls. improve ...&nbsp;<br/>
The proposal is based on a hint from John Denker for checking and correcting roundup / rounddown fails.&nbsp;<br/>
One of the checks I did against too big mistakes from my side: wolframalpha: &#39;log( 3, 8 )&#39; ( they use reverse operand ordering ) -&gt; 1.8927892607143723112985813430282825628987569203956412836119648315...&nbsp;</div>
</div></div></body></html>

--===============7742339317074401511==
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

_______________________________________________
gnumeric-list mailing list
[email protected]
https://mail.gnome.org/mailman/listinfo/gnumeric-list

--===============7742339317074401511==--