Gosper 的序列 循环检测
?
// 高效程序的奧秘hacker's delight Henry S. Warren, Jr. 著馮速譯
// 第5 章位計(jì)數(shù)5.4 后綴0 計(jì)數(shù)
// 假設(shè)序列X0, X1, X2, ... 是由Xn+1 = f(Xn) 定義的。如果f 的值域是有限的，那么序列必定是循環(huán)的。
// 給定函數(shù)f , 循環(huán)檢測(cè)問(wèn)題就是找到第一個(gè)重復(fù)出現(xiàn)的元素的下標(biāo)μ和周期λ.
// 參考:? 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 節(jié)的問(wèn)題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 函數(shù)是循環(huán)，參考knuth 的證明
// 1. 設(shè)前面存的T[] 最大位置p, 現(xiàn)在查詢最大位置kmax, 總有p = kmax
//??? 證：先證p 不小于且不大于kmax, p 存放的位置每次最大發(fā)生在2^i 處,
//?????? ? 如: 10, 100, 1000, 10000 ...
//?????? ? 設(shè)位長(zhǎng)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, 不會(huì)造成訪問(wèn)越界，也不會(huì)造成訪問(wèn)遺漏.
// 2. T[] 數(shù)組有些元素總會(huì)被后來(lái)的元素覆蓋, 但個(gè)別元素在一個(gè)周期內(nèi)不會(huì)被覆蓋.如果存在這種元素(稱為異元素)，
//?? ? 則總能在第二個(gè)周期λ中檢測(cè)到循環(huán)
//?? ? 證：如奇數(shù)，由于ntz(奇數(shù)) = 0, 前一個(gè)奇數(shù)總被后一個(gè)奇數(shù)覆蓋
//?????? ? 若偶數(shù)，是不是在一個(gè)周期λ內(nèi)有偶數(shù)不被覆蓋，且不覆蓋別的偶數(shù)呢？這一點(diǎn)很重要, 設(shè)這種偶數(shù)為異元素
//?????? ? 因?yàn)楸桓采w的元素在第二個(gè)周期都會(huì)反過(guò)來(lái)覆蓋先前覆蓋它的元素,
//?????? ? 而到第二個(gè)周期λ檢測(cè)中，這種異元素不會(huì)被反覆蓋，當(dāng)再一次走到這個(gè)值時(shí)，就可以檢測(cè)到相同,
//?????? ? 所以總能在第二個(gè)周期λ中檢測(cè)到循環(huán)
// 3. 這種異元素總是存在的
//?? ? 證：舉例λ= 10110101B, 從μ位置開始循環(huán), 如果μ= 0, 則10000000B 位置為異元素, 因?yàn)?span lang="EN-US">ntz(10000000B)
//?????? ? 是i < 11111111B 中最大的了, 如果 μ= 01001011B,? μ+ λ= 100000000B 這個(gè)元素獨(dú)一無(wú)二的，
//?????? ? 在下一個(gè)周期內(nèi)μ+ λ+ λ< 1000000000B, 不可能有更大的p, 又p 先于kmax 更新，所以在一個(gè)周期內(nèi)
//?????? ? 這種獨(dú)一無(wú)二的元素總是存在的
// 4. 對(duì)于周期λ和μ的求解也就迎刃而解了，可參考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;
}
?
?
總結(jié)
以上是生活随笔為你收集整理的Gosper 的序列 循环检测的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇： 种群计数 (pop_count)
- 下一篇： 最高标号预留与推进算法 --- 就是要比