標題看起來不太像是工程問題? 其實本質上是工程問題: '數'其實就是一堆石子數量的
符號,或代表石子本身. 電腦所處理的都是很實際的符號或'字串'. 沒有遐想空間 (其它
都是程式人員或數學家賦與的義意). 也因此,我主要聚焦工程上的方法: 如何處理'正確性'
這議題,但不多說明,因其它篇已談到.
前次文章 [C++ 類別到統一知識語言的想法 https://ithelp.ithome.com.tw/articles/10400528]
的考拉茲猜想部份有些錯誤,所以就有了這篇更新版.
文章不長,且含不少程式概念及技術展示, 加上考拉茲猜想易懂, 相信很快讀完及理解.
程式碼就是証明本身(Code is Proof 或 Code as Proof). 處理計算問題的第一步就是
要用程式(算法)表達問題,否則,你有可能正用電腦處理一個'不可計算'的問題.
+------+
| 前言 |
+------+
3x+1問題: 對於任意大於或等於1的整數n. 若n為奇數,則乘3加1. 若為偶數則除2.
問題是, 是否所有數經過這種計算一定會算至1? 以下嘗試証明此問題答案是肯定的.
皮亞諾公理系統難証此題(或會很長).
本文從算法模型(圖靈機與 C++ 狀態空間)出發探索 Collatz 猜想.透過建立
(a+b)/2^d 狀態轉移模型與窮舉樹狀搜尋，證明了「單靠奇偶運算與下降剪枝
(lt_org) 無法形成證明(ana 不終止，不存在 Bingo 狀態)」,從而揭示了純奇偶路徑
分析在解答 Collatz 猜想上的本質瓶頸。
考拉茲(Collatz)函數::=
wij2
int cop(int n) {
if(n<=1) {
if(n<1) {
throw Error;
}
return 1; // 1為疊代終點
}
if(n%2) {
return 3*n+1; // 奇數規則
} else {
return n/2; // 偶數規則
}
}
考拉茲數::= 若n, n∈N<1,+1>, cop疊代運算最後會算至1(即, cop(...cop(n))=1), 則n為
考拉茲數(記為C(n)),否則n不是考拉茲數.
考拉茲問題::= 是否對於所有整數n,n∈N<1,+1>, n為考拉茲數? 換句話說,問題等效於問
rcop程序會不會終止.
void rcop(int n) { // cop疊代
for(;n!=1;) {
n=cop(n);
}
}
命題: ∀n, n∈N<1,+1>, cop(n)疊代運算終止.
証: 程序rcop2 將rcop中的n等價分解成n=a+b形式表示.
void rcop2(int n) {
int pcnt=0,ocnt=0; // 僅用於証明解釋
Real D=1; // 同上
int a=n,b=0;
for(; a+b!=1;) { // 測量點A (a+b 等同cop疊代過程的n)
D*= 1.0/a;
if((a%2)!=0) {
--a; ++b; // a,b調整,使得a保持為偶數,以使以下算法能進行並保持對等
// 於cop(n)疊代.
}
D*= a;
// a/b測量點B
if((b%2)!=0) { // 等同(a+b)%2 (因a為偶)
a= 3*a;
b= 3*b+1; // 3*(a+b)+1= (3*a) +(3*b+1)
++ocnt;
}
a= a/2; // 每次疊代都會執行(a+b)/2
b= b/2;
++pcnt;
}
}
設n為奇數且無--a,++b過程並假設每次的奇運算搭配一個偶運算,則於測量點B可有:
a₁= n-1
aₖ= (3*aₖ₋₁)/2 =... = (n-1)*(3/2)ᵏ⁻¹
b₁= 1
bₖ= (3*bₖ₋₁+1)/2=... = 2*(3/2)ᵏ⁻¹ -1
aₖ/bₖ= (aₖ₋₁)/(bₖ₋₁) = ((n-1)*(3/2)ᵏ⁻¹)/(2*(3/2)ᵏ⁻¹ -1)
=... = (n-1)/(2- 1/(3/2)ᵏ⁻¹)
(1) aₖ/bₖ < aₖ₋₁/bₖ₋₁, 及 lim{k->∞} aₖ/bₖ= (n-1)/2.
因"約一半的極限值(n-1)/2"是遞迴成立的. 但此為估算值. aₖ,bₖ實際上是非連續
函數, 還會穿插a/2,b/2,--a,++b 運算使得aₖ更小,bₖ更大. 但無論如何,我們可
結論: 因每次(或幾乎)減半後都可視為另循環的開始, 因此,aₖ/bₖ值會遞減趨近0
(不是等於0,如同數列 1/10 -> 1/100 -> 1/1000,...).
(2) 若疊代循環, 則aₖ,bₖ可形成的組合數固定. 因除了當aₖ=0外, aₖ/bₖ 不可能無限地
(半)遞減. 因此,當aₖ≠0時,不會有循環.
(3) n的疊代過程可以乘某放大倍率的方式表示成 nₖ=n₀s₁s₂*...sₖ, sₖ= nₖ/nₖ₋₁
設S=s₁s₂*...sₖ, 則 nₖ=n₀S <=> S= nₖ/n₀
a疊代至1過程同様可表示成 aₖ=a₀*u₁*u₂*...*uₖ. 設U=u₁*u₂*...*uₖ, 則
U= aₖ/a₀ = 1/n₀. ... 公式1
設dₖ 表示 aₖ於調整運算段(--a)的放大率. dₖ= (aₖ₋₁-1)/aₖ₋₁, d<1
設u'ₖ 表示aₖ於奇偶運算段的放大率. uₖ=dₖ*u'ₖ, 公式1可寫成:
1/n₀= (d₁*u'₁)*(d₂*u₂')*...*(dₖ*uₖ') = D*U'
<=> D= 1/(n₀*U')
因U'可以a的奇偶運算比表達為 U'= 3^ocnt/2^pcnt, 故:
D= 1/(n₀*(3^ocnt/2^pcnt)) ... 公式2
....問題卡住了.
(4) 考慮到3x+1問題(命題)不太可能有第三個答案,即"不可判定". 因1. 問題的語義很
明確, 無抽象概念,沒出現自參問題. 2. 証"不可判定"需有一個無限迴圈. 若此能
証"不可判定",則等同証明了[命題]為假(存在n使得cop(n)疊代運算不終止).
因此猜想3x+1猜測的証明可能是種非線性程序証明. 因此構想由1開始逆向搜尋可能
疊代的數. 看由開始最大的連續數増長情形. 若能保証此連續區間[1,m]連續數増長,
則証明3x+1猜測成立.
VLInt narr(1u); // 設bit 0為1. VLInt 為一非負整數的大數(Bit Number)類別
int n_max= 1000000000; // 搜尋最大數值限制
int ffc_max= narr.ffc(); // 記錄narr第一個0的位置
// [Syn] 由n開始,逆向遞迴搜尋可能疊代的數並記錄所有產生的數. 最多level個
// 數.
//
void explore_coll(int n, int level) {
if((n<1)||(level<=0)||(n>n_max)) {
return;
}
narr._set_bit(n,true); // 記錄n對應的bit為1
size_t tfc=narr.ffc(); // 由bit-0位置開始找第一個0的位置
if(tfc>ffc_max) {
ffc_max=tfc;
cout << tfc-2 << ' '; // m= tfc-2. 每次刷新時顯示
}
explore_coll(2*n,level-1); // 偶分枝
int t= (n-1)/3;
if((3*t+1==n)&&(t%2)&&(n!=4)) {
explore_coll(t,level-1); // 奇分枝
}
};
int main()
{
explore_coll(1,400);
cout << "Ok\n";
return 0;
};
$ ./a.out
1 2 5 11 23 47 95 191 383 767 1535 3071 6143 12287 24575 49151 77670 Ok
這是有意思的(約)二倍增序列.
對於P= ∀n∈ℕ, p(n) 類的命題P必能以數學歸納法証明,但看起來這句話的義意
不直接: P不必依傳統Peano公理的形式的數學歸納法進行.
譬如,以'樹'証P的方式可先定義不同的'自然數'結構:
N<1,T>::= 1. 1 ∈ N<1,T>
2. 若n∈N<1,T>,則 2*n 及2*n+1 ∈N<1,T> // T的後繼有兩個數
N<0,S>::= 1. 0∈ N<0,S>
2. 若n∈N<0,S>,則 3*n+1, 3*n+2, 3*n+3 ∈N<0,S> // S有三個後繼
註: 這種'自然數'的看法有局限,不能取代正式的定義,除非還有更多條件.
(5) 設codd(n)為奇數的Collatz疊代(多對1)函數.
// [Syn] 計算 (3*n+1)/2^k (3n+1後除掉2因子)
int codd(int n) {
if(((n%2)==0)||(n<=1)) {
throw Error; // n 非大於1的奇數
}
n= 3*n+1;
for(;(n%2)==0;) { // 除去2因子
n=n/2;
}
return n;
}
F(n) ::= { m | m 為大於1的奇數，codd(m)=n }
性質1: n≠n' => F(n)與F(n')不相交. n=n' => F(n)=F(n').
証: 若n≠n',F(n)與F(n')相交, 則存在一個奇數x同時滿足 (3*x+1)/2^k=n <=>
3x+1= n*2^k 及 (3*x+1)/2^k'= n' <=> 3*x+1= n'*2^k'.
因此 <=> n*2^k = n'*2^k' => 不成立(不同奇數n,n'含不同的非2因子, 靠各増減
2因子不可能使等式成立). 性質的反向証法也成立(若n=n',則取k=k'即可).
性質2: 若n為3的倍數,則F(n)為空集.
証: F(3*n) => 3*n= (3*m+1)/2^k => n= (3*m+1)/(3*2^k) => 3不可能整除分子
定義Lev 為一個以奇數為輸入的決定性程序:
Lev(1)::= 0
Lev(n)::= Lev(codd(n))+1
codd疊代循環? 若某codd的循環疊代表示為 a -> b -> c -> a, 則有
Lev(a)=h -> Lev(b)=h-1 -> Lev(c)=h-2 -> Lev(a)=h-3.
Lev(a) 出現了矛盾的多值狀況(即h與h-3).
但, 若Lev 稱為'函數'有個漏洞: 是否存在某x, Lev(x)=無限大?
其實Lev值就是疊代次數pcnt值,無腦地把問題搞得過於'學究'是自找麻煩.
(6) 了解3x+1的運作方式後, 其整個過程可看作是一個有a,b兩個tape的圖靈機. 只是
a-tape的'讀取'方式較特殊一點:
void rcodd2(int n) {
if(((n&1)==0)||(n<3)) {
throw Error; // n非3以上的奇數
}
int a=n-1,b=1;
for(;a+b>1;) {
a= 3*a; // a-tape的(偶)值乘3
b=4;
do {
a= a/2; // a值右移(或a謮寫頭左移)
b= b/2; // b值右移
if((a&1)!=0) { // 讀a的bit-0
--a; ++b;
}
// Cond: a偶,b==1,2,3
} while((b&1)==0);
// Cond: a偶,b==1 或3
a+= b-1;
}
};
// 看看a的變化
void a3c(int n) {
int a=n;
for(;a>17;) { // 有循環 1-1, 7-5-7, 17-25-37-55-41-61-91-17
a= 3*a;
do {
a>>=1;
} while((a&1)==0); // 右移至a為奇數
}
};
小結: 等價3x-1(或原3x+1)問題,仍是原地打轉.
(7) 考慮'長步'方式: 每次疊代|n|次 (外加除去尾0)
void rcodd3(int n) {
int m=n;
for(;m>1;) {
int blen= hbs(m); // blen為m的有效bit數
for(int i=0; i<blen; ++i) { // 次數可調整,不影響結果
if(m<=1) {
break;
}
if(m&1) {
m= 3*m+1;
}
m= m/2;
}
for(;(m&1)==0;) { // 保持m為奇數以便循環分析(但不改變疊代可預測性)
m=m/2;
}
}
};
小結: 與codd類同. '長步'的次數可加倍使之更長,但已難解釋了.
由實測檢視(n<10^9),n最大的疊代值小於n^2.8, n=27時的長度倍率約2.8,
但隨n的增大,最大的疊代值的長度很快下降至約1.2倍. 以上這些顯示,拉茲猜想
可能可証,只是'証明'很難看懂(類同証明某亂數產生器或chaos序列對所有啓始的
自然數n最後都會出現1). 這種問題原則上可用電腦証明?
(8) 概念(concept)即class, 我們以C++分析一下這種想法.
/* CoNum.h */
#include <Wy.stdlib.h>
#include <CSCall/VLInt.h>
namespace Wy {
// [Syn] CoNum類別以 (a+b)/2^d 形式表示非負整數N<1,+1>, 其中a為3^y形式.
// CoNum類別也可表示集合{x| 2^d 整除 ax+b} 或{ minx+i2^d ∣ i∈N<0,+1>}
//
// 註: VLInt是個大(非負)整數類別
// 註: CoNum是以ClassGuidelines準則寫的.
//
class CoNum {
typedef VLInt::UIntType UIntType; // unsigned int (依VLInt的定義)
VLInt m_a, m_b;
UIntType m_d, m_y;
public:
WY_DECL_REPLY;
// [Syn] 建構Default(或稱原值)物件 (可表1或數集{1,2,3,...}).
//
CoNum() : m_a(1u), m_b(0u),m_d(0),m_y(0) {};
CoNum(const CoNum& s) : m_a(s.m_a), m_b(s.m_b), m_d(s.m_d), m_y(s.m_y) {};
CoNum(CoNum& s, ByMove_t) : m_a(s.m_a,ByMove), m_b(s.m_b,ByMove),
m_d(s.m_d), m_y(s.m_y) {
};
CoNum& operator=(const CoNum& s) {
CoNum(s).swap(*this);
return *this;
};
bool is_default() const {
return m_a.is_one()&&m_b.is_zero()&&(m_d==0)&&(m_y==0);
};
// [Syn] 計算 (a*n+b)/2^d 的整數值
// n必須是能產生此物件奇偶路徑的數,否則結果未定義
//
VLInt val(const VLInt& n) const {
VLInt v(m_a*n+m_b);
v>>=m_d;
return v;
};
const VLInt& a() const { return m_a; };
const VLInt& b() const { return m_b; };
UIntType d() const { return m_d; };
UIntType y() const { return m_y; };
// [Syn] 設定data成員.
//
// 註: 此成員為測試用途. 大部份的a,b,d,y數值組合不能由合法的odd(),even()
// 序列產生.
//
void _set_abdy(const VLInt& a, const VLInt& b, UIntType d, UIntType y) {
m_a=a;
m_b=b;
m_d=d;
m_y=y;
};
// [Syn] 執行 3n+1 奇運算
//
CoNum& odd() {
m_a*=3u; // m_a= 3^ocnt
m_b*=3u;
Errno r;
VLInt m;
if((r=m.set_bit(m_d,true))!=Ok) {
WY_THROW( Reply(r) );
} // m= 2^m_d (m_d=0時 m=1)
m_b+= m;
++m_y;
return *this;
};
// [Syn] 執行 (3n+1)/2 運算
//
CoNum& oddeven() {
odd();
even();
return *this;
};
// [Syn] 執行 n/2 偶運算
//
CoNum& even() {
if(m_d>=Limits<UIntType>::Max) {
WY_THROW( Reply(ERANGE) );
}
++m_d;
return *this;
};
// [Syn] 求最小的x, x∈N<1,+1>, 使得 (a*x + b)/ 2^d 為整數. 即, CoNum所表
// 集合中最小的數.
//
// 註: 此算法及實作是由AI產生
//
VLInt minx() const {
Errno r;
VLInt m;
if((r=m.set_bit(m_d,true))!=Ok) {
WY_THROW( Reply(r) );
} // m= 2^m_d (m_d=0時 m=1)
// 擴展歐幾里得(全非負): r0=a, r1=m ; s0=1, s1=0
// 遞迴: q = r0/r1, r2 = r0%r1, s2 = s0 + q*s1 <- 只有加、乘
// 正負號規律: 跑到第 i 輪(即 r1 變成 r_i)時, sign(s_{i-1}) = + 若 i 為
// 奇數, 否則 -
VLInt r0 = m_a, r1 = m;
VLInt s0(1u), s1(0u);
unsigned i = 1;
while (r1.is_zero()==false) {
VLInt q = r0 / r1;
VLInt r2 = r0 % r1;
VLInt s2 = s0 + q * s1;
r0 = r1; r1 = r2;
s0 = s1; s1 = s2;
++i;
}
// 迴圈結束: r0 = gcd(a,m) (恆為1,因a恆為3的冪、必為奇數)
// s0 = |反元素係數|
VLInt uMag = s0 % m;
VLInt u = (i % 2 == 1) ? uMag : (m - uMag); // a*u ≡ 1 (mod m)
// 解 a*x+b ≡0 (mod m) 的最小正整數 x
VLInt t = (m_b * u) % m;
VLInt x0(m);
if(t.is_zero()==false) { // t=0時取x=m(而非0),符合x∈{1,2,...}
x0 = m - t;
}
//return (m_a * x0 + m_b) / m;
return x0;
};
// [Syn] 判斷原值(實際疊代的值)經疊代後是否低於原值?
//
// [Ret] true: 低於原值
// false: 不低於原值
//
// (3^y*n+b)/2^d < n
// <=> 3^y*n+b < n*2^d
// <=> b < n*(2^d- 3^y) // 見說明
// <=> b/(2^d- 3^y) < n
//
// 設c=2^d-3^y:
// c>0, 命題總成立
// c=0, 命題不成立
// c<0, 命題不成立 (對所有n>=1, n*c<=c<0<=b, 故b<n*c恆假)
//
bool lt_org() const {
return m_d >= m_a.hbs(); // m_a.hbs() 為m_a的最高位bit為1的位置加1
// 等同判斷 2^d > m_a
};
void swap(CoNum& ano) {
this->m_a.swap(ano.m_a);
this->m_b.swap(ano.m_b);
Wy::swap(this->m_d, ano.m_d);
Wy::swap(this->m_y, ano.m_y);
};
};
String wrd(const CoNum& n) {
String str;
str << '(' << wrd(n.a()) << '+' << wrd(n.b())
<< ")/2^" << wrd(n.d()) << ", minx=" << wrd(n.minx())
<< " (" << (n.lt_org()? 'T':'F') << ')' ;
return str;
};
}; // end namespace Wy
/* t_idea.cpp
以奇偶運算作為探索域,處理Collatz猜測的証明.
ana() 預期兩種的終止方式：
path 下降，可以歸入較小的問題而証終止；
找到可覆蓋全域的 bingo 狀態。
但都不成立,因此,可能的答案是: 單靠奇偶運算不足以証Collatz猜測.
註: 証明0.999...(100%)的自然數都終結很容易,但0.999...!=1.
*/
#include <Wy.stdio.h>
#include "CoNum.h"
using namespace Wy;
bool hard_limit(const CoNum& n) { // 監測資源消耗
// 若存在Collatz猜測成立的証明,則存在這種限制.
// 目前可實測的n的疊代不超過n^2.8, 趨勢是n^1.2 (不過對此証法似乎沒用).
return false; // 目前不確定
};
bool is_bingo(const CoNum& n) {
// 若存在'bingo'物件n, 則 n.is_deault()==true 必成立才能判斷 ∀n∈N<1,+1>,
// n∈ Set(CoNum(n)). 但這是不可能的 (n.d()嚴格增加).
return false; // a,b,d含完整的路徑資訊. 純函式(pure function)方式無法判定.
};
// [Syn] ana(..) 探索可能的oddeven/even路徑,搜尋rcop必終結的証據.
//
// 目前結果: bing物件不存在,lt_org不可能全剪枝 (細節長), 因此結論: ana不終結.
//
// 註: 搜尋域封閉: 所有有限長度的oddeven/even路徑皆探索到.
//
void ana(const CoNum& n) try {
if(hard_limit(n)) {
cout << "Hard limit reached" WY_ENDL;
WY_THROW( Errno() ); // "n的疊代不超過n^2.8" 的猜測錯誤,或其它.
}
if(n.lt_org()) {
return; // 此路徑已下降，停止探索. lt_org 作為此用途有小瑕疵, 但不影響
// ana架構.
}
if(is_bingo(n)) {
throw "Bingo"; // ana全域終結
}
CoNum t(n);
t.even();
ana(t);
t=n;
t.oddeven();
ana(t);
}
catch(const char* msg) {
cout << msg << WY_ENDL;
};
int main(int argc, const char* argv[])
try {
CoNum n;
ana(n);
cout << "OK" WY_ENDL;
return 0;
}
catch(const Errno& e) {
cerr << wrd(e) << WY_ENDL;
return -1; // e.c_errno();
}
catch(...) {
cerr << "main() caught(...)" WY_ENDL;
throw;
};
iThome鐵人賽
看影片追技術