ルジャンドルの公式をJavaで実装する


よく見ているヨビノリさんの動画で、ハーバード&MIT数学トーナメントのこんな問題が紹介されていました。「n!の末尾に290個の0が並ぶ最小の自然数nを求めよ」。

https://www.youtube.com/watch?v=W3WruHbqjd8&t=66s

解法は動画の素晴らしい解説を見ていただくとして、今回はJavaで問題を解いてみます。類似の問題が毎年出題されているようですので、ある数（int）を与えると、その階乗（BigInteger）とその全ての素因数と冪指数の配列（int[][]）を保持できるクラスを作成し、より汎用的に使えるようにしてみましょう。数学クイズ以外の用途が思い浮かばないけど。

```java:LugendresFactorizer.java
package FactorialSolver;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.lang.Math;
import java.math.BigInteger;

public class LugendresFactorizer {
	int n;
	BigInteger factorial = BigInteger.ONE;	
	int[][] factorsAndExps;
	LugendresFactorizer(int inputnum){
		n = inputnum;
		if(n < 1) throw new ArithmeticException("Input value must be greater than 0.");

		for(int i=1; i<=n; i++) {
			factorial = factorial.multiply(BigInteger.valueOf(i));		
		}

		int factors[] = eratosthenesSeive(inputnum);
		int primesCount = factors.length;
		factorsAndExps = new int[primesCount][2];
		for(int i=0; i<primesCount; i++) {
			factorsAndExps[i][0]=factors[i];
			factorsAndExps[i][1] = 1;
		}	
		
		if(primesCount>1) {
			int powerCount;	
			for(int p=0; p<primesCount; p++) {
				powerCount=0;
				for(int i=1; Math.pow(factorsAndExps[p][0], i)<=n; i++){
					powerCount += (int)n/(Math.pow(factorsAndExps[p][0], i));
				}
				factorsAndExps[p][1]=powerCount;
			}
		}
	}
	
	private int[] eratosthenesSeive(int x) {
		List<Integer> primes = new ArrayList<Integer>();
	    boolean isPrime[] = new boolean[x+1];
	    Arrays.fill(isPrime, true);
	    isPrime[0] = false;
	    isPrime[1] = false;
	    
	    int sqrtX = (int) Math.sqrt(x);

	    for(int q=2; q<=sqrtX; q++) {
	        if (isPrime[q]) {
	            for (int powerQ =q*q; powerQ<isPrime.length; powerQ += q) {
	            	isPrime[powerQ] = false;
	            }
	        }
	    }
	    for (int i=2; i<isPrime.length; i++) {
	        if (isPrime[i]) {
	        	primes.add(i);
	        }
	    }
	    return primes.stream().mapToInt(i->i).toArray();
	}
	
	public String toString() {
		String s = n+"! = "+factorial;
		if(n>=3) {
			s+=" = ";
			for (int i = 0; i< factorsAndExps.length; i++) {
				if(factorsAndExps[i][1]==1) s += factorsAndExps[i][0];
				else s += factorsAndExps[i][0]+"^"+factorsAndExps[i][1];
				if(i != factorsAndExps.length-1) s += " * ";
			}
		}
		return s;
	}
}

```

入力値はintで良いとして、階乗は桁数が大きくなるので、BigIntergerクラスを用います。作るべき配列は、素因数の個数をlengthとして中に素因数とその冪指数が入る2次元配列で良いでしょう。

変なオブジェクトができないように入力値は1以上でなければ例外を返すようにします。なお、入力値1の場合、求めるべき素数の配列は2から始まるので、中身がnullになってしまいます。クラス設計のあるべき形としては変数がnullになるケースを避けるべきかもしれませんが、1!=1と言えるべきだし、このクラスは「最小の自然数」を求める問題のためのものなので1から使えるのが自然なのでは・・・と思い、そのままにしてしまいました。

コンストラクタにnが渡ると、まず、BigIntergerでn!を計算します。次に、nまでの素数をエラストテネスの篩を用いてリストアップしています。参照したWikipediaの記事を貼っておきます（なるべく変数名を記事の表記に近づけました）。なお、篩メソッドはクラス利用者側から見える必要が無いのでprivateです。

https://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%A9%E3%83%88%E3%82%B9%E3%83%86%E3%83%8D%E3%82%B9%E3%81%AE%E7%AF%A9

その後が今回のポイントになります。セットした素数に対して、ルジャンドルの公式を用います。任意の素数 p と正の整数 n があるとき、n!をpで何回割れるかは、nをpのi乗で割ったものの小数点以下切り捨て値（floor）の総計で求められます。
```math
v_p(n!)= \sum_{i=1}^{\infty} \biggl\lfloor\frac{n}{p^i}\biggl\rfloor
```
SUMが無限までありますが、実際には、pのi乗がnを上回った時点でfloor値は0なので、加算する必要がなくなります。プログラムとしてはfor文の条件部にそう書けばよいですね。今回nは必ず正のため、ガウス記号(floor)はintへのキャストのみで済ませています。これで、素因数とその冪指数の配列ができました。ルジャンドルの公式のWikipediaはこちら。

https://ja.wikipedia.org/wiki/%E3%83%AB%E3%82%B8%E3%83%A3%E3%83%B3%E3%83%89%E3%83%AB%E3%81%AE%E5%85%AC%E5%BC%8F

最後にtoString()をオーバーライドします。nに6を与えると、"6! = 720 = 2^4 * 3^2 * 5"という文字が読めるようにしておきます。nに1や2が入っても良いように、微妙に調整しました。

サンプルとして以下を実行してみます。

```java:Query.java
package FactorialSolver;
public class Query {
	public static void main(String[] args) {
		LugendresFactorizer sample = new LugendresFactorizer(6);
		System.out.println(sample.factorial);
		System.out.println(java.util.Arrays.deepToString(sample.factorsAndExps));
		System.out.println(sample);	
	}
}
```
結果は、以下の通り。
```
720
[[2, 4], [3, 2], [5, 1]]
6! = 720 = 2^4 * 3^2 * 5
```

では、冒頭の問題を解いてみましょう。まずはブルートフォース的に、1から順番にn!のBigIntergerの文字列末尾の0をカウントして、290個になるところで止めてみることにします。

```java:Query.java
package FactorialSolver;
public class Query {
	public static void main(String[] args) {
		int answer=0;
		for(int n=1;n<1200;n++) 
		{
			LugendresFactorizer l = new LugendresFactorizer(n);
			String str = l.factorial.toString();
			int zeroCount = 0;
			for(int i=1; i<str.length(); i++) {
				if(str.charAt(str.length()-(i))=='0') zeroCount++;
				else break;
			}
			if (zeroCount >= 290) {
				answer = n;
				break;
			}
		}
		System.out.println(answer);
	}
}
```
結果は、1170と分かりました。
```
1170
```
ちなみに1170の階乗は、以下らしいです。これを人間が計算するのは無理がありますね。
```
385592176713081905777902060537842937133444270763926883602505708504686828523855227736782513020240642477050060973716744485519220353584669286420027170808127302720699016954608238473674353770528409106749447927476477236411677646111332144137405703434344736931681418145778630072594748329716812084825177801594970266355143734157323691647373999306397991646320995132137373276929201067604618130494049409596259237229702387257422558840859528287825857910523170560588303635973317591255204546775135073285285878999917534483481914016611909634668946741353475159673090194800545290093801725368807301380867983212039014214265075015671756952790958574797599942904577418150219819188925707102238544747990821435098160386511692360160773478510741999591766968391163564106673059305143121362265015217415144282507739825314598458074771983211360774602775111642814771857044418793841119742035165824963835337868871088872381096259147487782688562909363841409964985060069231520990150218677187109335823954632872809344266246310724508740449261748597796411072248851088569241458572538194573316061120827086334791278876552180284013017811860384060797887094806367728668640229687480911270206188005954330023812714920892436721103616546630819076292692146948752666284442909800108685033454308083146705070678246512925168663800199686363976911540113404857996860080323014624931280392542362037914323938677568681743249721245815347866920743972468249027464230458639103789870896642945718045234361503690556177300957112100243651013552010598029223494728811615297077428597916371556113562640258973796712710956672106567303726858670156709778194593026920479341721876797300033531498046532317148488101819465339808819179802510491249935672952266513303607769862335587055757815071454115418799114339624558802961363266719780280640362520223947469571164185509943905318196725856928723203277687694649204992578641481148170808669974812844998900444515793233681341246883473172416321028987562625596545204666725057123072278950141098695714338879120868748166402995648937988676334898014515564994985745325407136426128517200366978204727138169309280193720398541518728985961305012276565415046078178323699570396804554663420103819595174809218301517464895569356696025580107615313075598468027313550716208814056419372495326703420023270078888865811090493524473060270050928442312899585585934941788822944497683818951238337611559226616340301545265720696567139051189275787730631377952371062791230226046999791565472869605081294204657336604974763107527377234870331566796421650534000892115906187724861971114775685373821195497869022738976176984608724133494327257024213321245889837196091331812225897490451092661718858474314925751977853723706885823318095811113816515263910506753838854997969254950504859692633315542537689338474630801472992662724878978439933475360492719262502466282099789861435709941557504665190400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
```

次に、人間にとっての現実的な解法である、ルジャンドルの公式を用いた方法でやってみます。前提として、末尾0の数は10(=2×5)で割れる数に等しく、2のi乗と5のj乗がある時、nを階乗した結果の今回の配列中では必ずi>jです。つまり、n!を5で何回割ることができるかが、n!の末尾の0の数に等しくなることになります。今回のクラスではfactorsAndExps配列に素数が格納され、その順序は2,3,5ですので、求める数はfactorsAndExps[2][1]ということになります。これが290以上になるnを計算すればよいわけです。配列範囲外を避けるためnは5以上を指定しています。

```java:Query.java
package FactorialSolver;
public class Query {
	public static void main(String[] args) {
		int anotherAnswer=0;
		for(int n=5;n<1200;n++) 
		{
			LugendresFactorizer l = new LugendresFactorizer(n);
			if (l.factorsAndExps[2][1] >= 290) {
				anotherAnswer = n;
				break;
			}
		}
		System.out.println(anotherAnswer);
	}
}
```
無事、等しい結果が出てきました。
```
1170
```
当記事のソース全体は以下に置いてあります。
https://github.com/KentAnak/FactorialSolver

また、元々の問題文はこちらですので、リファレンスを書いておきます。
"Find the smallest n such that n! ends in 290 zeros"(HMMT, 2003).
<b>Reference:</b>
HMMT.(2003) <i>Harvard-MIT Mathematics Tournament| March 15, 2003</i>.
https://hmmt-archive.s3.amazonaws.com/tournaments/2003/feb/alg/problems.pdf