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**

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:

**Code (Chapter 2)**

**Code (Pell equation)**

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.

Welcome to our group!

How long does it take to run?

ProjectEuler generally does not welcome sharing solutions to wide public. To quote from the About page:

I learned so much solving problem XXX, so is it okay to publish my solution elsewhere?It appears that you have answered your own question. There is nothing quite like that "Aha!" moment when you finally beat a problem which you have been working on for some time. It is often through the best of intentions in wishing to share our insights so that others can enjoy that moment too. Sadly, that will rarely be the case for your readers. Real learning is an active process and seeing how it is done is a long way from experiencing that epiphany of discovery. Please do not deny others what you have so richly valued yourself.

However, the rule about sharing solutions outside of Project Euler does not apply to the first one-hundred problems, as long as any discussion clearly aims to instruct methods, not just provide answers, and does not directly threaten to undermine the enjoyment of solving later problems. Problems 1 to 100 provide a wealth of helpful introductory teaching material and if you are able to respect our requirements, then we give permission for those problems and their solutions to be discussed elsewhere.

If you value ProjectEuler spirit as much as you appear to, I kindly suggest you remove the solution from the blog. At the very least, you should definitely remove the answer. Codeforces is not a shady answer leak page for cheaters to unfairly inflate their solve counts.

Oh bro, why so many downvotes? It is hard to imagine.