Cube Root

/*


      ___
  \3 /                           Input: integer x >= 0
   \/  x = y                    Output: integer y, such that y^3 <= x < (y+1)^3
   
  
Cube Roots are computed three to five times faster with the CR function,
compared to the Nth Root function.
For small values it's obvious to use: " y = (uint)Math.Pow(x, 1d / 3)) ".
It takes ~125 ns, but for example with x = 4503569204744003 (a 52 bits number),
it returns 165139, wrong, it should be 165140. So errors have to be corrected.
Another option: the "Integer Cube Root" algorithm from "Hacker's Delight" (see Refs).
For a uint it takes ~25 ns, not too bad for a C# version, for a ulong ~330 ns.
Up to ulong.MaxValue there are ~2.5 million perfect cubes, they can be generated like below.
  
3th powers:
 
    |-----|------|------|------|------|         dy=growth of y, ddy=growth of dy, ... 
    |  x  | y=x^3|  dy  | ddy  | dddy |          
    |-----|------|------|------|------|         y=x^3 , y'=3*x^2 , y''=6*x , y'''=6
    |  0  |   0  |      |   0  |      |               
    |     |      |   1  |      |   6  |         y=x^n ,,,,,,,,,,,,,, y''''...=n! A000142
    |  1  |   1  |      |   6  |      |         
    |     |      |   7  |      |   6  |             y: The cubes                 A000578
    |  2  |   8  |      |  12  |      |            dy: Central hexagonal numbers A003215
    |     |      |  19  |      |   6  |           ddy: Multiples of 6            A008588
    |  3  |  27  |      |  18  |      |          dddy: The  six sequence         A010722
    |     |      |  37  |      |   6  |         ddddy: The zero sequence         A000004    
    |  4  |  64  |      |  24  |      |     
    |-----|------|------|------|------|    y=x^4 : A000583 , A010863 , A101103 , A005914 , A005917    
 
 
    private static void abc()
    {
        uint n = 0, a = 0, b = 1, c = 6;
        Console.WriteLine("// {0,2}{1,4}{2,4}{3,3}", n, a, b, c);         //  0   0   1  6
        while (n < 7)                                                     //  1   1   7 12
        {                                                                 //  2   8  19 18
            n++;                                                          //  3  27  37 24
            a += b;                                                       //  4  64  61 30
            b += c;                                                       //  5 125  91 36
            c += 6;                                                       //  6 216 127 42
            Console.WriteLine("// {0,2}{1,4}{2,4}{3,3}", n, a, b, c);     //  7 343 169 48
        }
        Console.ReadLine();
    }

4th powers:

    private static void abcd()
    {
        uint n = 0;
        uint a = 0, b = 1, c = 14, d = 36;
        Console.WriteLine("// {0,2}{1,5}{2,5}{3,4}{4,4}", n, a, b, c, d);         //  0    0    1  14  36
        while (n < 8)                                                             //  1    1   15  50  60
        {                                                                         //  2   16   65 110  84
            n++;                                                                  //  3   81  175 194 108
            a += b;                                                               //  4  256  369 302 132
            b += c;                                                               //  5  625  671 434 156
            c += d;                                                               //  6 1296 1105 590 180
            d += 24;                                                              //  7 2401 1695 770 204
            Console.WriteLine("// {0,2}{1,5}{2,5}{3,4}{4,4}", n, a, b, c, d);     //  8 4096 2465 974 228
        }
        Console.ReadLine();
    }

The initialisation value's: "a = 0, b = 1, c = 14, d = 36, 24" : A019538
It's getting a bit off topic, hence another topic: Power Addition Only.
Above leads to cro12 / CR12, "Hacker's Delight" hardware algorithms, and
a new "cro12", while writing ... this. It's ~1 ns faster than the one used.
 
    private static uint cro12new(uint x)     //   |------|----|
    {                                        //   |    x | ns |
        uint y = 0, a = 0, b = 1, c = 0;     //   |------|----|
        while (a < x)                        //   |    0 |  4 |
        {                                    //   |    1 |  4 |
            y++;                             //   |    7 |  5 |
            b += c;                          //   |   63 |  7 |
            a += b;                          //   |  255 | 10 |
            c += 6;                          //   | 1000 | 13 |
        }                                    //   | 1023 | 15 |
        if (a != x) y--;                     //   | 4095 | 24 |
        return y;                            //   |------|----|
    }
                                       

    |----------------------------------|
    | Times in ns (Athlon X4, XP, 2GB) |
    |----------|---------|-------------|
    |          |         | output type |
    | function |  input  |------|------|
    |          |         | uint | Xint |
    |----------|---------|------|------|
    |   cro12  |  < 4096 | <=25 |  <60 |  
    |   cro32  |  <= ~0u |  ~25 |  ~80 |  
    |   cro64  | <= ~0uL | ~330 | ~360 |
    |          |         |      |      |
    |    CR12  |  < 4096 | <=35 | <130 |
    |    CR32  |  <= ~0u |  ~35 | ~135 |
    |    CR64  | <= ~0uL | ~340 | ~460 | 
    |----------------------------------|
 
  
Abbreviations: 
  ref  reference
   CR  Cube Root & Remainder CR(9)={2,1} 
  CRO  Cube Root Only       CRO(9)= 2
  RND  Random number
    ~  bitwise not operator, or depending 
       on context, approximately
  ~0u  uint.MaxValue
 ~0uL  ulong.MaxValue
   ns  nanosecond  (10^-9 s)
   us  microsecond (10^-6 s)
   ms  millisecond (10^-3 s)
  
Times:  
    |------------------------------------------------------------------------|
    |                             Athlon X4, XP, 2GB                         |
    |--------|--------|--------||----------|---------||------------|---------|               
    |        |   ns   |   ns   ||          |   us    ||            |   ms    |
    |      X | CRO(X) |  CR(X) || bits RND | CR(RND) ||  bits RND  | CR(RND) |
    |--------|--------|--------||----------|---------||------------|---------|
    |      0 |    35  |    97  ||      65  |      4  ||    100.000 |     13  | 
    |      1 |    36  |    98  ||     100  |      8  ||    200.000 |     31  | 
    |    100 |    42  |   103  ||     200  |     14  ||    500.000 |     92  | 
    |   1000 |    47  |   112  ||     500  |     27  ||  1.000.000 |    237  | 
    |   4095 |    59  |   124  ||    1000  |     38  ||  2.000.000 |    620  | 
    |   4096 |    66  |   129  ||    2000  |     59  ||  5.000.000 |   2330  | 
    |  ~0u   |    83  |   142  ||    5000  |    134  || 10.000.000 |   6250  | 
    |  ~0u+1 |   329  |   395  ||   10000  |    341  || 20.000.000 |  16500  | 
    | ~0uL   |   382  |   474  ||   20000  |   1060  ||            |         |
    | ~0uL+1 |  3870  |  3820  ||   50000  |   4360  ||            |         |
    |        |        |        ||  100000  |  12400  ||            |         |
    |--------|--------|--------||----------|---------||------------|---------|

Refs: Hacker's Delight   
      Henry S. Warren, Jr.   
      ISBN 0-201-91465-4   
      8th Printing February 2008   
      http://www.hackersdelight.org   
      (11-2 Integer Cube Root)
 
      Modern Computer Arithmetic    
      Richard Brent and Paul Zimmermann   
      version 0.1, October 2006   
      http://www.loria.fr/~zimmerma/mca/mca-0.1.pdf
      (1.5.2 k-th Root (Cube Root))   
 
*/
using System;
using System.Diagnostics;
using System.Threading.Tasks;
using Xint = System.Numerics.BigInteger;

class Cube_Root
{
    private static Xint[] CR(Xint D)
    {
        if (D <= ~0uL) return D < 4096 ? CR12((uint)D) : D <= ~0u ? CR32((uint)D) : CR64((ulong)D);
        int n = bL(D) / 6, m = 2 * n;
        Xint BA = D & ((Xint.One << m) - 1);                                  //  |-------D------|
        D >>= m;                                                              //  |              |
        Xint C = D & ((Xint.One << n) - 1);                                   //  |-D-|-C-|--BA--|
        D >>= n;
        Xint[] R = CR(D);                                                     //  R[0] = CubeRoot
        Xint R0R0 = SQ(R[0]);
        Xint[] Q = DQR(C + (R[1] << n), 3 * R0R0);                            //  Q[0] = Quotient
        Xint Q0Q0 = SQ(Q[0]);
        R[1] = BA + (Q[1] << m) - MTP(Q0Q0, Q[0] + ((3 * R[0]) << n));
        if (R[1] < 0)
        {
            R0R0 <<= m;
            R0R0 += Q0Q0 + (MTP(Q[0], R[0]) << (n + 1));
            R[0] <<= n;
            R[0] += Q[0];
            R[1] += 1 + 3 * (R0R0 - R[0]);
            R[0] -= 1;
            while (R[1] < 0)
            {
                R0R0 -= 1 + 2 * R[0];
                R[1] += 1 + 3 * (R0R0 - R[0]);
                R[0] -= 1;
            }
        }
        else
        {
            R[0] <<= n;
            R[0] += Q[0];
        }
        return R;
    }
    private static Xint[] CR12(uint x)
    {
        uint y = 0, z = 0, r = 0, s = 0;
        while (r < x)
        {
            s = r;
            z += y * 6;
            r += z + 1;
            y += 1;
        }
        return r == x ?
            new Xint[] { y, 0 } :          // no need to use s, replace "x - s" by "x - (r + z + 1)",
            new Xint[] { y - 1, x - s };   // but it's much slower, why?
    }
    private static Xint[] CR32(uint x)
    {
        uint y = 0, z = 0, b = 0;
        int s = 30;
        while (s >= 0)
        {
            y *= 2;
            z *= 4;
            b = 3 * y + 3 * z + 1 << s;
            s -= 3;
            if (x >= b)
            {
                x -= b;
                z += 2 * y + 1;
                y += 1;
            }
        }
        return new Xint[] { y, x };
    }
    private static Xint[] CR64(ulong x)
    {
        uint y = 0;
        ulong z = 0, b = 0, bs = 0;
        int s = 63;
        while (s >= 0)
        {
            y *= 2;
            z *= 4;
            b = 3 * y + 3 * z + 1;
            bs = b << s;
            if (x >= bs && b == bs >> s)
            {
                x -= bs;
                z += 2 * y + 1;
                y += 1;
            }
            s -= 3;
        }
        return new Xint[] { y, x };
    }
    private static Xint CRO(Xint X)
    {
        if (X <= ~0ul) return X < 4096 ? cro12((uint)X) : X <= ~0u ? cro32((uint)X) : cro64((ulong)X);
        Xint[] R = CR(X);
        return R[0];
    }
    private static uint cro12(uint x)
    {
        uint y = 0, z = 0, r = 0;
        while (r < x)
        {
            z += y * 6;
            r += z + 1;
            y += 1;                    //  y > 1625 causes overflow of r
        }
        return r == x ? y : y - 1;
    }
    private static uint cro32(uint x)
    {
        uint y = 0, z = 0, b = 0;
        int s = 30;
        while (s >= 0)
        {
            y *= 2;
            z *= 4;
            b = 3 * y + 3 * z + 1 << s;
            s -= 3;
            if (x >= b)
            {
                x -= b;
                z += 2 * y + 1;
                y += 1;
            }
        }
        return y;
    }
    private static uint cro64(ulong x)
    {
        uint y = 0;
        ulong z = 0, b = 0, bs = 0;
        int s = 63;
        while (s >= 0)
        {
            y *= 2;
            z *= 4;
            b = 3 * y + 3 * z + 1;
            bs = b << s;
            if (x >= bs && b == bs >> s)
            {
                x -= bs;
                z += 2 * y + 1;
                y += 1;
            }
            s -= 3;
        }
        return y;
    }

    // Use the faster versions from: Nth Root ////////////////////////////////////////////////////
    private static Xint SQ(Xint U) { return U * U; }                                            //
    private static Xint MTP(Xint U, Xint V) { return U * V; }                                   //
    private static Xint[] DQR(Xint U, Xint V)                                                   //
    {                                                                                           //
        Xint[] QR = new Xint[2];                                                                //
        QR[0] = Xint.DivRem(U, V, out QR[1]);                                                   //
        return QR;                                                                              //
    }                                                                                           //
    // Use the faster versions: http://www.bigintegers.blogspot.com/2013/07/nth-root-power.html //

    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        for (int i = 1; i < 100; i++)
        {
            Xint X = RND(8 * i);
            //X = Xint.Pow(X, 3);
            //X--;
            //X++;
            sw.Start(); Xint[] R = CR(X); sw.Stop();
            if (Xint.Pow(R[0], 3) > X) Console.WriteLine("WRONG1");
            if (Xint.Pow(R[0] + 1, 3) <= X) Console.WriteLine("WRONG2");
            if (Xint.Pow(R[0], 3) + R[1] != X) Console.WriteLine("WRONG3");
        }
        Console.WriteLine(sw.ElapsedMilliseconds + " ms");
        Console.ReadLine();
    }

    private static int bL(Xint U)
    {
        byte[] bytes = (U.Sign * U).ToByteArray();
        int i = bytes.Length - 1;
        return i * 8 | bitLengthMostSignificantByte(bytes[i]);
    }
    private static int bitLengthMostSignificantByte(byte b)
    {
        return b < 08 ? b < 02 ? b < 01 ? 0 : 1 :
                                 b < 04 ? 2 : 3 :
                        b < 32 ? b < 16 ? 4 : 5 :
                                 b < 64 ? 6 : 7;
    }
    private static int seed;
    private static Xint RND(int n)
    {
        if (n < 2) return n;
        if (seed == int.MaxValue) seed = 0; else seed++;
        Random rand = new Random(seed);
        byte[] bytes = new byte[(n + 15) >> 3];
        rand.NextBytes(bytes);
        int i = bytes.Length - 1;
        bytes[i] = 0;
        n = (i << 3) - n;
        i--;
        bytes[i] >>= n;
        bytes[i] |= (byte)(128 >> n);
        return new Xint(bytes);
    }
    private static Xint[] CRold(Xint D)
    {
        if (D < 64)
        {
            int d = (int)D;
            if (d >= 27) return new Xint[] { 3, d - 27 };
            if (d >= 08) return new Xint[] { 2, d - 08 };
            if (d >= 01) return new Xint[] { 1, d - 01 };
            return new Xint[] { 0, 0 };
        }
        int n = bL(D) / 6;
        Xint Mask = (Xint.One << n) - 1;
        Xint A = D & Mask; D >>= n;
        Xint B = D & Mask; D >>= n;
        Xint C = D & Mask; D >>= n;
        Xint[] R = CRold(D);                                  // R[0] = CubeRoot, R[1] = Remainder
        Xint R0R0 = SQ(R[0]);
        Xint[] Q = DQR(C + (R[1] << n), 3 * R0R0);            // Q[0] = Quotient, Q[1] = Remainder
        Xint Q0Q0 = SQ(Q[0]);
        R[1] = A + ((B + (Q[1] << n)) << n) - MTP(Q0Q0, Q[0] + ((3 * R[0]) << n));
        if (R[1] < 0)
        {
            R0R0 <<= n * 2;
            R0R0 += Q0Q0 + (MTP(Q[0], R[0]) << (n + 1));
            R[0] <<= n;
            R[0] += Q[0];
            R[1] += 1 + 3 * (R0R0 - R[0]);
            R[0] -= 1;
            while (R[1] < 0)
            {
                R0R0 -= 1 + 2 * R[0];
                R[1] += 1 + 3 * (R0R0 - R[0]);
                R[0] -= 1;
            }
        }
        else
        {
            R[0] <<= n;
            R[0] += Q[0];
        }
        return R;
    }
}

Power Addition Only

/*





Finding Powers using only additions?

    Multiplication is repeated addition: 
    3*2 = 3+3 or 2+2+2
    
    Powering is repeated multiplication: 
    3^3 = 3*3*3
        = 3*(3+3+3)
        = (3+3+3)+(3+3+3)+(3+3+3)
        
But here's something different. It came along while I wrote Cube Root code,
it went off topic, thus another topic, this one.

      x        x^4
      0            0       1     14     36   24   0
      1            1     15     50     60   24   0
      2          16     65   110     84   24   0
      3          81   175   194   108   24   0
      4        256   369   302   132   24   0

Above each black number is the sum of the numbers above and right above it,
because: "see Cube Root".
     _                     _       _               _ 
    |    0   1  14  36  24  |     |  A  B  C  D  E  |      F=A+B , G=B+C , ...
    |    1  15  50  60  24  |     |  F  G  H  I  .  |      
    |   16  65 110  84  24  |  =  |  J  K  L  .  .  |      
    |   81 175 194 108  24  |     |  M  N  .  .  .  |      
    |_ 256 369 302 132  24 _|     |_ P  .  .  .  . _|      
     
Each row depends on the row above it, depends on the first (magic?) row.
The other way round:     How to get the first row: " 0,1,14,36,24,0 ",
                       from the first colum (x^4): " 0,1,16,81,256  ".

    C=G-B G=J-F C=1J-2F+1A 
          B=F-A
          
    D=H-C H=K-G D=1K-2G+1B K=M-J D=M-J-2J+2F+F-A D=1M-3J+3F-1A 
          C=G-B            G=J-F   
                           B=F-A 
                           
    E=I-D I=L-H E=1L-2H+1C L=N-K E=N-K-2K+2G+G-B E=N-3K+3G-B N=P-M E=P-M-3M+3J+3J-3F-F+A E=1P-4M+6J-4F+1A 
          D=H-C            H=K-G                             K=M-J
                           C=G-B                             G=J-F
                                                             B=F-A
    A=      1.0      A=1.0^4=0                                                         A=           +2A+1.0
    B=    1F-1A      B=1.1^4-1.0^4                                                     B=         1F-2A+1A
    C=   1J-2F+1A    C=1.2^4-2.1^4+1.0^4             C=1J-2F+1A+1B-(1F-1A)             C=      1J-3F+2A+1B
    D=  1M-3J+3F-1A  D=1.3^4-3.2^4+3.1^4-1.0^4       D=1M-3J+3F-1A+1C-(1J-2F+1A)       D=   1M-4J+5F-2A+1C
    E=1P-4M+6J-4F+1A E=1.4^4-4.3^4+6.2^4-4.1^4+1.0^4 E=1P-4M+6J-4F+1A+1D-(1M-3J+3F-1A) E=1P-5M+9J-7F+2A+1D
     
    Pascal's Triangle

The first row is part of A131689 (~A019538): Triangle of numbers: T(n,k)=k!*Stirling2(n,k)
Next step: Use Pascal's triangle, use Binomial Coefficients ( C(.,.) ).

                                                 C(1,1).1^1 =  1   \
                                                                    \
                                                 C(1,1).1^2 =  1     \
                                    C(2,2).2^2 - C(2,1).1^2 =  2      \ 
                                                                       \
                                                 C(1,1).1^3 =  1        \
                                    C(2,2).2^3 - C(2,1).1^3 =  6         A019538
                       C(3,3).3^3 - C(3,2).2^3 + C(3,1).1^3 =  6        /
                                                                       /
                                                 C(1,1).1^4 =  1      /
                                    C(2,2).2^4 - C(2,1).1^4 = 14     /
                       C(3,3).3^4 - C(3,2).2^4 + C(3,1).1^4 = 36    /
          C(4,4).4^4 - C(4,3).3^4 + C(4,2).2^4 - C(4,1).1^4 = 24   /

*/
using System;
class Power_Addition_Only
{
    private static uint T(uint n, uint k)   // OEIS: A019538 =============================>>  //  1
    {                                                                                         //  1
        return F(k) * S2(n, k);                                                               //  2
    }                                                                                         //  1
    private static uint S2(uint n, uint k)  // Stirling numbers                               //  6
    {                                       //    2nd kind                                    //  6
        if (n == k) return 1;                                                                 //  1
        if (n == 0 || k == 0) return 0;                                                       //  14
        n--;                                                                                  //  36
        return k * S2(n, k--) + S2(n, k);                                                     //  24
    }                                                                                         //  1
    private static uint F(uint k)           // Factorial                                      //  30
    {                                                                                         //  150
        uint f = 1;                                                                           //  240
        while (k > 1) f *= k--;                                                               //  120
        return f;                                                                             //  1
    }                                                                                         //  62
                                                                                              //  540
    static void Main()                                                                        //  1560
    {                                                                                         //  1800
        Console.WriteLine();                                                                  //  720
        powers(4, 4);                                                                         //  1
        powers(3, 1);                                                                         //  126
        powers(21, 6);                                                                        //  1806
        powers(15, 7);                                                                        //  8400
        powers(9, 9);                                                                         //  16800
        Console.ReadLine();                                                                   //  15120
    }                                                                                         //  5040
    private static void powers(uint x, uint p) // up to x^p                                   //  1
    {                                                                                         //  254
        Console.WindowWidth = 105;                                                            //  5796
        Console.WindowHeight = Console.LargestWindowHeight - 10;                              //  40824
        Console.WriteLine("       x^" + p);                                                   //  126000
        Console.WriteLine();                                                                  //  191520
        int i = 0, j;                                                                         //  141120
        uint[] v = Init(p);                                                                   //  40320
        int vL = v.Length - 1;                                                                //  1
        for (j = 0; j <= vL; j++) Console.Write("{0,10}", v[j]);                              //  510
        Console.WriteLine();                                                                  //  18150
        for (; i < x; i++)                                                                    //  186480
        {                                                                                     //  834120
            for (j = 0; j < vL; ) v[j++] += v[j];                                             //  1905120
            for (j = 0; j <= vL; j++) Console.Write("{0,10}", v[j]);                          //  2328480
            Console.WriteLine();                                                              //  1451520
        }                                                                                     //  362880
        Console.WriteLine();                                                                  //  1
        Console.WriteLine();                                                                  //  1022
    }                                                                                         //  55980
    private static uint[] Init(uint i)      // get Initial value's                            //  818520
    {                                                                                         //  5103000
        uint[] v = new uint[i + 1];                                                           //  16435440
        v[0] = 0;                                                                             //  29635200
        for (uint j = 1; j <= i; j++)                                                         //  30240000
            v[j] = T(i, j);                                                                   //  16329600
        return v;                                                                             //  3628800
    }                                                                                         //  1      
}

Nth Root & Power

/*  Nth Root, Nth Power
  

     ____
 \n /                                 Input: integers x >= 0 , n > 0
  \/  x  = r                         Output: integer  r, such that r^n <= x < (r+1)^n
                               
 
Examples: 27^(1/3)=3, 26^(1/3)=2

           3th root 1234567890
                    1 234 567890
                    1 ... ......
           3th root 1=1
                   (1+1)*10=20
                    1 234 ......
                     (2*20 + 1234/20^2)/3
                     (2*20 + 12  /2 ^2)/3=14
                     (2*14 + 1234/14^2)/3=11
                     (2*11 + 1234/11^2)/3=10
                     (2*10 + 1234/10^2)/3=10
           3th root 1 234=10
                   (10+1)*100=1100
                    1 234 567890
                     (2*1100 + 1234567890/1100^2)/3
                     (2*1100 + 123456    /11  ^2)/3=1073
                     (2*1073 + 1234567890/1073^2)/3=1072
                     (2*1072 + 1234567890/1072^2)/3=1072
 
The remainder is one step more, one multiplication (and a subtraction),
in the last step 1072^2 is computed,
so the remainder is 1234567890-(1072^2)*1072
 
Finding roots involves finding powers.
The POW function is updated (the thirth power of an 1.000.000 bits number
is found ~20 times faster). It calls MTP_XP(Xint X, Xint P), which splits
a large P in pieces with X's bitlength, so fast multiplication can be used.
At the bottom a few versions of POW can be found, iterative and recursive, 
bitwise scanning the exponent from left to right or vice versa.
POWrRL, recursive, scans from R to L, fastest, POW is based on it.
POWiLR, iterative, marginally slower.
POWiRL, slow, why?
POWstd, standard System.Numerics.BigInteger.Pow method.
 
Abbreviations:
 NR      Nth Root
 PR      Partial Root, called by NR
 POW     Power
 MTP_XP  Multiply, see above, called by POW
 DIV     Division, returns quotient, calls DQR
 DQR     Division, returns quotient and remainder (Burnikel & Ziegler)
 bL      bitLength
 SR      Square Root & Remainder, SR(10)=3 R=1 (Zimmermann)
 SRO     Square Root Only, SRO(10)=3
 RND     Random number
 
 |---------------------------------------------------------------------------------------------------|
 |                                Times(in ms, Athlon X4, XP, 2GB)                                   |
 |-------------------------------------||------------------------------------------------------------|
 |                                     ||               RND of 1.000.000 bits                        |
 |------------|------------|-----------||----|------------|-----------|------------------------------|                           
 |   bits RND | POW(RND,3) | NR(RND,3) ||  n | POW(RND,n) | NR(RND,n) |  Using Square Root           |
 |------------|------------|-----------||----|------------|-----------|------------------------------|
 |     10.000 |       1    |      1    ||  2 |     160    |     840   |  187 SRO(RND)                |
 |     20.000 |       4    |      3    ||  3 |     550    |    1180   |                              |
 |     50.000 |       8    |     20    ||  4 |     618    |     777   |  261 SRO(SRO(RND))           |
 |    100.000 |      20    |     67    ||  5 |    1550    |     849   |                              |
 |    200.000 |      57    |    143    ||  6 |    1310    |    1130   |  602 NR(SRO(RND),3)          |
 |    500.000 |     213    |    414    ||  7 |    2440    |    1490   |                              |
 |  1.000.000 |     550    |   1200    ||  8 |    1790    |     908   |  291 SRO(SRO(SRO(RND)))      |
 |  2.000.000 |    1630    |   3120    ||  9 |    3300    |     868   |                              |
 |  5.000.000 |    6020    |   7860    || 10 |    3070    |    1050   |  537 NR(SRO(RND),5)          |
 | 10.000.000 |   15500    |  19700    || 11 |    4970    |    1090   |                              |
 |            |            |           || 12 |    3370    |    1360   |  410 NR(SRO(SRO(RND)),3)     |
 |            |            |           || 13 |    5700    |    1310   |                              |
 |            |            |           || 14 |    5150    |    1358   |  639 NR(SRO(RND),7)          |
 |            |            |           || 15 |    7870    |    1550   | 1010 NR(NR(RND,5),3)         |
 |            |            |           || 16 |    5020    |    1040   |  303 SRO(SRO(SRO(SRO(RND)))) |
 |            |            |           || 17 |   10400    |     750   |                              |
 |------------|------------|-----------||----|------------|-----------|------------------------------|
  
Reference: Modern Computer Arithmetic (Version 0.5.9 of 7 October 2010),
           Richard Brent and Paul Zimmermann, Cambridge University Press, 2010.
           (On-line book! (1.5.2 kth root))     
*/

using System;
using System.Diagnostics;
using System.Threading.Tasks;
using Xint = System.Numerics.BigInteger;
class Nth_Root_Power
{
    private static Xint NR(Xint X, int y)                 // Nth Root
    {
        int sr_X = (bL(X) - 1) / y * y - y;               // shift right
        Xint R0 = 1;
        int n = 1;
        while (sr_X >= 0)
        {
            R0++;
            R0 = PR(X >> sr_X, R0, n, y);
            n *= 2;
            sr_X -= n * y;
        }
        sr_X += n * y;
        if (sr_X == 0) return R0;
        R0++;
        R0 = PR(X, R0, sr_X / y, y);
        return R0;
    }
    private static Xint PR(Xint X, Xint R0, int n, int y) // Partial Root
    {
        Xint S, Z, T;
        S = R0 << n;
        Z = POW(R0, y - 1);
        T = (y - 1) * S + DIV((X >> (y - 1) * n), Z);
        R0 = T / y;
        do
        {
            S = R0;
            Z = POW(S, y - 1);
            T = (y - 1) * S + DIV(X, Z);
            R0 = T / y;
        }
        while (R0 < S);
        return S;
    }

    private static Xint POW(Xint X, int y)
    {
        if (y > 1) return ((y & 1) == 0) ? SQ(POW(X, y / 2)) : MTP_XP(X, SQ(POW(X, y / 2)));
        return y == 0 ? 1 : X;
    }
    private static Xint MTP_XP(Xint X, Xint P)
    {
        int n = bL(X);
        if (n <= 7500) return X * P;
        Xint Mask = (Xint.One << n) - 1;
        Xint T = MTP(X, P & Mask, n);
        P >>= n;
        int s = n;
        while (P > 0)
        {
            T += MTP(X, P & Mask, n) << s;
            P >>= n;
            s += n;
        }
        return T;
    }

    #region
    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        Xint X = RND(1000000);
        if (NR(X, 15) != NR(NR(X, 5), 3)) Console.WriteLine("WRONG");
        sw.Restart();
        NR(NR(X, 5), 3);
        sw.Stop();
        Console.WriteLine(sw.ElapsedMilliseconds);
        sw.Restart();
        NR(X, 15);
        sw.Stop();
        Console.WriteLine(sw.ElapsedMilliseconds);
        X = RND(500000);
        X = SQ(X); //X++; //X--;
        if (NR(X, 16) != SRO(SRO(SRO(SRO(X))))) Console.WriteLine("WRONG");
        int i, j = 1;
        sw.Restart();
        for (i = 0; i < j; i++)
            SRO(SRO(SRO(SRO(X))));
        sw.Stop();
        Console.WriteLine(sw.ElapsedMilliseconds);
        sw.Restart();
        for (i = 0; i < j; i++)
            NR(X, 16);
        sw.Stop();
        Console.WriteLine(sw.ElapsedMilliseconds);
        Console.ReadLine();
    }

    private static Xint[] SR(Xint A)
    {
        return SR(A, bL(A));
    }
    private static Xint[] SR(Xint A, int n)
    {
        if (n < 53) return SR52(A);
        int m = n >> 2;
        Xint Mask = (Xint.One << m) - 1;
        Xint A0 = A & Mask; A >>= m;
        Xint A1 = A & Mask; A >>= m;
        Xint[] R = SR(A, n - (m << 1));
        Xint[] D = DQR(R[1] << m | A1, R[0] << 1);
        R[0] = (R[0] << m) + D[0];
        R[1] = (D[1] << m | A0) - SQ(D[0], m);
        if (R[1] < 0)
        {
            R[0] -= 1;
            R[1] += (R[0] << 1) | 1;
        }
        return R;
    }
    private static Xint[] SR52(Xint A)
    {
        double a = (double)A;
        long q = (long)Math.Sqrt(a);
        long r = (long)(a) - q * q;
        Xint[] QR = { q, r };
        return QR;
    }

    private static Xint SRO(Xint A)
    {
        return SRO(A, bL(A));
    }
    private static Xint SRO(Xint A, int n)
    {
        if (n < 53) return (int)Math.Sqrt((double)A);
        Xint[] R = SROr(A, n, 1);
        return R[0];
    }
    private static Xint[] SROr(Xint A, int n, int rc) // rc=recursion counter
    {
        if (n < 53) return SR52(A);
        int m = n >> 2;
        Xint Mask = (Xint.One << m) - 1;
        Xint A0 = A & Mask; A >>= m;
        Xint A1 = A & Mask; A >>= m;
        Xint[] R = SROr(A, n - (m << 1), rc + 1);
        Xint[] D = DQR((R[1] << m) | A1, R[0] << 1);
        R[0] = (R[0] << m) + D[0];
        rc--;
        if (rc != 0)
        {
            R[1] = (D[1] << m | A0) - SQ(D[0], m);
            if (R[1] < 0)
            {
                R[0] -= 1;
                R[1] += (R[0] << 1) | 1;
            }
            return R;
        }
        n = (bL(D[0]) << 1) - bL(D[1] << m | A0);
        if (n < 0) return R;
        if (n > 1)
        {
            R[0] -= 1;
            return R;
        }
        int shift = (bL(D[0]) - 31) << 1;
        long d0 = (int)(D[0] >> (shift >> 1));
        long d = (long)((D[1] >> (shift - m)) | (A0 >> shift)) - d0 * d0;
        if (d < 0)
        {
            R[0] -= 1;
            return R;
        }
        if (d > d0 << 1) return R;
        R[0] -= (1 - (((D[1] << m) | A0) - SQ(D[0], m)).Sign) >> 1;
        return R;
    }

    private static int bL(Xint U)
    {
        byte[] bytes = (U.Sign * U).ToByteArray();
        int i = bytes.Length - 1;
        return i * 8 | bitLengthMostSignificantByte(bytes[i]);
    }
    private static int bitLengthMostSignificantByte(byte b)
    {
        return b < 08 ? b < 02 ? b < 01 ? 0 : 1 :
                                 b < 04 ? 2 : 3 :
                        b < 32 ? b < 16 ? 4 : 5 :
                                 b < 64 ? 6 : 7;
    }

    private static int fL2(int i)
    {
        return
        i < 1 << 15 ? i < 1 << 07 ? i < 1 << 03 ? i < 1 << 01 ? i < 1 << 00 ? -1 : 00 :
                                                                i < 1 << 02 ? 01 : 02 :
                                                  i < 1 << 05 ? i < 1 << 04 ? 03 : 04 :
                                                                i < 1 << 06 ? 05 : 06 :
                                    i < 1 << 11 ? i < 1 << 09 ? i < 1 << 08 ? 07 : 08 :
                                                                i < 1 << 10 ? 09 : 10 :
                                                  i < 1 << 13 ? i < 1 << 12 ? 11 : 12 :
                                                                i < 1 << 14 ? 13 : 14 :
                      i < 1 << 23 ? i < 1 << 19 ? i < 1 << 17 ? i < 1 << 16 ? 15 : 16 :
                                                                i < 1 << 18 ? 17 : 18 :
                                                  i < 1 << 21 ? i < 1 << 20 ? 19 : 20 :
                                                                i < 1 << 22 ? 21 : 22 :
                                    i < 1 << 27 ? i < 1 << 25 ? i < 1 << 24 ? 23 : 24 :
                                                                i < 1 << 26 ? 25 : 26 :
                                                  i < 1 << 29 ? i < 1 << 28 ? 27 : 28 :
                                                                i < 1 << 30 ? 29 : 30;
    }

    private static int seed;
    private static Xint RND(int n)
    {
        if (n < 2) return n;
        if (seed == int.MaxValue) seed = 0; else seed++;
        Random rand = new Random(seed);
        byte[] bytes = new byte[(n + 15) >> 3];
        rand.NextBytes(bytes);
        int i = bytes.Length - 1;
        bytes[i] = 0;
        n = (i << 3) - n;
        i--;
        bytes[i] >>= n;
        bytes[i] |= (byte)(128 >> n);
        return new Xint(bytes);
    }

    private static Xint SQ(Xint U)
    {
        return SQ(U, U.Sign * U.ToByteArray().Length << 3);
    }
    private static Xint SQ(Xint U, int n)
    {
        if (n <= 700) return U * U;
        if (n <= 3000) return Xint.Pow(U, 2);
        if (n <= 6000) return SQ2(U, n);
        if (n <= 10000) return SQ3(U, n);
        if (n <= 40000) return SQ4(U, n);
        return SQ2P(U, n);
    }
    private static Xint SQr(Xint U, int n)
    {
        if (n <= 3000) return Xint.Pow(U, 2);
        if (n <= 6000) return SQ2(U, n);
        if (n <= 10000) return SQ3(U, n);
        return SQ4(U, n);
    }
    private static Xint SQ2(Xint U1, int n)
    {
        n >>= 1;
        Xint U0 = U1 & ((Xint.One << n) - 1); U1 >>= n;
        Xint P0 = SQr(U0, n);
        Xint P2 = SQr(U1, n);
        return ((P2 << n) + (SQr(U0 + U1, n) - (P0 + P2)) << n) + P0;
    }
    private static Xint SQ3(Xint U2, int n)
    {
        n = (int)((long)(n) * 0x55555556 >> 32);
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U2 & Mask; U2 >>= n;
        Xint U1 = U2 & Mask; U2 >>= n;
        Xint W0 = SQr(U0, n);
        Xint W4 = SQr(U2, n);
        Xint P3 = SQr((((U2 << 1) + U1) << 1) + U0, n);
        U2 += U0;
        Xint P2 = SQr(U2 - U1, n);
        Xint P1 = SQr(U2 + U1, n);
        Xint W2 = (P1 + P2 >> 1) - (W0 + W4);
        Xint W3 = W0 - P1;
        W3 = ((W3 + P3 - P2 >> 1) + W3) / 3 - (W4 << 1);
        Xint W1 = P1 - (W4 + W3 + W2 + W0);
        return ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint SQ4(Xint U3, int n)
    {
        n >>= 2;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U3 & Mask; U3 >>= n;
        Xint U1 = U3 & Mask; U3 >>= n;
        Xint U2 = U3 & Mask; U3 >>= n;
        Xint W0 = SQr(U0, n);                                   //  0
        U0 += U2;
        U1 += U3;
        Xint P1 = SQr(U0 + U1, n);                              //  1
        Xint P2 = SQr(U0 - U1, n);                              // -1
        U0 += 3 * U2;
        U1 += 3 * U3;
        Xint P3 = SQr(U0 + (U1 << 1), n);                       //  2
        Xint P4 = SQr(U0 - (U1 << 1), n);                       // -2
        Xint P5 = SQr(U0 + 12 * U2 + ((U1 + 12 * U3) << 2), n); //  4
        Xint W6 = SQr(U3, n);                                   //  inf
        Xint W1 = P1 + P2;
        Xint W4 = (((((P3 + P4) >> 1) - (W1 << 1)) / 3 + W0) >> 2) - 5 * W6;
        Xint W2 = (W1 >> 1) - (W6 + W4 + W0);
        P1 = P1 - P2;
        P4 = P4 - P3;
        Xint W5 = ((P1 >> 1) + (5 * P4 + P5 - W0 >> 4) - ((((W6 << 4) + W4) << 4) + W2)) / 45;
        W1 = ((P4 >> 2) + (P1 << 1)) / 3 + (W5 << 2);
        Xint W3 = (P1 >> 1) - (W1 + W5);
        return ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint SQ2P(Xint A, int n)
    {
        n >>= 1;
        Xint[] U = new Xint[3];
        U[0] = A & ((Xint.One << n) - 1); A >>= n; U[2] = A; U[1] = U[0] + A;
        Xint[] P = new Xint[3];
        Parallel.For(0, 3, (int i) => P[i] = SQr(U[i], n));
        return ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0];
    }

    private static Xint MTP(Xint U, Xint V)
    {
        return MTP(U, V, Xint.Max(U.Sign * U, V.Sign * V).ToByteArray().Length << 3);
    }
    private static Xint MTP(Xint U, Xint V, int n)
    {
        if (n <= 3000) return U * V;
        if (n <= 6000) return TC2(U, V, n);
        if (n <= 10000) return TC3(U, V, n);
        if (n <= 40000) return TC4(U, V, n);
        return TC2P(U, V, n);
    }
    private static Xint MTPr(Xint U, Xint V, int n)
    {
        if (n <= 3000) return U * V;
        if (n <= 6000) return TC2(U, V, n);
        if (n <= 10000) return TC3(U, V, n);
        return TC4(U, V, n);
    }
    private static Xint TC2(Xint U1, Xint V1, int n)
    {
        n >>= 1;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U1 & Mask; U1 >>= n;
        Xint V0 = V1 & Mask; V1 >>= n;
        Xint P0 = MTPr(U0, V0, n);
        Xint P2 = MTPr(U1, V1, n);
        return ((P2 << n) + (MTPr(U0 + U1, V0 + V1, n) - (P0 + P2)) << n) + P0;
    }
    private static Xint TC3(Xint U2, Xint V2, int n)
    {
        n = (int)((long)(n) * 0x55555556 >> 32); // n /= 3;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U2 & Mask; U2 >>= n;
        Xint U1 = U2 & Mask; U2 >>= n;
        Xint V0 = V2 & Mask; V2 >>= n;
        Xint V1 = V2 & Mask; V2 >>= n;
        Xint W0 = MTPr(U0, V0, n);
        Xint W4 = MTPr(U2, V2, n);
        Xint P3 = MTPr((((U2 << 1) + U1) << 1) + U0, (((V2 << 1) + V1 << 1)) + V0, n);
        U2 += U0;
        V2 += V0;
        Xint P2 = MTPr(U2 - U1, V2 - V1, n);
        Xint P1 = MTPr(U2 + U1, V2 + V1, n);
        Xint W2 = (P1 + P2 >> 1) - (W0 + W4);
        Xint W3 = W0 - P1;
        W3 = ((W3 + P3 - P2 >> 1) + W3) / 3 - (W4 << 1);
        Xint W1 = P1 - (W4 + W3 + W2 + W0);
        return ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint TC4(Xint U3, Xint V3, int n)
    {
        n >>= 2;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U3 & Mask; U3 >>= n;
        Xint U1 = U3 & Mask; U3 >>= n;
        Xint U2 = U3 & Mask; U3 >>= n;
        Xint V0 = V3 & Mask; V3 >>= n;
        Xint V1 = V3 & Mask; V3 >>= n;
        Xint V2 = V3 & Mask; V3 >>= n;

        Xint W0 = MTPr(U0, V0, n);                               //  0
        U0 += U2; U1 += U3;
        V0 += V2; V1 += V3;
        Xint P1 = MTPr(U0 + U1, V0 + V1, n);                     //  1
        Xint P2 = MTPr(U0 - U1, V0 - V1, n);                     // -1
        U0 += 3 * U2; U1 += 3 * U3;
        V0 += 3 * V2; V1 += 3 * V3;
        Xint P3 = MTPr(U0 + (U1 << 1), V0 + (V1 << 1), n);       //  2
        Xint P4 = MTPr(U0 - (U1 << 1), V0 - (V1 << 1), n);       // -2
        Xint P5 = MTPr(U0 + 12 * U2 + ((U1 + 12 * U3) << 2),
                       V0 + 12 * V2 + ((V1 + 12 * V3) << 2), n); //  4
        Xint W6 = MTPr(U3, V3, n);                               //  inf

        Xint W1 = P1 + P2;
        Xint W4 = (((((P3 + P4) >> 1) - (W1 << 1)) / 3 + W0) >> 2) - 5 * W6;
        Xint W2 = (W1 >> 1) - (W6 + W4 + W0);
        P1 = P1 - P2;
        P4 = P4 - P3;
        Xint W5 = ((P1 >> 1) + (5 * P4 + P5 - W0 >> 4) - ((((W6 << 4) + W4) << 4) + W2)) / 45;
        W1 = ((P4 >> 2) + (P1 << 1)) / 3 + (W5 << 2);
        Xint W3 = (P1 >> 1) - (W1 + W5);
        return ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint TC2P(Xint A, Xint B, int n)
    {
        n >>= 1;
        Xint Mask = (Xint.One << n) - 1;
        Xint[] U = new Xint[3];
        U[0] = A & Mask; A >>= n; U[2] = A; U[1] = U[0] + A;
        Xint[] V = new Xint[3];
        V[0] = B & Mask; B >>= n; V[2] = B; V[1] = V[0] + B;
        Xint[] P = new Xint[3];
        Parallel.For(0, 3, (int i) => P[i] = MTPr(U[i], V[i], n));
        return ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0];
    }

    private static Xint DIV(Xint A, Xint B)
    {
        Xint[] QR = DQR(A, B);
        return QR[0];
    }
    private static Xint[] DQR(Xint A, Xint B)
    {
        int n = bL(B);
        int m = bL(A) - n;
        if (m <= 6000) return SmallDivRem(A, B);
        int signA = A.Sign; A *= signA;
        int signB = B.Sign; B *= signB;
        Xint[] QR = new Xint[2];
        if (m <= n) QR = D21(A, B, n);
        else
        {
            Xint Mask = (Xint.One << n) - 1;
            m /= n;
            Xint[] U = new Xint[m];
            int i = 0;
            for (; i < m; i++)
            {
                U[i] = A & Mask;
                A >>= n;
            }
            QR = D21(A, B, n);
            A = QR[0];
            for (i--; i >= 0; i--)
            {
                QR = D21(QR[1] << n | U[i], B, n);
                A = A << n | QR[0];
            }
            QR[0] = A;
        }
        QR[0] *= signA * signB;
        QR[1] *= signA;
        return QR;
    }
    private static Xint[] SmallDivRem(Xint A, Xint B)
    {
        Xint[] QR = new Xint[2];
        QR[0] = Xint.DivRem(A, B, out QR[1]);
        return QR;
    }
    private static Xint[] D21(Xint A, Xint B, int n)
    {
        if (n <= 6000) return SmallDivRem(A, B);
        int s = n & 1;
        A <<= s;
        B <<= s;
        n++;
        n >>= 1;
        Xint Mask = (Xint.One << n) - 1;
        Xint B1 = B >> n;
        Xint B2 = B & Mask;
        Xint[] QR1 = D32(A >> (n << 1), A >> n & Mask, B, B1, B2, n);
        Xint[] QR2 = D32(QR1[1], A & Mask, B, B1, B2, n);
        QR2[0] |= QR1[0] << n;
        QR2[1] >>= s;
        return QR2;
    }
    private static Xint[] D32(Xint A12, Xint A3, Xint B, Xint B1, Xint B2, int n)
    {
        Xint[] QR = new Xint[2];
        if ((A12 >> n) != B1) QR = D21(A12, B1, n);
        else
        {
            QR[0] = (Xint.One << n) - 1;
            QR[1] = A12 + B1 - (B1 << n);
        }
        QR[1] = (QR[1] << n | A3) - MTP(QR[0], B2, n);
        while (QR[1] < 0)
        {
            QR[0] -= 1;
            QR[1] += B;
        }
        return QR;
    }

    private static Xint POWstd(Xint X, int y) // standard
    {
        return Xint.Pow(X, y);
    }
    private static Xint POWrRL(Xint X, int y) // recursive Right to Left 
    {
        if (y > 1) return ((y & 1) == 0) ? SQ(POWrRL(X, y / 2)) : SQ(POWrRL(X, y / 2)) * X;
        return y == 0 ? 1 : X;
    }
    private static Xint POWiRL(Xint X, int y) // iterative Right to Left
    {
        Xint P = 1;
    L0: if ((y & 1) == 1) P *= X;
        y /= 2;
        if (y == 0) return P;
        X = SQ(X);
        goto L0;
    }
    private static Xint POWiLR(Xint X, int y) // iterative Left to Right
    {
        switch (y)
        {
            case 0: return 1;
            case 1: return X;
            case 2: return SQ(X);
            default:
                Xint P = SQ(X);
                for (int i = fL2(y) - 1; i > 0; i--)
                {
                    if (((y >> i) & 1) != 0) P *= X;
                    P = SQ(P);
                }
                return (y & 1) == 0 ? P : P * X;
        }
    }
    #endregion
}

Pascal's triangle

// Horizontal row in Pascal's triangle

using System;
class Pascal_Triangle_Horizontal_Row
{
    private static int[] fr(int n)            // full row
    {                                         // 0=>{1} 1=>{1,1} 2=>{1,2,1} 3=>{1,3,3,1} ...
        int i = 0, j = n / 2, k = 1;
        int[] c = new int[n + 1];
        c[0] = 1;
        for (; i < j; i++, k++)
            c[k] = c[i] * (n - i) / k;
        for (i = 0; n > i; i++, n--)
            c[n] = c[i];
        return c;
    }

    private static int[] pfr(int[] d)         // previous full row
    {                                         // ...=>{1,3,3,1}=>{1,2,1}=>{1,1}=>{1}=>{}=>{} 
        int n = d.Length - 2;
        if (n < 0)
            return new int[] { };
        int i = 1, m = n / 2;
        int[] c = new int[n + 1];
        c[0] = 1;
        for (; i <= m; i++)
            c[i] = d[i] - c[i - 1];
        for (i = 0; n > i; i++, n--)
            c[n] = c[i];
        return c;
    }

    private static int[] nfr(int[] b)         // next full row
    {                                         // {}=>{1}=>{1,1}=>{1,2,1}=>{1,3,3,1}=>...
        int i = 1, n = b.Length, m = n / 2;
        int[] c = new int[n + 1];
        c[0] = 1;
        for (; i <= m; i++)
            c[i] = b[i - 1] + b[i];
        for (i = 0; n > i; i++, n--)
            c[n] = c[i];
        return c;
    }

    private static int[] hr(int n)            // half row
    {                                         // 0=>{1} 1=>{1} 2=>{1,2} 3=>{1,3} 4=>{1,4,6} ...
        int i = 0, j = n / 2, k = 1;
        int[] c = new int[j + 1];
        c[0] = 1;
        for (; i < j; i++, k++, n--)
            c[k] = c[i] * n / k;
        return c;
    }

    private static int[] phr(int[] d)         // previous half row
    {                                         // ...=>{1,4,6}=>{1,3}=>{1,2}=>{1}=>{}=>{}
        if (d.Length < 2)
            return new int[] { };
        int i = 1, n = (d[1] + 1) / 2;
        int[] c = new int[n];
        c[0] = 1;
        for (; i < n; i++)
            c[i] = d[i] - c[i - 1];
        return c;
    }

    private static int[] nhr(int[] b)         // next half row
    {                                         // {}=>{1}=>{1,2}=>{1,3}=>{1,4,6}=>...
        switch (b.Length)
        {
            case 0: return new int[] { 1 };
            case 1: return new int[] { 1, 2 };
            default:
                int i = 1, b1 = b[1], n = (b1 + 3) / 2;
                b1 &= 1;
                int[] c = new int[n--];
                c[0] = 1;
                for (; i < n; i++)
                    c[i] = b[i] + b[i - 1];
                c[i] = (1 - b1) * b[i - 1] + (1 + b1) * b[i - b1];
                return c;
        }
    }

    static void Main()
    {
        int[] c = new int[] { };              // 1
        for (int i = 0; i < 9; i++)           // 1 1
        {                                     // 1 2 1
            c = nfr(c);                       // 1 3 3 1
            foreach (int j in c)              // 1 4 6 4 1
                Console.Write(j + " ");       // 1 5 10 10 5 1
            Console.WriteLine();              // 1 6 15 20 15 6 1      
        }                                     // 1 7 21 35 35 21 7 1   
        Console.ReadLine();                   // 1 8 28 56 70 56 28 8 1
    }
}

A218076 three more with C on 64 bits Linux

/*
The next terms are:
F(177416979) with 54 consecutive zeros,
  198423144       56
  275354607       58
Next one will be larger than 3*10^8, presumable 18*10^8.
 
Computed with a new silent cheap PC.


 

 


















CPU: AMD Phenom II X4 965, 3.4 GHz (fastest in GMP benchmarks)
RAM: 8 GB
 OS: Debian 6.0.6 (Squeeze) Kernel Linux 2.6.32-5-amd64
GMP: 5.1.0
GCC: 4.4.5
IDE: MonoDevelop 2.4

Results: Fibonacci(10^9) in 12 seconds, factorial(32*10^6) 30 seconds,
both using 1 core, what will happen with a parallel version of Karatsuba? 
The C code below is about two times faster than the C++ code before,
probably because a 64 bits system is used instead of a 32 bits one,
most time is spent adding big numbers, with 64 bits it should be 
about two times faster, shouldn't it?

Is it possible to approximate the chance which Fibonacci number, F(n),
has a record number of consecutive zeros, z?
For larger values of n it seems to be.
Suppose we write the numbers as a long binary string.
F1=1, F2=1, F3=10, F4=11, F5=101, F6=1000, ...
The long binary string ends with: ...1000 101 11 10 1 1

How many bits are there from F1 up to Fn?
                             F1       F1?  0.7 * 1
                             F1       F2?  0.7 * (1+2)
                             F1       F3?  0.7 * (1+2+3)
                             F1       F4?  0.7 * (1+2+3+4)
                             F1       Fn?  0.7 * n * (n+1) / 2

How many coin flips on average does it take to get z consecutive heads?

Answer: 2^(z+1) - 2

So:     2^(z+1) - 2 = 0.7 * n * (n+1) / 2

For larger values of z (and n):

 n ~ 2.40 * 2^(z/2)

 z ~ 2.89 * ln(n) - 2.53

          n  cz(F(n))  2.89*ln(n)-2.53

          3        1    1 
          6        3    3
         12        4    5
         19        5    6
         38        6    8
         42        7    8
         68       12   10
        243       13   13
        384       14   15
        515       16   16
        740       19   17
       1709       24   19
       5151       27   22
      11049       28   24
      45641       30   28
      94729       31   31
     185610       34   33
     644593       35   36
     726681       39   36
    2296396       40   40
    3098358       42   41
    6178778       43   43
   15743325       45   45 
   22436908       48   46
   80141430       49   50
   84300971       51   50
  127495932       53   51
  177416979       54   52
  198423144       56   53
  275354607       58   54
 
Conclusion: 
A218076 seems to have properties of a random sequence.
Of course it isn't a random sequence, the terms are the: 
Lowest positive indices of Fibonacci numbers whose binary
expansions have record numbers of consecutive zeros.
And...it's getting harder and harder to find new ones.
The next term, with 59 consecutive zeros, will be 
~6.6 times larger than the last term found (50% chance).
*/

#include <stdio.h>
#include <gmp.h>
#include <time.h>

int lz(unsigned char b)
{
    return b < 1 << 4 ? b < 1 << 2 ? b < 1 << 1 ? 7 : 6 :
                                     b < 1 << 3 ? 5 : 4 :
                        b < 1 << 6 ? b < 1 << 5 ? 3 : 2 :
                                     b < 1 << 7 ? 1 : 0;
} 

int tz(unsigned char b)
{
    return 7 - lz((unsigned char)(b & -(char)b));
} 

time_t rawtime;
struct tm * timeinfo;
char buffer[80];
void printTime()
{
    time( &rawtime);
    timeinfo = localtime(&rawtime);
    strftime (buffer, 80, "%c", timeinfo);
    puts(buffer);    
}

void main ()
{
    int i, i0, i1, j, k, z, zmax, size; 
    printf(" low: "); scanf("%d", &i0);
    printf("high: "); scanf("%d", &i1);
    printf("zmax: "); scanf("%d", &zmax);
    k = (zmax - 6) / 8;
    mpz_t F0, F1;
    mpz_init(F0); mpz_init(F1);
    printTime();    
    printf(" compute F(%d)\n", i1);
    mpz_fib_ui(F0, i1);
    mpz_set(F1, F0); 
    printTime();    
    printf(" compute F's(%d & %d)\n", i0 - 2 , i0 - 1);
    mpz_fib2_ui(F1, F0, i0 - 1); 
    printTime();
    printf("enter a number: "); scanf("%d", &i);
    printTime();
    printf("   S T A R T \n");
    unsigned char *p , *q;     
    for(i = i0; i <= i1; i++)
    {
        if(i % 100000 == 0) 
        {
            printTime();
            printf("F(%d)\n", i);
        }
        mpz_add(F0,F0,F1);
        p = (unsigned char *)F0[0]._mp_d;
        j = 1;
        if (*p == 0)
        {
            for (z = 8, p++, j++; *p == 0; p++, j++, z += 8);
            if (zmax < z + 7)
            {
                z += tz(*p);
                if (zmax < z)
                {
                    zmax = z;
                    k = (z - 6) / 8;
                    printTime();
                    printf("  LOWZS: %d F(%d) \n", z, i);                   
                }
            }
        }
        size = mpz_sizeinbase(F0, 256);
        for (j += k, p += k; j < size; j += k, p += k)
        {
            if (*p == 0)
            {
                for (z = 8, q = p - 1; *q == 0; q--, z += 8);
                for (p++, j++; *p == 0; p++, j++, z += 8);
                if (zmax < z + 14)
                {
                    z += lz(*q);
                    if (zmax < z + 7)
                    {
                        z += tz(*p);
                        if (zmax < z)
                        {
                            zmax = z;
                            k = (z - 6) / 8;
                            printTime();
                            printf("  ZEROS: %d F(%d)  ZEROS \n", z, i);
                        }
                    }
                 }
            }
        }
        mpz_swap(F0, F1);     
    } 
    printTime();
    printf("   R E A D Y\n");
    printf("enter a number: "); scanf("%d", &i);
    printTime();
    printf(" compute F(%d)\n", i1);
    mpz_fib_ui(F0, i1); 
    printTime();
    if (mpz_cmp(F0, F1) != 0) 
    {
        printf("   E R R O R\n");
    }
    else
    {
        printf("  ZMAX: %d \n", zmax);
    }
    mpz_clear(F0); mpz_clear(F1);
}