1
$\begingroup$

I have a rather long expression (https://pastebin.com/jUsxdCCs) that is an analytical solution of a set of differential equations generated symbolically from Maple. I need to solve a set of equations like this in C, in a 32-bit system.

The long expression has a variable omega2, which when substituted leads to a number that I am interested in. When this operation is carried out in Maple, it results in a rather decent floating point number, but when I try to do it in C, (try with omega2 = 100), I get nan, the reason being that some of the sub-expressions in the long expression go out of range of double i.e. 2.3E-308 to 1.7E+308. I have been able to resolve this using 'long double` datatype available on some systems, but this isn't platform independent and is definitely not available on micro-controller that I plan to target which has a 32-bit ARM core.

One solution that I can imagine is to use the GMP library, that can handle numbers of infinite precision. Just want to know if there are any other hacks as GMP'ing the whole code can be fairly cumbersome and I don't have enough memory for GMP on the target machine.

The expression:

0.137197706359762e-3 * (-0.708981540362208e2 * cosh(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * (-0.727331795502673e6 * cosh(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) - 0.209935630488671e8 * cosh(0.354490770181104e2 * sqrt(omega2)) * cos(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) - 0.727331795502673e6 * cos(0.354490770181104e2 * sqrt(omega2)) * sinh(0.354490770181104e2 * sqrt(omega2))) * sqrt(omega2) / (0.348294222590750e3 * cosh(0.354490770181104e2 * sqrt(omega2)) * sinh(0.354490770181104e2 * sqrt(omega2)) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sqrt(omega2) - 0.348294222590750e3 * cos(0.354490770181104e2 * sqrt(omega2)) * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * cos(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * cosh(0.354490770181104e2 * sqrt(omega2)) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sinh(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.120668159897616e2 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.120668159897616e2 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.502654824574368e4 * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2 + 0.502654824574368e4 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2) - 0.708981540362208e2 * cosh(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) / (0.348294222590750e3 * cosh(0.354490770181104e2 * sqrt(omega2)) * sinh(0.354490770181104e2 * sqrt(omega2)) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sqrt(omega2) - 0.348294222590750e3 * cos(0.354490770181104e2 * sqrt(omega2)) * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * cos(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * cosh(0.354490770181104e2 * sqrt(omega2)) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sinh(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.120668159897616e2 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.120668159897616e2 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.502654824574368e4 * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2 + 0.502654824574368e4 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2) * sqrt(omega2) * (-0.209935630488671e8 * sinh(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) - 0.727331795502673e6 * cosh(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) + 0.727331795502673e6 * cos(0.354490770181104e2 * sqrt(omega2)) * sinh(0.354490770181104e2 * sqrt(omega2))) - 0.708981540362208e2 * cos(0.354490770181104e2 * sqrt(omega2)) * sinh(0.354490770181104e2 * sqrt(omega2)) * (-0.727331795502673e6 * cosh(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) - 0.209935630488671e8 * cosh(0.354490770181104e2 * sqrt(omega2)) * cos(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) - 0.727331795502673e6 * cos(0.354490770181104e2 * sqrt(omega2)) * sinh(0.354490770181104e2 * sqrt(omega2))) * sqrt(omega2) / (0.348294222590750e3 * cosh(0.354490770181104e2 * sqrt(omega2)) * sinh(0.354490770181104e2 * sqrt(omega2)) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sqrt(omega2) - 0.348294222590750e3 * cos(0.354490770181104e2 * sqrt(omega2)) * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * cos(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * cosh(0.354490770181104e2 * sqrt(omega2)) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sinh(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.120668159897616e2 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.120668159897616e2 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.502654824574368e4 * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2 + 0.502654824574368e4 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2) + 0.708981540362208e2 * cos(0.354490770181104e2 * sqrt(omega2)) * sinh(0.354490770181104e2 * sqrt(omega2)) / (0.348294222590750e3 * cosh(0.354490770181104e2 * sqrt(omega2)) * sinh(0.354490770181104e2 * sqrt(omega2)) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sqrt(omega2) - 0.348294222590750e3 * cos(0.354490770181104e2 * sqrt(omega2)) * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * cos(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * cosh(0.354490770181104e2 * sqrt(omega2)) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sinh(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.120668159897616e2 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.120668159897616e2 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.502654824574368e4 * pow(sinh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2 + 0.502654824574368e4 * pow(cosh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2) * sqrt(omega2) * (-0.209935630488671e8 * sinh(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) - 0.727331795502673e6 * cosh(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) + 0.727331795502673e6 * cos(0.354490770181104e2 * sqrt(omega2)) * sinh(0.354490770181104e2 * sqrt(omega2)))) * pow(omega2, -0.1e1 / 0.2e1) - 0.833338826540431e3;

Simplifying the expression in Maple, gives:

.414741582884880e315/(.702517310158453e311*pow(.487478610966388,2.)+.702517310158453e311*pow(-.873134928776922,2.)+502666.891390358*pow(.449112744897583e154,2.)*pow(.487478610966388,2.)+502666.891390358*pow(-.873134928776922,2.)*pow(.449112744897583e154,2.))*pow(.1e3,-.500000000000000)-833.338826540431

As can be seen, the sinh and cosh terms are blowing up inside, but cancel out ultimately.

$\endgroup$
  • 4
    $\begingroup$ One curious thing about your expression is that I checked some cases with MPFR, and it seems to suffer only from overflows, it's otherwise numerically stable. Can you manipulate the expression directly (it seems to be autogenerated)? I think you could try cancelling leading powers of cosh/sinh so that you only have non-exponential terms and negative powers of cosh/sinh. Then only harmless underflow would remain. $\endgroup$ – Kirill Jul 8 '17 at 23:15
  • $\begingroup$ You should simplify the expression to to an MWE that still produces the nan; this makes it easier to come up with an answer. $\endgroup$ – Nico Schlömer Jul 9 '17 at 13:45
  • $\begingroup$ The general case of this problem, given an arbitrary auto-generated expression that seems to fail with floating-point arithmetic, come up with an equivalent that can be evaluated directly, is actually very difficult. I believe there is no general solution, other than using MPFR like Henri Menke's answer, mostly just tricks that you apply by hand. CASes like Maple usually fall back on arbitrary-precision arithmetic. $\endgroup$ – Kirill Jul 9 '17 at 19:31
  • 1
    $\begingroup$ Would solving the original set of differential equations numerically instead of analytically be a solution? This avoids the problem Kirill mentions that "simplifying" such expressions (and coming up with numerically stable forms) is extremely difficult.... $\endgroup$ – GertVdE Jul 9 '17 at 20:40
  • $\begingroup$ @NicoSchlömer : I am sorry I cannot provide a smaller example. I have however, included a simplified version of the expression which brings out the problem. Hope this helps. $\endgroup$ – Chintan Pathak Jul 11 '17 at 22:32
5
$\begingroup$

I think there's a simple way to do this. You have a rational function of identical cosh/sinh terms, where every expression is a homogeneous polynomial in cosh/sinh, and the only problem is that these exponential terms overflow. The function does not diverge as these terms approach infinity, so if you divide every numerator and denominator by the same power of cosh, you'll get rid of overflow and be left with harmless underflow only.

Since the powers of cosh/sinh must match (the expression is bounded), we don't need to determine them by hand: we can (literally) replace every $\cosh(\cdots)$ by $1$, and every $\sinh(\cdots)$ by $\tanh(\cdots)$, getting the following expression:

0.137197706359762e-3 * (-0.708981540362208e2 * 1 * sin(0.354490770181104e2 * sqrt(omega2)) * (-0.727331795502673e6 * 1 * sin(0.354490770181104e2 * sqrt(omega2)) - 0.209935630488671e8 * 1 * cos(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) - 0.727331795502673e6 * cos(0.354490770181104e2 * sqrt(omega2)) * tanh(0.354490770181104e2 * sqrt(omega2))) * sqrt(omega2) / (0.348294222590750e3 * 1 * tanh(0.354490770181104e2 * sqrt(omega2)) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sqrt(omega2) - 0.348294222590750e3 * cos(0.354490770181104e2 * sqrt(omega2)) * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * pow(1, 0.2e1) * cos(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * 1 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * tanh(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.120668159897616e2 * pow(1, 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.120668159897616e2 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.502654824574368e4 * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2 + 0.502654824574368e4 * pow(1, 0.2e1) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2) - 0.708981540362208e2 * 1 * sin(0.354490770181104e2 * sqrt(omega2)) / (0.348294222590750e3 * 1 * tanh(0.354490770181104e2 * sqrt(omega2)) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sqrt(omega2) - 0.348294222590750e3 * cos(0.354490770181104e2 * sqrt(omega2)) * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * pow(1, 0.2e1) * cos(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * 1 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * tanh(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.120668159897616e2 * pow(1, 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.120668159897616e2 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.502654824574368e4 * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2 + 0.502654824574368e4 * pow(1, 0.2e1) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2) * sqrt(omega2) * (-0.209935630488671e8 * tanh(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) - 0.727331795502673e6 * 1 * sin(0.354490770181104e2 * sqrt(omega2)) + 0.727331795502673e6 * cos(0.354490770181104e2 * sqrt(omega2)) * tanh(0.354490770181104e2 * sqrt(omega2))) - 0.708981540362208e2 * cos(0.354490770181104e2 * sqrt(omega2)) * tanh(0.354490770181104e2 * sqrt(omega2)) * (-0.727331795502673e6 * 1 * sin(0.354490770181104e2 * sqrt(omega2)) - 0.209935630488671e8 * 1 * cos(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) - 0.727331795502673e6 * cos(0.354490770181104e2 * sqrt(omega2)) * tanh(0.354490770181104e2 * sqrt(omega2))) * sqrt(omega2) / (0.348294222590750e3 * 1 * tanh(0.354490770181104e2 * sqrt(omega2)) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sqrt(omega2) - 0.348294222590750e3 * cos(0.354490770181104e2 * sqrt(omega2)) * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * pow(1, 0.2e1) * cos(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * 1 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * tanh(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.120668159897616e2 * pow(1, 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.120668159897616e2 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.502654824574368e4 * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2 + 0.502654824574368e4 * pow(1, 0.2e1) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2) + 0.708981540362208e2 * cos(0.354490770181104e2 * sqrt(omega2)) * tanh(0.354490770181104e2 * sqrt(omega2)) / (0.348294222590750e3 * 1 * tanh(0.354490770181104e2 * sqrt(omega2)) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sqrt(omega2) - 0.348294222590750e3 * cos(0.354490770181104e2 * sqrt(omega2)) * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * pow(1, 0.2e1) * cos(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.348294222590750e3 * 1 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * tanh(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) + 0.120668159897616e2 * pow(1, 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.120668159897616e2 * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) + 0.502654824574368e4 * pow(tanh(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * pow(sin(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2 + 0.502654824574368e4 * pow(1, 0.2e1) * pow(cos(0.354490770181104e2 * sqrt(omega2)), 0.2e1) * omega2) * sqrt(omega2) * (-0.209935630488671e8 * tanh(0.354490770181104e2 * sqrt(omega2)) * sin(0.354490770181104e2 * sqrt(omega2)) * sqrt(omega2) - 0.727331795502673e6 * 1 * sin(0.354490770181104e2 * sqrt(omega2)) + 0.727331795502673e6 * cos(0.354490770181104e2 * sqrt(omega2)) * tanh(0.354490770181104e2 * sqrt(omega2)))) * pow(omega2, -0.1e1 / 0.2e1) - 0.833338826540431e3;

This only works when every expression is homogeneous in cosh/sinh and the degrees of numerators and denominators are the same, but it looks to me like they are, and numerically it checks out also. I think the same approach would work if the powers didn't match, but you'd probably have to cancel the powers manually or with a proper CAS.

$\endgroup$
  • $\begingroup$ Thanks a lot for your explanation and implementation. This approach is manual and may or may not work for me, and so I will wait one more day to see if something comes up or I will accept it. $\endgroup$ – Chintan Pathak Jul 11 '17 at 22:36
  • $\begingroup$ @ChintanPathak You could maybe try automating the process by introducing a dummy variable $t=\cosh(\cdots)$, replacing $\cosh\to t$, $\sinh \to t \tanh$, and using something like collect/factor/degree to extract and cancel powers of $t$. I'm not sure, but it doesn't sound too difficult. It doesn't generalize to other possible outputs of dsolve, I guess. $\endgroup$ – Kirill Jul 11 '17 at 23:23
  • $\begingroup$ The method you suggested works. The expression is of the form, (ax^2 + bx + c) / (dx^2 + ex + f), where x is represented by some cosh or sinh term that overflows. Taking the limit of this quantity for x->inf, leads to a/d and I applied this approximation to all terms and the answer is pretty close to Maple's solution. Thanks again for spending the time to make sense of the expression. $\endgroup$ – Chintan Pathak Jul 24 '17 at 20:52
  • $\begingroup$ @ChintanPathak Why not the exact equivalent $(a+b/x+c/x^2)/(d+e/x+f/x^2)$? $\endgroup$ – Kirill Jul 24 '17 at 22:54
  • 1
    $\begingroup$ @ChintanPathak If $1/x$ underflows, the underflow is harmless and would not affect accuracy (because you add zero to the answer instead of a tiny number, which is different from what happened with overflow). Or do you mean the case $x\approx 0$? Then the standard approach is to pick one of the two expressions depending on whether $x$ is above/below some threshold. (All I want to point out is that the approximation is unnecessary if you want the exact result.) $\endgroup$ – Kirill Jul 25 '17 at 0:26
2
$\begingroup$

If you do not want to (or just can't) juggle with the expression, just use the multiprecision library of your choice. Below an example with Boost.Multiprecision in C++. This has the invaluable advantage that you do not have to touch the original expression at all, thanks to operator and function overloading.

#include <iostream>
#include <cmath>

#include <boost/multiprecision/cpp_dec_float.hpp>

using boost::multiprecision::cpp_dec_float_50;

// Forward declaration for brevity
cpp_dec_float_50 expr(cpp_dec_float_50 const &omega2);


int main ()
{
  std::cout << expr (100) << '\n';
}

// Indentation by GNU indent
cpp_dec_float_50 expr(cpp_dec_float_50 const &omega2)
{
  return
    0.137197706359762e-3 * (-0.708981540362208e2 *
                cosh (0.354490770181104e2 * sqrt (omega2)) *
                sin (0.354490770181104e2 * sqrt (omega2)) *
                (-0.727331795502673e6 *
                 cosh (0.354490770181104e2 * sqrt (omega2)) *
                 sin (0.354490770181104e2 * sqrt (omega2)) -
                 0.209935630488671e8 * cosh (0.354490770181104e2 *
                             sqrt (omega2)) *
                 cos (0.354490770181104e2 * sqrt (omega2)) *
                 sqrt (omega2) -
                 0.727331795502673e6 * cos (0.354490770181104e2 *
                            sqrt (omega2)) *
                 sinh (0.354490770181104e2 * sqrt (omega2))) *
                sqrt (omega2) / (0.348294222590750e3 *
                         cosh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         sinh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         pow (sin
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) * sqrt (omega2) -
                         0.348294222590750e3 *
                         cos (0.354490770181104e2 *
                          sqrt (omega2)) *
                         pow (sinh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         sin (0.354490770181104e2 *
                          sqrt (omega2)) *
                         sqrt (omega2) +
                         0.348294222590750e3 *
                         pow (cosh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         cos (0.354490770181104e2 *
                          sqrt (omega2)) *
                         sin (0.354490770181104e2 *
                          sqrt (omega2)) *
                         sqrt (omega2) +
                         0.348294222590750e3 *
                         cosh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         pow (cos
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         sinh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         sqrt (omega2) +
                         0.120668159897616e2 *
                         pow (cosh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         pow (sin
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) +
                         0.120668159897616e2 *
                         pow (cos
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         pow (sinh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) +
                         0.502654824574368e4 *
                         pow (sinh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         pow (sin
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) * omega2 +
                         0.502654824574368e4 *
                         pow (cosh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         pow (cos
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) * omega2) -
                0.708981540362208e2 * cosh (0.354490770181104e2 *
                            sqrt (omega2)) *
                sin (0.354490770181104e2 * sqrt (omega2)) /
                (0.348294222590750e3 *
                 cosh (0.354490770181104e2 * sqrt (omega2)) *
                 sinh (0.354490770181104e2 * sqrt (omega2)) *
                 pow (sin (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * sqrt (omega2) -
                 0.348294222590750e3 * cos (0.354490770181104e2 *
                            sqrt (omega2)) *
                 pow (sinh (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * sin (0.354490770181104e2 *
                        sqrt (omega2)) *
                 sqrt (omega2) +
                 0.348294222590750e3 *
                 pow (cosh (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * cos (0.354490770181104e2 *
                        sqrt (omega2)) *
                 sin (0.354490770181104e2 * sqrt (omega2)) *
                 sqrt (omega2) +
                 0.348294222590750e3 * cosh (0.354490770181104e2 *
                             sqrt (omega2)) *
                 pow (cos (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * sinh (0.354490770181104e2 *
                         sqrt (omega2)) *
                 sqrt (omega2) +
                 0.120668159897616e2 *
                 pow (cosh (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * pow (sin (0.354490770181104e2 *
                             sqrt (omega2)),
                        0.2e1) +
                 0.120668159897616e2 *
                 pow (cos (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * pow (sinh (0.354490770181104e2 *
                              sqrt (omega2)),
                        0.2e1) +
                 0.502654824574368e4 *
                 pow (sinh (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * pow (sin (0.354490770181104e2 *
                             sqrt (omega2)),
                        0.2e1) * omega2 +
                 0.502654824574368e4 *
                 pow (cosh (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * pow (cos (0.354490770181104e2 *
                             sqrt (omega2)),
                        0.2e1) * omega2) *
                sqrt (omega2) * (-0.209935630488671e8 *
                         sinh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         sin (0.354490770181104e2 *
                          sqrt (omega2)) *
                         sqrt (omega2) -
                         0.727331795502673e6 *
                         cosh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         sin (0.354490770181104e2 *
                          sqrt (omega2)) +
                         0.727331795502673e6 *
                         cos (0.354490770181104e2 *
                          sqrt (omega2)) *
                         sinh (0.354490770181104e2 *
                           sqrt (omega2))) -
                0.708981540362208e2 * cos (0.354490770181104e2 *
                               sqrt (omega2)) *
                sinh (0.354490770181104e2 * sqrt (omega2)) *
                (-0.727331795502673e6 *
                 cosh (0.354490770181104e2 * sqrt (omega2)) *
                 sin (0.354490770181104e2 * sqrt (omega2)) -
                 0.209935630488671e8 * cosh (0.354490770181104e2 *
                             sqrt (omega2)) *
                 cos (0.354490770181104e2 * sqrt (omega2)) *
                 sqrt (omega2) -
                 0.727331795502673e6 * cos (0.354490770181104e2 *
                            sqrt (omega2)) *
                 sinh (0.354490770181104e2 * sqrt (omega2))) *
                sqrt (omega2) / (0.348294222590750e3 *
                         cosh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         sinh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         pow (sin
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) * sqrt (omega2) -
                         0.348294222590750e3 *
                         cos (0.354490770181104e2 *
                          sqrt (omega2)) *
                         pow (sinh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         sin (0.354490770181104e2 *
                          sqrt (omega2)) *
                         sqrt (omega2) +
                         0.348294222590750e3 *
                         pow (cosh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         cos (0.354490770181104e2 *
                          sqrt (omega2)) *
                         sin (0.354490770181104e2 *
                          sqrt (omega2)) *
                         sqrt (omega2) +
                         0.348294222590750e3 *
                         cosh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         pow (cos
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         sinh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         sqrt (omega2) +
                         0.120668159897616e2 *
                         pow (cosh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         pow (sin
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) +
                         0.120668159897616e2 *
                         pow (cos
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         pow (sinh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) +
                         0.502654824574368e4 *
                         pow (sinh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         pow (sin
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) * omega2 +
                         0.502654824574368e4 *
                         pow (cosh
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) *
                         pow (cos
                          (0.354490770181104e2 *
                           sqrt (omega2)),
                          0.2e1) * omega2) +
                0.708981540362208e2 * cos (0.354490770181104e2 *
                               sqrt (omega2)) *
                sinh (0.354490770181104e2 * sqrt (omega2)) /
                (0.348294222590750e3 *
                 cosh (0.354490770181104e2 * sqrt (omega2)) *
                 sinh (0.354490770181104e2 * sqrt (omega2)) *
                 pow (sin (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * sqrt (omega2) -
                 0.348294222590750e3 * cos (0.354490770181104e2 *
                            sqrt (omega2)) *
                 pow (sinh (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * sin (0.354490770181104e2 *
                        sqrt (omega2)) *
                 sqrt (omega2) +
                 0.348294222590750e3 *
                 pow (cosh (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * cos (0.354490770181104e2 *
                        sqrt (omega2)) *
                 sin (0.354490770181104e2 * sqrt (omega2)) *
                 sqrt (omega2) +
                 0.348294222590750e3 * cosh (0.354490770181104e2 *
                             sqrt (omega2)) *
                 pow (cos (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * sinh (0.354490770181104e2 *
                         sqrt (omega2)) *
                 sqrt (omega2) +
                 0.120668159897616e2 *
                 pow (cosh (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * pow (sin (0.354490770181104e2 *
                             sqrt (omega2)),
                        0.2e1) +
                 0.120668159897616e2 *
                 pow (cos (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * pow (sinh (0.354490770181104e2 *
                              sqrt (omega2)),
                        0.2e1) +
                 0.502654824574368e4 *
                 pow (sinh (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * pow (sin (0.354490770181104e2 *
                             sqrt (omega2)),
                        0.2e1) * omega2 +
                 0.502654824574368e4 *
                 pow (cosh (0.354490770181104e2 * sqrt (omega2)),
                  0.2e1) * pow (cos (0.354490770181104e2 *
                             sqrt (omega2)),
                        0.2e1) * omega2) *
                sqrt (omega2) * (-0.209935630488671e8 *
                         sinh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         sin (0.354490770181104e2 *
                          sqrt (omega2)) *
                         sqrt (omega2) -
                         0.727331795502673e6 *
                         cosh (0.354490770181104e2 *
                           sqrt (omega2)) *
                         sin (0.354490770181104e2 *
                          sqrt (omega2)) +
                         0.727331795502673e6 *
                         cos (0.354490770181104e2 *
                          sqrt (omega2)) *
                         sinh (0.354490770181104e2 *
                           sqrt (omega2)))) *
    pow (omega2, -0.1e1 / 0.2e1) - 0.833338826540431e3;
}
$\endgroup$
  • $\begingroup$ Thank you for the code sample. However, this approach wont work for me due to limited memory on my target CPU. As per my reading, boost.multiprecision is a wrapper around GMP, MPFR, MPIR etc. and downloading from here was >8 Mb after unzipping. $\endgroup$ – Chintan Pathak Jul 11 '17 at 22:27
  • $\begingroup$ @ChintanPathak Do you have to compile on the device? $\endgroup$ – Henri Menke Jul 11 '17 at 22:39
  • $\begingroup$ No, I can compile on my system and deploy there. Does the target then not need anything else but the machine code ? $\endgroup$ – Chintan Pathak Jul 11 '17 at 22:41
  • $\begingroup$ @ChintanPathak You only need the binary and the dynamically linked libraries. If I compile the above example with -Os for file-size optimisation I end up with a 307Kb binary. $\endgroup$ – Henri Menke Jul 11 '17 at 22:43
  • $\begingroup$ Conflicting answere here: quora.com/… $\endgroup$ – Chintan Pathak Jul 11 '17 at 23:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.