前回の漸化式の強引な解法(※失敗談)

このブール代数の漸化式を誰か解いてください! - Plus Le Toolの漸化式を数学的な方法で解くのは(少なくとも自力では)無理っぽいので、別の強引な方法で無理矢理解くことにした。

なお、以下で使う記号の定義については[id:plusletool:20130427:p2]を参照。
また、以下では独自に定義した関数 \mathfrak{C}_{r} を使うが、これについては[id:plusletool:20130612:p1]を参照。


「漸化式を解け」という場合、通常は一般解を求めることを要求している。
しかし前回の問題の場合は 1\le i\le 32 の範囲について c_{i} を求めれば充分だった。
つまり極端な話、 c_{1} から c_{32} を手計算で求めてしまっても構わないのだ。
……といっても、実際にやってみると c_{6}c_{7} あたりで余裕で挫折できる程度の項数になるので、さすがに手計算ではやってられない。
こういう単純計算はコンピュータに任せるべきだ。

というわけで今回は、前回の漸化式を解くプログラムを書いて強引に解くことにする。


さて、前回解きたかったのは次の漸化式だった。

[id:plusletool:20130616:p1](★1)再掲
\begin{array}{lclcl}c_{1}&=&0\\c_{2}&=&\mathfrak{C}_{2}\left(\Phi_{1}\right)\\c_{3}&=&\mathfrak{C}_{2}\left(\Phi_{2}\right)\,\oplus\,\mathfrak{C}_{1}\left(\Phi_{2}\right)\mathfrak{C}_{2}\left(\Phi_{1}\right)\\c_{i+1}&=&\mathfrak{C}_{2}\left(\Phi_{i}\right)\,\oplus\,\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\,\oplus\,\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-2}\right)\\&&\,\oplus\,\left(\mathfrak{C}_{1}\left(\Phi_{i}\right)\,\oplus\,\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-2}\right)\right)c_{i}\\&&\,\oplus\,\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i-1}\;\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i}c_{i-1}&\;&\left(\,\mathrm{as}\;i\ge 3\,\right)\end{array}

このまま安直に実装して計算させると項数が多すぎて c_{12} くらいで System.OutOfMemoryException吐いて落ちた(経験談)ので、一工夫加える必要がある。

まずは補題:32bit符号なし整数変数間の算術和演算をブール関数で表す - Plus Le Toolの【失敗例】と同様に、 c_{i}' を次のように定義する。

(*1)
\begin{array}{lcl}c_{1}'&=&0\\c_{i+1}'&=&c_{i+1}\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i}\right)\end{array}

これを(★1)に代入して整理すると次のようになる。

(*2)
\begin{array}{lclcl}c_{1}'&=&0\\c_{2}'&=&0\\c_{3}'&=&\mathfrak{C}_{1}\left(\Phi_{2}\right)\mathfrak{C}_{2}\left(\Phi_{1}\right)\\c_{i+1}'&=&\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\,\oplus\,\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\mathfrak{C}_{1}\left(\Phi_{i}\right)\\&&\,\oplus\,\left(\mathfrak{C}_{1}\left(\Phi_{i}\right)\,\oplus\,\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\right)c_{i}'\\&&\,\oplus\,\left(\mathfrak{C}_{1}\left(\Phi_{i}\right)\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\right)\mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i-1}'\;\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i}'c_{i-1}'&\;&\left(\,\mathrm{as}\;i\ge 3\,\right)\end{array}

補題:32bit符号なし整数変数間の算術和演算をブール関数で表す - Plus Le Toolの最後にも書いたように、これを解く際、 \mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i}'c_{i-1}' の項が非常に邪魔だった。
なのでこの項をやや強引に消してしまおう。
そのためにとりあえず、 b_{i} を次のように定義する。

(*3)
\begin{array}{lclcl}b_{1}&=&0\\b_{2}&=&0\\b_{i}&=&\mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i}'c_{i-1}'&\;&\left(\,\mathrm{as}\;i\ge 3\,\right)\end{array}

このとき、(*2)の漸化式部分は次のようになる。

(*4)
\begin{array}{lcl}c_{i+1}'&=&\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\,\oplus\,\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\mathfrak{C}_{1}\left(\Phi_{i}\right)\\&&\,\oplus\,\left(\mathfrak{C}_{1}\left(\Phi_{i}\right)\,\oplus\,\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\right)c_{i}'\;\,\oplus\,\left(\mathfrak{C}_{1}\left(\Phi_{i}\right)\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\right)\mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i-1}'\;\,\oplus\,b_{i}\end{array}

また、(*2)(*3)より b_{i+1} は次のようになる。

(*5)
\begin{eqnarray}b_{i+1}&=&\mathfrak{C}_{2}\left(\Phi_{i-1}\right)c_{i+1}'c_{i}'\\&=&\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\left({\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\,\oplus\,\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\mathfrak{C}_{1}\left(\Phi_{i}\right)\\\;\;\,\oplus\,\left(\mathfrak{C}_{1}\left(\Phi_{i}\right)\,\oplus\,\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\right)c_{i}'\\\;\;\,\oplus\,\left(\mathfrak{C}_{1}\left(\Phi_{i}\right)\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\right)\mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i-1}'\;\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i}'c_{i-1}'}\right)c_{i}'\\&=&\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\mathfrak{C}_{4}\left(\Phi_{i-2}\right)c_{i}'\,\oplus\,\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i}'c_{i-1}'\\&=&\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\mathfrak{C}_{4}\left(\Phi_{i-2}\right)c_{i}'\,\oplus\,\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-1}\right)b_{i}\end{eqnarray}

これで、邪魔だった \mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i}'c_{i-1}' の項を(見た目上は)消去できた。
また、 c_{*}うしの論理積がなくなったことで計算量が大幅に減り、その分だけメモリ使用量も減った(はず)。

さて、ここまでの結果をまとめると次のようになる。

(*6)
\begin{array}{rclcrcl}b_{1}&=&0&\,,\,&c_{1}'&=&0\\b_{2}&=&0&\,,\,&c_{2}'&=&0\\b_{3}&=&0&\,,\,&c_{3}'&=&\mathfrak{C}_{1}\left(\Phi_{2}\right)\mathfrak{C}_{2}\left(\Phi_{1}\right)\end{array}
\begin{array}{lclcl}b_{i+1}&=&\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\mathfrak{C}_{4}\left(\Phi_{i-2}\right)c_{i}'\,\oplus\,\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-1}\right)b_{i}&\;&\left(\,\mathrm{as}\;i\ge 3\,\right)\\c_{i+1}'&=&\mathfrak{C}_{1}\left(\Phi_{i}\right)\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\,\oplus\,\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\mathfrak{C}_{1}\left(\Phi_{i}\right)\\&&\,\oplus\,\left(\mathfrak{C}_{1}\left(\Phi_{i}\right)\,\oplus\,\mathfrak{C}_{4}\left(\Phi_{i-2}\right)\right)c_{i}'\;\,\oplus\,\left(\mathfrak{C}_{1}\left(\Phi_{i}\right)\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i-1}\right)\right)\mathfrak{C}_{2}\left(\Phi_{i-2}\right)c_{i-1}'\;\,\oplus\,b_{i}&\;&\left(\,\mathrm{as}\;i\ge 3\,\right)\end{array}

\begin{array}{lcl}c_{1}&=&0\\c_{i+1}&=&c_{i+1}'\,\oplus\,\mathfrak{C}_{2}\left(\Phi_{i}\right)\end{array}

この(*6)の連立漸化式を実装して解くことにする。
(もっと効率がいい式に変形できるかもしれないけど、それはまた後で考える。)
とりあえずスパゲッティソース貼ってみる。

【ソース】

using System;
using System.Collections.Generic;
using System.Text;

namespace Test {
	public class Calculator {
		#region static constructor
		static Calculator() {
			try {
				encode = Encoding.GetEncoding("shift_jis");
			} catch( ArgumentException ) {
				encode = Encoding.Default;
			}
		}
		#endregion

		// static field
		private const string dirPath = @"result\";
		private static readonly Encoding encode;
		// field
		private readonly bool doWriteFile = true;
		private readonly bool doOpenFile = false;
		private readonly bool doDisplay = false;
		private readonly bool doGC = true;
		private readonly bool doNormalizeTerm = true;

		// method
		public void Run() {
			Term[] cPrev = new Term[] { Term.Zero };	// c'[i-1]
			Term[] cTemp = new Term[] { Term.Zero };	// c'[i]
			Term[] bTemp = new Term[] { Term.Zero };	// b[i]

			for( int i=-1;i<32;i++ ) {
				List<Term> cNext = new List<Term>(3 + 2 * cTemp.Length + 2 * cPrev.Length + bTemp.Length);	// c'[i+1]
				List<Term> bNext = new List<Term>(cTemp.Length + bTemp.Length);		// b[i+1]

				// c'[i+1]
				if( i >= 1 ) {
					Term t = Term.GetTerm(0u,1u<<( i-1 ),1u<<i);	// C1[i]C2[i-1]
					if( !t.Equals(Term.Zero) ) {
						cNext.Add(t);
					}
				}
				if( i >= 2 ) {
					{
						Term t = Term.GetTerm(1u<<( i-2 ),1u<<( i-1 ),0u);	// C2[i-1]C4[i-2]
						if( !t.Equals(Term.Zero) ) {
							cNext.Add(t);
						}
					}
					{
						Term t = Term.GetTerm(1u<<( i-2 ),0u,1u<<i);	// C1[i]C4[i-2]
						if( !t.Equals(Term.Zero) ) {
							cNext.Add(t);
						}
					}
				}
				foreach( Term tT in cTemp ) {
					if( i >= 0 ) {
						Term t = Term.GetTerm(tT.C4,tT.C2,tT.C1|( 1u<<i ));			// C1[i]c'[i]
						if( !t.Equals(Term.Zero) ) {
							cNext.Add(t);
						}
					}
					if( i >= 2 ) {
						Term t = Term.GetTerm(tT.C4|( 1u<<( i-2 ) ),tT.C2,tT.C1);	// C4[i-2]c'[i]
						if( !t.Equals(Term.Zero) ) {
							cNext.Add(t);
						}
					}
				}
				if( i >= 2 ) {
					foreach( Term tP in cPrev ) {
						{
							Term t = Term.GetTerm(tP.C4,tP.C2|( 1u<<( i-2 ) ),tP.C1|( 1u<<i ));		// C1[i]C2[i-2]c'[i-1]
							if( !t.Equals(Term.Zero) ) {
								cNext.Add(t);
							}
						}
						{
							Term t = Term.GetTerm(tP.C4,tP.C2|( 1u<<( i-1 ) )|( 1u<<( i-2 ) ),tP.C1);	// C2[i-1]C2[i-2]c'[i-1]
							if( !t.Equals(Term.Zero) ) {
								cNext.Add(t);
							}
						}
					}
				}
				foreach( Term tB in bTemp ) {
					Term t = Term.GetTerm(tB.C4,tB.C2,tB.C1);	// b[i]
					if( !t.Equals(Term.Zero) ) {
						cNext.Add(t);
					}
				}

				// b[i+1]
				foreach( Term tT in cTemp ) {
					Term t = Term.GetTerm(tT.C4|( 1u<<( i-2 ) ),tT.C2|( 1u<<( i-1 ) ),tT.C1|( 1u<<i ));	// C1[i]C2[i-1]C4[i-2]c'[i]
					if( !t.Equals(Term.Zero) ) {
						bNext.Add(t);
					}
				}
				foreach( Term tB in bTemp ) {
					Term t = Term.GetTerm(tB.C4,tB.C2|( 1u<<( i-1 ) ),tB.C1|( 1u<<i ));	// C1[i]C2[i-1]b[i]
					if( !t.Equals(Term.Zero) ) {
						bNext.Add(t);
					}
				}

				Normalize(ref cNext);
				Normalize(ref bNext);
				cPrev = cTemp;
				cTemp = cNext.ToArray();
				bTemp = bNext.ToArray();
				if( i >= 0 ) {
					Term t = Term.GetTerm(0u,1u<<i,0u);	// C2[i]
					if( !t.Equals(Term.Zero) ) {
						cNext.Add(t);
					}
				}
				Normalize(ref cNext);

				Console.WriteLine("######## c[" + ( i+2 ) + "] : (項数" + cNext.Count + ") ########");
				string filePath = string.Format("c{0}.txt",i+2);
				SaveAndOpenFile(filePath,cNext);
				Console.WriteLine();
				Console.WriteLine();

				if( doGC ) {
					GC.Collect();
				}
			}
			return;
		}

		private void Normalize(ref List<Term> exp) {
			if( exp == null || exp.Count == 0 ) {
				exp = new List<Term>(new Term[] { Term.Zero });
				return;
			}

			// exp をソート
			exp.Sort();

			// 不要な term を削除
			if( doNormalizeTerm ) {
				for( int i=0;i<exp.Count; ) {
					Term termTemp = exp[i];
					if( termTemp == null || termTemp.Equals(Term.Zero) ) {
						exp.RemoveAt(i);
						continue;
					} else if( i + 1 >= exp.Count ) {
						break;
					}

					Term termNext = exp[i+1];
					if( termTemp.Equals(termNext) ) {
						exp.RemoveAt(i+1);
						exp.RemoveAt(i);
					} else {
						i++;
					}
				}

				if( exp.Count == 0 ) {
					exp.Add(Term.Zero);
				}
			}

			return;
		}

		private void SaveAndOpenFile(string fileName,ICollection<Term> exp) {
			// ディレクトリ作成
			if( !System.IO.Directory.Exists(dirPath) ) {
				System.IO.Directory.CreateDirectory(dirPath);
			}

			// 保存
			string filePath = dirPath + fileName;
			if( doWriteFile || doDisplay ) {
				try {
					if( doWriteFile ) {
						System.IO.File.WriteAllText(filePath,string.Empty,encode);	// ファイル内容消去
					}

					StringBuilder buf = new StringBuilder("   ");
					bool isFirstTerm = true;
					foreach( Term term in exp ) {
						if( isFirstTerm ) {
							isFirstTerm = false;
						} else {
							try {
								buf.Append(" ^ ");
							} catch( ArgumentOutOfRangeException ) {
								if( doWriteFile ) {
									System.IO.File.AppendAllText(filePath,buf.ToString(),encode);
								}
								if( doDisplay ) {
									Console.Write(buf.ToString());
								}
								buf = new StringBuilder(" ^ ");
							}
						}
						try {
							buf.Append(term.ToString());
						} catch( ArgumentOutOfRangeException ) {
							if( doWriteFile ) {
								System.IO.File.AppendAllText(filePath,buf.ToString(),encode);
							}
							if( doDisplay ) {
								Console.Write(buf.ToString());
							}
							buf = new StringBuilder(term.ToString());
						}
						try {
							buf.AppendLine();
						} catch( ArgumentOutOfRangeException ) {
							if( doWriteFile ) {
								System.IO.File.AppendAllText(filePath,buf.ToString(),encode);
							}
							if( doDisplay ) {
								Console.Write(buf.ToString());
							}
							buf = new StringBuilder(Environment.NewLine);
						}
					}
					if( doWriteFile ) {
						System.IO.File.AppendAllText(filePath,buf.ToString(),encode);
					}
					if( doDisplay ) {
						Console.Write(buf.ToString());
					}
				} catch( Exception ) {
					if( doWriteFile ) {
						Console.WriteLine(string.Format("ファイルの保存に失敗しました。 : {0}",filePath));
					}
					return;
				}
			}
			// 表示
			if( doWriteFile && doOpenFile ) {
				try {
					System.Diagnostics.Process.Start(filePath);
				} catch( Exception ) {
					Console.WriteLine(string.Format("ファイルを開けませんでした。 : {0}",filePath));
					return;
				}
			}

			return;
		}

	}
}
using System;
using System.Collections.Generic;
using System.Text;

namespace Test {
	public class Term : IComparable<Term>,IEquatable<Term> {
		// constructor
		private Term(uint c4,uint c2,uint c1) {
			this.c4 = c4;
			this.c2 = c2;
			this.c1 = c1;
		}
		// construct method
		public static Term GetTerm(uint c4,uint c2,uint c1) {
			if( ( c4 & c2 ) != 0u ) {
				return zero;
			} else if( ( c4 | c2 | c1 ) == 0u ) {
				return one;
			} else {
				return new Term(c4,c2,c1);
			}
		}

		// singleton
		private static readonly Term zero = new Term(~0u,~0u,~0u);
		public static Term Zero { get { return zero; } }
		private static readonly Term one = new Term(0u,0u,0u);
		public static Term One { get { return one; } }

		// field property
		private readonly uint c1;
		public uint C1 { get { return c1; } }
		private readonly uint c2;
		public uint C2 { get { return c2; } }
		private readonly uint c4;
		public uint C4 { get { return c4; } }

		// implement method
		public int CompareTo(Term other) {
			if( (object)other == null ) {
				return -1;
			}

			for( int i=31;i>=0;i-- ) {
				int r1 = this.R(i);
				int r2 = other.R(i);
				if( r1 == 0 && r2 == 0 ) {
					continue;
				} else if( r1 == 0 || r2 == 0 ) {
					return r2 - r1;
				}
			}
			for( int i=31;i>=0;i-- ) {
				int r1 = this.R(i);
				int r2 = other.R(i);
				if( r1 != r2 ) {
					return r1 - r2;
				}
			}
			return 0;
		}
		public bool Equals(Term other) {
			if( (object)other == null ) {
				return false;
			}

			return this.c1 == other.c1
				&& this.c2 == other.c2
				&& this.c4 == other.c4;
		}

		// override method
		public override string ToString() {
			if( this.Equals(zero) ) {
				return "<0>";
			} else if( this.Equals(one) ) {
				return "<1>";
			}

			StringBuilder buf = new StringBuilder();
			bool isFirstFactor = true;
			for( int i=31;i>=0;i-- ) {
				int r = this.R(i);
				if( r > 5 ) {
					return "<0>";
				} else if( r > 0 ) {
					if( isFirstFactor ) {
						isFirstFactor = false;
						buf.Append(string.Format("C{0}[{1}]",r,i+1));
					} else {
						buf.Append(string.Format("&C{0}[{1}]",r,i+1));
					}
				}
			}
			return buf.ToString();
		}

		// method
		/// <summary>
		/// この項に含まれる C_r(Ω_i) の r を求める。
		/// </summary>
		/// <param name="omegaIndex">Ω_i の添字 i (0≦i≦31)</param>
		/// <returns>C_r の添字 r (0≦r≦7)</returns>
		public int R(int omegaIndex) {
			uint mask = 1u;
			uint r = ( ( c1 >> omegaIndex ) & mask )
				 | ( ( ( c2 >> omegaIndex ) & mask ) << 1 )
				 | ( ( ( c4 >> omegaIndex ) & mask ) << 2 );
			return (int)r;
		}

	}
}


ここまで書いておいて残念なお知らせ。
このプログラムを実行しても、 c_{1} から c_{16} までしか求められなかった。
(厳密に言うと c_{17} までは求められたけど、ファイル出力時に System.OutOfMemoryException 吐いて出力できなかった。)
原因は明らかに、項数が爆発的に増加したことによるメモリ不足。
実行してみて分かったのだが、 c_{i+1} の項数は c_{i} の項数のほぼ3倍になるのだ。
c_{16} の項数が(リード・マラー標準形の最簡易形で)1,793,621項だったから、 c_{32} の項数はその 3^{16} 倍の約77兆項になると思われる。
……ああ、そりゃ無理だわ。
ちなみに、 c_{16} を記録したファイルのファイルサイズが133MBだった。
c_{i}i が増加するにつれて1項あたりの(平均)因子数も増えるので、ファイルサイズはさらに爆発的に増大することが予想できるわけで……
もちろん投稿文字数制限を余裕で越えるので、このプログラムで求めた c_{i} の式を貼るのは諦めた。


ちなみに、上のプログラムの途中で項数を出力させて確認したところ、(少なくとも c_{17} までの範囲では) Normalize() メソッドで正規化しても項数が減らないことが分かった。
これはつまり、(*6)の連立漸化式の計算途中の式に同一の項が含まれないことを示している。
言い換えれば、 Normalize() メソッドを使う必要がないということになる(表示用に項のソートはする必要があるが)。
(★1)をそのまま実装したプログラムの場合は Normalize() メソッドで正規化すると項数が(√項数)くらいまで減るので、(★1)は無駄な項の計算にメモリと時間を割いていることになる。
メモリ効率・時間効率ともに(★1)より(*6)の漸化式のほうが優秀ということだ。
あとはもう少しコードを最適化できれば、 c_{18} くらいまでなら求められるかも……?
もっとも、どんな方法を使っても推定約77兆項もある c_{32} まで求めるのは無理っぽいけど。

結局のところ、 c_{i} の一般項を求めないとダメっぽいことか。
……残念だけど、今回はそんな結論で〆ておくことにする。