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;
}
?
?
總結
以上是生活随笔為你收集整理的Gosper 的序列 循环检测的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 种群计数 (pop_count)
- 下一篇： 最高标号预留与推进算法 --- 就是要比