---
tags: 数学, 整数問題
---
今年の整数問題（2017/01/10）
==================

## 問題

【問題1】
$6,5,6,9,29,4,29,1,61,24,\dots$ という数列の、$n$番目は$209$であり、$50$番目は$n$である。$n$はいくつか。ただし今日(2017/01/10)は私の$n$歳の誕生日です。

【問題2】
この数列の$k$番目の値を$a_k$と書くと、$k=3$は$a_k=2k$となる最小の値である（$a_3=6$）。$a_k=k$となる最小の$k$を求めよ。

## 解答

【問題1】
$n=41$

【問題2】
$k=52892561984$

### 解説

この数列 $\{a_i\}$ は、以下のような法則で成り立っている：

+ $b_i$ を $i$ と $a_i$ を10進表記で続けて記述してできる数としたとき、$b_i$ は平方数である。
+ （補足）$a_i$ を10進表記で$m$桁とした場合、$b_i = 10^mi+a_i$ である。
+ $a_i$ の（10進表記の）上1桁は$0$ではない。
+ $a_i$ は、上記2つの条件を満たす数のうちの最小値である。

例：

| $i$ | $a_i\ (b_i)$ |
| --: | :-- |
| $1$ | $6\ (16=4^2)$ |
| $2$ | $5\ (25=5^2)$ |
| $3$ | $6\ (36=6^2)$ |
| $4$ | $9\ (49=7^2)$ |
| $5$ | $29\ (529=23^2)$ |
| $6$ | $4\ (64=8^2)$ |
| $7$ | $29\ (729=27^2)$ |
| $8$ | $1\ (81=9^2)$ |
| $9$ | $61\ (961=31^2)$ |
| $10$ | $24\ (1024=32^2)$ |
| $11$ | $56\ (1156=34^2)$ |
| $12$ | $1\ (121=11^2)$ |
| : | : |
| $24$ | $336\ (24336=156^2)$ |
| : | : |

例えば $2401=49^2$ だが、$a_{24}$ は $1$ ではない（$241$ は平方数ではない）。一方で $a_{240}=1$ である。

【問題1】は、$a_{50}$ を求めればOK。その値（$m$と置く）で $a_m=209$ であることを確認（検算）すれば良い。

一般の $n$ について $a_n$ を求めるアルゴリズム（概要）は、以下の通りとなる：

1. $d := 10$, $m := 10n + 1$ と置く
2. $r := \lceil \sqrt{m} \rceil^2 - dn$ と置く
3. $r < d$ ならば、$r$ を返して終了。
4. $d \leftarrow 10d$, $m \leftarrow 10m$ として、2. に戻る

[Julia](https://julialang.org) で記述したコード例がこちら：

function an(n::T) where {T<:Integer}
    d = n10 = T(10)
    m = n * n10 + one(T)
    while true
        r = ceil(T, sqrt(m)) ^ 2 - n * d
        if r < d
            return r
        end
        d *= n10
        m *= n10
    end
end

これにより、`an(50) # => 41`（$a_{50}=41$）が得られる。また `an(41) # => 209`（$a_{41} = 209$）も確かめられる。

【問題2】は、このアルゴリズムを用いて $a_k=k$ となる $k$ を探すのだが、それを愚直に $k=1$ から順に実施すると全く終了しない。

$a_k=k$ ということは、$k$ を $m$桁とすると、$(10^m+1)k$ が平方数になる、ということである。よって少なくとも $10^m+1$ が平方因子をもたなければならない（$\because$ $10^m+1$ が平方因子を持たない（$\Leftrightarrow$素因数分解したときの各素数の冪数がすべて$1$ の）場合、$(10^m+1)k$ が平方数であるためには $k$ が $10^m+1$ の倍数である必要があるが、この時 $k$ は $m+1$ 桁以上になる）。
これを利用して $m$ の候補を絞り込んで $k$ を探せば良い。

詳細は省略するが、例えば [Julia](https://julialang.org) なら以下のようなコードを記述すれば求解出来る（要：`Primes.jl`）。

# Pkg.add("Primes")
using Primes

function solve(verbose=false)
    n1 = big"1"
    n10 = big"10"
    i = n1
    while true
        i *= n10    # = 10^m
        q = i + n1  # = 10^m + 1
        pd = factor(q)  # 素因数分解（要：`Primes.jl`）
        if any(v > 1 for v ∈ values(pd))
            r = prod(v == 1 ? p : p ^ ((v + 1) ÷ 2) for (p, v) ∈ pd)
            l = r * r % i
            d = isqrt(i ÷ n10 ÷ l)
            r *= d + n1
            k_cand = r * r ÷ q  # ← 10^(m-1) ≤ k < 10^m かつ (10^m + 1)k が平方数となるような k の最小値
            for k = [k_cand, 4k_cand, 9k_cand]
                if k >= i
                    break
                elseif an(k) == k
                    if verbose
                        println("a(", k, ") = ", k, " : O")
                    end
                    return k
                elseif verbose
                    println("a(", k, ") = ", an(k), " : X")
                end
            end
        end
    end
end

`solve(true)` とすると、以下のように途中経過も含めて求解された結果が表示される。

julia> solve(true)
a(13223140496) = 9685186064 : X
a(52892561984) = 52892561984 : O
52892561984