どうすればうまく解けるんでしょうね。
D = Pk - Pj = Pi
として、
k(3k - 1) / 2 - j(3j - 1) / 2 = i(3i - 1) / 2
(36k2 - 12k) - (36j2 - 12j) = (36i2 - 12j)
(6k - 1)2 - (6j - 1)2 = (6i - 1)2 - 1
ここで、
a = 6i - 1
b = 6j - 1
c = 6k - 1
と置くと、
c2 - b2 = a2 - 1
(c - b)(c + b) = (a - 1)(a + 1)
積の形に持ち込んだら素因数分解です。(a - 1)(a + 1)を素因数分解して、それを2つの数に分解してそれぞれをc - bとc + bにあてはめます。ここで、(a - 1)(a + 1)の素因数分解を求めるとき、a - 1とa + 1の素因数分解を求めたほうがエラトステネスのふるい的アルゴリズムが使えて速いです。
c - b ≡ 0(mod 6)
c + b ≡ 4(mod 6)
なので、(a - 1)(a + 1)の素因数分解のうち、3の因子は全てc - bに割り当てます。それから2は両方に一つずつ割り当てる必要があります。例えば、a = 11のとき、
a - 1 = 2・5
a + 1 = 22・3
(a - 1)(a + 1) = 23・3・5
3と2を2つ割り当てると、
c - b ← 2・3
c + b ← 2
(a - 1)(a + 1) ← 2・5
なので、残りは(1, 2・5), (2, 5), (5, 2), (2・5, 1)という分配が考えられます。それぞれについて(c - b, c + b)は、
(2・3, 22・5), (22・3, 2・5), (2・3・5, 22), (22・3・5, 2)
このうち、c - b < c + bでc + b ≡ 4(mod 6)を満たすものを選び、Pk + Pjが五角数ならその最初のものが答えです。
import scala.collection.mutable.ArrayBuffer def int_sqrt(n :Long) :Long = { def f(m :Long) :Long = { val s = (m + n / m) / 2 if(s >= m) m else f(s) } f(n) } def is_pentagonal(n :Long) = { val m = int_sqrt(1 + 24 * n) m % 3 == 2 && m * m == 1 + 24 * n } type Facts = List[(Int,Int)] object Primes { private val primes = new ArrayBuffer[Int] private val L = 1e5.toInt private var offs = 0 def apply(n :Int) = { if(n >= primes.size) sieve primes(n) } private def sieve() = { if(offs == 0) sieve_1st else sieve_2nd_or_later offs += L } private def sieve_1st() = { val a = new Array[Boolean](L) for(i <- 2 to (L - 1)) a(i) = true for(p <- 2 to (L - 1) if a(p); m <- (p * 2) to (L - 1) by p) a(m) = false primes ++= (2 to (L - 1)).filter(a(_)).toArray } private def sieve_2nd_or_later() = { val N = offs + L - 1 val a = new Array[Boolean](L) for(i <- 0 to (L - 1)) a(i) = true for(p <- primes.takeWhile(p => p * p <= N); m <- ((offs + p - 1) / p * p) to N by p) a(m - offs) = false primes ++= (offs to N).filter(k => a(k - offs)).toArray } } class genFactors extends Iterator[Facts] { private val L = 1e5.toInt private val b = new Array[Facts](L) private var offs = 1 private var k = 0 def hasNext() = true def next() = { k += 1 if(k >= offs) sieve b(k - offs + L) } private def sieve() = { val N = offs + L - 1 val a = (offs to N).toArray for(i <- 0 to (L - 1)) b(i) = Nil for(p <- (0 to N).map(Primes(_)).takeWhile(p => p * p <= N); m <- ((offs + p - 1) / p * p) to N by p) { val (q, e) = div_pow(a(m - offs), p) a(m - offs) = q b(m - offs) = b(m - offs) ++ List((p, e)) } for(i <- 0 to (L - 1) if a(i) > 1) b(i) = b(i) ++ List((a(i), 1)) offs += L } private def div_pow(n :Int, p :Int, e :Int = 0) :(Int,Int) = if(n % p != 0) (n, e) else div_pow(n / p, p, e + 1) } def pow(n :Int, e :Int) = (1 to e).foldLeft(1)((x, y) => x * n) def num_divs(f :Facts) = f.foldLeft(1)((x, y) => x * (y._2 + 1)) def divs(fs :Facts) :Iterator[Facts] = fs match { case Nil => Iterator(Nil) case (p, e) :: fst => for(fs2 <- divs(fst); e1 <- Iterator.range(0, e + 1)) yield if(e1 == 0) fs2 else (p, e1) :: fs2 } def value(f :Facts) = f.foldLeft(1L)((x, y) => x * pow(y._1, y._2)) def mul(f1 :Facts, f2 :Facts) :Facts = (f1, f2) match { case (f1, Nil) => f1 case (Nil, f2) => f2 case ((p1, e1) :: t1, (p2, e2) :: t2) if p1 < p2 => (p1, e1) :: mul(t1, f2) case ((p1, e1) :: t1, (p2, e2) :: t2) if p1 > p2 => (p2, e2) :: mul(f1, t2) case ((p1, e1) :: t1, (p2, e2) :: t2) => (p1, e1 + e2) :: mul(t1, t2) } def div(f1 :Facts, f2 :Facts) :Facts = (f1, f2) match { case (f1, Nil) => f1 case (Nil, f2) => Nil case ((p1, e1) :: t1, (p2, e2) :: t2) if p1 < p2 => (p1, e1) :: div(t1, f2) case ((p1, e1) :: t1, (p2, e2) :: t2) if p1 > p2 => div(f1, t2) case ((p1, e1) :: t1, (p2, e2) :: t2) if e1 == e2 => div(t1, t2) case ((p1, e1) :: t1, (p2, e2) :: t2) => (p1, e1 - e2) :: div(t1, t2) } def P(n :Long) = n * (3 * n - 1) / 2 def gen_solutions() = { // (6k - 1)^2 - (6j - 1)^2 = (6i - 1)^2 - 1 // c^2 - b^2 = a^2 - 1 // (c - b)(c + b) = (a - 1)(a + 1) val it_fac = new genFactors def next(fs :Facts) = { val a = it_fac.take(6).toArray mul(a(3), a(5)) } // remove all of 3 and 2 of 2 def remove3(fs :Facts) :Facts = fs match { case (2, 2) :: fst => remove3(fst) case (2, e) :: fst => (2, e - 2) :: remove3(fst) case (3, e) :: fst => fst case _ => fs } // number of 3 def number3(fs :Facts) :Int = fs match { case (3, e) :: fst => e case f :: fst => number3(fst) case _ => 0 } // 2^3 3, 2^3 3 5, ... val it_a2 = Iterator.iterate(Nil :Facts)(next).drop(1) // c - b = 6m c + b = 6n + 4 def gen_PPP(fs :Facts) = { val e3 = number3(fs) val fs2 = remove3(fs) for(d1 <- divs(fs2); d2 = div(fs2, d1); d2v = value(d2) if d2v * 2 % 6 == 4; d1v = value(d1); p = d1v * 2 * pow(3, e3); // c - b q = d2v * 2 if p < q; // c + b (c, b) = ((q + p) / 2, (q - p) / 2); (k, j) = ((c + 1) / 6, (b + 1) / 6); i = (int_sqrt(c * c - b * b + 1) + 1) / 6) yield (P(i), P(j), P(k)) } for(fs <- it_a2; (pi, pj, pk) <- gen_PPP(fs); pl = pk + pj if is_pentagonal(pl)) yield pi } println (gen_solutions().next)
