Project Euler 399

发表于 分类于 阅读次数： 本文字数： 3.7k 阅读时长 ≈ 3 分钟

Project Euler 399

题目

Squarefree Fibonacci Numbers

The first \(15\) fibonacci numbers are:

\(1,1,2,3,5,8,13,21,34,55,89,144,233,377,610.\)

It can be seen that \(8\) and \(144\) are not squarefree: \(8\) is divisible by \(4\) and \(144\) is divisible by \(4\) and by \(9\).

So the first \(13\) squarefree fibonacci numbers are:

\(1,1,2,3,5,13,21,34,55,89,233,377\) and \(610.\)

The \(200\text{th}\) squarefree fibonacci number is:

\(971183874599339129547649988289594072811608739584170445.\)

The last sixteen digits of this number are: \(1608739584170445\) and in scientific notation this number can be written as \(9.7e53\).

Find the \(100 000 000\text{th}\) squarefree fibonacci number.

Give as your answer its last sixteen digits followed by a comma followed by the number in scientific notation (rounded to one digit after the decimal point).

For the \(200\text{th}\) squarefree number the answer would have been: \(1608739584170445,9.7e53\)

Note:

For this problem, assume that for every prime \(p\), the first fibonacci number divisible by \(p\) is not divisible by \(p^2\) (this is part of Wall’s conjecture). This has been verified for primes \(\le 3\cdot 10^{15}\), but has not been proven in general.

If it happens that the conjecture is false, then the accepted answer to this problem isn’t guaranteed to be the \(100 000 000\text{th}\) squarefree fibonacci number, rather it represents only a lower bound for that number.

解决方案

令\(N=10^8\)。我们给题目额外假设（Wall 猜想相关）：对任意素数 \(p\)，第一次出现 \(p\mid F_{a_p}\) 时不会发生 \(p^2\mid F_{a_p}\)（也就是 \(v_p(F_{a_p})=1\)）。

论文里定义了 rank of apparition：\(\kappa(n)=\min\{m\ge 1:n\mid F_m\}.\)并在引理1给出关键性质：对任意正整数 \(m,n\)，有\(n\mid F_m \iff \kappa(n)\mid m.\)

因此，对素数 \(p\)：

  • 令 \(a_p = \kappa(p)\)（最小使得 \(p\mid F_{a_p}\) 的下标）。
  • 那么 \(p\mid F_n\) 当且仅当 \(a_p\mid n\)。

论文的引理4给出：设 \(p^l\) 是 \(F_{\kappa(p)}\) 中 \(p\) 的最高幂（即 \(p^l\mid F_{\kappa(p)}\) 且 \(p^{l+1}\nmid F_{\kappa(p)}\)），则\(\kappa(p^{l+1})=p\kappa(p).\)

由于 \(\kappa(p)\) 的定义保证 \(p\mid F_{\kappa(p)}\)，因此 \(F_{\kappa(p)}\) 中 \(p\) 的指数至少为 \(1\)。而题目给定的 Wall 假设进一步要求 \(p^2\nmid F_{\kappa(p)}\)，这就意味着该指数不可能达到 \(2\)。因此，引理4中的\(l\)只能是\(1\)。

把 \(l=1\) 代入引理4，得到关键结论：\(\kappa(p^2)=p\kappa(p)=pa_p.\)

由引理 1 对 \(n=p^2\) 的情形，得到\(p^2\mid F_n \iff \kappa(p^2)\mid n.\)再代入上一步得到的 \(\kappa(p^2)=pa_p\)，可化为 \(p^2\mid F_n \iff (pa_p)\mid n.\)

因此，\(F_n\) 不是 squarefree 当且仅当存在某个素数 \(p\) 使得 \(p^2\mid F_n\)，等价于存在 \(p\) 使得\((pa_p)\mid n.\)换句话说，令\(b_p=pa_p,\)则 squarefree 的下标 \(n\) 恰好是所有不被任何 \(b_p\) 整除的正整数。

于是原问题等价为：在 \(n=1,2,3,\dots\) 中筛去所有满足 \(b_p\mid n\) 的下标，剩下的下标按从小到大第 \(N\) 个记为 \(n^{\ast}\)，答案即为 \(F_{n^{\ast}}\)

直接模拟斐波那契模 \(p\) 直到出现 \(0\) 虽然可行，但总体步数很大。更快的做法是利用论文证明定理1到定理2时引用的经典结论：对任意素数 \(p\ne 5\)，有 \(p\mid F_{p-\left(\frac{p}{5}\right)}.\)其中\(\left(\dfrac{a}{b}\right)\)是指勒让德符号。

因此 \(a_p=\kappa(p)\)（最小使得 \(p\mid F_{a_p}\) 的下标）一定整除

\[ B= \begin{cases} p-1,&\text{if}\quad \left(\frac{p}{5}\right)=1\quad(p\equiv \pm 1\pmod 5),\\ p+1,&\text{if}\quad \left(\frac{p}{5}\right)=-1\quad(p\equiv \pm 2\pmod 5). \end{cases} \]

于是可以用按素因子缩小阶的方法求最小 \(a_p\)：

  1. 先令 \(k=B\)（由上面的结论可知 \(p\mid F_k\)，因此 \(a_p\mid k\)）。
  2. 分解 \(k\)。对每个素因子 \(q\)，尝试把 \(k\) 除以 \(q\)：若 \(p\mid F_{k/q}\) 仍成立，就令 \(k\leftarrow k/q\) 并继续；否则停止对这个 \(q\) 的缩小。
  3. 最终的 \(k\) 就是 \(a_p\)（因为经过逐个素因子的“能除就除”，\(k\) 已经被削到仍满足 \(p\mid F_k\) 的最小可能值）。

判断 \(p\mid F_t\) 用斐波那契快速倍增法在 \(O(\log t)\) 时间算出 \(F_t\bmod p\) 即可。

注意，斐波那契快速倍增是通过把下标从 \(k\) 一次性翻到 \(2k\) 或 \(2k+1\)，用 \(O(\log n)\) 步就能算出 \(F_n\)（常同时得到 \(F_{n+1}\)）。基于下面两条倍增恒等式：设 \(F_0=0,F_1=1\)，对任意 \(k\ge 0\)，\(F_{2k}=F_k(2F_{k+1}-F_k),F_{2k+1}=F_{k+1}^2+F_k^2.\)因此用 \(O(\log n)\) 步就能算出 \(F_n\)。

因此，有了所有 \(b_p=p\cdot a_p\)（且只需保留 \(b_p\) 不超过上界 \(L\) 的那些），就可以像筛合数一样筛下标：

  • 对每个 \(b\)，把 \(b,2b,3b,\dots\) 这些下标标记为“坏”（对应 \(F_n\) 含某个 \(p^2\) 因子）。
  • 未被标记的下标就是 squarefree Fibonacci 的下标。

因此，实际目标下标在 \(1.31 N\) 左右，取 \(L=2N\) 足够覆盖；最终找到第 \(N\) 个未被标记的下标 \(n^{\ast}\)。

对于\(F_{n^{\ast}}\)的低16位，直接使用斐波那契快速倍增计算\(F_{n^{\ast}}\bmod M\)即可，其中\(M=10^{16}\)。

对于\(F_{n^{\ast}}\)的科学计数法，用 Binet 近似的对数形式即可：

\[ \log_{10}F_n \approx n\log_{10}\varphi-\log_{10}\sqrt 5, \varphi=\frac{1+\sqrt5}{2}. \]

令 \(x=\log_{10}(F_n)\)，则指数 \(e=\lfloor x\rfloor\)，尾数 \(m=10^{x-e}\)（\(1\le m < 10\)）。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
import math
from typing import List, Tuple
from tools import *

N = 100_000_000
L = 200_000_000
P_LIMIT = 2_200_000
MOD = 10 ** 16


# ---------------------------
# Fast doubling: Fibonacci mod
# F0 = 0, F1 = 1
# ---------------------------
def fib_mod(n: int, mod: int) -> int:
    def rec(k: int) -> Tuple[int, int]:
        if k == 0:
            return 0, 1
        a, b = rec(k >> 1)          # a = F_k, b = F_{k+1} (for k = floor(k/2))
        c = (a * ((2 * b - a) % mod)) % mod
        d = (a * a + b * b) % mod
        if k & 1:
            return d, (c + d) % mod
        else:
            return c, d
    return rec(n)[0]

# ---------------------------
# Legendre symbol (5/p) via Euler criterion
# returns +1 or -1 for odd prime p != 5
# ---------------------------
def legendre_5(p: int) -> int:
    t = pow(5, (p - 1) // 2, p)
    if t == 1:
        return 1
    if t == p - 1:
        return -1
    return 0  # happens only if p == 5

# ---------------------------
# Rank of apparition a_p = kappa(p):
# smallest n>=1 with F_n == 0 (mod p)
# using that a_p | (p - (p/5))
# ---------------------------
def rank_ap(p: int, factorizer) -> int:
    if p == 2:
        return 3
    if p == 5:
        return 5

    s = legendre_5(p)
    base = p - 1 if s == 1 else p + 1

    k = base
    primes = [p for p,e in factorizer.factor_small(base)]
    # unique prime factors (order doesn't matter)
    last = 0
    for q in sorted(primes):
        if q == last:
            continue
        last = q
        while k % q == 0:
            cand = k // q
            if fib_mod(cand, p) == 0:
                k = cand
            else:
                break
    return k

# ---------------------------
# Build moduli b_p = p * a_p (<= L)
# ---------------------------
def build_bad_moduli(L: int, P_LIMIT: int) -> List[int]:
    factorizer = Factorizer(P_LIMIT + 1)
    primes, spf = factorizer.primes, factorizer.spf

    moduli = set()

    # handle primes up to P_LIMIT
    for p in primes:
        a = rank_ap(p, factorizer)
        b = p * a
        if b <= L:
            moduli.add(b)

    # (Optional completeness) If p > P_LIMIT and p*a_p <= L then a_p <= L/(P_LIMIT+1).
    # Choosing P_LIMIT=2_200_000 and L=200_000_000 gives a_p <= 90, and F_90 < 2^64.
    # In practice for this problem, primes up to P_LIMIT already cover all needed b.
    return sorted(moduli)

# ---------------------------
# Segmented sieve over indices to find the N-th "good" index
# good means: not divisible by any b_p => F_n is squarefree under the assumption
# ---------------------------
def nth_squarefree_fib_index(N: int, L: int, moduli: List[int], SEG: int = 5_000_000) -> int:
    count = 0
    start = 1
    while start <= L:
        end = min(L + 1, start + SEG)
        seg = bytearray(end - start)  # 0 = good, 1 = bad

        for b in moduli:
            # find first multiple of b in [start, end)
            m = ((start + b - 1) // b) * b
            if m >= end:
                continue
            # mark all multiples
            for x in range(m, end, b):
                seg[x - start] = 1

        good_in_seg = seg.count(0)
        if count + good_in_seg < N:
            count += good_in_seg
            start = end
            continue

        # locate exact position in this segment
        pos = 0
        while True:
            i = seg.find(0, pos)
            if i == -1:
                break
            count += 1
            if count == N:
                return start + i
            pos = i + 1

        start = end

    raise RuntimeError("L is too small; N-th squarefree index not found within [1..L].")

# ---------------------------
# Scientific notation for F_n: mantissa with 1 decimal + exponent
# ---------------------------
def fib_scientific(n: int) -> Tuple[float, int]:
    # Use logs: log10(F_n) ~ n*log10(phi) - log10(sqrt(5))
    phi = (1.0 + math.sqrt(5.0)) / 2.0
    x = n * math.log10(phi) - 0.5 * math.log10(5.0)
    e = int(math.floor(x))
    m = 10.0 ** (x - e)

    # round to 1 decimal
    m = round(m, 1)
    if m >= 10.0:
        m = 1.0
        e += 1
    return m, e



moduli = build_bad_moduli(L, P_LIMIT)
idx = nth_squarefree_fib_index(N, L, moduli)

last16 = fib_mod(idx, MOD)
mant, exp = fib_scientific(idx)

print(f"{last16:016d},{mant:.1f}e{exp}")