-4

I am trying to invert the exponent of a long double.

Suppose x = 3.5e1356. I want x to be 3.5e-1356.

I have this code:

long double x = 3.5e1356L;
int exponent;
long double fraction = frexpl(x, &exponent);

// recreate number with an inverted exponent
long double newNumber =  ldexpl(fraction, -exponent);

after this code, newNumber is 1.14732677619641872902e-1357

that has nothing to do with the original number.

What am I missing?

Duck
  • 34,902
  • 47
  • 248
  • 470
  • 5
    The exponent is in binary, not decimal. – user2357112 Jul 15 '17 at 04:57
  • Try looking at the (implementation defined) values of `LDBL_MAX` and `LDBL_MIN` (largest and smallest positive values a `long double` can represent, respectively). Even "typical" values of `LDBL_MAX` are no more than about `1E308`, and your input value is a lot larger than that - so is probably overflowing. Once that happens, the behaviour is undefined. Basic thing to remember: floating point representations cannot represent arbitrary values of any magnitude. – Peter Jul 15 '17 at 04:57
  • 1
    @peter: on intel hardware, long double has a 16-bit exponent, which allows for much larger numbers than 1E308. – rici Jul 15 '17 at 05:00
  • I am running on MacOS. There, `LDBL_MAX` is reported to be `1.18973149535723176502E+4932` and `LDBL_MIN` to be `3.36210314311209350626E-4932`, so I am far from overflowing... if my code is wrong how do I do that? thanks – Duck Jul 15 '17 at 05:01
  • @rici - before my comment, the OP had given no information whatsoever about the target platform or compiler. The values being played with are WAY outside the range that is required of that type. – Peter Jul 15 '17 at 05:18
  • 1
    @peter: floating point requirements are very conservative. The values in the OP are certainly well beyond the minimum requirements. But the 80-bit long double is a very common implementation (which dates back to the 80387 FP coprocessor), and it is a reasonable guess that someone who doesn't feel the need to explain that they are using unusual hardware is in fact on a platform with 80-bit long doubles. – rici Jul 15 '17 at 05:38
  • @rici - I'd hardly agree with that. Ignorance of the need to provide information when asking questions is hardware agnostic. In any event, I simply asked him to check what the range of values were, and explained why I asked that. – Peter Jul 15 '17 at 05:46
  • 1
    I have to ask - why? Is your lab, is 2712 powers of 10 within the bounds of experimental error? – Martin James Jul 15 '17 at 06:01
  • I'll echo @MartinJames's question: why? This seems like a rather unnatural operation. What are you using it for? Maybe there's a better way to achieve whatever it is you're trying to achieve. – Mark Dickinson Jul 15 '17 at 08:50

3 Answers3

3

You've inverted the exponent, but the exponent was never 1356 in the first place. The exponent of a long double is its binary exponent.

You may have written 3.5 * 10^1356, but internally, the computer stores it as something * 2^something else. What you've done is produce something * 2^-something else, not 3.5 * 10^-1356.*

If you want to get 3.5 * 10^-1356, you'll probably need to implement it in terms of base 10 logarithms.

*frexpl uses a different normalization convention than what your computer is probably using, so it's not quite inverting the machine-level exponent, but it is inverting a binary exponent instead of the decimal one you wrote.

user2357112
  • 260,549
  • 28
  • 431
  • 505
  • Ok, I understand it now but please use the answer to show the solution to the problem, not to point what is wrong. Thanks. I can imagine a lame solution where I create a loop dividing x by 10 and counting the number of divisions until the remainder is >= 1 and < 10. if there is something better, please explain. – Duck Jul 15 '17 at 05:05
  • 1
    @SpaceDog: Are you aware of [`log10l`](http://en.cppreference.com/w/c/numeric/math/log10)? – user2357112 Jul 15 '17 at 05:16
  • ok I can use that to discover the exponent but see my answer... I am not seeing how this can help there. Sorry, my brain is not seeing... – Duck Jul 15 '17 at 05:32
1

Use the fact that

a * (10^n) * (10^m) = a * 10^(n+m)

If you calculate m so that it is -2n you'll get:

a * (10^n) * (10^-2n) = a * 10^(n+(-2n)) =  a * 10^-n

In other words - just multiply the original number with 10^-2n

I would try something like this:

#include <stdio.h>
#include <math.h>

int main(void) {
    long double x = 8.9e-100;
    long double y = 0.0;
    int exp10;

    if (x != 0.0)
    {
        exp10 = log10l(fabsl(x));
        exp10 = (exp10 >= 0) ? -2*exp10 : -2*(exp10-1);
        y = x * powl(10, exp10);
    }

    printf("%.36Le\n", x);
    printf("%.36Le\n", y);

    return 0;
}

Example for 3.5 * 10^12 :

3.5 * 10^12 * 10^-24 --> 3.5 * 10^-12
Support Ukraine
  • 42,271
  • 4
  • 38
  • 63
  • 1
    PERFECT, my loop solution works well too but probably loses precision as the loop goes, as pointed by Mateo. The approach you have used on this line `exp10 = (exp10 >= 0) ? -2*exp10 : -2*(exp10-1);`, multiplying the whole thing by the opposite, is interesting. THANKS – Duck Jul 16 '17 at 01:50
  • Maybe not perfect. Corner case: `exp10 = log10l(fabsl(x));` has trouble when `x == 0.0`. – chux - Reinstate Monica Jul 16 '17 at 20:56
  • `exp10 = (exp10 >= 0) ? -2*exp10 : -2*(exp10-1); y = x * powl(10, exp10);` invites unnecessary overflow.. Could use 2 steps: `exp10 = (exp10 >= 0) ? -exp10 : -(exp10-1); y = x * powl(10, exp10) * powl(10, exp10);` – chux - Reinstate Monica Jul 16 '17 at 20:59
0

Based on the help from you guys, I have created this code that is working but probably can be improved...

  NSInteger exponent = 0;
  long double remainder = x;

  // stop the loop if remainder is >= 1 and less than 10
  while (!((fabsl(remainder) >= 1.0L) && (fabsl(remainder) < 10.0L))) {
    remainder /= 10.0L;
    exponent ++;
  }

  long double finalNumber =  remainder * powl(10.0L, -exponent);
Duck
  • 34,902
  • 47
  • 248
  • 470