ScalaでProject Euler（106）

Scala

Problem 70

この問題は、

Aを満たす数でBが最小のものを求めよ

というパターンなので、

Bの昇順に並べてAを満たす最初の数を求めよ

と読み替えればよいです。すなわち、n / φ(n)（以下、phi比と呼びます）の小さい順に並べればよいです。幸いこれはそんなに難しくないです。

簡単のために100より小さい自然数で考えましょう。

phi比が最小のnは97ですね。しかし、φ(97) = 96で97の並び替えにはなっていません。一般に素数pに対してφ(p) = p - 1は下1桁が違うだけなので並び替えには成り得ません。

素数以外でphi比が最小なのは72 = 49です。phi比は7/6になります。その次に小さい候補は、73, 7 * 13, 52です。このうち73 = 243は100以上なので捨てられます。
phi比が自分以上の整数を選ぶルールを決めて、それに従うと漏れなく、重複なく整数が列挙できる、そんなルールが必要です。このルールで列挙してPriorityQueueに入れていけばphi比の昇順に並べることができます。今回使ったルールは、素因数分解して、

  1. 最大の素数の指数を1増す ex) 72 -> 73
  2. 最大の素数の指数が2以上ならその素数の一つを可能な最大の素数に変える ex) 72 -> 7 * 13
  3. 最大の素数の指数が1なら次に大きい素数に変えることができる。ただし、最大の素数の次に大きい素数より大きくなければならない ex) 7 * 13 -> 7 * 11
  4. 素数の平方なら次に大きい素数べきに ex) 72 -> 52

これで、phi比で並べられました。
かかった時間は、

limit time
10^7 0.7s
10^8 0.04s
10^9 9s
10^10 28s
10^11 169s

意外と遅いんですよね。実はPythonより遅いです。きっと何かが間違っています。

（追記）
間違っていたのは2か所です。1か所はmax_primeのところでRangeを使っているところですね。Iteratorを使えばまっとうな速度になります。Rangeを安易に使ってはいけませんね。もう1か所は、expandで毎回sieveしているところです。凡ミスですね。
これで速くなりました。

limit time
10^7 0.03s
10^8 0.03s
10^9 0.04s
10^10 0.04s
10^11 0.14s
10^12 0.4s
10^13 2s
10^14 4s
10^15 13s
10^16 43s 
import scala.math.sqrt
import scala.collection.mutable.PriorityQueue
import scala.testing.Benchmark

object Primes {
    private var L = 0
    private var N = 0
    private var bits = new Array[Int](0)
    
    def resize(n :Int) = {
        L = n / 96 + 1
        N = L * 96
        sieve()
    }
    
    def enumerate() :Iterator[Int] =
        Iterator(2, 3) ++ enumerate2()
    
    // maximal p < n
    def max_prime(n :Int) :Int =
        if(n <= 5)
            if(n == 5) 3 else if(n >= 3) n - 1 else 0
        else {
            val k0 = if(n % 6 <= 1) n / 3 - 1 else n / 6 * 2
            expand(k0 >> 5)
            val it = for(k <- Iterator.range(k0, -1, -1);
                        m = k >> 5; l = k & 31;
                        if ((bits(m) >> l) & 1) == 1)
                            yield m * 96 + l * 3 + 1 + (l & 1)
            it.next
        }
    
    private def sieve() = {
        bits = new Array[Int](L)
        for(k <- 0 to L - 1) bits(k) = -1
        bits(0) &= ~1
        for(p <- enumerate2().takeWhile(p => p * p < N)) {
            val k0 = p / 3
            for(k <- k0 + p * 2 to N / 3 - 1 by p * 2)
                erase_bit(k)
            for(k <- -k0 - 1 + p * 2 to N / 3 - 1 by p * 2)
                erase_bit(k)
        }
    }
    
    private def enumerate2() :Iterator[Int] =
        for((b, m) <- bits.toIterator.zipWithIndex;
            k <- Iterator.range(0, 32)
            if ((b >> k) & 1) == 1)
                yield m * 96 + k * 3 + 1 + (k & 1)
    
    private def erase_bit(k :Int) = {
        val m = k >> 5
        val l = k & 31
        bits(m) &= ~(if(l != 31) 1 << l else -1 << 31)
    }
    
    private def expand(n :Int) :Unit = {
        if(n >= L) {
            N <<= 1
            L <<= 1
            sieve()
            expand(n)
        }
    }
}

def pow(n :Int, e :Int) :Long =
    if(e == 0) 1L else n * pow(n, e - 1)

type Facts = List[(Int,Int)]

def value(fs :Facts) :Long =
    fs.foldLeft(1L)((x, y) => x * pow(y._1, y._2))

def phi(fs :Facts) :Long =
    fs.foldLeft(1L)((x, y) => x * pow(y._1, y._2 - 1) * (y._1 - 1))

def phi_ratio(fs :Facts) :(Long,Long) = {
    val a = fs.foldLeft(1L)((x, y) => x * y._1)
    val b = fs.foldLeft(1L)((x, y) => x * (y._1 - 1))
    (a, b)
}

class numbers_sorted_by_phi_ratio(N :Long) extends Iterator[Facts] {
    private val o = new Ordering[((Long,Long),Facts)]() {
        def compare(x :((Long,Long),Facts), y :((Long,Long),Facts)) =
            BigInt(x._1._2) * y._1._1 compare BigInt(x._1._1) * y._1._2
    }
    private val pq = new PriorityQueue[((Long,Long),Facts)]()(o)
    init()
    
    def hasNext() = pq.nonEmpty
    
    def next() = {
        val (_, fs) = pq.dequeue
        val (p, e) = fs.head
        
        // [ (5, 2) ] -> [ (5, 3) ]
        if(value(fs) < N / p)
            put((p, e + 1) :: fs.tail)
        
        // [ (5, 2) ] -> [ (19, 1), (5, 1) ]
        if(e >= 2) {
            val n = (N - 1) / value((p, e - 1) :: fs.tail) + 1
            val p_max = Primes.max_prime(n.toInt)
            if(p_max > p) fs match {
                case (p, e) :: tail => put((p_max, 1) :: (p, e - 1) :: tail)
                case Nil => ()
            }
        }
        // [ (19, 1), (5, 1) ] -> [ (17, 1), (5, 1) ]
        else {
            val p_max = Primes.max_prime(p)
            if(p_max > fs.tail.head._1)
                put((p_max, 1) :: fs.tail)
        }
        
        // [ (5, 2) ] -> [ (3, 2) ]
        if(fs.size == 1 && e == 2 && p >= 3)
            put(List((Primes.max_prime(p), e)))
        
        fs
    }
    
    private def init() = {
        val p_max = Primes.max_prime(sqrt(N.toDouble).toInt + 1)
        val fs0 = List((p_max, 2))
        put(fs0)
    }
    
    private def put(fs :Facts) :Unit = pq.enqueue((phi_ratio(fs), fs))
}

def digits(n :Long) :List[Int] =
    if(n == 0) Nil else (n % 10).toInt :: digits(n / 10)

def sorted_digits(n :Long) :List[Int] =
    digits(n).sortWith(_ < _)

def is_permutation(m :Long, n :Long) :Boolean =
    sorted_digits(m) == sorted_digits(n)

def solve() = {
    def is_valid(fs :Facts) = is_permutation(value(fs), phi(fs))
    
    val N = 1e7.toLong
    Primes.resize(sqrt(N.toDouble).toInt * 3)
    val it = new numbers_sorted_by_phi_ratio(N)
    println (it.filter(is_valid).map(value).next)
}

object Test extends Benchmark {
    def run() = solve
}

println (Test.runBenchmark(10))