Re: A numerical issue related to fp limitations

--===============3630931622161281194==
Content-Type: multipart/alternative; boundary="0000000000003009070652f69bae"

--0000000000003009070652f69bae
Content-Type: text/plain; charset="UTF-8"

Am Do., 28. Mai 2026 um 20:13 Uhr schrieb Friedrich Romstedt <
[email protected]>:

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

So I discussed the topic with ChatGPT, and we found a viable solution.

In essence, the problem boils down two the following situation:

   1. Let there be a basis *b*, either 1000 for SI prefixes, 1024 for IEC
   prefixes, or 10 for the remaining SI prefixes corresponding to 10^{-2,-1,
   1, 2}.
   2. Let there further be a *prefix exponent, p*, denoting the exponent of
   *b* for a given prefix.
   3. An finally, let there be a *unit exponent* *u*, indicating the
   exponent the entire prefix has in the unit (e.g. -2 for a unit of the
   form '1 / _m^2', where the underscore '_' indicates the position of the
   prefix).

Then two numbers can be deduced:

   1. The *numerical value* N = 1000 ** (p * u). This is the *numerical
   value* handed over by the user, the value for which the best-fitting
   prefix shall be determined. The name *numerical value* is a little
   misleading, it refers to the *numerical value* combined with the *unit* to
   the *quantity*.
   2. The overall *prefix value*, P = (1000 ** p) ** u. This is the number
   which is compared to N.

Mathematically, both are equal, but technically P and N can be a little
different, e.g. in the SI case for p = -2 and u = -2. These corner cases I
am interested in here.

To decide about the applicability of a prefix, P is compared to N as N >= P.
In the example p = -2 & u = -2, however N is slighly smaller than P.

ChatGPT proposed to iterate about all possible combinations of p and u. For
each pair, 1000.0 ** (p * u) can be compared to (1000.0 ** p) ** u. I am
comparing the difference of P and N, P - N, to the sum S of the two ULPs
deduced for P and N (like X - numpy.nextafter(X, 0)). More specifically, I
am calculating an "error coefficient" r = (P - N) / S. I did so both for SI
as well as for IEC prefixes.

The unit exponent u is always in [-3.0, -2.5, ..., 2.0. 2.5, 3.0].

For SI, p is in [-10.0, -9.0, ..., 9.0, 10.0]. Here, the maximum r is 1.0
(seldom), more often 0.5, most often 0; sometimes -0.5, rarely -1.0, one
time -1.5. S is the sum of the two ULPs, and apparently the ULPs are equal
for P and N, leading to the results for r being an integer multiple of 0.5.

For 10-based SI units (d, c, ...), the results are similar, but
bounded by [-0.5,
0.5].

With less regular exponents, e.g. with p = 3.0 and u = -2.1, r reaches
8.00. With even more irregular values, I am reaching r = 38.50 (p approx.
-8.95, u = 3.0).

For the IEC prefixes, I obtained no deviations between N and P at all
(p in [0.0,
1.0, ..., 7.0, 8.0]).

So my issue is sorted. I will install 2 * S as the safety margin when
comparing a *numerical value* to an *overall prefix value*, and this safety
margin will suffice in all regular situations tested.

Ideally, I would put this report on my web site and just link to it here,
but this web site is not yet there, so I am posting my results here, to
provide the solution of my problem.

Friedrich

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

<div dir=3D"ltr"><div dir=3D"ltr"><div dir=3D"ltr"><div class=3D"gmail_quot=
e"><div dir=3D"ltr" class=3D"gmail_attr">Am Do., 28. Mai 2026 um 20:13=C2=
=A0Uhr schrieb Friedrich Romstedt &lt;<a href=3D"mailto:friedrichromstedt@g=
mail.com" target=3D"_blank">[email protected]</a>&gt;:<br></div><=
blockquote class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8ex;border-l=
eft:1px solid rgb(204,204,204);padding-left:1ex"><div dir=3D"ltr"><div>Any =
idea, comment or pointer would be welcome very much. I&#39;d like to sort t=
his problem thoroughly.</div></div></blockquote><div><br></div><div>So I di=
scussed the topic with ChatGPT, and we found a viable solution.</div><div><=
br></div><div>In essence, the problem boils down two the following situatio=
n:</div><div><ol><li>Let there be a basis=C2=A0<i>b</i>, either 1000 for SI=
 prefixes, 1024 for IEC prefixes, or 10 for the remaining SI prefixes corre=
sponding to <span style=3D"font-family:monospace">10^{-2,-1, 1, 2}</span>.<=
/li><li>Let there further be a=C2=A0<i>prefix exponent, p</i>, denoting the=
 exponent of=C2=A0<i>b</i>=C2=A0for a given prefix.</li><li>An finally, let=
 there be a=C2=A0<i>unit exponent</i>=C2=A0<i>u</i>, indicating the exponen=
t the entire prefix has in the unit (e.g. <span style=3D"font-family:monosp=
ace">-2</span> for a unit of the form <span style=3D"font-family:monospace"=
>&#39;1 / _m^2&#39;</span>, where the underscore &#39;_&#39; indicates the =
position of the prefix).</li></ol><div>Then two numbers can be deduced:</di=
v><div><ol><li>The <i>numerical=C2=A0value</i>=C2=A0<span style=3D"font-fam=
ily:monospace">N =3D 1000 ** (p * u)</span>. This is the=C2=A0<i>numerical =
value</i>=C2=A0handed over by the user, the value for which the best-fittin=
g prefix shall be determined. The name=C2=A0<i>numerical value</i>=C2=A0is =
a little misleading, it refers to the=C2=A0<i>numerical value</i>=C2=A0comb=
ined with the=C2=A0<i>unit</i>=C2=A0to the=C2=A0<i>quantity</i>.</li><li>Th=
e overall=C2=A0<i>prefix value</i>, <span style=3D"font-family:monospace">P=
 =3D (1000 ** p) ** u</span>. This is the number which is compared to <span=
 style=3D"font-family:monospace">N</span>.</li></ol><div>Mathematically, bo=
th are equal, but technically P and N can be a little different, e.g. in th=
e SI case for <span style=3D"font-family:monospace">p =3D -2</span> and <sp=
an style=3D"font-family:monospace">u =3D -2</span>. These corner cases I am=
 interested in here.</div><div><br></div><div>To decide about the applicabi=
lity of a prefix, P is compared to N as <span style=3D"font-family:monospac=
e">N &gt;=3D P</span>. In the example <span style=3D"font-family:monospace"=
>p =3D -2</span> &amp; <span style=3D"font-family:monospace">u =3D -2</span=
>, however N is slighly smaller than P.</div><div><br></div><div>ChatGPT pr=
oposed to iterate about all possible combinations of p and u. For each pair=
, <span style=3D"font-family:monospace">1000.0 ** (p * u)</span> can be com=
pared to <span style=3D"font-family:monospace">(1000.0 ** p) ** u</span>. I=
 am comparing the difference of P and N, <span style=3D"font-family:monospa=
ce">P - N</span>, to the sum S of the two ULPs deduced for P and N (like <s=
pan style=3D"font-family:monospace">X - numpy.nextafter(X, 0)</span>). More=
 specifically, I am calculating an &quot;error coefficient&quot;=C2=A0<span=
 style=3D"font-family:monospace">r =3D (P - N) / S</span>. I did so both fo=
r SI as well as for IEC prefixes.</div><div><br></div><div>The unit exponen=
t <span style=3D"font-family:monospace">u</span> is always in <span style=
=3D"font-family:monospace">[-3.0, -2.5, ..., 2.0. 2.5, 3.0]</span>.</div><d=
iv><br></div><div>For SI, <span style=3D"font-family:monospace">p</span> is=
 in <span style=3D"font-family:monospace">[-10.0, -9.0, ..., 9.0, 10.0]</sp=
an>. Here, the maximum r is 1.0 (seldom), more often 0.5, most often 0; som=
etimes -0.5, rarely -1.0, one time -1.5. S is the sum of the two ULPs, and =
apparently the ULPs are equal for P and N, leading to the results for <span=
 style=3D"font-family:monospace">r</span> being an integer multiple of 0.5.=
</div><div><br></div><div>For 10-based SI units (d, c, ...), the results ar=
e similar, but bounded by <span style=3D"font-family:monospace">[-0.5, 0.5]=
</span>.</div><div><br></div><div>With less regular exponents, e.g. with <s=
pan style=3D"font-family:monospace">p =3D 3.0</span> and <span style=3D"fon=
t-family:monospace">u =3D -2.1</span>, r reaches 8.00. With even more irreg=
ular values, I am reaching <span style=3D"font-family:monospace">r =3D 38.5=
0</span> (<span style=3D"font-family:monospace">p</span> approx. <span styl=
e=3D"font-family:monospace">-8.95</span>, <span style=3D"font-family:monosp=
ace">u =3D 3.0</span>).</div><div><br></div><div>For the IEC prefixes, I ob=
tained no deviations between N and P at all (p in <span style=3D"font-famil=
y:monospace">[0.0, 1.0, ..., 7.0, 8.0]</span>).</div><div><br></div><div>So=
 my issue is sorted. I will install <span style=3D"font-family:monospace">2=
 * S</span> as the safety margin when comparing a <i>numerical value</i> to=
 an <i>overall prefix value</i>, and this safety margin will suffice in all=
 regular situations tested.</div><div><br></div><div>Ideally, I would put t=
his report on my web site and just link to it here, but this web site is no=
t yet there, so I am posting my results here, to provide the solution of my=
 problem.</div><div><br></div><div>Friedrich</div></div></div></div></div>
</div>
</div>

--0000000000003009070652f69bae--

--===============3630931622161281194==
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]

--===============3630931622161281194==--