ScalaでProject Euler(74)

Problem 44

どうすればうまく解けるんでしょうね。

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 - bc + 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 + bc + 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)