Gosper 的序列 循环检测

 

// 高效程序的奥秘hacker's delight Henry S. Warren, Jr. 著冯速译

// 第5 章位计数5.4 后缀0 计数

// 假设序列X0, X1, X2, ... 是由Xn+1 = f(Xn) 定义的。如果f 的值域是有限的，那么序列必定是循环的。

// 给定函数f , 循环检测问题就是找到第一个重复出现的元素的下标μ和周期λ.

// 参考:  http://home.pipeline.com/~hbaker1/hakmem/flows.html#item132  LOOP DETECTOR

// Knuth, Donald E. The Art of Computer Programming, Volume 2, Third Edition 3.1 节的问题6

 

 

#include "iostream"

#include "bitset"

#include "limits"

#include "time.h"

 

using namespace std;

 

#define out2(T) cout<<bitset<numeric_limits<unsigned int>::digits>(T)<<endl;

#define SIZE(T) (numeric_limits<unsigned T>::digits)

 

bool SHR31(int x)           //逻辑右移位

{

     return (unsigned int)x>>31;

}

 

#define HOUT(x,y,z)    hex_out(x,y,z)

 

inline void hex(int x)

{

     int w = cout.width();

     cout << "0x";

     cout.width(8); cout.fill('0');

     cout << hex << x << "  ";

     cout.width(w);

}

inline void hex(unsigned __int64 x)

{

     int w = cout.width();

     cout << "0x";

     cout.width(16); cout.fill('0');

     printf("%I64x   ",x);

     cout.width(w);

}

void hex_out(int x, int y, int z)

{

     hex(x); hex(y); hex(z); cout << endl;

}

 

int nlz(unsigned int x) {

     int n;

     if (x == 0) return 32;

     n = 1;

     if ((x >> 16) == 0) {n = n +16; x = x <<16;}

     if ((x >> 24) == 0) {n = n + 8; x = x << 8;}

     if ((x >> 28) == 0) {n = n + 4; x = x << 4;}

     if ((x >> 30) == 0) {n = n + 2; x = x << 2;}

     n = n - (x>>31);

     return n;

}

 

int ntz(unsigned int x)

{

     int n;

 

     if (x == 0) return 32;

     n = 1;

     if ((x & 0x0000FFFF) == 0) {n +=16; x >>= 16;}

     if ((x & 0x000000FF) == 0) {n += 8; x >>= 8;}

     if ((x & 0x0000000F) == 0) {n += 4; x >>= 4;}

     if ((x & 0x00000003) == 0) {n += 2; x >>= 2;}

     return n - (x & 1);

}

 

// 平方取中

int f(int a)

{

     return (a*a / 120) % 14400;

}

 

// 这里分五步论证算法的正确性：

// 0. f值域有限，则f 函数是循环，参考knuth 的证明

// 1. 设前面存的T[] 最大位置p, 现在查询最大位置kmax, 总有p = kmax

//    证：先证p 不小于且不大于kmax, p 存放的位置每次最大发生在2^i 处,

//         如: 10, 100, 1000, 10000 ...

//         设位长W = 8, 则 1000B 处是交界, ntz(1000B) = 3; kmax = 8-1 - nlz(1000B) = 7 - 4 = 3; 若i < 1000B, kmax <= p

//         若i > 1111B , 如i = 10000B, 则p = 4, kmax = 8-1 - nlz(10000B) = 4, kmax = p, 即kmax 随p 更新

//         又p 取所有位置的最大位置所以总有p = kmax, 不会造成访问越界，也不会造成访问遗漏.

// 2. T[] 数组有些元素总会被后来的元素覆盖, 但个别元素在一个周期内不会被覆盖.如果存在这种元素(称为异元素)，

//     则总能在第二个周期λ中检测到循环

//     证：如奇数，由于ntz(奇数) = 0, 前一个奇数总被后一个奇数覆盖

//         若偶数，是不是在一个周期λ内有偶数不被覆盖，且不覆盖别的偶数呢？这一点很重要, 设这种偶数为异元素

//         因为被覆盖的元素在第二个周期都会反过来覆盖先前覆盖它的元素,

//         而到第二个周期λ检测中，这种异元素不会被反覆盖，当再一次走到这个值时，就可以检测到相同,

//         所以总能在第二个周期λ中检测到循环

// 3. 这种异元素总是存在的

//     证：举例λ= 10110101B, 从μ位置开始循环, 如果μ= 0, 则10000000B 位置为异元素, 因为ntz(10000000B)

//         是i < 11111111B 中最大的了, 如果 μ= 01001011B,  μ+ λ= 100000000B 这个元素独一无二的，

//         在下一个周期内μ+ λ+ λ< 1000000000B, 不可能有更大的p, 又p 先于kmax 更新，所以在一个周期内

//         这种独一无二的元素总是存在的

// 4. 对于周期λ和μ的求解也就迎刃而解了，可参考ITEM 132 和下面的程序

 

 

// 找出序列 的周期 λ (lambda) 和 开始位置 μ 的上下界 mu_l , mu_u

void ld_Gosper(int (*f)(int), int X0, int *mu_l, int *mu_u, int *lambda)

{

     int Xn, k, m, kmax, n, lgl;

     int T[33];

 

     T[0] = X0;

     Xn = X0;

     for (n = 1; ; n++) {

         Xn = f(Xn);

         kmax = 31 - nlz(n);                       // Floor(log2 n).

         for (k = 0; k <= kmax; k++) {

              if (Xn == T[k]) goto L;

         }

         T[ntz(n+1)] = Xn;                    // No match.

     }

L:

     // Compute m = max{i | i < n and ntz(i+1) = k}.

 

     m = ((((n >> k) - 1) | 1) << k) - 1;

     *lambda = n - m;

     lgl = 31 - nlz(*lambda - 1);              // Ceil(log2 lambda) - 1

     *mu_u = m;                                     // Upper bound on mu.

     *mu_l = m - __max(1, 1 << lgl) + 1;       // Lower bound on mu.

}

 

 

 

void test1();

void main()

{

     test1();

}

 

void test2()

{

     int (*p)(int) = f;

     int X0 = 32690;

     int *mu_l = new int; int *mu_u = new int; int *lambda = new int;

 

     ld_Gosper(p, X0, mu_l, mu_u, lambda);

 

     cout << (*mu_l) << " " << (*mu_u) << " " << (*lambda) << endl;

 

     delete mu_l; delete mu_u; delete lambda;

}

 

void test1()

{   

     int (*p)(int) = f;

     int X0 = 286;//123456;//2375;

     int *mu_l = new int; int *mu_u = new int; int *lambda = new int;

 

//   ld_Gosper(p, X0, mu_l, mu_u, lambda);

 

//   cout << (*mu_l) << " " << (*mu_u) << " " << (*lambda) << endl;

 

     int max = 0; int q = X0;

 

     for (int j = 1; ;j++) {

         X0 = j;

         ld_Gosper(p, X0, mu_l, mu_u, lambda);

         if (*lambda >= max) { max = *lambda; q = X0;}

         if (j >= 0x7FFF) break;

     }

 

     ld_Gosper(p, q, mu_l, mu_u, lambda);

     cout << (*mu_l) << " " << (*mu_u) << " " << (*lambda) << endl << endl;

 

     X0 = q;

     for (int i = 0; i < 50; i++) {

         cout << X0 << endl;

         X0 = p(X0);

     }

 

     delete mu_l; delete mu_u; delete lambda;

}