This blog targets at persons who are interested in (computational) number theory, or struggling to solve this problem. It contains many spoilers (even the final answer). Therefore, if you want to solve it on yourself, please close immediately (or add an upvote). The idea of this blog is quite similar to ecnerwala's comments on the Project Euler thread, but the thread is not open and only visible to persons who solved this problem. In fact, I have finished this problem by myself.

1. Problem Statements

Link to problem.

Let $$$C(n)$$$ be the number of square-free integers of the form $$$x^2+1$$$ ($$$1 \leq x \leq n$$$). For example, $$$C(10)=9$$$ (only $$$7 \times 7+1 = 2 \times 5 \times 5$$$ is not square-free) and $$$C(1000)=895$$$. Find $$$C(123567101113)$$$.

2. The basic idea, and obstacle

The first basic idea is the principle of inclusion-exclusion (PIE). For $$$d$$$ not necessarily prime, the final answer $$$C(n)$$$ is $$$C(n) = \sum\limits_{d=1}^n\mu(d) \#\{x \text{ such that } d \mid (x^2+1) \} \tag{1}$$$.

Here, $$$\mu (\mathbb{N}^* \rightarrow \{-1, 0, 1\})$$$ is the Mobius function, i.e., If $$$n$$$ is not square free, $$$\mu(n) = 0$$$. If $$$n$$$ is square free, then $$$\mu(n)$$$ is $$$1$$$ ($$$-1$$$) if $$$n$$$ has even (odd) number of distinct prime factors. Specially, $$$\mu(1) = 1$$$, as $$$1$$$ is square-free and has zero prime factor. Here is a simple example, $$$268^2 + 1 = 71825 \equiv 0 (\mod 65 ^ 2 = 4225)$$$, so $$$268$$$ should be discounted once. However, $$$5$$$ and $$$13$$$ discounts two times, so $$$5 \times 13 = 65$$$ should add once. $$$\#$$$ is the cardinality.

There are some basic number theory facts. First, the Legendre symbol $$$\left(\frac{a}{p}\right)$$$, where $$$p$$$ has to be a prime, is defined as:

. There are many interesting facts of $$$\left(\frac{a}{p}\right)$$$, e.g., the Gauss's Lemma， and Quadratic Reciprocity. However, we are only interested in one important lemma: $$$\left(\frac{-1}{p}\right)$$$ is $$$1$$$ iff $$$p=2$$$ or $$$p=4k+1$$$. If $$$p = 4k+1$$$, $$$x^2 + 1 \equiv 0 (\mod p)$$$ has two distinct solutions modulo $$$p$$$, which are $$$(\frac{p-1}{2})!$$$ and $$$-(\frac{p-1}{2})!$$$ respectively. If you are not familiar with such lemma, see tutorial Chapter 9.6. For $$$x^2+1 \equiv 0 (\mod p^2)$$$, obviously $$$p \neq 2$$$, and we can use Hensel lifting to uniquely lift a solution of $$$x^2+1 \equiv 0 (\mod p)$$$ to $$$x^2+1 \equiv 0 (\mod p^2)$$$, so the latter equation also has two solutions modulo $$$p^2$$$. For example, when $$$p=29$$$, the two solutions are $$$41$$$ and $$$800$$$ ($$$41^2+1 = 1682, 29^2=841$$$). By CRT, if $$$d$$$ is an odd number with no $$$4k+3$$$ type prime factor, then there are $$$2^{\omega(d)}$$$ ($$$\omega(d)$$$ is the number of distinct prime factors) solutions of $$$x^2+1 \equiv 0 (\mod n)$$$. After we solve $$$x^2 \equiv -1 (\mod d)$$$, $$$\#\{x \text{ such that } d \mid (x^2+1) \} = \lfloor \frac{n}{d} \rfloor 2^{\omega(d)}+ \text{Some Round Up}$$$. For example, when $$$d=25$$$, $$$7$$$ and $$$18$$$ are solutions to $$$x^2 \equiv -1 (\mod d)$$$. If $$$n=32$$$, $$$32$$$ will round up. If $$$n=31$$$, no such round up. I find the round up really annoying, it seems that the best way to deal with the round up that I can come up with is bisecting the whole solution list of length $$$2^{\omega(d)}$$$.

Such a process could be shorten as: Factor integer -> Find quadratic residue (e.g., the Cipolla algorithm) -> Hensel Lifting -> CRT -> bisecting to calculate RoundUps. However, when $$$n$$$ is large (e.g., $$$\sim 10^{11}$$$), every step is so difficult.

3. Balancing for large d, Negative Pell equations

If $$$d$$$ is large, and $$$x^2+1=kd^2$$$, then $$$k$$$ is small. Such equation is called the Negative Pell's equation, also known as Pell equation of the second type, if $$$k$$$ is square free. The key idea is to use the direct method for small $$$d$$$, and the Pell equation method for large $$$d$$$ (here, $$$k$$$ is small). The relation is:

(1)Small $$$d$$$ are only dealt using the method in chapter 2;

(2)The pell equation generates both small $$$d$$$ and large $$$d$$$, so we need to do some de-duplication. However, the pell equation does not generate "too many solutions".

I choose the SymPy library, which uses the LMM algorithm to get a fundamental solution ($$$x_0, d_0$$$). Here, fundamental means $$$x_0 + \sqrt{k}d_0$$$ is the smallest among all solutions. For a negative pell equation, the fundamental solution does not necessarily exist, but as long as it exists, the equation has infinitely many solutions, all of which are of the form $$$x + \sqrt{k}d = (x_0 + \sqrt{k}d_0)^{2m+1}$$$. Hence, each $$$k$$$ only generates $$$O(log n)$$$ solutions. Here, we need to enumerate all solutions, therefore the binary exponentiation technique is useless here.

Here we need to pay attention that $$$k$$$ is required to be square-free. Hence, some solutions are omitted. For example, if $$$d=65$$$, $$$268^2 - 17 \times (65^2) = -1$$$, the solution $$$(268, 65)$$$ is ok, not omitted. However, for $$$d=13, k=17 \times 25=425$$$, $$$268^2 - 17 \times 25 \times (13^2) = -1$$$, $$$(268, 13)$$$ is omitted as $$$17 \times 25$$$ is not square free. Be careful!

4. Implementation

I set the upper bound of $$$k$$$ to $$$160000$$$, hence the method in chapter 2 only deals $$$d \leq \lfloor \sqrt{ \frac{123567101113^2+1}{160000} }\rfloor = 308917752$$$. Large $$$d > 308917752$$$ are dealt via Pell equations in chapter 3.

I use the SymPy library, Chinese Zhihu, as there are three very powerful functions:

(1)fast sympy.factorint;

(2)from sympy.ntheory import sqrt_mod to do all the steps in Chapter. 2 except bisecting (for example, sqrt_mod(-1, 65**2, all_roots=True));

(3)The most important, diop_DN to find fundamental solutions or report no solution. diop_DN returns either a singleton list containing a fundamental solution $$$(x_0, d_0)$$$ represented by a Python tuple, or returns an empty list.

The algorithm in Chapter. 2 and Chapter.3 can be run in parallel, you might organize them into two Python files.

Code:

import sympy
import sys
from sympy.solvers.diophantine.diophantine import diop_DN
from sympy.ntheory import sqrt_mod 
from math import isqrt
import bisect

debug_pipeline = False
THRES = (int(sys.argv[1]) if len(sys.argv) > 1 else 20000) if debug_pipeline else 123567101113 
BALANCE_PARAM = 10**2 if debug_pipeline else 400 ** 2
ENUMERATE_THRESHOLD = THRES
if BALANCE_PARAM > 1:
    ENUMERATE_THRESHOLD = isqrt(((THRES**2 + 1) + (BALANCE_PARAM - 1))//BALANCE_PARAM)
    if (ENUMERATE_THRESHOLD ** 2) * BALANCE_PARAM == THRES**2 + BALANCE_PARAM:
        ENUMERATE_THRESHOLD -= 1
print("ENUMERATE_THRESHOLD", ENUMERATE_THRESHOLD)
FILE_LARGE = 'dump_large.txt'
FILE_DEBUG = 'dump_debug.txt'
global_set = set()
debug_set = set()
filtered = set()

class SmallNumber:
    def __init__(self, n):
        assert n % 4 == 1
        self.n = n
        self.error_msg = ""
        if n in filtered:
            self.ok = False
            return
        self.ok = True
        self.d = sympy.factorint(n)
        for k, v in self.d.items():
            if k % 4 == 3:
                self.ok = False
                self.error_msg = "Find prime %s of 4k+3"%k
                return
            if v >= 2:
                self.ok = False
                self.error_msg = "Find prime square %s^%s for n%s"%(k, v, n)
                return

    def need_computation(self):
        return self.ok, self.error_msg

    def compute(self, thres=THRES):
        num = 0
        need_computation, msg = self.need_computation()
        if not need_computation:
            return num, msg
        sol = sorted(list(sqrt_mod(self.n**2 - 1, self.n**2, True)))
        num = (thres // self.n**2) * len(sol) + bisect.bisect_right(sol, thres % self.n**2)            
        return num * (1 if len(self.d) % 2 == 1 else -1), msg


class NegativePell:
    def __init__(self, k):
        self.x = -1
        self.y = -1
        d = sympy.factorint(k)
        self.error_msg = ""
        for k1, v in d.items():
            if v >= 2:
                self.state = 2
                #self.error_msg = "Find prime square %s^%s for k %s"%(k1, v, k) #Accelerate
                return
        
        l = diop_DN(k, -1)
        if not l:
            self.state = 1
            #self.error_msg = "No solution found for pell equation x^2 - %sy^2 = -1!"%(k)
            return
        x, y = l[0]
        assert(x**2 - k*(y**2) == -1)
        self.state = 0
        self.x = x
        self.y = y
        self.x2 = x**2 + k * (y**2)
        self.y2 = 2*x*y
        self.k = k

    def meta(self):
        return self.state, self.error_msg, self.x, self.y

    def next(self, curx, cury, check=False):
        nextx = self.x2 * curx + self.y2 * self.k * cury
        nexty = self.x2 * cury + self.y2 * curx
        if check:
            assert nextx**2 - self.k * (nexty**2) == -1
        return nextx, nexty


def solveNegativePell(k, check=True):
    pellsolver = NegativePell(k)
    state, msg, x, y = pellsolver.meta()
    if state != 0:
        return False, msg
    while x <= THRES:
        if y != 1:
            global_set.add(x)
            debug_set.add((x, y, k))
        x, y = pellsolver.next(x, y, check)
        if check:
            assert x not in global_set, "%s[n] is unexpectly duplicated!"%x 
    return True, msg


if __name__ == '__main__':
    ans = 0
    i = 0
    print("FILTERING...")
    for i in range(1, ENUMERATE_THRESHOLD+1, 2):
        j = 3
        while i * j <= ENUMERATE_THRESHOLD:
            filtered.add(i*j)
            j += 4
        k = i**2
        j = 1
        while k != 1 and k % 4 == 1 and k * j <= ENUMERATE_THRESHOLD:
            filtered.add(k*j)
            j += 4
    print("AFTER FILTERING, len(filtered)==%s"%len(filtered))
    print(sorted(list(filtered)))
    i = 0
    for n in range(5, ENUMERATE_THRESHOLD+1, 4):
        sn = SmallNumber(n)
        contrib, msg = sn.compute()
        if i % 1000 == 0:
            print(i, n, ans)
            with open("stage2.log", "w") as f:
                f.write("%s %s %s\n"%(i, n, ans))
        ans += contrib
        i += 1
    print(ans) #12994164947

#The below code generates pell solutions
for i in range(2, BALANCE_PARAM+1):
    has_solution, _ = solveNegativePell(i)
    if i % 1000 == 0: print(i, len(global_set))
large_solutions = sorted(list(global_set))
debug_solutions = sorted(list(debug_set))
with open(FILE_LARGE, 'w') as f:
    for sol in large_solutions:
        f.write("%s\n" % sol)
with open(FILE_DEBUG, 'w') as f:
    for x, y, n in debug_solutions:
        f.write("%s %s %s\n" % (x, y, n))

#The below code get statistics
from sympy import factorint

gt = 308917752+1
cnt = 0
with open("dump_debug.txt") as f:
    l = f.readlines()
    for line in l:
        line = list(map(int, line.strip().split()))
        if line[1] >= gt:
            print(line)
            d = factorint(line[1])
            ok = True
            l = list(d.keys())
            print(l, 'l')
            contrib = 0
            for msk in range(1<<len(l)):
                omitted = 1
                popcount = 0
                for k in range(len(l)):
                    if msk & (1<<k):
                        popcount += 1
                        omitted *= l[k]
                if omitted >= gt:
                    contrib += (1 if popcount % 2 == 1 else -1)
                print(line, d, msk, 'msk', omitted, 'omitted', popcount, 'popcnt', contrib, 'contrib')
                
            cnt += contrib
            print(line, contrib, cnt, d)
print(cnt) #-11

Sorry for the extremely poor code quality, I get insomnia after every CF round (都是网瘾害的)!

5. Answer

Not shown.

6. (For Chinese Readers) An invitation to our QQ group:

My grandma Aveiro_quanyue and me are co-organizing a QQ chat group. If you are interested, please add my grandma （QQ number $$$3381896043$$$, nickname "全月"). It focuses on three aspects: MATH, DS (Data Structure) and CP (Competitive Programming). Here are the reasons why you should join:

(1) The CF ratings of our group members are between $$$1600-2800$$$. Therefore, I believe you can almost always find a member with similar rating to compete and/or share ideas. Although CF scores vary widely among group members, we communicate with each other in a very friendly and equal manner.

(2) Our group is informative. We are sharing brilliant ideas and useful learning materials (e.g., PDF e-book or learning notes) with others, and we hold reading seminars regularly. Currently we are reading Donald Knuth's Concrete Mathematics and some number theory stuff. I believe our group is much better than some other XCPC groups that actually focus on some sexy stuff. Our group is very small, currently only 32 people, so it’s relatively easy to manage (filter useless information).

(3) Everybody in our group has her (or his) strength, so never look down upon anyone (e.g., low-rated like me) in our group. Some people have outstanding CF ratings, some constantly won XCPC gold medals and entered ICPC World Final, some have incredibly high GPA rankings, some are data structure masters, and some have extraordinary business talents. As for me, I am almost the most low-rated in our group, but I think I am a slow thinker and good at solve hard problems (especially math). This group offers a good chance for you to work with outstanding partners.

(4) The group leader is a kind old lady who gives each of her members a nickname. In addition, she will give award to students who solve difficult problems.