dfsof's blog

By dfsof, 7 months ago, In English

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:

$$$ \left(\frac{a}{p}\right) = \begin{cases} \\{1\\}, x^2 \equiv a (\mod p) \text{ has a solution} \\ \\{-1\\}, x^2 \equiv a (\mod p) \text{ has no solution} \end{cases} $$$

. 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 (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.

Full text and comments »

  • Vote: I like it
  • -20
  • Vote: I do not like it

By dfsof, history, 8 months ago, In English

Orz.hardstone gives me four problems. However, I am so dumb that can only solve the easiest one:

First, we define $$$A$$$ as a non-negative integer array $$$A:=\{a_i\}$$$. We call $$$A$$$ is valid if $$$A$$$ is a [degree sequence](https://en.wikipedia.org/wiki/Degree_(graph_theory)) of a simple undirected graph.

P1: Give $$$A$$$, decide whether $$$A$$$ is valid. (Solved using Erdos-Gallai, $$$O(n)$$$);

P2: Give $$$A$$$ and $$$q$$$ independent queries $$$(l, r)$$$, decide whether $$$A[l...r]$$$ is valid.

P3: Give $$$A$$$, count how many continuous subarrays of $$$A$$$ are valid.

P4: Solve $$$P2$$$ if modification on $$$A$$$ is allowed.

Full text and comments »

  • Vote: I like it
  • +22
  • Vote: I do not like it

By dfsof, history, 9 months ago, In English

Hello lady/bros, I am struggling with 100551E.Disconnected Graph.

I have considered: (1) Online fully dynamic graph connectivity. I copied a piece of code here: https://www.luogu.com.cn/problem/solution/P5247. I passed LuoguP5247 and SPOJDYNACON2, however I could not pass this problem;


Codeforces submission is similar, but it always gets TLE on test15, due to a relatively large constant:

Codeforces Disconnected Graph Submission

(2)Retractable DSU, however it seems that DSU only supports rolling back the add options, it cannot roll back delete operations...

Full text and comments »

  • Vote: I like it
  • 0
  • Vote: I do not like it

By dfsof, 10 months ago, In English

Problem Link: https://mirror.codeforces.com/contest/1864/problem/H

Submission Link: https://mirror.codeforces.com/contest/1864/submission/220919821

Let $$$f(x), x \in \mathbb{Z}$$$ be the expected times starting from $$$x$$$. There are three basic facts:

(1)$$$f(x) = 1 + 1/2f(x+1) + 1/2f(2x)$$$.

(2)$$$f(x) = 0$$$ if and only if $$$x \geq n$$$.

(3)The final answer is $$$f(1)$$$, and $$$f(1)$$$ could be found in $$$O(n)$$$ time naively.

This blog is an extension of CristianoPenaldo's blog. CristianoPenaldo, also known as CP, is one of my best friends besides bfsof. First, similar to CP's idea, I process $$$f$$$ in a reversed and coarse-to-fine manner. I calculate $$$f$$$ from $$$n$$$ to $$$1$$$, and divide the interval $$$[1,n]$$$ into scales as what CP did. When $$$scale$$$ is small, $$$S(scale)$$$ is a "coarse" scale. Otherwise, $$$S(scale)$$$ is a "fine" scale, that is why we call it "coarse-to-fine". Formally, let $$$S(scale)$$$ be $$$\{x \in \mathbb{Z}| x \times scale \geq n, x \times scale/2 < n\}$$$. For example, if $$$n=7$$$, scale $$$1$$$ is $$$[7, 7]$$$, scale $$$2$$$ is $$$[4, 6]$$$, scale $$$4$$$ is $$$[2, 3]$$$, scale $$$8$$$ is $$$[1, 1]$$$. It is guaranteed that the scale is a power of $$$2$$$ in my algorithm.

$$$ S(scale)= \begin{cases} \\{n\\}, \text{scale==1} \\ [\lceil \frac{n}{scale} \rceil, \lceil \frac{n}{scale/2} \rceil - 1] \cap \mathbb{Z},\text{Otherwise} \end{cases} $$$

After some brute force computation, I find that the closed-form formula for scale $$$1$$$ is $$$f(x) = 0$$$ (because there is only one element $$$n$$$ on scale $$$1$$$, and $$$f(n) = 0$$$). The closed-form formula for scale $$$2$$$ is $$$2 - 2(1/2)^{(n-x)}$$$, by the fact that $$$f(2x) = 0$$$ for $$$x$$$ on scale $$$2$$$. The closed-form formula for scale $$$4$$$ is $$$4 - ?(1/2)^{(n-x)} - ?(1/2)^{(n-2x)}$$$, I failed to compute it due to my poor computation ability.

The key obstacle lies in the difficulty handling $$$f(2x)$$$. For $$$x \in S(scale), scale \neq 1$$$, $$$2x$$$ belongs to $$$S(scale/2)$$$. The key idea is to determine the closed-form solutions for each scale recursively. When we determine the closed-form solutions for $$$S(scale)$$$, we can utilize the closed-form solution from $$$S(scale/2)$$$ by simply expanding $$$f(2x)$$$. CP considered using the Berlekamp-Massey algorithm, but it involves matrix binary exponentiation. We will handle the closed-form solution in a more explicit manner.

Theorem: The closed form solution for $$$S(scale)$$$ is

$$$f(x) = C_0+\sum\limits_{i=1}^{\log_2{scale}}C_i (1/2)^{(n-2^{(i-1)}x)}$$$. $$$C_i$$$ are undetermined coefficients.

We can prove via induction. Suppose $$$scale > 1$$$, then $$$f(x) = 1 + 1/2f(x+1) + 1/2f(2x) = 1 + 1/2f(x+1) + 1/2(C_0+\sum\limits_{i=1}^{\log_2{scale}-1}C_i (1/2)^{(n-2^{i}x)})$$$

$$$f(x)-C_0-2 = 1/2(f(x+1)-C_0-2) + \sum\limits_{i=2}^{\log_2{scale}}C_i (1/2)^{(n-2^{i-1}x)}$$$

Suppose $$$f(x) - C_0 - 2 + \sum\limits_{i=2}^{\log_2{scale}}D_i (1/2)^{(n-2^{i-1}x)} = 1/2(f(x+1) - C_0 - 2 + \sum\limits_{i=2}^{\log_2{scale}}D_i (1/2)^{(n-2^{i-1}(x+1))})$$$.

for each $$$D_i$$$, $$$(2^{(2^{i-1}-1)}-1)D_i = 1/2C_i$$$, and $$$D_i = \frac{1}{2^{(2^{i-1})} - 2}C_i$$$.

Let $$$mx$$$ (short for maximum) be the maximum element from this scale. For example, when $$$n=7$$$, the $$$mx$$$ for scale $$$1, 2, 4, 8$$$ are $$$7, 6, 3, 1$$$ respectively.

$$$f(x) - C_0 - 2 + \sum\limits_{i=2}^{\log_2{scale}}D_i (1/2)^{(n-2^{i-1}x)} = (1/2)^{(mx - x)}(f(mx) - C_0 - 2 + \sum\limits_{i=2}^{\log_2{scale}}D_i (1/2)^{(n-2^{i-1}mx)})$$$. And if we calculate $$$f(mx)$$$ in advance, $$$f(mx) - C_0 - 2 + \sum\limits_{i=2}^{\log_2{scale}}D_i (1/2)^{(n-2^{i-1}mx)}$$$ would be a constant, and $$$(1/2)^{(mx-x)}$$$ could be transformed into $$$Constant*(1/2)^{(n-x)}$$$, where $$$Constant$$$ is $$$(1/2)^{(mx-n)}$$$.

By the proof of the above theorem, we can almost get the closed-form solution of $$$S(scale)$$$ from $$$S(scale/2)$$$ except one term $$$(1/2)^{m-x}$$$. To handle this issue, we just fetch the $$$mx$$$ element from that scale and use the method of undeterminated coefficients to calculate $$$C_1$$$, i.e., the coefficient of $$$(1/2)^{m-x}$$$. The closed-form solution of each scale has length $$$O(log scale)$$$, and calculate $$$f(mx)$$$ takes $$$O(log scale \times log n)$$$ time (because the length of closed-form is $$$O(log scale)$$$, and calculate each item, for example, $$$(1/2)^{(n-x)}$$$, involves binary exponentiation, therefore the overall time is $$$\sum\limits_{log scale=1}^{log n} O(log scale \times log n) = O((log n)^3)$$$ per test case.

Full text and comments »

  • Vote: I like it
  • +44
  • Vote: I do not like it

By dfsof, history, 10 months ago, In English

hardstone.Orz gives you an integer array $$$A$$$. The length of $$$A$$$ is $$$n$$$ and there are $$$m$$$ distinct numbers in $$$A$$$. Count the number of tuple $$$(l, r)$$$, $$$1 \leq l \leq r \leq n$$$, such that:

Numbers that appear in the interval $$$a[l...r]$$$ appear the same number of times.

For example, $$$A=[1,2,1,2]$$$, then there are $$$8$$$ legal tuples: $$$(1, 1), (2, 2), (3, 3), (4, 4), (1, 2), (2, 3), (3, 4), (1, 4)$$$.

This is an open problem with brain storm. $$$O(n^2m)$$$ brute force using the prefix sum and $$$O(n2^m)$$$ brute force using bitmasks and hashtable are easy to come up with. I am looking for a $$$O(nmlog^k)$$$ solution. Are there any smart data structures?

Note that when $$$A = [1,2,3]$$$, all intervals are legal. For example, $$$[1, 2]$$$ is legal, as both $$$1$$$ and $$$2$$$ appear once. We do not care about $$$3$$$ because $$$3$$$ does not appear.

amenotiomoi proposes a genius randomized idea, which could make my yesterday's idea work: Similar to the Zobrist hashing, we assign a random value to each distinct integer. We record the prefix sum of the hash values in a hashtable (let $$$h_r$$$ be the prefix sum of hash values of $$$a[1...r]$$$). Then, we fix $$$l$$$ and count the number of $$$r$$$ with respect to this $$$l$$$. For each $$$l$$$, we denote $$$p(l, j)$$$ be the first place $$$j$$$ appears after $$$l$$$ (inclusive), somehow like std::string.find(j, l). If $$$j$$$ never appears after $$$l$$$, $$$p(l, j) = \infty$$$. For example, if $$$A=[4,1,2,3]$$$, then $$$p(2, 1)=2, p(2,2)=3, p(2,3)=4, p(2,4) = \infty$$$. The array $$$p$$$ could be fould via binary search in $$$O(mlogn)$$$. Note that $$$p(l, j) \neq p(l, k)$$$ if $$$j \neq k$$$. Then, we sort the pair $$${p(l, j), j}$$$ in the ascending order of $$$p(l, j)$$$, and let $$$q$$$ be the sorted list. The complexity of sorting is $$$O(mlogm)$$$. For two adjacent elements of $$$q$$$, the present and absent numbers could be uniquely determined. For example, $$$A=[1,2,2,2,3]$$$, $$$l=1$$$, $$$2 \leq r \leq 4$$$, then $$$1, 2$$$ appear and $$$3$$$ is absent. Therefore we need to find the number of $$$r$$$, $$$2 \leq r \leq 4$$$, such that $$$1$$$ and $$$2$$$ appear the same number of times with in $$$a[l...r]$$$. Yesterday I was stuck here. But with the genius hashtable, we only need to count $$$r$$$ that $$$(hashvalue(1) + hashvalue(2)) \mid h_r - h_{l-1}$$$. By the pigeon hole principle, the number appear the least number of times appear at most $$$\frac{n}{i}$$$ times, then we only need to enumerate $$$\frac{n}{i}$$$ items for each adjacent pair of $$$q$$$, there are $$$m$$$ adjacent pairs, and querying the hash table is $$$O(logn)$$$ (using std::map) or amortized $$$O(1)$$$ (using std::unordered_map), therefore the overall complexity could be reduced to $$$O(n(mlogn + mlogm + \sum\limits_{i=1}^m\frac{n}{i}logn)) = O(n(mlogn + mlogm+nlogmlogn))$$$. But this is not deterministic, and the error probability is hard to estimate, heavily depending on implementation.

Full text and comments »

  • Vote: I like it
  • +38
  • Vote: I do not like it

By dfsof, history, 11 months ago, In English

Where are CodeChef contests? Are these contests going to stop due to financial shortage? CodeChef_admin

Full text and comments »

  • Vote: I like it
  • +8
  • Vote: I do not like it

By dfsof, history, 12 months ago, In English

I think the official editorial of problem G is a little bit hard to understand for me, therefore I write a learning note, with an example, in English:

Google Drive (Typo corrected): https://drive.google.com/file/d/1imUKYXcxQNw8wC28YFeTzEAkBV2wswMo/view?usp=sharing

Tencent Docs (Typo corrected): https://docs.qq.com/pdf/DU2VOa09uYUJHYU9E

The above pdf files do not contain code. My code:


and my submission: 211097938.

Be careful when handling indices! Here is a wrong submission with Runtime Error: 211089726.

Full text and comments »

  • Vote: I like it
  • +29
  • Vote: I do not like it

By dfsof, history, 12 months ago, In English

These three problems are almost the same.

CF618F [Double Knapsack]: https://mirror.codeforces.com/problemset/problem/618/F

CF1836E [Twin Cluster]: https://mirror.codeforces.com/contest/1836/problem/E

Beijing College Entrance Exam:

Given two positive integer arrays $$$A$$$ and $$$B$$$, such that:

$$$len(A) = len(B) = n$$$


$$$\forall 1 \leq i \leq n$$$, $$$1 \leq a_i, b_i \leq n$$$.

Prove there are subsegments $$$[x, y] \subseteq [1, n]$$$, $$$[z, w] \subseteq [1, n]$$$ such that $$$\sum\limits_{i=x}^y A[i] = \sum\limits_{j=z}^w B[j]$$$. For example, $$$n=4, A=[1,2,3,4], B=[4,4,4,4]$$$, then $$$x=4, y=4, z=1, w=1$$$ is a solution.

Full text and comments »

  • Vote: I like it
  • -7
  • Vote: I do not like it