Pagina's

2011/12/29

Factorial by binary splitting, part 2

using Xint = System.Numerics.BigInteger;
using System;
class Factorial
{
    public static Xint F(uint n)
    {
        uint a = 0;
        uint s = 0;
        Xint P = 1;
        Xint Q = 1;
        uint b = 1;
        for (int i = fL2((int)(n / 2)); i >= 0; i--)
        {
            a = n >> i;
            s = s + a / 2;
            a = a - 1 | 1;
            P = Q * P;
            Q = OddP(a, b) * Q;
            b = a + 2;
        }
        return Q * P << (int)s;
    }

    private static Xint OddP(uint a, uint b)
    {
        if (a == b) return a;
        uint m = (a + b) / 2;
        m += m & 1;
        return OddP(a, m + 1) * OddP(m - 1, b);
    }

    private static int fL2(int n) { int i = -1; for (; n > 0; n /= 2) i++; return i; }

    // 42!=42*41*..
    //    =41*39*...*3*1 * 42*40*...*4*2 
    //    =41,1? * 42,2?
    //             42,2?=21! *                                                                  2^(42/2)
    //                   21!=21*19*...*3*1 * 20*18*...*4*2                                               
    //                      =21,1? * 20,2?                                                               
    //                               20,2?=10! *                                                2^(20/2)
    //                                     10!=9*7*5*3*1 * 10*8*6*4*2                                
    //                                        =9,1? * 10,2?                                          
    //                                                10,2?=5! *                                2^(10/2)
    //                                                      5!=5*3*1 * 4*2                            
    //                                                        =5,1?  * 4,2?                                  
    //                                                                 4,2?=2! *                2^( 4/2)
    //                                                                      2!=1*1  * 2                
    //                                                                        =1,1? * 2,2?                          
    //                                                                                2,2?=1! * 2^( 2/2)
    //        
    //    = 41,1?  * 21,1?  * 9,1?     * 5,1?   * 2^(42/2+20/2+10/2+4/2+2/2)
    //    = 5,1?   * 9,1?   * 21,1?    * 41,1?  * 2^(21+10+5+2+1)
    //    = 5,1?^4 * 9,7?^3 * 21,11?^2 * 41,23? * 2^39
    //    =    a^4 *    b^3 *      c^2 *    d^1 << 39
    //
    //    = 1                    *  a      
    //    = a^1                  *  a*b    
    //    = a^2 * b^1            *  a*b*c  
    //    = a^3 * b^2 * c^1      *  a*b*c*d
    //    = a^4 * b^3 * c^2* d^1                << 39

    static void Main()
    {
        Console.WriteLine(F(0));
        Console.WriteLine(F(100));
        Console.ReadLine();
    }
}