A numerical issue related to fp limitations

--===============2412736579754537316==
Content-Type: multipart/alternative; boundary="00000000000058bd8e0652e4b0b0"

--00000000000058bd8e0652e4b0b0
Content-Type: text/plain; charset="UTF-8"

Hi,

This is about an fp issue I stumbled upon while determining the order of
magnitude of some fp number. It occured to me while determining the
appropriate *prefix* of a unit, given some *numerical value* carrying the
unit.

The prefix to determine (denoted by an underscore '_') can arise in the
unit (based on e.g. 'm' for meter) in any arbitrary power. In the most
simple case, as in '_m', the power of the prefix is 1, but e.g. in '1 /
_m^2', the power is -2. I am explaining the problem here as it occurs with:

   - A prefix "u" (micro), carrying 1000 ** -2 as the coefficient;
   - A power of the prefix of -2;
   - A numerical value of 10 ** 12, for which the prefix shall be
   determined.

I am comparing the numerical value with the prefix coefficient, raised to
the prefix power. The power of the prefix coefficient is called the "prefix
value". For *negative* prefix powers, the correct prefix is the *smallest*
one, where it holds:

numerical value >= prefix value

Here, for a prefix power of -2, and the prefix 'u' (1000 ** -2 = 10 ** -6),
mathematically the prefix value is (10 ** -6) ** -2 = 10 ** 12, which is
just the numerical value, so the abovementioned comparison should be
fulfilled as numerical value == prefix value.

Technically I am obtaining:

   1. The numerical value (10 ** 12) is 1000000000000.
   2. The prefix value however, (10 ** -6) ** -2, is 10000000000.0001.

Here, where the numerical value and the prefix value should coincide,
the numerical
value is smaller than the prefix value:

numerical value - prefix value = -0.000122...

so the 'u' prefix is mistakenly rejected, as "numerical value >= prefix
value" is not fulfilled.

I kind-of-fixed it utilising numpy.nextafter:

stepnumerical = numericalvalue - numpy.nextafter(numericalvalue, 0)
stepprefix = prefixvalue - numpy.nextafter(prefixvalue, 0)

and replacing the condition "numerical value >= prefix value" by:

numerical value - prefix value >= -(stepnumerical + stepprefix)

In the singleton example discussed here, this remedies the flaw.

It *seems*, that the error of the final figures "numerical value" and "prefix
value" is not larger than the error resulting from their limited number of
fp digits. However, I am unsure if this can be used *always* as the
decisive threshold. I could speculate that the fp uncertainty of
intermediate numbers might carry through, so that the overall error margin
is even larger than the one implied by the limitedness of the *final*
number of binary digits. This suspection applies most prominently to
the "prefix
value".

One idea would be to use a multiple of the numpy.nextafter result. However,
2x (2 bits), 4x (3 bits), ...? I think that the approach using
numpy.nextafter is valid, but my use is currently a little heuristic.

Any idea, comment or pointer would be welcome very much. I'd like to sort
this problem thoroughly.

Best,
Friedrich

--00000000000058bd8e0652e4b0b0
Content-Type: text/html; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable

<div dir=3D"ltr"><div>Hi,</div><div><br></div><div>This is about an fp issu=
e I stumbled upon while determining the order of magnitude of some fp numbe=
r. It occured to me while determining the appropriate=C2=A0<i>prefix</i>=C2=
=A0of a unit, given some <i>numerical value</i> carrying the unit.</div><di=
v><br></div><div>The prefix to determine (denoted by an underscore <span st=
yle=3D"font-family:monospace">&#39;_&#39;</span>) can arise in the unit (ba=
sed on e.g. <span style=3D"font-family:monospace">&#39;m&#39;</span> for me=
ter) in any arbitrary power. In the most simple case,=C2=A0as in <span styl=
e=3D"font-family:monospace">&#39;_m&#39;</span>, the power of the prefix is=
 1, but e.g. in <span style=3D"font-family:monospace">&#39;1 / _m^2&#39;</s=
pan>, the power is <span style=3D"font-family:monospace">-2</span>. I am ex=
plaining the problem here as it occurs with:</div><div><ul><li>A prefix &qu=
ot;u&quot; (micro), carrying <span style=3D"font-family:monospace">1000 ** =
-2</span> as the coefficient;</li><li>A power of the prefix of <span style=
=3D"font-family:monospace">-2</span>;</li><li>A numerical value of <span st=
yle=3D"font-family:monospace">10 ** 12</span>, for which the prefix shall b=
e determined.</li></ul><div>I am comparing the numerical value with the pre=
fix coefficient, raised to the prefix power. The power of the prefix coeffi=
cient is called the &quot;<span style=3D"font-family:monospace">prefix valu=
e</span>&quot;. For <i>negative</i> prefix powers, the correct prefix is th=
e <i>smallest</i> one, where it holds:<br><br></div><div style=3D"margin-le=
ft:40px"><span style=3D"font-family:monospace">numerical value &gt;=3D pref=
ix value</span></div><br></div><div>Here, for a prefix power of <span style=
=3D"font-family:monospace">-2</span>, and the prefix <span style=3D"font-fa=
mily:monospace">&#39;u&#39;</span> <span style=3D"font-family:monospace">(1=
000 ** -2 =3D 10 ** -6</span>), mathematically the <span style=3D"font-fami=
ly:monospace">prefix value</span> is <span style=3D"font-family:monospace">=
(10 ** -6) ** -2 =3D 10 ** 12</span>, which is just the <span style=3D"font=
-family:monospace">numerical value</span>, so the abovementioned comparison=
 should be fulfilled as <span style=3D"font-family:monospace">numerical val=
ue =3D=3D prefix value</span><font face=3D"arial,sans-serif">.</font></div>=
<div><font face=3D"arial,sans-serif"><br></font></div><div><font face=3D"ar=
ial,sans-serif">Technically I am obtaining:</font></div><div><ol><li><font =
face=3D"arial,sans-serif">The </font><span style=3D"font-family:monospace">=
numerical value</span><font face=3D"arial,sans-serif"> (</font><span style=
=3D"font-family:monospace">10 ** 12</span><font face=3D"arial,sans-serif">)=
 is </font><span style=3D"font-family:monospace">1000000000000<font face=3D=
"arial,sans-serif">.</font></span></li><li><span style=3D"font-family:monos=
pace"><font face=3D"arial,sans-serif">The </font>prefix value<font face=3D"=
arial,sans-serif"> however, </font>(10 ** -6) ** -2<font face=3D"arial,sans=
-serif">, is </font>10000000000.0001<font face=3D"arial,sans-serif">.</font=
></span></li></ol><div><font face=3D"arial,sans-serif">Here, where the </fo=
nt><span style=3D"font-family:monospace">numerical value</span><font face=
=3D"arial,sans-serif"> and the </font><span style=3D"font-family:monospace"=
>prefix value</span><font face=3D"arial,sans-serif"> should coincide, the <=
/font><span style=3D"font-family:monospace">numerical value</span><font fac=
e=3D"arial,sans-serif"> is smaller than the </font><span style=3D"font-fami=
ly:monospace">prefix value</span><font face=3D"arial,sans-serif">:</font></=
div><div><font face=3D"arial,sans-serif"><br></font></div><div style=3D"mar=
gin-left:40px"><span style=3D"font-family:monospace">numerical value - pref=
ix value =3D -0.000122...</span></div><font face=3D"arial,sans-serif"><br><=
/font></div><div><font face=3D"arial,sans-serif">so the </font><span style=
=3D"font-family:monospace">&#39;u&#39;</span><font face=3D"arial,sans-serif=
"> prefix is mistakenly rejected, as &quot;</font><span style=3D"font-famil=
y:monospace">numerical value &gt;=3D prefix value</span><font face=3D"arial=
,sans-serif">&quot; is not fulfilled.</font></div><div><font face=3D"arial,=
sans-serif"><br></font></div><div><font face=3D"arial,sans-serif">I kind-of=
-fixed it utilising </font><span style=3D"font-family:monospace">numpy.next=
after</span><font face=3D"arial,sans-serif">:</font></div><div><font face=
=3D"arial,sans-serif"><br></font></div><div style=3D"margin-left:40px"><spa=
n style=3D"font-family:monospace">stepnumerical =3D numericalvalue - numpy.=
nextafter(numericalvalue, 0)</span></div><div style=3D"margin-left:40px"><s=
pan style=3D"font-family:monospace">stepprefix =3D prefixvalue - numpy.next=
after(prefixvalue, 0)</span></div><div><br></div><div><font face=3D"arial,s=
ans-serif">and replacing the condition &quot;</font><span style=3D"font-fam=
ily:monospace">numerical value &gt;=3D prefix value</span><font face=3D"ari=
al,sans-serif">&quot; by:</font></div><div><font face=3D"arial,sans-serif">=
<br></font></div><div style=3D"margin-left:40px"><span style=3D"font-family=
:monospace">numerical value - prefix value &gt;=3D -(stepnumerical=C2=A0+ s=
tepprefix)</span></div><div><br></div><div>In the singleton example discuss=
ed here, this remedies the flaw.</div><div><br></div><div>It=C2=A0<i>seems<=
/i>, that the error of the final figures &quot;<span style=3D"font-family:m=
onospace">numerical value</span>&quot; and &quot;<span style=3D"font-family=
:monospace">prefix value</span>&quot; is not larger than the error resultin=
g from their limited number of fp digits. However, I am unsure if this can =
be used=C2=A0<i>always</i>=C2=A0as the decisive threshold. I could speculat=
e that the fp uncertainty of intermediate numbers might carry through, so t=
hat the overall error margin is even larger than the one implied by the lim=
itedness of the <i>final</i> number of binary digits. This suspection appli=
es most prominently to the <span style=3D"font-family:monospace">&quot;pref=
ix value</span>&quot;.</div><div><br></div><div>One idea would be to use a =
multiple of the <span style=3D"font-family:monospace">numpy.nextafter</span=
> result. However, 2x (2 bits), 4x (3 bits), ...?
 I think that the approach using <span style=3D"font-family:monospace">nump=
y.nextafter</span> is valid, but my use is currently a little heuristic.</d=
iv><div><br></div><div>Any idea, comment or pointer would be welcome very m=
uch. I&#39;d like to sort this problem thoroughly.</div><div><br></div><div=
>Best,</div><div>Friedrich</div></div>

--00000000000058bd8e0652e4b0b0--

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

_______________________________________________
NumPy-Discussion mailing list -- [email protected]
To unsubscribe send an email to [email protected]
https://mail.python.org/mailman3//lists/numpy-discussion.python.org
Member address: [email protected]

--===============2412736579754537316==--