1
$\begingroup$

First of all, apologies for the large numbers, I couldn't find a working example for a far lower order and the zeros are necessary... my problem is this: I have a large polynomial (from high-order FEM) I want to evaluate in a compiled function:

fc=Compile[{x},-820.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+447720.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 x-(9.1446810000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(7)) x^2+(9.9250271120000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(9)) x^3-(6.6993933006000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(11)) x^4+(3.0817209182760000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(13)) x^5-(1.0259562557093850000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(15)) x^6+(2.5795471572121680000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(16)) x^7-(5.0559124281358492800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(17)) x^8+(7.9158224884955216000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(18)) x^9-(1.0092673672831790040000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(20)) x^10+(1.0643183145895342224000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(21)) x^11-(9.4014784455408856312000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(21)) x^12+(7.0294131146659544873280000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(22)) x^13-(4.4875271223090691593210000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(23)) x^14+(2.4637403808755673815880000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(24)) x^15-(1.1702766809158945062543000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(25)) x^16+(4.8332789236588336574094000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(25)) x^17-(1.7426655674747683575881670000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(26)) x^18+(5.5031544236045316555415800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(26)) x^19-(1.5258746356358019590365290000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(27)) x^20+(3.7214913473684776557868140000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(27)) x^21-(7.9927484619618440561784982500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(27)) x^22+(1.5123670133242584927169054080000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^23-(2.5206116888737641545281756800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^24+(3.6968971436815207599746576640000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^25-(4.7633097812819594407365781440000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^26+(5.3777341591995557415212453120000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^27-(5.3009093854967049452137989504000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^28+(4.5402672156089686405057009920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^29-(3.3579059615441330570406746920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^30+(2.1270020166965769510873481920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^31-(1.1416996119033096869807089560000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^32+(5.1203497745966616264589371360000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(27)) x^33-(1.8824815347781844214922563000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(27)) x^34+(5.5238840402757535148035320000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(26)) x^35-(1.2436814944480497679528420000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(26)) x^36+(2.0167808018076482723559600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(25)) x^37-(2.0963905703000554410015900000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(24)) x^38+(1.0488508169105968435280000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(23)) x^39];

If I want to evaluate this, I get the numerical error message: "CompiledFunction::cfn: Numerical error encountered at instruction 1; proceeding with uncompiled evaluation."

If I use the uncompiled version it gets strange:

f=-820.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000+447720.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 x-(9.1446810000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(7)) x^2+(9.9250271120000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(9)) x^3-(6.6993933006000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(11)) x^4+(3.0817209182760000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(13)) x^5-(1.0259562557093850000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(15)) x^6+(2.5795471572121680000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(16)) x^7-(5.0559124281358492800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(17)) x^8+(7.9158224884955216000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(18)) x^9-(1.0092673672831790040000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(20)) x^10+(1.0643183145895342224000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(21)) x^11-(9.4014784455408856312000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(21)) x^12+(7.0294131146659544873280000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(22)) x^13-(4.4875271223090691593210000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(23)) x^14+(2.4637403808755673815880000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(24)) x^15-(1.1702766809158945062543000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(25)) x^16+(4.8332789236588336574094000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(25)) x^17-(1.7426655674747683575881670000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(26)) x^18+(5.5031544236045316555415800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(26)) x^19-(1.5258746356358019590365290000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(27)) x^20+(3.7214913473684776557868140000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(27)) x^21-(7.9927484619618440561784982500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(27)) x^22+(1.5123670133242584927169054080000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^23-(2.5206116888737641545281756800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^24+(3.6968971436815207599746576640000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^25-(4.7633097812819594407365781440000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^26+(5.3777341591995557415212453120000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^27-(5.3009093854967049452137989504000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^28+(4.5402672156089686405057009920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^29-(3.3579059615441330570406746920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^30+(2.1270020166965769510873481920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^31-(1.1416996119033096869807089560000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(28)) x^32+(5.1203497745966616264589371360000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(27)) x^33-(1.8824815347781844214922563000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(27)) x^34+(5.5238840402757535148035320000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(26)) x^35-(1.2436814944480497679528420000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(26)) x^36+(2.0167808018076482723559600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(25)) x^37-(2.0963905703000554410015900000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(24)) x^38+(1.0488508169105968435280000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^(23)) x^39; 

f/.x->N[7/10,16]
f/.x->0.7

results in: output

Weirdly enough both evaluate to the same value if I delete all the zeros in the uncompiled function...

The function is highly oscilliating: plot of function I think what Mathematica does is to first try to evaluate the uncompiled function symbolically and then finds that the numerical error is very high which results in the "No significant digits to display" error...

My questions are the following:

1) I want the compiled function to evaluate no matter what, I prefer a false value to a numerical error message. How can I enforce that?

2) How can I force Mathematica to always evaluate the uncompiled function?

3) Are there some "rules" for compiled functions to prevent this. I found that when I use 

N[f,MachinePrecision]

and then copy/paste the result into the compiled function it works. But I would need to do that for thousands of functions...What does compile do if there is a body with numbers containing more digits than MachinePrecision can handly? Simply cut it off?

I really love Mathematica but sometimes it gets really confusing. If I plug the function into MATLAB it always evaluates for all values...its a different value though but at least I always get a value.

I really can not think of any better title than this, but I am trying to understand this and fix it for hours now and would really appreciate some help. Thanks a lot in advance!

€dit: So, the problem was that Compile[] in this case needs a body with MachinePrecision numbers. While something like

f=Compile[{x},
2934893242394239492349.239482934892394239422342342343*x
];
Precision[2934893242394239492349.239482934892394239422342342343]
f[0.333]

works perfectly fine even though the body contains numbers far higher than MachinePrecision in my case it does not. ...sometimes...for higher orders it would only work if the body has MachinePrecision. And also the code has to be copy-pasted so it would always contains the "`" after each number...If it didnt, it would also not work... If copypaste is not possible wrapping the body with Evaluate[N[...,MachinePrecision]] should do the trick. At least I hope so after some more testing... Thanks everyone for the help!

Improve this question
$\endgroup$
5
  • $\begingroup$ With these many digits, this number cannot be processed in machine precision. But compiled functions cannot handle arbitrary precision; they work internally with the C-type double (and its complex counterpart). So, N[f,MachinePrecision] rounds all numbers in the expression f to machine precision and this is why it works afterwards. $\endgroup$
    – Henrik Schumacher
    May 21, 2018 at 14:21
  • $\begingroup$ Btw.: Polynomials of degree 39 in FEM? No wonder that you experience problems. Already Runge found out that this would be a bad idea... $\endgroup$
    – Henrik Schumacher
    May 21, 2018 at 14:28
  • $\begingroup$ Related/duplicate: (3152) $\endgroup$
    – Michael E2
    May 21, 2018 at 14:48
  • $\begingroup$ Thank you, I know that Compile only works with MachinePrecision, but I wonder why it doesn't "accept" such input and simply cuts of or something like that... It btw does not work when I round to 16... I guess I have to rewrite all my functions, so that all input is in MachinePrecision. FEM works (although still limited) with such high polynomials if you use (among other methods) a special kind of node distribution (such as Gauss-Legendre-Lobatto nodes or Chebychev-Gauss-Legendre nodes). $\endgroup$
    – John P. Wiener
    May 21, 2018 at 16:16
  • 1
    $\begingroup$ For item (1), could use fc = Compile[{x}, Evaluate[N[f]]]; $\endgroup$
    – Daniel Lichtblau
    May 21, 2018 at 17:59

1 Answer 1

Reset to default
6
$\begingroup$

As @Henrik says, those numbers with all the zeros are not machine numbers. They are arbitrary precision numbers with high precision. For example:

number = -820.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000;
Precision[number]
number //InputForm

127.914

-820.`127.91381385238373

So, the first number has a precision of almost 128. Compile does not work with arbitrary precision numbers.

Regarding the uncompiled f. If you want a correct answer, you need to use arbitrary precision input with high precision, since the subtractive cancellation that occurs is substantial. Compare:

f /. x->.7
f /. x->.7`100

-3.22123*10^9

-0.1048536426633009587564239155085961687000000000000000000000000000000000000

The MachinePrecision calculation is wrong by 10 orders of magnitude. Using high precision will also give you a better plot:

Plot[f, {x, 0.6, .7}, WorkingPrecision->100]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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