Splitting a floating point number as sums of floating point of fixed precision

Question

Suppose i have an algorithm by which i can compute an infinitely precise floating point number (depending from a parameter N) lets say in pseudocode:

arbitrary_precision_float f = computeValue(n); //it could be a function which compute a specific value, like PI for instance.

I guess i can implement computeValue(int) with the library mpf of the gnump library for example...

Anyway how can i split such number in sums of floating point number where each number has L Mantissa digits?

//example
f = x1 + x2 + ... + xn;
/*
for i = 1:n
  xi = 2^ei * Mi
 Mi has exactly p digits.
*/

I don't know if i'm clear but i'm looking for something "simple".

Answer 1

You can use a very simple algorithm. Assume without loss of generality that the exponent of your original number is zero; if it's not, then you just add that exponent to all the exponents of the answer.

Split your number f into groups of L digits and treat each group as a separate xi. Any such group can be represented in the form you need: the mantissa will be exactly that group, and the exponent will be negated start position of the group in the original number (that is, i*L, where i is the group number).

If any of the resulting xis starts from zero, you just shift its mantissa correcting the exponent correspondingly.

For example, for L=4

f = 10010011100
    1001     
        0011 
            100
-> x1=1.001 *2^0
   x2=0.011 *2^{-4} = 1.1*2^{-6}
   x3=1.00  *2^{-8}

---

Another question arises if you want to minimize the amount of numbers you get. In the example above, two numbers are sufficient: 1.001*2^0+1.11*2^{-6}. This is a separate question, and in fact is a simple problem for dynamic programming.