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);
}

A218076 extended with C++

/*
The next term is: F(127495932) with 53 consecutive zeros.
I used Visual C++ 2010 Express, GMP, and MPIR,
see: GMP on windows
The C++ code below is about five times faster than the C# code before.
And... the next term will be larger than 150 * 10^6. 
*/

#include "stdafx.h"
#include "C:\projects\mpir-2.6.0\lib\Win32\Release\gmpxx.h"
#include "ctime"
#include "iostream"
using namespace std;

int lz(unsigned __int8 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 __int8 b)
{
    return 7 - lz((unsigned __int8)(b & -(__int8)b));
}

void main()
{
    int i, i0, i1, j, k, z, zmax, size;
    cout << " low: "; cin >> i0;
    cout << "high: "; cin >> i1;
    cout << "zmax: "; cin >> zmax;
    k = (zmax - 6) / 8;
    mpz_t F0, F1;
    mpz_inits(F0, F1, NULL);
    time_t now = time(0); cout << ctime(&now);
    cout << "compute F(" << dec << i1 << ")" << endl ;
    mpz_fib_ui(F0, i1); // mpz_out_str(stdout,16,F0)
    now = time(0); cout << ctime(&now);
    mpz_set(F1, F0);
    cout << "compute F's(" << dec << i0 - 2 << " & " << dec << i0 - 1 << ")" << endl;
    mpz_fib2_ui(F1, F0, i0 - 1);
    now = time(0); cout << ctime(&now) << endl;
    cout << "enter a number: ";
    cin >> i;
    now = time(0); cout << ctime(&now);
    cout << "S T A R T" << endl << endl;
    unsigned __int8 *p , *q;
    for(i = i0; i <= i1; i++)
    {
        if(i % 100000 == 0) 
        {
            now = time(0); cout << ctime(&now);
            cout << "F(" << dec << i << ")" << endl << endl;
        }
        mpz_add(F0,F0,F1);
        p = (unsigned __int8 *)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;
                    now = time(0); cout << "  " << ctime(&now);
                    cout << "  LOWZS: " << dec << z <<"  F(" << dec << i << ")" << endl << endl;
                }
            }
        }
        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;
                            now = time(0); cout << "  "  << ctime(&now);
                            cout << "  ZEROS: " << dec << z << "  F(" << dec << i << ") :ZEROS" << endl << endl;
                        }
                    }
                 }
            }
        }
        mpz_swap(F0,F1);
    }
    now = time(0); cout << ctime(&now);
    cout << "F(" << dec << i1 << ")" << endl;
    cout << "enter a number: ";
    cin >> i;
    now = time(0); cout << ctime(&now);
    cout << "compute F(" << dec << i1 << ")" << endl ;
    mpz_fib_ui(F0, i1);
    now = time(0); cout << ctime(&now) << endl;
    if (mpz_cmp(F0, F1) != 0) cout << endl << "E R R O R ! ! ! !" << endl;
    mpz_clears(F0, F1, NULL);
    cout << "zmax: " << dec << zmax << endl;
    cout << endl << "R E A D Y";
    cin >> i;
}

A218076 extended

// The next two terms are:
// F(80141430) 49 consecutive zeros, and
// F(84300971) 51 consecutive zeros
// I used Emil Stefanov's wrapper for GMP 
// (The GNU Multiple Precision Arithmetic Library).
// www.emilstefanov.net/Projects/GnuMpDotNet/
// Copied "libgmp-3.dll" to "C:\WINDOWS\System32\"
// This general compilation computes F(10^9) in 42 seconds.
// Two copies of the program below, for example one
// scanning the range 80000000-82000000, the other 
// scanning the range 82000000-84000000, take ~32 hours.
// And... the next term will be larger than 95999999.

using System;
using System.Threading.Tasks;
using Emil.GMP;  // Project>>Add Reference>>Emil.GMP.dll
class F_cz_gmp
{
    static void Main()
    {
        int i = 80000000;                          // 82000000
        int zmax = 48;
        int[] z = { 0, 0 };
        int k = (zmax - 6) / 8;
        Console.WriteLine("first i: " + i);
        BigInt[] F = new BigInt[2];
        F[1] = BigInt.Fibonacci(i - 1, out F[0]);
        for (; i < 82000000; i += 2)               // 84000000
        {
            if (i % 500000 == 0) Console.WriteLine(i + " " + DateTime.Now);
            F[0] += F[1];
            F[1] += F[0];
            Parallel.For(0, 2, (int j) => z[j] = zeros(F[j], zmax, k));
            GC.Collect();
            if (z[0] > zmax)
            {
                zmax = z[0];
                k = (zmax - 6) / 8;
                Console.WriteLine(i + " " + zmax);
            }
            if (z[1] > zmax)
            {
                zmax = z[1];
                k = (zmax - 6) / 8;
                Console.WriteLine((i + 1) + " " + zmax);
            }
        }
        Console.WriteLine("last i: " + (i - 1) + " " + DateTime.Now);
        Console.ReadLine();
    }
    private static int zeros(BigInt F, int zmax, int k)
    {
        byte[] b = F.ToByteArray(); // F.ToUintArray ~4 times faster
        int bL = b.Length;
        int j = 0;
        int z;
        if (b[0] == 0)
        {
            z = 8;
            while (b[++j] == 0) z += 8; // add zerobits
            if (zmax < z + 7)
            {
                z += tz(b[j]);
                if (zmax < z)
                {
                    zmax = z;
                    k = (z - 6) / 8;
                }
            }
        }
        j += k;
        while (j < bL)
        {
            if (b[j] == 0)
            {
                z = 8;
                int j0 = j;
                while (b[--j0] == 0) z += 8; // add zerobits
                while (b[++j] == 0) z += 8;
                if (zmax < z + 14)
                {
                    z += lz(b[j0]);
                    if (zmax < z + 7)
                    {
                        z += tz(b[j]);
                        if (zmax < z)
                        {
                            zmax = z;
                            k = (z - 6) / 8;
                        }
                    }
                }
            }
            j += k;
        }
        return zmax;
    }
    private static int lz(byte b) // leading zeros, b > 0
    {
        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;
    }
    private static int tz(byte b) // trailing zeros
    {
        return 7 - lz((byte)(b & -(sbyte)b));
    }
}

Consecutive zeros in Fibonacci numbers

// When I looked at the Fibonacci Sequence Binary Plot
// http://mathworld.wolfram.com/FibonacciNumber.html
// and the statement: F(2^n+2^(n+1)) ends in n+2 zeros.
// I wondered: Are there n+2 consecutive zeros for smaller indexes?
// Or: Indexes of Fibonacci numbers whose binary expansions have 
//     record numbers of consecutive zeros.
// It seems there are.
// Example: The first four records occur at 3, 6, 12, and 19:
// F(3)=10 2 (one zero)
// F(6)=1000 2 (three zeros)
// F(12)=10010000 2 (four zeros)
// F(19)=1000001010101 2 (five zeros)
//
//      Index  Zeros 
//          3    1  
//          6    3  
//         12    4  
//         19    5
//         38    6
//         42    7
//         68   12
//        243   13
//        384   14  
//        515   16
//        740   19
//       1709   24
//       5151   27
//      11049   28
//      45641   30
//      94729   31
//     185610   34
//     644593   35
//     726681   39
//    2296396   40
//    3098358   42
//    6178778   43
//   15743325   45
//   22436908   48
//
// In the program below, 
// the first loop "for (; zmax < 15; i++)" 
// searches for consecutive zerobits, bit by bit.
// As soon as 15 zeros are found, 
// the second loop "for (; i < int.MaxValue; i++)"
// searches for consecutive zerobytes, byte by byte,
// adding leading/trailing zerobits of the bytes before/after those zerobytes.
// It's much faster than searching bit by bit.

using System;
using Xint = System.Numerics.BigInteger; // Project>>Add Reference>>System.Numerics
class fibonacciRepeatingZeros
{
    static void Main()
    {
        Xint F0 = 1, F1 = 0, F2 = 0;
        int i = 1, z, zi = 0, zmax = 0; // z=zerobits,zi=maxzerobitsF(i)
        for (; zmax < 15; i++)
        {
            F2 = F1 + F0; F0 = F1; F1 = F2;
            z = 0;
            for (; !F2.IsOne; F2 /= 2)
            {
                if (F2.IsEven)
                {
                    if (zi < ++z) zi = z;
                }
                else z = 0;
            }
            if (zmax < zi)
            {
                zmax = zi;
                Console.WriteLine(i + " " + zi);
            }
        }
        byte[] b; int j, j0, bL; z = 0;
        for (; i < int.MaxValue; i++)
        {
            F2 = F1 + F0; F0 = F1; F1 = F2;
            b = F2.ToByteArray();
            bL = b.Length - 1;
            if (b[bL] == 0) bL--;
            j = 0;
            if (b[0] == 0)
            {
                z = 8;
                while (b[++j] == 0) z += 8; // add zerobits
                if (zi < z + 7)
                {
                    z += tz(b[j]);
                    if (zi < z) zi = z;
                }
                z = 0;
            }
            while (++j < bL)
            {
                if (b[j] == 0)
                {
                    j0 = j - 1;
                    z = 8;
                    while (b[++j] == 0) z += 8; // add zerobits
                    if (zi < z + 14)
                    {
                        z += lz(b[j0]);
                        if (zi < z + 7)
                        {
                            z += tz(b[j]);
                            if (zi < z) zi = z;
                        }
                    }
                    z = 0;
                }
            }
            if (zmax < zi)
            {
                zmax = zi;
                Console.WriteLine(i + " " + zi);
            }
        }
        Console.ReadLine();
    }
    private static int lz(byte b) // leading zeros, b > 0
    {
        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;
    }
    private static int tz(byte b) // trailing zeros
    {
        return 7 - lz((byte)(b & -(sbyte)b));
    }
}

Number of bits Fibonacci number

// Binet's formula: 
//
//         Fib(n) = ( Phi^n - (-1)^n / Phi^n ) / 5^(1/2)
//
// Where Phi = (1+5^(1/2))/2 ~ 1.61803398874989484820458683436563811772030917980576....(Knott)
//                           ~ 1.6180339887498948482045868343656381177203+             (Knuth)
//
// For larger values of n:
//
//         Fib(n) ~ Phi^n / 5^(1/2)
//
// For the number of bits, x, of Fib(n):
//
//            2^x ~ Fib(n) ~ Phi^n / 5^(1/2)
//
// Taking natural logarithms:
//
//        ln(2^x) ~ ln(Phi^n / 5^(1/2))
//
// So:
//      x * ln(2) ~ n * ln(Phi) - 1/2 * ln(5)
//              x ~ n * ln(Phi)/ln(2) - ln(5)/ln(2)/2
//
// Using "windows calculator" which seems to have a precision of ~105 bits ~31 decimal digits:
//
//              x ~ n * 0,6942419136306173017387902668986 - 1,1609640474436811739351597147447
//
// A partial bit (digit) is a bit (digit):
//
//        x ~ (int)(n * 0,6942419136306173017387902668986 - 1,1609640474436811739351597147447 + 1)
//
// The mantissa of a double has a precision of 52 bits ~ 16 decimal digits:
//
//        x = (int)(n * 0,6942419136306173 - 0,1609640474436812)
//
// Below it is correct for 0 <= n <= 24 * 10^6, after two days the output is: "R E A D Y"
// Also for n = 5 * 10^8 it's correct
//      for n = 10^9 it seems to be correct too.

using System;
using Xint = System.Numerics.BigInteger;
class bitLength_Fibonacci
{
    private static int bL_F(int n)
    {
        return n < 6 ? ++n / 2 : (int)(0.6942419136306173 * n - 0.1609640474436812);
    }

    static void Main()
    {
        Xint[] F = { 0, 1 };
        for (int n = 1; n <= 24000000; n++)
        {
            F = new Xint[2] { F[1], F[1] + F[0] };
            if (bL(F[0]) != bL_F(n))
            {
                Console.WriteLine(n + "   " + bL(F[0]) + "   " + bL_F(n));
            }
        }
        Console.WriteLine("R E A D Y");
        Console.ReadLine();
    }

    private static int bL(Xint X)
    {
        byte[] bytes = (X.Sign * X).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;
    }
}

Fibonacci sequence binary plot, Ed Pegg Jr.(2003)

At Wolfram Mathworld you can find an intriguing plot of the
Fibonacci sequence represented in binary.

The one below shows it for small value's, 

Fibonacci(5)=5,  101 binary,        black, white, black    

Fibonacci(6)=8, 1000 binary, black, white, white, white

  

A larger one:

Fibonacci(1000000000) Challenge in 0.1 ms

// Find the first 10 bytes and the last 10 bytes of f(1_000_000_000).
// That's right, fibonacci of one billion.

// I suppose the first 80 bits and the last 80 bits are more interesting,
// when you try to (dis)prove .... things.
// And, if you try, start with small numbers, first bit, last bit, first/last two bits, etc.
// By the way, the first bit is always one,
// though there is an exception, do you know it? (hint: start with small numbers)

// You can find the challenge at:
// weblogs.java.net/blog/kabutz/archive/2012/02/24/fibonacci-1000000000-challenge
// Or google "fibonacci 1000000000".

// My result(using 1 core):
//
//         F(1000000000)
//         bits: 694241913
//          MSB: 01 62 80 B8 2D 8C BE 0E DC 1B
//         time: 46 us
//          LSB: A9 53 2D F4 D2 D2 5B 5D B6 3B
//         time: 46 us

// Slightly less than 0.1 ms!!!!
// (the conversion of an 80 bits number to a string takes about 2 us, eternal question:
// can this conversion be done faster? I'm not going to answer it;)

// When the Most Significant Bytes are calculated with a precision of 80 bits,
// the lowest 2 or 3 bytes of these MSB's are wrong, if it's doubled to 160 bits,
// things look fine, are the same as those of Rex Young, which can be found at:
// weblogs.java.net/blog/kabutz/archive/2012/02/24/fibonacci-1000000000-challenge
// 
// Big question: "how to (dis)prove it?"
// It probably boils down to "more bits, higher accuracy, analyse the error". 
// Amazing: MSB's (multiplication of 160 bits numbers)
//      and LSB's (multiplication of  80 bits numbers)
//      take nearly the same time.
// I started with LowFib, returning the LSB's,
//   and thought HighFib, returning the MSB's,
//   would be 4 times slower.

// A few ways to speed it up:
// -upto a bitlength of 80 bits LowFib and HighFib do the same (iterations).
// -in the last iteration F(n) and F(n+1) are computed, get rid of F(n+1)
// -don't use a biginteger library, use arrays of uints,
//  and (Knuth's) classical multiplication algorithm.
// -translate it to c (or assembler)

// LSB's, 80 bits result from 80 bits multiplied by 80 bits.
// Looks like some time can be saved within the classical * algo.

// The final digits of fibonacci numbers repeat with a certain cycle length, see: 
// www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/fibmaths.html#cycles
// Looks like a modulo (remainder) operation.

// An underbound for the number of bits of fibonacci(n) is (2*n+1)/3
// Are there better bounds? 

// FibonacciHighLow(10^12)? The trillionth one.
// Use a few longs, it would surprise me if it takes more than 0.2 ms.

//       !!! Project >> Add Reference >> System.Numerics !!!
using Xint = System.Numerics.BigInteger;
using System.Diagnostics;
using System;

class FibonacciHighLow
{
    private static int nob;                           // number of bits
    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        for (int n = 1953125; n <= 1000000000; n *= 2)
        {
            Console.WriteLine("F(" + n + ")");
            string s = hexXint(HighFib(n));
            Console.WriteLine("bits: " + nob);
            Console.WriteLine(" MSB: " + s);
            HighFib(n);
            sw.Restart();
            for (int i = 0; i < 1000; i++)
            {
                HighFib(n);
                //hexXint(HighFib(n));               
            }
            sw.Stop();
            Console.WriteLine("time: " + sw.ElapsedMilliseconds + " us");
            s = hexXint(LowFib(n));
            Console.WriteLine(" LSB: " + s);
            LowFib(n);
            sw.Restart();
            for (int i = 0; i < 1000; i++)
            {
                LowFib(n);
                //hexXint(LowFib(n));
            }
            sw.Stop();
            Console.WriteLine("time: " + sw.ElapsedMilliseconds + " us");
            Console.WriteLine();
        }
        Console.ReadLine();
    }

    public static Xint LowFib(int n)
    {
        Xint Mask = (Xint.One << 80) - 1;
        Xint[] F = { 0, 1 };
        for (int i = fL2(n); i >= 0; i--)             // floorLog2
        {
            Xint L = 2 * F[1] - F[0];                 // L is a Lucas Number(2,1,3,4,7,..)
            switch ((3 << i & n) >> i)
            {
                case 0: F = new Xint[2] { L * F[0], L * F[1] - 1 }; break;
                case 2: F = new Xint[2] { L * F[0], L * F[1] + 1 }; break;
                case 1: F = new Xint[2] { L * F[1] - 1, L * (F[1] + F[0]) - 1 }; break;
                case 3: F = new Xint[2] { L * F[1] + 1, L * (F[1] + F[0]) + 1 }; break;
            }
            F[0] &= Mask;
            F[1] &= Mask;
        }
        return F[0];
    }
    // To do: 
    // case 0: F[1] = (L * F[1] - 1) & Mask;
    // if (F[1] < 0) F[1] = Mask; // ????
    // Looks like LowFib can be exact.
    // Test it with small Mask values (decimal 1,3,7,15,.../binary 1,11,111,1111,....)
    //    A few years ago I wrote Knuth's classical division algorithm in TWITS (TWo bITS)
    //    It's quite difficult to test it when bytes, ushorts or uints are used.
    //    With twits it's a lot easier.
    // Or, even better, don't, because:
    // mathworld.wolfram.com/FibonacciNumber.html
    // F(2^n + 2^(n+1)) ends in n+2 zeroes binary.
    // So 80 zeroes, n=78, 2^78+2^79, F(906.694.364.710.971.881.029.632) ends in 80 zeroes.

    public static Xint HighFib(int n)
    {
        Xint Mask = (Xint.One << 160) - 1;
        Xint[] F = { 0, 1 };
        int i = fL2(n);
        for (; i >= 0 && bL(F[0]) <= 160; i--)        // bitLength
        {
            Xint L = 2 * F[1] - F[0];
            switch ((3 << i & n) >> i)
            {
                case 0: F = new Xint[2] { L * F[0], L * F[1] - 1 }; break;
                case 2: F = new Xint[2] { L * F[0], L * F[1] + 1 }; break;
                case 1: F = new Xint[2] { L * F[1] - 1, L * (F[1] + F[0]) - 1 }; break;
                case 3: F = new Xint[2] { L * F[1] + 1, L * (F[1] + F[0]) + 1 }; break;
            }

        }
        int s = bL(F[0]) - 160;                       // shift
        int t = s;                                    // total shift
        F[0] >>= s;
        F[1] >>= s;
        for (; i >= 0; i--)
        {
            Xint L = 2 * F[1] - F[0];
            if ((1 << i & n) >> i == 0)
            {
                F[0] *= L; F[1] *= L;
            }
            else
            {
                F = new Xint[2] { L * F[1], L * (F[1] + F[0]) };
            }
            t *= 2;
            s = bL(F[0]) - 160;
            t += s;
            F[0] >>= s;
            F[1] >>= s;

        }
        t += bL(F[0]);
        nob = t;
        return F[0] >>= 80 + (8 - t & 7);
    }
    // mathworld.wolfram.com/FibonacciNumber.html
    // Edd Pegg Jr.'s plot(2003) might indicate that
    // repeating zeroes (binary) in Fibonacci Numbers,
    // in our case 160 zeroes, arise for F(VERY VERY LARGE)

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

    private static int bL(Xint X)
    {
        byte[] bytes = X.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 string hexXint(Xint X)
    {
        byte[] bytes = X.ToByteArray();
        byte[] zeros = new byte[10];
        int i = bytes.Length;
        if (i == 11) i--;
        for (--i; i >= 0; i--)
        {
            zeros[i] = bytes[i];
        }
        string s = zeros[9].ToString("X2");
        for (i = 8; i >= 0; i--)
        {
            s += " " + zeros[i].ToString("X2");
        }
        return s;
    }
}