# serious scientific calculating with small numbers, using different methods

I got stuck in a problem with pgfmath module. When trying to calculate the following functions values the result that gets printed out is always zero. I guess the problem is an underflow of pgfmaths co-domain, and the mathematical engine cuts off the small values.

This is an minimal example of what i was trying to explain in the introduction:

\documentclass[tikz]{standalone}
\usepackage[fleqn]{amsmath}
% physical constants:
\pgfmathdeclarefunction{m0}{0}{%
\pgfmathparse{4*pi*1e-7}%
}

\pgfmathdeclarefunction{K}{0}{%
\pgfmathparse{m0*pi/4}%
}
\pgfmathdeclarefunction{c1}{0}{% c1 = K/1.45
\pgfmathparse{K/1.45}%
}
\pgfmathdeclarefunction{c2}{1}{% c2(gamma) = c1/gamma^2
\pgfmathparse{c1/(#1)^2}%
}
\pgfmathdeclarefunction{c3}{2}{% c3(gamma,lambda)
\pgfmathparse{K/((#1)*((#2)+0.45))}%
}
\pgfmathdeclarefunction{DL_rel}{2}{% DL_rel(lambda, gamma)
\pgfmathparse{(2*sqrt((580*#1+261)*#2^3)+40*#1+18)/(29*#2^3-20*#1-9)}%
}

% the problematic pgfmath-function
\pgfmathdeclarefunction{N1}{5}{% N1(soluition, Lmin, gamma, lambda, d1)
\pgfmathsetmacro\numa{(DL_rel(#4,#3)*#2)/2+#2}%
\pgfmathsetmacro\numb{sqrt(\numa^2-c1*c3(#3,#4))/(4*c2(#3)^2)}%
\pgfmathsetmacro\denom{2*c1*#5}%
\pgfmathparse{(#1 -1) ?%
(sqrt((\numa + \numb)/\denom))%
:%
(sqrt((\numa - \numb)/\denom))%
}%
}

\begin{document}
This is the equation im trying to solve:
$$\label{eq:N1} N_{1_{1,2}} =\sqrt{ \frac{\frac{\Delta L}{2} +L_{min} \pm \sqrt{\left(\frac{\Delta L}{2} + L_{min}\right)^2 - \frac{c_1 c_3}{4\cdot c_2^2}\cdot \Delta L^2}} {2\cdot c_1 d_1} }$$
This should print out its solution:
\begin{equation*}
N_{1_{1,2}} =
\left\{ \begin{array}{l}
\pgfmathparse{N1(1, 4e-4, 1.2, 0.466, 0.115)}\pgfmathresult\\
\pgfmathparse{N1(2, 4e-4, 1.2, 0.466, 0.115)}\pgfmathresult
\end{array}\right.
\end{equation*}
\end{document}


It took me already way to long to get so far and it would be very frustrating for me to stop here, and get all the math done in a different application.

Maybe someone out there has a solution to my problem.

Anyway i want use this post to say thank you! this has been a wonderful place to learn Latex so far, good job everybody!

• Trying to compile your code, it misses the definition of \stellbereich – Alain Remillard Feb 16 '16 at 15:25
• sorry was an translation error..i've fixed it above – M Arpogaus Feb 16 '16 at 15:27
• The default PGF math cannot do this calculation because one of the inputs is too small (4e-6). You could try putting \usepackage{fp} \usetikzlibrary{fixedpointarithmetic} \tikzset{fixed point arithmetic} in your preamble to get round this. It still won't work however because (a) the \pgfmathparse expression in function c3 appears to missing an operator in the denominator, and (b) possibly the final pgfmathparse in function N1 is missing an operator in the conditional. There maybe other things, but I guess these two need to be fixed first. – Mark Wibrow Feb 16 '16 at 16:21
• ...an alternative to get round the small numbers problem is \usetikzlibrary{fpu} \pgfkeys{/pgf/fpu=true} in the preamble. But some extra work will be required at the end to transform the output of the parser in this case. – Mark Wibrow Feb 16 '16 at 16:35
• Hi Mark! Thanks for your advice. I've added the missing * at c3. The condition at N1 is right by the way. I subtract 1 from the input so im able to decide more natural between the first solution (N1(1,..)) and second solution (N1(2,..) in the source later on. trying to add \usetikzlibrary{fpu} \pgfkeys{/pgf/fpu=true} now. – M Arpogaus Feb 16 '16 at 16:41

### update to re-insist on numerical instability

Here is from the log output, using gamma=1.2, lambda=0.466 and Lmin=4e-6, which were the values from the original OP. Similar, but differing result are observed with gamma=1.2, lambda=0.44.

4: -1.00000e-13
5: 0
6: 0
7: 0
8: 0
9: 0
10: 2.00000e-19
11: 1.00000e-20
12: 2.00000e-21
13: 1.00000e-22
14: 1.00000e-23
15: 1.00000e-24
16: 0
17: 0
18: 0
19: 0
20: 1.00000e-29
21: 1.00000e-30
22: 0
23: 0
24: 1.00000e-33
25: 1.00000e-34
26: 0
27: 2.00000e-36
28: 0
29: 0
30: 0
31: -1.00000e-40
32: 2.00000e-41
33: 0
34: 2.00000e-43
35: 0
36: 1.00000e-45
37: 1.00000e-46
38: 0
39: -1.00000e-48
40: 0
41: 0
42: -1.00000e-51
43: 1.00000e-52
44: 1.00000e-53
45: 1.00000e-54
46: 0
47: 1.00000e-56
48: 1.00000e-57
49: 0
50: 1.00000e-59
51: 2.00000e-60
52: 0
53: 0
54: 2.00000e-63
55: 1.00000e-64
56: 2.00000e-65
57: 0
58: 0
59: 1.00000e-68
60: 1.00000e-69
61: 0
62: -1.00000e-71
63: 0
64: -1.00000e-73
65: 0
66: 0
67: 1.00000e-76
68: 0
69: 0
70: -1.00000e-79
71: 1.00000e-80
72: 0
73: 0
74: 0
75: 0
76: 0
77: 0
78: -1.00000e-87
79: 1.00000e-88
80: 1.00000e-89
81: 1.00000e-90
82: 1.00000e-91
83: 1.00000e-92
84: 1.00000e-93
85: 0
86: 0
87: 0
88: 1.00000e-97
89: 1.00000e-98
90: 1.00000e-99
91: 0
92: 1.00000e-101


Note: I have set constant K to value 1 to skip the computations with Pi. This modifies the result here as the rounding float operations are not exactly the same.

Source code to generate the above.

\documentclass{article}
\usepackage{xintexpr}% tested with 1.2e release
%\xintverbosetrue

\usepackage[fleqn]{amsmath}

\begin{document}
This is the equation im trying to solve:
$$\label{eq:N1} N_{1_{1,2}} =\sqrt{ \frac{\frac{\Delta L}{2} +L_{min} \pm \sqrt{\left(\frac{\Delta L}{2} + L_{min}\right)^2 - \frac{c_1 c_3}{4\cdot c_2^2}\cdot \Delta L^2}} {2\cdot c_1 d_1} }$$

The $\Delta L/L_{min}$, $c_1$, $c_2$, $c_3$ are functions of some parameters
($L_{min}$ is only an overall scaling) and these functions are set-up in such
a way that actually the square root at the numerator of the fraction vanishes
exactly, independently of the values of the parameters. Numerically though it
may be found non zero and this will then induce a catastrophic drop in
precision in the final result. And it may even happen that it will be found
negative, thus potentially raising an error in the complete evaluation.

See the compilation log for this problematic $\left(\frac{\Delta L}{2} + L_{min}\right)^2 - \frac{c_1 c_3}{4\cdot c_2^2}\cdot \Delta L^2$ evaluated
with the float precision set to various values. You will see it turns negative
at times.

% \xintdeffloatvar pi:=
% 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342;

\xintFor* #1 in {\xintSeq{4}{92}}\do
{
\xintDigits := #1;

% constants
% \xintdeffloatvar m0:= 4pi*1e-7;
% \xintdeffloatvar K := m0*pi/4;
\xintdeffloatvar K := 1;
\xintdeffloatvar c1:= K/1.45;

% functions
\xintdeffloatfunc c2(u)  := c1/u^2;
\xintdeffloatfunc c3(u,v):= K/(u(v+0.45));

% attention arguments like in OP-update, permuted compared to the OP-original
\xintdeffloatfunc DL_rel(u,v):= (2sqrt((580v+261)u^3)+40v+18)/(29u^3-20v-9);

\xintdeffloatfunc DL(t,u,v)  := DL_rel(u,v)*t;

% Notice that t=Lmin acts only as an overall scaling factor.
\xintdeffloatfunc numbsquared(t,u,v):=
subs((Z/2+t)^2-c1*c3(u,v)/(4*c2(u)^2)*Z^2, Z=DL(t,u,v));

\typeout{#1: \xintthefloatexpr [6] numbsquared (4e-6, 1.2, 0.466)\relax }

}

\end{document}


I have also tested with Maple. Again setting K to 1.

numbsquared := proc (N)
local m0, K, c1, c2, c3, DL_rel, DL, localnumbsquared;
Digits:=N;
# m0 := 4*Pi*1e-7;
# K  :=  m0*Pi/4;
K  := 1;
c1 := K/1.45;
c2 := u->c1/u^2;
c3 := (u,v)->K/(u*(v+0.45));
DL_rel := (u,v)->(2*sqrt((580*v+261)*u^3)+40*v+18)/(29*u^3-20*v-9);
DL := (t,u,v)->DL_rel(u,v)*t;
localnumbsquared := (t,u,v)->subs(Z=DL(t,u,v),(Z/2+t)^2-c1*c3(u,v)/(4*c2(u)^2)*Z^2);
return localnumbsquared(4e-6, 1.2, 0.466)
end proc:
for N from 4 to 92 do printf("%2d, %e\n", N, numbsquared(N)) end do;


The results are of the same type but with both coincidences and differences.

 4, 0.000000e+00
5, 0.000000e+00
6, 0.000000e+00
7, 0.000000e+00
8, 0.000000e+00
9, 1.000000e-18
10, 2.000000e-19
11, 1.000000e-20
12, -1.000000e-21
13, 1.000000e-22
14, -1.000000e-23
15, 1.000000e-24
16, 0.000000e+00
17, 0.000000e+00
18, 0.000000e+00
19, 0.000000e+00
20, 1.000000e-29
21, 1.000000e-30
22, 0.000000e+00
23, 1.000000e-32
24, -1.000000e-33
25, 1.000000e-34
26, 0.000000e+00
27, 2.000000e-36
28, 1.000000e-37
29, 0.000000e+00
30, 0.000000e+00
31, -1.000000e-40
32, 2.000000e-41
33, 0.000000e+00
34, 2.000000e-43
35, -1.000000e-44
36, 1.000000e-45
37, 1.000000e-46
38, 0.000000e+00
39, -1.000000e-48
40, 0.000000e+00
41, 0.000000e+00
42, -1.000000e-51
43, 1.000000e-52
44, -1.000000e-53
45, 1.000000e-54
46, 0.000000e+00
47, -1.000000e-56
48, 1.000000e-57
49, 0.000000e+00
50, 1.000000e-59
51, 0.000000e+00
52, 0.000000e+00
53, 0.000000e+00
54, 2.000000e-63
55, 1.000000e-64
56, 0.000000e+00
57, 0.000000e+00
58, 0.000000e+00
59, 1.000000e-68
60, 1.000000e-69
61, 0.000000e+00
62, -1.000000e-71
63, 0.000000e+00
64, -1.000000e-73
65, 0.000000e+00
66, 0.000000e+00
67, 1.000000e-76
68, -2.000000e-77
69, -1.000000e-78
70, -1.000000e-79
71, 1.000000e-80
72, 0.000000e+00
73, 0.000000e+00
74, 0.000000e+00
75, 0.000000e+00
76, 0.000000e+00
77, 0.000000e+00
78, 0.000000e+00
79, 1.000000e-88
80, -1.000000e-89
81, -1.000000e-90
82, 1.000000e-91
83, 0.000000e+00
84, -1.000000e-93
85, 0.000000e+00
86, 0.000000e+00
87, 0.000000e+00
88, -1.000000e-97
89, -2.000000e-98
90, 1.000000e-99
91, -2.000000e-100
92, 1.000000e-101


### update to comment an intriguing numerical curiosity/instability

First version of OP asked for N1(1, 4e-6, 1.2, 0.466, 0.115), which is treated below. Turns out apparently that the \numb in that case is zero. But depending on the float precision it may be found or not to be zero. The \numa being only about 1e-5 the precision of final result may be drastically reduced by small but non zero \numb.

I compared xint with maple and obtained similar results (differ only in last digit) for 16, 20, 24 digits of precision (takes time to use maple, more testing on the way). For Pi, I did my tests in xint starting with 94 decimals and multiplying first by 1.0 to reduce to stated \xintDigits. Thus I ran the code below with

\xintDigits := 16; % or 20, 24, 28, ...
\xintdeffloatvar pi:= 1.*3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342;


On the Maple side, Pi is treated as a symbol until the final evalf.

Here is how the results look like on the xint side with higher and higher float precision:

8.038962683509860
8.0389626847662074875
8.03896268352242165100730
8.038962683509858157707745289
8.0389626835098594140570752438847
8.03896268350985817027123858823346261
8.038962683509858157707745288681429203705
8.0389626835098581577090016380113844070447859
8.03896268350985815770775785217472875573691848736


(this is for the solution with \numa+\numb). Notice how the result with 16 digits is much better than the one with 20 digits!!!!!!! and it is even quite better than the one for 24 digits! (but not as good as 28 digits). This is due to the fact that at 16 digits, both xint and maple find \numb to be zero, but about 3e-15 at 20 digits which induces a big error in the sum, as \numa is about 1e-5.

With 92 digits of precision one finds \numb about 3e-51 numerically. It the exact value is zero, this means that it spoils the digits of the result after about 46 of them...

With 92 digits of precision, Maple finds \numb to be

> evalf(Q(4e-6, 1.2, 0.466, 0.115));
0.3162277660168379331998893544432718533719555139325216826857504852792594438\

-50
6392382213442481084 10


And xint obtains

3.1622776601683793319988935444327185337195551393252168268575048527925944386392382213442481084e-51  you can see it matches up to the end. ;-)

update what sort of fool I am! How many years before recognizing the square root of 10 ???? notice that the above is essentially the square root of 1e-101 ... the reason appears simple enough that \numb is a square root of a difference, and somehow this difference rather than being found to be zero turns out be 1e-101 due to a rounding error in the last, 92th digit of each term, which probably are of the order of 1e-10 !!! Yes, this should explain it for all levels N of float precision. I guess sometimes the difference is zero, sometimes it gives 1e-(9+N). For example with N=20 one can expect a difference of 1e-29, hence about 3e-15 on square root which is exactly what is observed. Strangely the numerics never seem to give a difference of -1e-(9+N) which would raise on error on square root.

As the returned values decrease with increasing precision perhaps I can trust the exact one is zero (I have not done the algebra). If the exact value is really zero, adding or subtracting the above to something which is 0.000010117182975... will corrupt it after about 46 significant digits, destroying the 92 digits float evaluations to get it...

Very surprising ! (but do read quoted block above)

This should be taken into account when comparing with any other math engine: the formula is numerically instable due to possible catastrophic cancellation in a subtraction.

Here is an approach using another math engine. It knows only square root, but this is enough here. Note that in the given example \numb turns out to be exactly zero.

\documentclass[tikz]{standalone}
\usepackage{xintexpr}% tested with 1.2e release
\usepackage[fleqn]{amsmath}

% constants
\xintdeffloatvar pi:= 3.14159265358979323846;
\xintdeffloatvar m0:= 4pi*1e-7;
\xintdeffloatvar K := m0*pi/4;
\xintdeffloatvar c1:= K/1.45;

% functions
\xintdeffloatfunc c2(x)  := c1/x^2;

\xintdeffloatfunc c3(x,y):= K/(x(y+0.45));

\xintdeffloatfunc DL_rel(u,v):=
(2*sqrt((580*u+261)*v^3)+40*u+18)/(29*v^3-20*u-9);

% This is allowed by xint parser also (tacit multiplications):
% \xintdeffloatfunc DL_rel(u,v):= (2sqrt((580u+261)v^3)+40u+18)/(29v^3-20u-9);

% Of course we could simplify here by defining more intermediate functions.
% We could define "numa" and "numb" functions, and set them up as functions
% of an already computed "DL_rel" which serves in both.
% It is possible to use the "subs(expression, x=...)" syntax.
% Limitation is that the dummy parameter must be a single letter.
% Also, the inner-most subs will have the last defined thing, and the
% outer-most subs the first defined thing.
\xintdeffloatfunc N1(a,t,u,v,w):=
subs(subs(subs(subs(
if(a=1, sqrt((P+Q)/D), sqrt((P-Q)/D)),
% debugging because something is strange with Q = \numb which is zero
% (P, sqrt(c1*c3(u,v))/c2(u)*X ),
% well after all it was CORRECT that Q was zero with these numerics
Q = sqrt(P^2-c1*c3(u,v)/(c2(u)^2)*X^2)% =\numb,
),
P = X+t % P=\numa, and I think t is Lmin
),
X = DL_rel(v,u)*t/2 % X= DeltaL/2
),
D = 2c1*w % D=\denom
)% must use single letters in subs
;%

\begin{document}
This is the equation im trying to solve:
$$\label{eq:N1} N_{1_{1,2}} =\sqrt{ \frac{\frac{\Delta L}{2} +L_{min} \pm \sqrt{\left(\frac{\Delta L}{2} + L_{min}\right)^2 - \frac{c_1 c_3}{4\cdot c_2^2}\cdot \Delta L^2}} {2\cdot c_1 d_1} }$$

This should print out its solution:
\begin{equation*}
N_{1_{1,2}} =
\left\{ \begin{array}{l}
\xintthefloatexpr N1(1, 4e-6, 1.2, 0.466, 0.115)\relax\\
\xintthefloatexpr N1(2, 4e-6, 1.2, 0.466, 0.115)\relax
\end{array}\right.
\end{equation*}
\end{document}


In this example, the two solutions are the same, as the \numb vanishes ...

Notice that the solutions are computed expandably. That matters to some (fools...).

\xintDigits := 32;
\xintdeffloatvar pi:= 3.141592653589793238462643383279503;

• Thank you so much! That's exactly what i need. And much more readable than pgfmath! need to get deeper into xintexpr. Do you know a good place to learn the basics? – M Arpogaus Feb 16 '16 at 17:07
• not really. The doc is not as famed as TikZ's one ;-). Besides, so far the sole non elementary function implemented is square root. It works in arbitrary precision (say up to one hundred digits, beyond that, time cost). – user4686 Feb 16 '16 at 17:10
• defined constants may be multi-letters with digits underscore ..., but variables for functions must be single-letters lowercase or uppercase (and they should not naturally have been used as constants). Also, symbols for the "substitution" on-the-fly intermediate computations must be single letters, I could not use gamma=... for example. There should be a better syntax than subs(subs(subs,...). Of course you can do the sub-evaluations first in macros (using \edef\numa{\xintthefloatexpr ....\relax}) or constants/functions. – user4686 Feb 16 '16 at 17:13
• if(a=1, sqrt((P+Q)/D), sqrt((P-Q)/D)) syntax a priori evaluates both branches and the syntax for lazy evaluation is rather (a=1)?{YES}{NO}, but there is an exception when it is used to define a function like here. Then, the (a=1)?{YES}{NO} is not allowed and one must use if(a=1,YES,NO) but when the computation will really be done only the correct branch will be evaluated. This info is buried somewhere in the xintexpr doc... – user4686 Feb 16 '16 at 17:22
• All right! Thanks again for your help. Im doing some research tomorrow, the doc isn't to bad. Do you know if its possible to plot \xintdeffloatfunc using pgfplot? or is there a better way? Sorry for all the questions but its my first week using latex! – M Arpogaus Feb 16 '16 at 17:24

You could try going the lualatex route and then have all the precision of the math in lua, so no pgf at all. The following may show the wrong answer, but that isn't due to the mathematical capabilities of lua, more likely just a rather hasty translation of the required formula:

\documentclass[preview,border=5]{standalone}
\usepackage[fleqn]{amsmath}
\usepackage{luacode}
\begin{luacode*}
pi = math.pi
sqrt = math.sqrt
m0 = 4 * pi * 1e-7
K = m0 * pi / 4
c1 = K / 1.45
c2 = function (x) return c1 / (x^2); end
c3 = function (x, y) return K / (x * (y + 0.45)); end
DLrel = function (x, y)
return (2 * sqrt((580 * x + 261) * y^3) + 40 * x + 18) /
(29 * y^3 - 20 * x - 9)
end
N1 = function (s, Lmin, g, l, d1)
dL = DLrel(l, g) * Lmin
nm = dL / 2 + Lmin
rt = (nm^2 - c1 * c3(g, l) * 0.25 * c2(g)^-2 * dL^2)
dn = 2 * c1 * d1
s = -(s % 2) * 2 + 1
return sqrt((nm + s * sqrt(rt)) / dn)
end
\end{luacode*}
\def\luaprint#1{\directlua{tex.print(#1)}}
\begin{document}
\begin{equation*}
N_{1_{1,2}} = \left\{
\begin{array}{l}
\luaprint{N1(1, 4e-4, 1.2, 0.44, 0.115)}
\\
\luaprint{N1(2, 4e-4, 1.2, 0.44, 0.115)}
\end{array}\right.
\end{equation*}
\end{document}


• +1 does it use the new math libraries in LuaTeX or was it already there earlier ? – user4686 Feb 17 '16 at 16:20
• @jfbu the math used here is described in the lua documentation and (in principle) doesn't need any special libraries, although I used the luacode package to make writing the code easier. – Mark Wibrow Feb 17 '16 at 16:25
• I will probably delete my comments, as I incorporated the lesson learned in an update to my answer. – user4686 Feb 17 '16 at 18:24
• i expected the square root to be almost zero, that's why i chose those values to verify my functions. i think the reason for the small failure on higher precisions is caused by used arguments.Those are rounded and therefore not the exact algebraic value to get the \numb exactly zero. Anyways what ive learned so far is that pgfmath isn't giving me the results i expect, there is some strange inaccuracy wich causes a negative root and pgfmath prints out NaN. – M Arpogaus Feb 18 '16 at 12:49
• Thats why i favor Marks solution. but i am till stuck in some basic problems, i need to get used to lua now. going to update the post later the day. – M Arpogaus Feb 18 '16 at 12:49

## Update 1: pgfmath + fpu

I've played a bit around with XINT today but due to its lack of mathematical functions i returned to my pgfmath approach. There are a few more formulas i need to embed later on, for whom i need trigonometric functions etc. and as far as i can see there is no implementation for those in xint available. With Marks hint on using \usetikzlibrary{fpu} for higher precision i got it finally working, after some tweaking.

That's my code so far:

\documentclass[tikz]{standalone}
\usepackage[fleqn]{amsmath}
\usetikzlibrary{fpu}

% physical constants:
\pgfmathdeclarefunction{m0}{0}{%
\pgfmathparse{4*pi*1e-7}%
}

\pgfmathdeclarefunction{K}{0}{%
\pgfmathparse{m0*pi/4}%
}
\pgfmathdeclarefunction{c1}{0}{% c1 = K/1.45
\pgfmathparse{K/1.45}%
}
\pgfmathdeclarefunction{c2}{1}{% c2(gamma) = c1/gamma^2
\pgfmathparse{c1/(#1)^2}%
}
\pgfmathdeclarefunction{c3}{2}{% c3(gamma,lambda)
\pgfmathparse{K/((#1)*((#2)+0.45))}%
}
\pgfmathdeclarefunction{DL_rel}{2}{% DL_rel(gamma, lambda)
\pgfmathparse{(2*sqrt((580*#2+261)*#1^3)+40*#2+18)/(29*#1^3-20*#2-9)}%
}

% the problematic pgfmath-function
%                                   #1         #2    #3     #4      #5
\pgfmathdeclarefunction{N1}{5}{% N1(soluition, Lmin, gamma, lambda, d1)
\pgfmathsetmacro\DL{(DL_rel(#3,#4)*#2)}%
\pgfmathsetmacro\numa{\DL/2+#2}%
\pgfmathsetmacro\numb{sqrt(\numa^2 - c1*c3(#3,#4)/(4*c2(#3)^2)*\DL^2)}%
\pgfmathsetmacro\denom{2*c1*#5}%
\pgfmathfloatparse{(#1 == 1) ?%
(sqrt((\numa + \numb)/\denom))%
:%
(sqrt((\numa - \numb)/\denom))%
}%
}

\begin{document}
\pgfkeys{/pgf/fpu=true}
\pgfmathparse{N1(1, 4e-6, 1.2, 0.44, 0.115)}\pgfmathprintnumber[sci, precision=2]{\pgfmathresult}\\
\pgfkeys{/pgf/fpu=false}
This is the equation im trying to solve:
$$\label{eq:N1} N_{1_{1,2}} =\sqrt{ \frac{\frac{\Delta L}{2} +L_{min} \pm \sqrt{\left(\frac{\Delta L}{2} + L_{min}\right)^2 - \frac{c_1 c_3}{4\cdot c_2^2}\cdot \Delta L^2}} {2\cdot c_1 d_1} }$$
This should print out its solution:
\pgfkeys{/pgf/fpu=true}
\begin{equation*}
N_{1_{1,2}} =
\left\{ \begin{array}{l}
\pgfmathparse{N1(1, 4e-6, 1.2, 0.44, 0.115)}\pgfmathprintnumber[fixed, precision=2]{\pgfmathresult}\\
\pgfmathparse{N1(2, 4e-6, 1.2, 0.44, 0.115)}\pgfmathprintnumber[fixed, precision=2]{\pgfmathresult}
\end{array}\right.
\end{equation*}
\pgfkeys{/pgf/fpu=false}
\end{document}


but i am unsure about the results precision. compared with jfbus solution and my output from wxMaxima the squareroot in the numerator now isn't zero anymore, although its very small (e-18). Even if its a very small difference, i'm interested in pgfmaths overall precision. because i'am planing to get all the math done in Latex later on, when i'm compiling, and i'am seeing my self questioning every computation afterwards.

So all in all guys, do you think its a good idea to use pgfmath for this kind of application? Has someone done something similar, or is it just a bad idea to trust a word processor on math?

By the way, i've now used \pgfmathsetmacro for saving results i need later on. Is there a better solution, or is this the common way?

## Update 2: LuaTeX + luacode

So far so good. I know gave up my pgfmath approach, due to its inaccuracy and concentrate on Marks luaTeX version now.

This solution offers me a nice way to embeed the calculated solutions in formulas and text. but to make this stuff really useful i need to get it working in pgfplot enviroment. Here is an Example:

\documentclass[preview,border=5]{standalone}
\usepackage[fleqn]{amsmath}
\usepackage{luacode}
\usepackage{tikz}
\usepackage{pgfplots}
\begin{luacode*}
-- test
pi = math.pi
sqrt = math.sqrt
m0 = 4 * pi * 1e-7
K = m0 * pi / 4
c1 = K / 1.45
c2 = function (g) return c1 / (g^2); end
c3 = function (g, l) return K / (g * (l + 0.45)); end
DLrel = function (g, l)
return  (2 * sqrt((580 * l + 261) * g^3) + 40 * l + 18) /
(29 * g^3 - 20 * l - 9)
end
DLrel_lmd = function (l) return DLrel(sqrt(l^2+1),l); end

N1 = function (s, Lmin, g, l, d1)
dL = DLrel(g, l) * Lmin
nm = dL / 2 + Lmin
rt = (nm^2 - c1 * c3(g, l) * 0.25 * c2(g)^-2 * dL^2)
dn = 2 * c1 * d1
s = -(s % 2) * 2 + 1
return sqrt((nm + s * sqrt(rt)) / dn)
end

print = function (d,s)
if d == 0 then
format = "%d"
else
format = "%." .. d .. "f"
end
tex.sprint(string.format(format,s));
end

\end{luacode*}
\newcommand{\lp}[2][16]{\directlua{print(#1,#2)}}

\begin{document}
now i can do calculations in lua and print them out in formulas:
\begin{equation*}
N_{1_{1,2}} = \left\{
\begin{array}{l}
\lp[4]{N1(1, 4e-6, 1.2, 0.446, 0.115)}
\\
\lp[4]{N1(2, 4e-6, 1.2, 0.446, 0.115)}
\end{array}\right.
\end{equation*}

Some serious research done here..\\
..oh some new values i came across, i can embed them in my text:\\

\begin{luacode}
now  = 42.42424242424242
\end{luacode}

$$now = \lp[4]{now}$$
as you can see in this fancy plot:

\begin{tikzpicture}
\begin{axis}[
width=0.5\linewidth,
height=10cm,
xmin=0, xmax=180,
ymin=0, ymax=50,
xlabel=$time$,
ylabel=Questions accouring while learning LaTeX,
]
thick,
blue,
domain=0:180,
samples=200,
]{42*sin(x)}; % This value should be taken from lua
node[pos=0] {$now = 42$}; % This value should be taken from lua