Project Euler 686 / 2の累乗数

問題概略
をとって評価
だけで評価
の候補をしぼりこむ

問題概略

$2^7=128$ は 12 からはじまる最初の 2 の累乗数であり，2 つ目は $2^{80}$ である。
$L$ からはじまる 2 の累乗数のうち小さい方から $n$ 番目の数を $p(L, n)$ とする。
$p(12,1)=7,\, p(12,2)=80,\, p(123, 45)=12710$
$p(123,678910)$ を求めよ。
https://projecteuler.net/problem=686

解説の pdf も作りました。きれいなレイアウトで読みたい方はこちらをどうぞ。

drive.google.com

$\log_2$ をとって評価

$2^n$ の首位 3 桁が 123 のとき $123\cdot 10^k < 2^n <124\cdot 10^k$ のようにおけます。 $k$ は整数です。

\begin{align*}
&k\log_2 10+\log_2 123 < n< k\log_2 10+\log_2 124\mbox{　……(1)}\\
&\therefore \lceil k\log_2 10+\log_2 123\rceil \leqq n
\leqq \lfloor k\log_2 10+\log_2 124\rfloor\mbox{　……(2)}
\end{align*}

(1)の左辺を $l(k)$ ，右辺を $r(k)$ とおいて，(2)の左辺を $L(k)$ ，右辺を $R(k)$ とおきます。

$L(k)=\lceil l(k)\rceil,\, R(k)=\lfloor r(k)\rfloor,\, L(k)\leqq n\leqq R(k)\mbox{　……(3)}$

$l(k)$ と $r(k)$ は無理数で，差は非常に小さい値です。

$r(k)-l(k)=\log_2 124-\log_2 123\approx 0.0116818$

$l(k)$ と $r(k)$ の整数部分が等しいとき， $l(k)\to L(k)$ の変換では小数部分の切り上げがおきて， $r(k)\to R(k)$ では切り捨てがおきます。
$L(k)=R(k)+1$ なので(3)をみたす $n$ はありません。

$l(k)$ と $r(k)$ の整数部分が異なるときは $L(k)=R(k)$ になって $n=L(k)=R(k)$ です。

$L(k)=R(k)$ になる 678910 番目の $k$ に対する $L(k)$ が答えです。約 36 秒かかります。

    import time
    from math import ceil, floor, log2
    
    def left(k):
        return ceil(log2(123) + k * log2(10))
    
    def right(k):
        return floor(log2(124) + k * log2(10))
    
    def solve(n):
        k, cnt = 0, 0
        while cnt < n:
            if left(k) == right(k):
                cnt += 1
            k += 1
        return ceil(left(k - 1))
    
    if __name__ == "__main__":
        start = time.time()
        n = 678910
        print(solve(n))
        elapsed_time = time.time() - start
        print("elapsed_time:{0} [sec]".format(elapsed_time))

$l(k)$ だけで評価

$r(k)=l(k)+\log_2 124-\log_2 123$ なので $l(k)$ と $r(k)$ の整数部分が異なるのは $l(k)$ の小数部分と $\log_2 124-\log_2 123$ の和が 1 以上のときです。

$\left(l(k)-\lfloor l(k)\rfloor\right) +\log_2 124-\log_2 123\geqq 1$

整理すると次のようになります。等号は成立しないので外しました。

$k \log_2 10-\lfloor \log_2 123+k\log_2 10\rfloor> 1-\log_2 124$

この解法は約 18 秒でした。

    import time
    from math import ceil, log2
    
    A, B, C = log2(10), log2(123), log2(124)
    
    def left(k):
        return B + k * A
    
    def isok(k):
        return k * A - int(B + k * A) > 1 - C
    
    def solve(n):
        k, cnt = 0, 0
        while cnt < n:
            if isok(k):
                cnt += 1
            k += 1
        return ceil(left(k - 1))
    
    if __name__ == "__main__":
        start = time.time()
        n = 678910
        print(solve(n))
        elapsed_time = time.time() - start
        print("elapsed_time:{0} [sec]".format(elapsed_time))

両側からの評価をやめたことに加えて，対数を定数扱いしたことが効いた感じです。
ためしに pe686a.py を対数をすべて定数扱いするように書き直したところ，計算時間が 36 秒から 23 秒に減りました。

$k$ の候補をしぼりこむ

階差の候補は3つだけ

20000 以下の $k$ で $L(k)=R(k)$ をみたすものをすべて求めて，その階差をとってみました（mathematica）。
条件をみたす $k$ は 234 個あって，階差は 59, 87, 146 の 3 つの値しかありません。

    In[]:= AbsoluteTiming[
        left[k_] := Ceiling[Log2[123] + k*Log2[10]];
        right[k_] := Floor[Log2[124] + k*Log2[10]]; 
        lst = Select[Range[0, 20000], left@# == right@# &];
        {Length@lst, Sort@DeleteDuplicates@Differences@lst}]
       
    Out[]= {1.12065, {234, {59, 87, 146}}}

なぜ $k$ の階差が限られた値しかとらないのか考えましょう。

ある $k$ で $L(k)=R(k)$ になったとします。
$l(k)$ と $r(k)$ の整数部分は異なるので次のようにおけます。 $m$ は整数で $\alpha$ と $\beta$ は 0 と 1 の間の小数です。

\begin{align*}
&l(k)=m-\alpha,\, r(k)=m+\beta\\
&\qquad (\alpha+\beta=r(k)-l(k)=\log_2 124-\log_2 123\approx 0.0116818)
\end{align*}

整数 $m$ をはさんでそのちょっと下に $l(k)$ があって，ちょっと上に $r(k)$ があるイメージです。

この状態で $k$ を 1 増やすと $l(k)$ , $r(k)$ は $\log_2 10\approx 3.32193$ 増えます。
この数の小数部分は約 $0.32193$ です。
$\alpha\approx 0$ , $\beta\approx 0$ なので $l(k+1)$ と $r(k+1)$ の整数部分は $m$ になって一致してしまいます。

このように「整数のちょっと下」「ちょっと上」を同時にみたすのは難しく， $k$ の階差が制約を受けるわけです。
一方， $\log_2 10$ を 59, 87, 146 回足すときはうまくいきます。以下， $x$ の小数部分を $\langle x\rangle$ であらわします。

\begin{align*}
&\langle 59\cdot \log_2 10\rangle \approx 0.993758=1-0.0062424\\
&\langle 87\cdot \log_2 10\rangle \approx 0.00774426\\
&\langle 146\cdot \log_2 10\rangle \approx 0.00150185
\end{align*}

整数から高々 $\alpha+\beta\approx 0.0116818$ しか離れていない $l(k)$ , $r(k)$ にこれらの値を足しても「整数のちょっと下」「ちょっと上」という状態は保たれます。

コードに落とす

「 $k$ の階差は 59, 87, 146 だけ」を示すのは難しそうですが，「58 回以下足しても無駄」は示せます。
$l(k)$ , $r(k)$ で $k$ を $n$ 増やしたときの小数部分の変化量の絶対値は次の式で与えられます。

$f(n)=\min\left\{\langle n \log_2 10\rangle,\, 1-\langle n \log_2 10\rangle\right\}$

$k=1,\, 2,\, \cdots,\, 58$ に対する最小値は $n=28$ のときの $0.0139867$ です。mathematica で計算しました。

    In[]:= AbsoluteTiming[
        L123 = Log2@123; L124 = Log2@124; L10 = Log2@10;
        f[n_] := N[Min[FractionalPart[L10*n], 1 - FractionalPart[L10*n]]];
        m = Min[f /@ Range@58];
        {m, m < L124 - L123}]
       
    Out[]= {0.0011467, {0.0139867, False}}

$l(k)$ が整数から目一杯離れたとしてその距離は $\alpha+\beta\approx 0.0116818$ で， $k=1,\, 2,\, \cdots,\, 58$ のときの移動量はこれを超えてしまいます。
条件をみたす $k$ の階差として 1～58 はすべて不適です。

まとめると，「ある $k$ が条件をみたすとき，次の $k$ は $k+59$ 以降」です。
この性質を利用することで計算時間を $5.7$ 秒に縮めることができました。

    import time
    from math import ceil, log2
    
    L10, L123, L124 = log2(10), log2(123), log2(124)
    
    def left(k):
        return L123 + k * L10
        
    def isok(k):
        return k * L10 - int(L123 + k * L10) > 1 - L124
    
    def solve(n):
        k, cnt = 0, 0
        while True:
            if not (isok(k)):
                k += 1
            cnt += 1
            if cnt == n:
                return ceil(left(k))
            else:
                k += 59
    
    if __name__ == "__main__":
        start = time.time()
        n = 678910
        print(solve(n))
        elapsed_time = time.time() - start
        print("elapsed_time:{0} [sec]".format(elapsed_time))

未証明の解法＝嘘解法

「 $k$ の階差は 59, 87, 146 だけ」は証明できていませんが，これでコードを書いてみました。
ある $k$ が OK だったら，次の $k$ として $k+59$ , $k+87$ , $k+146$ を順次試すコードです。
これだと $0.45$ 秒で終わります。最初にあげたコードの約 80 倍のスピードです。

    import time
    from math import ceil, log2
    
    L10, L123, L124 = log2(10), log2(123), log2(124)
    
    def left(k):
        return L123 + k * L10
    
    def isok(k):
        return k * L10 - int(L123 + k * L10) > 1 - L124
    
    def solve(n):
        k = 0
        while not (isok(k)):
            k += 1
        cnt = 1
    
        lst = [59, 87, 146]
        while cnt < n:
            flag = False
            for x in lst:
                if isok(k + x):
                    k += x
                    cnt += 1
                    flag = True
                    break
            if flag is False:
                return -1
        return ceil(left(k))
    
    if __name__ == "__main__":
        start = time.time()
        print(solve(678910))
        elapsed_time = time.time() - start
        print("elapsed_time:{0} [sec]".format(elapsed_time))

$146=59+87$ に注目して $k+59$ , $k+87$ だけ試すコードも書いてみたのですが，こちらはうまくいきません。

前に見たように $k+59$ は小数部分を減らす操作で， $k+87$ は増やす操作です。これらだけではうまくいかなくて，適当なタイミングで $k+146$ して，ほんのちょっとだけ小数部分を増やすことが必要なようでした。

砂場で遊ぼう

c++/python/mathematica などの練習帳