# Inaccuracy in floating number computation with fpu library

So I am trying to do some a bit involved computations with LaTeX, and it kept spitting out a nonsense answer. I am trying to compute number of layers you could cover a ball with, given some conditions, and LaTeX keeps giving me a negative answer! After pulling out my hair for hours, I was able to track down the error, which is shown in the MWE below

\documentclass[border=1mm]{article}
\usepackage[utf8]{inputenc}

\usepackage{mathtools}
\usepackage{pgfplots}

\begin{document}

\pgfmathsetmacro{\coinHeightM}{1.7 / 1000}

\pgfkeys{/pgf/fpu, /pgf/fpu/output format=fixed}

\pgfmathsetmacro{\coinsTotalHeight}{3.27*10^17}

\pgfmathsetmacro{\R}{
}

\pgfmathsetmacro{\layers}{
}

\pgfkeys{/pgf/fpu=false}

$\sqrt{(R_\oplus^3)^{1/3} - R_\oplus}$ equals $0$ not \radiusCoinsLayerCubedMtest !

\begin{align*}
R = \sqrt[3]{R_\oplus^3 + \frac{3}{2}r_m h_c}
\approx
\R
\end{align*}
%
Which means that the total number of layers are
%
\begin{align*}
n &= \frac{R - R_\oplus}{h_m} \\
\approx \layers
\end{align*}
\end{document}


The problem is that

(something^3)^(1/3) - something


does not equal zero, presumably because of rounding errors. It is clear that the expression above should evaluate to zero, however it does not. Instead I get -1400.0 which is complete nonsense. How can I get the fpu library too accurately calculate square roots?

My actual example is a little more involved, but it boils down to calculating the same thing.

• FWIW, running \directlua{earthRadiusM=6371*1000; tex.sprint((earthRadiusM^(1/3))^3 - earthRadiusM)} under LuaLaTeX returns -6.5192580223083e-09. Not exactly equal to zero either, but nevertheless about 12 orders of magnitude closer...
– Mico
Commented Sep 17, 2020 at 11:52
• I dont need hyper precision, but it would be nice that my answer was correct to at least two digits. Any suggestions on how to make it work without invoking the powers of lua? Commented Sep 17, 2020 at 11:57
• While pgfmath uses the FPU, it still converts the results back to lengths and text each step of the way. So far, only \pgfmathparseFPU retains full FPU accuracy through the whole computation. Commented Sep 17, 2020 at 15:02
• Read "What every programmer should know about floating point arithmetic". Long version: docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html. TL:DR version: floating-point-gui.de/basic Commented Sep 18, 2020 at 0:24

With xfp I get a more accurate result:

\documentclass{article}
\usepackage{xfp}

\begin{document}

\fpeval{((6371*1000)^(1/3))^3 - 6371*1000}

\end{document}


• This seems much more reasonable. Is there a better way to define variables with this method than \newcommand{\R}{\fpeval{((\earthRadiusM)^3 + 1.5 * (\coinRadiusM) * (\coinsTotalHeight))^(1/3)}}? Also, do you have any idea why pfu has so low accuracy? It is not even close to being correct Commented Sep 17, 2020 at 12:04
• @N3buchadnezzar \def\fpset#1#2{\edef#1{\fpeval{#2}}} then \fpset\whatever{pi^2}. Commented Sep 17, 2020 at 12:07
• @PhelypeOleinik perhaps with a \newcommand\whatever first (one of the problems with the pgf commands is that they overwrite commands without warnings ...) Commented Sep 17, 2020 at 12:09

Use the fp module of expl3 along with some syntactic sugar for variables that also ensures that we're not redefining existing commands.

However, you can't expect that (x3)1/3 = x.

\documentclass{article}

\usepackage{mathtools,xfp}

\ExplSyntaxOn

\NewDocumentCommand{\setfpvar}{mm}
{
\fp_zero_new:c { nebu_var_#1_fp }
\fp_set:cn { nebu_var_#1_fp } { #2 }
}
\NewExpandableDocumentCommand{\fpvar}{m}
{
\fp_use:c { nebu_var_#1_fp }
}

\ExplSyntaxOff

\begin{document}

\setfpvar{coinHeightM}{1.7 / 1000}
\setfpvar{coinsTotalHeight}{3.27*10^17}

}

\setfpvar{R}{
}

\setfpvar{layers}{
}

$\sqrt{(R_\oplus^3)^{1/3} - R_\oplus}$ equals
$\fpvar{radiusCoinsLayerCubedMtest}$

\bigskip

\begin{align*}
R = \sqrt[3]{R_\oplus^3 + \frac{3}{2}r_m h_c}
\approx
\fpvar{R}
\end{align*}
which means that the total number of layers is
\begin{align*}
n &= \frac{R - R_\oplus}{h_m} \\
`