3

I found myself with the need to compute the exponential of a large number, e. g.exp(709). Such a number would be represented, in floating point precision, as 8.2184074615549724e+307.

It seems that numbers with exponents larger than that would be simply translated into Inf, which creates problems in my code. I can only guess that things can be fixed using more bits to represent the exponent, but I am not aware of a pragmatical way to proceed.

Here is a code snippet:

double expon = exp(500); /*here I also tried `long double`, with no effect */
printf("%e\n", expon ); /*gives INF*/
double Wa = LambertW<0>( expon); /*gives error, as it can't handle inf*/

Is there a way to compute this?

This problem has been debated in general, but I did not find an useful answer. Also, it seems that GCC supports multiple-precision floating-point arithmetics since version 4.3. How does it help?

Edit: The suggested possible-duplicate questions turned out irrelevant because as I need huge decimals, not exact decimals. This is not a duplicate.

altroware
  • 940
  • 3
  • 13
  • 22
  • 3
    Whats with big/unlimited integer arithmetic? See https://gmplib.org/ – Superlokkus Oct 29 '15 at 14:39
  • See also: [List of arbitrary-precision arithmetic software](https://en.wikipedia.org/wiki/List_of_arbitrary-precision_arithmetic_software) – Ivan Aksamentov - Drop Oct 29 '15 at 14:42
  • There is no language C/C++. Your code is not C. – too honest for this site Oct 29 '15 at 14:56
  • 1
    Use logarithmic. If you're working on something with this big numbers, logs are your friend. For manipulation with exp, log, ln, etc. see wiki https://en.wikipedia.org/wiki/Natural_logarithm – Anders Schou Oct 29 '15 at 14:58
  • 1
    I looked at the definition of LambertW (never heard of it before) so maybe I'm missing the point, but: You seem to want `y=W(exp(x))` where `exp(x)` is too big. IIUC, that is the y which solves `x==y+ln(y)`. Wouldn't it be easier to do a crude successive approximation computation of that y directly, rather than try to extend LambertW to support larger numbers? – JSF Oct 29 '15 at 16:00
  • @Olaf I used the `c` tag as I would accept a solution in both `c` or `c++` . Also, I edited the question to point out that I don't really need to improve the significang, but only the exponent of the floating-point representation. – altroware Oct 29 '15 at 16:01
  • Assuming I understood LambertW correctly, W(exp(x)) is super easy by crude successive approximation when x is large. You need a better method only when X is small. Start at `y=x;` then repeat `oldy=y; y=x-log(y);` until `y==oldy` – JSF Oct 29 '15 at 16:14
  • @JSF Very interesting! Actually I need something like `W[A*exp(x)]`, but your argument seems equally good. I guess that @Anders Schous meant the same thing. Let us also wait for a solution in terms of floating-point representation – altroware Oct 29 '15 at 16:17
  • Then how big is `x+log(A)`? If that is big, use it instead of x in my crude successive approximation plan. If it is not big, then you have a way to use LambertW (if ln(A) is negative, you would want `W(exp(x+log(A)))` but if ln(A) were positive either you have my solution or you wouldn't have a problem with the original exp). – JSF Oct 29 '15 at 16:20
  • @JSF well, `log(A)` is (very much) positive... I'm looking at you previous solution :) . – altroware Oct 29 '15 at 16:35
  • As a quick comparison, I found a simpler version of LambertW online using the newton method. For W(exp(100)) my version converges much faster. For W(exp(3)) my version converges much slower. (In both cases the two methods converge to the same answer). For W(exp(2)) my version oscillates near convergence but never lands. For W(exp(500)) my version works :) – JSF Oct 29 '15 at 17:11
  • The precision tags were misleading and may have contributed to the (I think incorrect) designation of your question as a duplicate. Even before your edit, it was clear you thought you wanted a bigger exponent, not more precision. If you dig past the popular answers in the linked question (all of which is about precision) you could find the bigger exponent answer (rather well hidden). As it turns out, I don't think you really wanted any of that. – JSF Oct 29 '15 at 17:31
  • @JSP I tested your solution with the Newton method, and I found it robust enough for my purposes. Thanks! Anyway, I'm also trying to learn how to handle large numbers. Is it possible to remove the tag [duplicate]? – altroware Oct 30 '15 at 14:34
  • @Drop, The issue of removing the duplicate tag (because the linked question barely relates to the topic of larger exponents and covers only precision) should be addressed to those who added the tag. In my opinion a clearer new question (without the sidetrack of LambertW) would be appropriate: **Which libraries extending floating point support bigger exponents?** Following each link in answers to the other question, I found a hint at "yes" only for http://www.boost.org/doc/libs/1_57_0/libs/multiprecision/doc/html/boost_multiprecision/tut/floats/cpp_dec_float.html – JSF Oct 31 '15 at 12:15
  • @JSF Libraries on linked wiki provide arbitrary significand and exponent sizes (provided you have enough storage space). Example for GMP: [link](https://gmplib.org/list-archives/gmp-discuss/2003-December/000906.html). I don't understand what you are trying to say and why you are addressing your last message to me. See also: [What Every Computer Scientist Should Know About Floating-Point Arithmetic](https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) – Ivan Aksamentov - Drop Oct 31 '15 at 17:57
  • @Drop, I addressed you in an attempt to help altroware with questioning the duplicate status. The meta posts I found on that topic say to "ping" the users who marked it as duplicate, but don't say how. I assumed altroware knows less how to do that than I do. Unless I'm badly misunderstanding, your link (like many in the original question) shows how to get a very high pow using extreme precision. But the question here was how to get a floating point format whose exponent part is larger, while still keeping reasonable precision – JSF Oct 31 '15 at 20:28
  • @Drop, your "see also" was an enormous amount of information **zero** of which was relevant to the question. That kind of link is very disappointing to someone looking for an answer and wasting lots of time reading something not so far away from topic that it can be rejected out of hand, but ultimately useless. – JSF Oct 31 '15 at 20:33
  • @Drop, I saw that at least one of the libraries in the links from the original question **does** provide for larger exponents (but that one looks like a bad choice for other reasons). At least two provide just unlimited precision and **not** extended exponent size. Those details are not at all evident from the "answers" nor from the pages directly linked. So considering that an answer is less helpful than "It exists somewhere, google it". It does exist, but an answer would be a clearer list of choices. – JSF Oct 31 '15 at 20:40
  • @JSF I still don't understand what you want to say. Link with GMP example shows how to get a number as a result of raising any base to **any** exponent. Size of the exponent is hidden there, you don't need to care about it. Second link explains how IEEE-754 floating point numbers work, just in case. If you still have something to add here I can't help much other than to vote for reopening the question. – Ivan Aksamentov - Drop Oct 31 '15 at 23:16

1 Answers1

4

You should be able to perform your computation with adequate precision using long double arithmetic:

The maximum value for 80 bit long double is 1.18×10^4932, much larger than e^709.

In order for the computation to be performed as long double, your must use expl instead if exp:

long double expon = expl(500);
printf("%Le\n", expon);

The LambertW function will handle the long double if it is properly overloaded for this type, otherwise expon will be converted to double and produce inf and the computation will fail as you mentioned.

I don't know which implementation of Lambert W function you use, Darko Veberic's does not support long double arguments, but it might be feasible to extend the implementation to type long double as it is available in source form: https://github.com/DarkoVeberic/LambertW . You might want to contact him directly.

Another approach is to consider that exp(709) is just too close to the maximum precision of the double type, 10^308. If you can alter your computation using just smaller exponents and a different formula, the computation might be done with regular double types.

chqrlie
  • 131,814
  • 10
  • 121
  • 189