Euler number with tasks

Question

I want to calculate Euler number with multiple threads using this formula =∑((3k) ^ 2 + 1)/(3k)! , k =0,... ,∞ , but I am not getting the right results so far and once of the problems is that when I use fairly big number I am going out of range for decimal for factorial function here is what I've done so far

static void Main(string[] args)
{
   Console.WriteLine(Program.Calculate(5, 1));
}
public static decimal Calculate(int x, byte taskNumber)
{
   var tasks = new List<Task<decimal>>();

   for (int i = 0; i < x; i += (x / taskNumber))
   {
        int step = i;
        tasks.Add(Task.Run(() =>
        {
            int right = (step + x / taskNumber) > x ? x : (step + x / taskNumber);
            return ChunkE(step + 1, right);
        }));
    }

    Task.WaitAll(tasks.ToArray());

    return tasks.Select(t => t.Result).Aggregate(((i, next) => i + next));
    }

then I have simple factorial and euler functions

public static decimal ChunkFactorial(int left, int right)
{
    //Console.WriteLine("ChunkFactorial Thread ID :" + Thread.CurrentThread.ManagedThreadId);
    if (left == right)
    {
        return left == 0 ? 1 : left;
    }
    else
    {
        return right * ChunkFactorial(left, right - 1);
    }
}

public static decimal ChunkE(int left, int right)
{
    if(left == right)
    {
        return left == 0 ? 1 : left;
    }
    else
    {
        return ((3 * right) * (3 * right) + 1) / ChunkFactorial(left, right) + ChunkE(left, right - 1);
    }
}

What I want to achieve is calculating Euler number up until x precision using different amount of tasks.

What I am getting with this call is 41.01666..7 if I increase x decimal will eventually overflow. How can I fix this problem I tried using BigInteger but then it begin to be mess and I loose precision of the result.. Any ideas?

Also when I start the program with 1 task I get one result and when I start program with 4(or different of 1) I get different result I don't know what I am missing..

Answer 1

If you are allowed to transform the terms before implementation, consider that

(3k)^2/(3k)! = (3k)/(3k-1)! = 1/(3k-2)!+1/(3k-1)!

which then proves that your formula really computes the Euler number.

But you need to include the factor 3 in (3k)! in the factorial computations.

It should be perfectly fine to use the floating point type for the factorial computation, since the required precision should be reached long before overflow happens. Note that an error bound is twice the next term not in the partial sum.

Answer 2

As the main point of this question was already answered I just add some stuff:

You can do few upgrades to your code like while computing factorial each time you can divide both divident and divider by 2 or 5 do it to lower the bits needed. This will also boost the speed a lot if done right.

But in the end Your formula will still take forever to converge to e! Not to mention the overflows due to factorial. I am using something better suited for computers (not sure where the formula comes from though it was ages ago)...

e = (1+1/x)^x

where x -> +inf This has many advantages when using binary representations of numbers. I refined the process by using powers of 2 for x which simplifies things a lot... My computation code (based on my arbnum class) is as this:

    arbnum c,x;
    int bit=512;        // min(int_bits,fract_bits)/2 ... this is just remnant from fixed point code where bitwidth matters
    // e=(1+1/x)^x  ... x -> +inf
    c.one(); c>>=bit; c++;  // c = 1.000...0001b = (1+1/x)          = 2^-bits + 1
//  x.one(); x<<=bit;       // x = 1000...000.0b =    x   = 1/(c-1) = 2^+bits
        for (;bit;bit--)        // c = c^x = c^(2^bits) = e
            {
            c*=c;
            c._normalize(2048); // this just cut off the result to specific number of fractional bits only to speed up the computation instead you should cut of only last zeros !!!
            }

As you can see I am computing target precision from the start (bits) and truncating the result (to manageable bit width for fractional part). If you got fixed point arithmetics then you do not need to do this at all (this was just my fast attempt to port it from my ancient fixed point code to my new arbnum class).

So set the bits constant to what you want and also the truncating size. Both should be based on your target precision. As you can see this is not iterative process... The only thing that is in loop is the power. It is optimized a bit to understand it you need to realize what are you computing:

(1.000...0001b) ^ (1<<bits)

so I just square the c until first signed bit in x is hit. Beware each squaring doubles the needed fractional bit width ... that is the reason for truncation (it lowers accuracy bat boosts performance much much more)

As you can see this approach is quite nice as it does not need any division ... only bit operations and multiplication.

Here comparison:

[e]
reference          2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274
my e               2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274
00000200 method 0: 2.7182818280709941626033162692782994930797533824790948079298224728031965377356362865659816488677135209

The reference is first 100 digits of e. The my e is result from this code and method 0 is result from your equation after 200 iterations in arbitrary precision with last zero truncating and the %2,%5 optimization to speed it up.

Mine approach took just few milliseconds and yours ~20 seconds to get to that point...

Here the same computed in decadic base directly on strings:

• Calculating the Euler's Number gives a limited number of decimals