NPCA Advent Calendar 2015 12/24


今日はクリスマス・イヴですね。
だからといって特にクリスマス関係のことはありません。
他のプロ達と違って私は
実力皆無人間なので

文が稚拙なうえに内容も薄いものとなるでしょう。

eに探す素数



問題:
r進表記のeの連続する桁で、
最初に出てくるn桁の素数

(Σ{k=0}ƒ(k)とはkの値を0より始めたƒ(k)の無限和のこととする)
私の解法の主な点としては、
・整数をm進表記の各桁の配列で表す
・eをテイラー展開した式を、下記の条件を満たすƒ(k)におけるΣƒ(k)となるように変形
kが偶数ならƒ(k) > 0, kが奇数ならƒ(k) < 0, 絶対収束する

が挙げられる

1 整数を配列で

整数を result[k]*(m^k)の総和 (0 <= result[k] < m)の形で表し、resultを求める配列とする。
但し、任意の配列とそれに0を0個以上pushして得られる配列は同値とする。
通常の記法とは逆に桁が並んでいることに注意。
例:
10進法で495 ⇔ [5, 9, 4]
[0, 0, 5] ≠ [0, 0, 0, 5] (∵ 500 ≠ 5000)
[0, 0, 5] = [0, 0, 5, 0] (∵ 500 = 0500)

1.1 大小関係

先ず、末尾に0をいくらか付け足し、2数の桁数を一致させる。
その上で、末尾から各桁の大小関係を確認する。

1.2 四則演算

筆算の方式で計算。
(2数a,bの積を求める時、bの任意の桁kにおいてdp[k] = a * kとすれば、
a * kを求める回数を各々のkにおいて1回以内に抑えられる。)

1.3 インクリメント

配列の末尾から(r-1)(rは基数)と等しくないものを探し、
そこまでの(r-1)を0に変えた上で探索結果をインクリメント。
すべて(r-1)であった場合は全て0にして末尾に1を付加。

1.4 デクリメント

1.3と同様に、配列の末尾から0と等しくないものを探し、
そこまでの0を(r-1)に変えた上で探索結果をデクリメント。
但し、全て0であった場合は即ち0なので、NaNを返す。

2 eを式に表す

e = Σ{k=0} 1 / k!
これを変形し
e = Σ{k=0}ƒ(k)
但し、非負整数xにおいて
ƒ(2x) = 1 / 2x-1
ƒ(2x+1) = 1/x! – ƒ(2x)

また この式より
e = 4 – Σ{k=0} (1 / 2k-1 – 1/k!)

3 効率化

3.1

私は1で”末尾に0を付加しても同値”といいました。
しかし仕様上、そのままでは効率が悪くなるので, 適宜短くします。

3.2

私の実装では配列によって疑似的な分数を表現しています。
そのため、適切なタイミングで約分しなければなりません。
そこで、分子と分母の最大公約数を求める必要があります。
そのアルゴリズムとして、有名なユークリッドの互換法を使っています。

(桁が求まったらその分eから減じるのも手ですが、実装していません)

実装例

​
function solve(digits, radix){
	radix = parseInt(radix);
	// Invalid radix
	if(!radix){
		// default radix
		radix = 10;
	}
	digits = parseInt(digits);
	// Invalid
	if(isNaN(digits)){
		return NaN;
	}
	var forceNum = function (n){
		// Force be a number.
		if(!n)
			return 0;
		else
			return Number(n);
	};
	var bigInt = {
		wash: function (n){
			while(n.length){
				var i = n.pop();
				if(i){
					n.push(i);
					return n;
				}
			}
			n.push(0);
			return n;
		},
		incf: function (result){
			var r = result.slice();
			var s = [];
			while(r.length){
				var i = r.shift();
				if(i !== radix - 1){
					return s.concat(i + 1, r);
				}else{
					s.push(0);
				}
			}
			return s.concat(1);
		},
		decf: function (result){
			var r = result.slice();
			var s = [];
			while(r.length){
				var i = r.shift();
				if(i){
					return s.concat(i - 1, r);
				}else{
					s.push(radix - 1);
				}
			}
			return NaN;
		},
		pure: function (radix, big){
			if(!big.length){
				return 0;
			}else{
				return big[0] + radix * big.slice(1);
			}
		},
		conversion: function (radix, molecule, denomi){
			var single = function (m){
				var result = [];
				while(m){
					result.push(m % radix);
					m = Math.floor(m / radix);
				}
				if(!result.length){
					result.push(0);
				}
				return result;
			};
			var result = [];
			if(denomi){
				// denomi isn't 0, undefined, etc...
				return [single(molecule), single(denomi)];
			}else{
				// denomi is 0, undefined, etc...
				return single(molecule);
			}
		},
		comparison: function (a,b){
			var A = a.slice();
			var B = b.slice();
			var result = 0;
			single = function (x,y){
				if(x>y)
					return +1;
				else if(x<y)
					return -1;
				else
					return 0;
			};
			while(A.length || B.length){
				var i = single(
					forceNum(A.shift()),
					forceNum(B.shift())
				);
				if(i)
					result = i;
			}
			return result;
		},
		add: function (a,b){
			// a + b
			var result = [];
			var A = a.slice();
			var B = b.slice();
			var increase = 0;
			while(A.length || B.length || increase){
				var tiny = forceNum(A.shift()) + forceNum(B.shift()) + increase;
				increase = 0;
				while(tiny >= radix){
					tiny -= radix;
					increase++;
				}
				result.push(tiny);
			}
			return result;
		},
		sub: function (a,b){
			// a - b (a >= b)
			// Use comparison(a,b) before use me
			var result = [];
			var A = a.slice();
			var B = b.slice();
			var decrease = 0;
			while(A.length || B.length || decrease){
				var tiny = forceNum(A.shift());
				var parameter = forceNum(B.shift()) + decrease;
				decrease = 0;
				while(tiny < parameter){
					decrease++;
					tiny += radix;
				}
				result.push(tiny - parameter);
			}
			return result;
		},
		multi: function (a,b,m){
			// a * b
			var result = [0];
			var A = a.slice();
			var B = b.slice();
			var height = [];
			var map = (m ? m.slice() : []);
			var single = function (a,b){
				var memo = map[b];
				if(memo){
					// defined
					return memo;
				}
				var result = [];
				var A = a.slice();
				var increase = 0;
				var simple = function (a,b,i){
					var r = 0;
					var s = 0;
					while(b--){
						r += a;
						if(r >= radix){
							r -= radix;
							s++;
						}
					}
					r += i;
					if(r >= radix){
						r -= radix;
						s++;
					}
					return [r,s];
				};
				while(A.length || increase){
					var i = simple(forceNum(A.shift()), b, increase);
					result.push(i[0]);
					increase = i[1];
				}
				return map[b] = result;
			};
			while(B.length){
				result = bigInt.add(
					result, height.concat(single(A, B.shift()))
				);
				height.push(0);
			}
			return [result, map];
		},
		div: function (a,b){
			// a / b, a % b
			var A = a.slice();
			var B = b.slice();
			var result = [];
			var height = [0];
			while(bigInt.comparison(A,B) >= 0){
				height.unshift(0);
				B.unshift(0);
			}
			height.shift();
			B = height.concat(b);
			while(height.length){
				var i = bigInt.conversion(radix, 0);
				var lower, gauss = bigInt.conversion(radix, 0);
				B.shift();
				for(;;){
					lower = bigInt.multi(B,i)[0];
					if(bigInt.comparison(A, lower) < 0){
						break;
					}
					gauss = lower;
					i = bigInt.incf(i);
				}
				i = bigInt.decf(i);
				A = bigInt.sub(A, gauss);
				result.unshift(i[0]);
				height.shift();
			}
			return [result,A];
		},
		gcd: function (a,b){
			if(bigInt.conversion(a,b) < 0){
				return arguments.callee(b,a);
			}else if(!b.some(function (e){return e})){
				return a;
			}else{
				return arguments.callee(b, bigInt.wash(bigInt.div(a,b)[1]));
			}
		},
		isPrime: function (n){
			// simple test
			var i = bigInt.conversion(radix, 2);
			while(bigInt.comparison(bigInt.multi(i,i)[0], n) < 0){
				var divided = bigInt.div(n,i)[1];
				if(!bigInt.comparison(divided, bigInt.conversion(radix,0))){
					return false;
				}
				i = bigInt.incf(i);
			}
			return true;
		}
	};
	var fraction = {
		floor: function (q,h){
			var molecule = q[0].slice();
			var height = h.slice();
			while(height.some(function (e){return e})){
				molecule.unshift(0);
				height = bigInt.decf(height);
			}
			return bigInt.div(molecule, q[1])[0];
		},
		add: function (p,q){
			// a/b + c/d = (ad+bc)/bd
			var i = bigInt.multi(q[1], p[0]);
			var j = bigInt.multi(p[1], q[0]);
			var m = bigInt.add(i[0], j[0]);
			var d;
			if(i[1].length > j[1].length){
				d = bigInt.multi(q[1], p[1], i[1]);
			}else{
				d = bigInt.multi(p[1], q[1] ,j[1]);
			}
			return [m,d[0]];
		},
		sub: function (p,q){
			// a - b (a >= b)
			// a/b - c/d = (ad - bc)/bd
			var i = bigInt.multi(q[1], p[0]);
			var j = bigInt.multi(p[1], q[0]);
			var m = bigInt.sub(i[0], j[0]);
			var d;
			if(i[1].length > j[1].length){
				d = bigInt.multi(q[1], p[1], i[1]);
			}else{
				d = bigInt.multi(p[1], q[1] ,j[1]);
			}
			return [m,d[0]];
		},
		wash: function (q){
			var divisor = bigInt.gcd(q[0], q[1]);
			return [
				bigInt.wash(bigInt.div(q[0], divisor)[0]),
				bigInt.wash(bigInt.div(q[1], divisor)[0])
			];
		}
	};
	var factiorial = {
		results: (function (){
			var x = {};
			x[bigInt.conversion(radix,0)] = bigInt.conversion(radix,1);
			return x;
		})(),
		solve: function (k){
			bigInt.wash(k);
			var n = factiorial.results[k];
			console.log('fact', k);
			if(n){
				// defined
				return n;
			}else{
				var lower = arguments.callee(bigInt.decf(k));
				return factiorial.results[k] = bigInt.multi(lower, k)[0];
			}
		}
	};
	var power = {
		results: (function (){
			var x = {};
			x[bigInt.conversion(radix,0)] = bigInt.conversion(radix,1);
			return x;
		})(),
		solve: function (k){
			bigInt.wash(k);
			var n = power.results[k];
			console.log('pow', k);
			if(n){
				//defined
				return n;
			}else{
				var two = bigInt.conversion(radix, 2);
				var lower = arguments.callee(bigInt.decf(k));
				return power.results[k] = bigInt.multi(lower, two)[0];
			}
		}
	};
	var terms = {
		min_sum: bigInt.conversion(radix, 0,1), // the first term = 0
		max_sum: bigInt.conversion(radix, 4,1), // the first term = 4
		min_general: function (k){
			// 1/(k!)
			return [bigInt.conversion(radix, 1), factiorial.solve(k)];
		},
		max_general: function (k){
			// 1/(2^(k-1)) - 1/(k!)
			var i = [bigInt.conversion(radix, 2), power.solve(k)];
			var j = [bigInt.conversion(radix, 1), factiorial.solve(k)];
			return fraction.sub(i,j);
		}
	};
	function main(){
		var height = bigInt.conversion(radix,0);
		var e = [];
		var i = bigInt.conversion(radix,0);
		for(;;){
			console.log('min start');
			terms.min_sum = fraction.add(terms.min_sum, terms.min_general(i));
			console.log(terms.min_sum);
			console.log('min end');
			console.log('max start');
			terms.max_sum = fraction.sub(terms.max_sum, terms.max_general(i));
			console.log('max end');
			console.log(i);
			var tiny = fraction.floor(terms.min_sum, height);
			var ceil = fraction.floor(terms.max_sum, height);
			console.log(tiny, ceil);
			console.log(height);
			if(!bigInt.comparison(tiny, ceil)){
				terms.min_sum = fraction.wash(terms.min_sum);
				terms.max_sum = fraction.wash(terms.max_sum);
				e.unshift(tiny[0]);
				height = bigInt.incf(height);
				var E = e.slice(0,digits);
				console.log(e,E);
				if(E.length >= digits && bigInt.isPrime(E)){
					return E.reverse();
				}
			}
			i = bigInt.incf(i);
		}
	}
	return main();
}
​
var begin = new Date();
console.log(solve(10,10));
var end = new Date();
​
console.log(end - begin); // ms
​

r = 10, n = 10とし、手元の環境(Chrome 47.0.2526.106)で試した所、
8188184ms(約2時間15分)かかりました。
原因は素数判定の関数が遅いからと思われます。

素数がメインだったはずだが。

さいごに

こんな駄文におつきあい頂きまして、ありがとうございました。
取り掛かるのが遅かったせいで実装が雑で文も提出前日に書き始める始末。
こんなことにならないように、何事も余裕をもって早めに始めましょう(戒め)

GitHubは挫折しました。

さて、明日は誰が何と言おうとクリスマスですが、
日付が日付なので、こんな記事よりもよほど
秀逸な記事が現れるでしょう。

明日は最終日、担当はhideo54さんです。ご期待ください。




コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です