grhkm's blog

By grhkm, history, 7 months ago, In English

Rejoice.

https://arxiv.org/pdf/2404.05177.pdf

From what I heard, this seems practical, so maybe there will be uses in OI?

Full text and comments »

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

By grhkm, history, 13 months ago, In English

Hi, UKIEPC 2023 just concluded yesterday, the problems can be found here and the solutions can be found here.

I don't understand the solution to H. You should read the problems documentation linked above, but here is a short description. The problem asks for a data structure that supports two types of queries: (1) range add / sub (2) after removing duplicate consecutive elements of the array, check whether all "local minimas" of the array is strictly increasing, where "local minima" is defined as a[i] such that a[i — 1] > a[i] && a[i] < a[i + 1], ignoring the condition when out of bounds.

The solution says we can maintain a segment tree on the ranges of "increasing" and "non-increasing" sequences for the array. However, I am not sure what the segment tree should contain at all. Can someone provide a more detailed explanation? Any similar problems in the past (hopefully with writeup) will be great. I know the basic stuff about segment tree (range add range max etc.) and I want to understand the solution of this.

Thanks!

Full text and comments »

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

By grhkm, history, 17 months ago, In English

I was reading adamant's great blog posts when I found multiple comments, some by LGM's, just shitting on him and his blogs for no reason, e.g. by ko_osaga at https://codeforces.me/blog/entry/115164?#comment-1023672, early-morning-dreams under the same post, Xaxixin again under https://codeforces.me/blog/entry/117195?#comment-1036456 and many other examples. I don't understand what these people's intentions are, if it's not your type of posts then why comment? Yes, you might not understand the content, but codeforces is such a big community now that there's definitely many people who are still interested, even if it's remotely related to OI. Back when MiFaFa (jqdai0815) posted his blog on differential equations a few years ago, the community at least respected the post and had (some) constructive discussion below, whereas nowadays it seems that people just shit on others for not much reason. Am I the only one feeling this?

Even worse, there's no block option on CF as far as I know, which means every time I scroll and read I see these people just going "what's the point of learning this".

By the way, if adamant is reading this, I really enjoy reading your posts. As a math enjoyer myself, I love reading your posts since they are usually explained very well. You should try PE830 and PE831 if you haven't, seems like your kind of problem :D

Full text and comments »

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

By grhkm, history, 3 years ago, In English

As promised, here's "interpolating multiple values" in $$$O(n\log^2 n)$$$.

Pre-req: Polynomial multiplication and modulos, basic lagrange interpolation knowledge (you can read my previous post here)

As a reminder, polynomial mod is $$$O(n\log n)$$$.

Main results:

Given $$$f(x)$$$ and $$$x_1, x_2, \cdots, x_n$$$, we can find $$$f(x_1), f(x_2), \cdots, f(x_n)$$$ in $$$O(n\log^2 n)$$$.

Given a bunch of points $$${(x_0, y_0), (x_2, y_2), \cdots, (x_n, y_n)}$$$, we can find a degree $$$n$$$ polynomial $$$f(x)$$$ such that

$$$\forall 0\leq i\leq n, f(x_i)=y_i$$$

Why useful?

It helps you AC problems that you can't with $$$O(n^2)$$$.

It also gives you contribution on CF


Idea 1:

Since the pre-req says polynomial modulo, let's see if we can mod something to calculate some of $$$f(x_1), f(x_2), \cdots$$$ quickly.

It can be seen that we can split the set of $$$x$$$ into two halves — $$$X_0={x_1, x_2, \cdots, x_{\lfloor\frac{n}{2}\rfloor}}$$$ and the rest $$$X_1={x_{\lfloor\frac{n}{2}\rfloor+1}, x_{\lfloor\frac{n}{2}\rfloor+2}, \cdots, x_n}$$$

Now the crucial observation is that if we consider the polynomial $$$g_0(x)=\prod_{i=1}^{\lfloor\frac{n}{2}\rfloor} (x-x_i)$$$, then $$$\forall x\in X_0, g_0(x)=0$$$.

We can then see that polynomial mod helps us, as writing

$$$f(x)=g_0(x)\cdot Q(x)+f_0(x)$$$
$$$f(x)\equiv f_0(x)\mod g_0(x)$$$

Moreover, $$$\forall x\in X_0, f(x)=f_0(x)$$$ and the RHS can be calculated by polynomial mod. It is similar for $$$X_1$$$.

Now just recur and time complexity is $$$T(n)=2T(\frac{n}{2})+O(n\log n)=O(n\log^2 n)$$$.


Idea 2:

Given $$$(x_1,y_1),\cdots,(x_{n+1},y_{n+1})$$$, find a $$$(n-1)$$$-degree polynomial $$$f(x)$$$.

From my first post, we know that

$$$f(x)=\sum_{i=1}^{n+1} \prod_{j\neq i} \frac{x-x_j}{x_i-x_j}y_i$$$

Now look at that denominator $$$\prod_{j\neq i}x_i - x_j$$$. We can think of it as the whole product divided by $$$(x_i-x_i)$$$. We can't divide by zero but limits are cool:

$$$\prod_{j\neq i}(x_i-x_j)=\lim_{x\to x_i}\frac{\prod_{j=1}^{n+1} (x-x_j)}{x-x_i}=\lim_{x\to x_i}\frac{d}{dx}(\prod_{i=1}^{n+1} x-x_j)$$$

We can write $$$P(x)=\prod_{i=1}^{n+1} x-x_j$$$, which can be calculated in $$$O(n\log^2 n)$$$ via divide and conquer. Then $$$P'(x)$$$ can also be found.

Now going back to the original expression, $$$f(x)=\sum_{i=1}^{n+1} \frac{y_i}{M'(x_i)}\prod_{j\neq i}(x-x_j)$$$. Simplifying a bit we shall write $$$v_i=\frac{y_i}{M'(x_i)}$$$ and that $$$f(x)=\sum_{i=1}^{n+1} v_i\prod_{j\neq i}(x-x_j)$$$. I swear this is not becoming a graph problem.

Now let's consider what happens if we try to divide and conquer on this term. By divide and conquering we will get the terms

$$$f_0(x)=\sum_{i=1}^{\lfloor\frac{n+1}{2}\rfloor} v_i\prod_{j\neq i,1\leq j\leq \lfloor\frac{n+1}{2}\rfloor}(x-x_j)$$$

,

$$$f_1(x)=\sum_{i=\lfloor\frac{n+1}{2}+1\rfloor}^{n+1} v_i\prod_{j\neq i,\lfloor\frac{n+1}{2}\rfloor\leq j\leq n+1}(x-x_j)$$$

Comparing those with the $$$f(x)$$$ we want we can notice that the product part in $$$f_0(x)$$$ is missing the later half's terms, so we can "put them back" via multiplying $$$f_0(x)$$$ by $$$M_1(x)=\prod_{\lfloor\frac{n+1}{2}\rfloor\leq j\leq n+1}(x-x_j)$$$. Doing the same for $$$f_1(x)$$$, we get that

$$$f(x)=f_0(x)M_1(x)+f_1(x)M_0(x)$$$

which can then be calculated by divide and conquer.

Also note that $$$M_0(x)$$$ and $$$M_1(x)$$$ is calculated in the process of calculating the original $$$M(x)=\prod (x-x_i)$$$, so those don't require extra time.

Time complexity is $$$O(n\log^2 n)$$$ for similar reasons.


Thank you for reading the post. I'm too lazy to code it up but I still hope this blog helps someone.

I might start working through the "maths blog post recommendation" topics and post a few more blogs~

HKOI

Full text and comments »

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

By grhkm, history, 4 years ago, In English

Hello as title, any math topic suggestions that have few resources on, and are not too hard (i.e. it's explainable in one blog, perhaps with high-school maths knowledge)? IMO my previous blogs on maths (lagrange interpolation and linear modulo inversions) are pretty pog and received decent feedback, and I think (not to boast or anything) I am a good explainer? So I want to contribute with my grade 10 maths knowledge and also practice my english and writing skills in general.

It's free too for the lazy people out there, so you save your time and get your credit while I get my contribution points. Win-win condition (win for the cf readers too).

tldr please read blog title and reply thx

Full text and comments »

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

By grhkm, history, 4 years ago, In English

I am just wondering, how does your place's (country?) OI team selection work?

I am from Hong Kong and to join the HKOI selection process, we must register through our school. Unfortunately, my teacher did not confirm through the three confirmation emails and thus did not sign me up for the heat event, therefore I miss an annual chance to join the HKOI. Therefore, I am wondering if it is the same anywhere else, that you must register through your school? That sounds stupid to me but HKOI (hi tony) simply replied with "rules are rules".

In addition, I have joined the mirror of China OI first round selection and did decent. I asked HKOI if they can make an exception for me, especially since it wasn't even my fault (again, the teacher did not check the confirmation emails), but they won't let me join the second round of HKOI with those results.

Hope this makes sense, since I am probably leaving out details of the HKOI selection process.

Full text and comments »

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

By grhkm, history, 4 years ago, In English

PROBLEM: (Source)

Given $$$n$$$ $$$d$$$-dimensional vectors $$$v_1, v_2, \cdots, v_n$$$, output $$$i > j$$$ such that $$$v_i \cdot v_j \equiv 0 \mod k$$$, or determine that such indexes do not exist.

Constraints: For 50% cases $$$n \leq 10^5, d \leq 30, 2 \leq k \leq 3$$$, for all others $$$n \leq 10^4, d \leq 100, 2 \leq k \leq 3$$$.


$$$k = 2$$$

The main idea is check for every $$$i$$$ if there exists a $$$j$$$ satisfying the requirement. Keep in mind that here, $$$i$$$ can be treated as a constant. Therefore, requirement is equivalent to negating

$$$\forall 1\leq j\leq i-1, v_i \cdot v_j \not\equiv 0\mod 2\\\\$$$
$$$\forall 1\leq j\leq i-1, \sum_{k=1}^d v_{ik}v_{jk} \equiv 1\mod 2\\\\$$$
$$$\sum_{j=1}^{i-1}\sum_{k=1}^d v_{ik}v_{jk} \equiv i - 1\mod 2\\\\$$$
$$$\sum_{k=1}^d v_{ik}\sum_{j=1}^{i-1}v_{jk} \equiv i - 1\mod 2\\\\$$$

The inner summation can be computed quickly using prefix sum. Therefore complexity is O($$$nd$$$) where the $$$n$$$ comes from looping through every $$$i$$$.


$$$k = 3$$$

Now we deal with $$$k = 3$$$. Same idea, but now we don't have the nice property of $$$\not\equiv 0 \implies \equiv 1$$$. However, we can do something clever:

$$$\forall 1\leq j\leq i-1, v_i \cdot v_j \not\equiv 0\mod 3\\\\$$$
$$$\forall 1\leq j\leq i-1, (\sum_{k=1}^d v_{ik}v_{jk})^2 \equiv 1\mod 3\\\\$$$
$$$\sum_{j=1}^{i-1}\Big(\sum_{k=1}^d v_{ik}v_{jk}\Big)^2 \equiv i - 1\mod 3\\\\$$$
$$$\sum_{j=1}^{i-1}\sum_{k=1}^d \sum_{l=1}^d v_{ik}v_{jk}v_{il}v_{jl} \equiv i - 1\mod 3\\\\$$$
$$$\sum_{k=1}^d \sum_{l=1}^d v_{ik}v_{il} \Big(\sum_{j=1}^{i-1}v_{jk}v_{jl}\Big) \equiv i - 1\mod 3\\\\$$$

The inner summation can again be stored from previous queries. Therefore the time complexity is $$$O(nd^2)$$$.


Details:

1) The logic above implies that if $$$\sum \not\equiv i-1 \mod 3$$$, there will be a pair of vectors with inner product divisible by $$$3$$$. HOWEVER, it is possible that such pair exists but still bypasses all tests. For example taking $$$d = 1, k = 2, v_1 = 0, v_i = 1 \forall 2\leq i\leq n$$$, the vectors will satisfy that test above. Therefore, it is necessary to shuffle the vectors before applying the check.

2) In the original problem, we also have to restore the index. If we know a given $$$i$$$ fails, we can simply check all $$$j<i$$$ in $$$O(nd)$$$ time.

My implementation can be found here

Full text and comments »

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

By grhkm, history, 4 years ago, In English

I am messing around in godbolt.org (compiler), and I tried to compile this code with -O2

unsigned int func(unsigned int x) {
    return x / 7;
}

Surprisingly, this is what I get

func:
        mov     eax, edi             # edi is argument (x). Move x to eax (register)
        imul    rax, rax, 613566757  # eax *= 613566757 (unsigned 64-bit multiplication, rax is 64-bit version of eax)
        shr     rax, 32              # eax >>= 32 (gets upper 32 bit of product)
        sub     edi, eax             # x -= eax
        shr     edi                  # x >>= 1
        add     eax, edi             # eax += x
        shr     eax, 2               # x >>= 2
        ret                          # return x

What is imul rax, rax, 613566757 doing? Where does the constant come from?

Also, this translated code works as intended:

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

int main() {
        int i = 137;

        long long edi = i;
        long long eax = i;
        eax *= 613566757LL;
        eax >>= 32;
        edi -= eax;
        edi >>= 1;
        eax += edi;
        eax >>= 2;
        cout << i << " " << i / 7 << " " << eax << endl;
}

Thank you very much.

(Originally I had int instead of unsigned int, but it's more complicated)

Full text and comments »

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

By grhkm, history, 4 years ago, In English
Thank you

Hello!

We're going to learn how to find inverses mod p today (efficiently).

Pre-req: Know how to find inverses

Main results:

As we know, finding the inverse of n numbers is $$$O(n\log{p})$$$. That is too slow, especially when time limit is tight. Therefore, we want a faster way. I present:

  • Find inverse of all numbers between 1 and n in $$$O(n)$$$

  • Find inverse of n numbers in $$$O(n + \log{p})$$$

Idea 1:

Notice that $$$p = i \cdot \lfloor \frac{p}{i} \rfloor + (p \% i)$$$. (If you don't know why, compare it to $$$13=3\cdot 4+1$$$)

$$$p = i \cdot \lfloor \frac{p}{i} \rfloor + (p \% i)$$$
$$$\implies i \cdot \lfloor \frac{p}{i} \rfloor + (p \% i) \equiv 0 \mod p$$$
$$$p \% i \equiv -i \cdot \lfloor \frac{p}{i} \rfloor \mod p$$$

Now NOTICE that $$$p \% i$$$ is actually less than $$$i$$$. (Big brain mode)

So if we've calculated the inverses of $$$1-(i-1)$$$ already, we can find inverse of $$$(p \% i)$$$. This is what we're going to use.

$$$1 \equiv i \cdot (-\lfloor\frac{p}{i}\rfloor \cdot (p\% i)^{-1}) \mod p$$$

And therefore, by definition of inverse:

$$$i^{-1} \equiv (-\lfloor\frac{p}{i}\rfloor \cdot (p\% i)^{-1}) \mod p$$$

Or to make it easier to read: (Integer division)

$$$inv[i] \equiv (-(p/i) \cdot inv[p\% i]) \mod p$$$

So we can calculate inverse of $$$1~n$$$ in $$$O(n)$$$.

Code: (Be aware of overflow) (Also, -(p/i) mod p = (p-p/i) mod p)

int n = 10, p = 1000000007;
int inv[n + 1];
inv[1] = 1;
for (int i = 2; i <= n; i ++) inv[i] = 1LL * (p - p / i) * inv[p % i] % p;

Idea 2:

Our target is $$$O(n+\log p)$$$, so we shall only take inverse once. How?

Let's look at what we want to find:

$$$\frac{1}{a_1}, \frac{1}{a_2}, \cdots, \frac{1}{a_n} \mod p$$$

Let's try to make them have the same denominator!

$$$\frac{a_2a_3\cdots a_n}{a_1a_2\cdots a_n}, \frac{a_1a_3a_4\cdots a_n}{a_1a_2\cdots a_n}, \cdots, \frac{a_1a_2\cdots a_{n-1}}{a_1a_2\cdots a_n} \mod p$$$

As we can see, the numerator is a prefix times a suffix, and the denominator is constant. Therefore, we can precompute those and then take inverse of denominator once.

Code:

#include "bits/stdc++.h"
using namespace std;

const long long P = 1000000007;
long long qpow(long long a, long long b) {
	long long ans = 1;
	while (b) {
		if (b & 1) ans = ans * a % P;
		a = a * a % P;
		b >>= 1;
	}
	return ans;
}

long long n, a[1000010], pre[1000010], suf[1000010], pr = 1; // prefix, suffix, product
int main() {
	cin >> n; pre[0] = suf[n + 1] = 1;
	for (int i = 1; i <= n; i ++) cin >> a[i];
	for (int i = 1; i <= n; i ++) {
		pre[i] = pre[i - 1] * a[i] % P;
		suf[n + 1 - i] = suf[n + 2 - i] * a[n + 1 - i] % P;
		pr = pr * a[i] % P;
	}
	pr = qpow(pr, P - 2);
	for (int i = 1; i <= n; i ++) {
		cout << (pre[i - 1] * suf[i + 1] % P) * pr % P << endl;
	}
}

Thank you for reading. If you have any suggestions or any topic you want to learn about I guess I can try writing! If this is helpful vote pls thx

Спойлер

Full text and comments »

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

By grhkm, history, 4 years ago, In English

Lagrange Interpolation by grhkm

No polynomial stuff because I don't know $$$O(n\log{n})$$$ polynomial multiplication :(

Main results:

Given $$$f(x_1)=y_1, f(x_2)=y_2, ..., f(x_n)=y_n$$$ we can calculate $$$f(X)$$$ where $$$f$$$ is the unique $$$(n-1)$$$-degree polynomial

For general {$$$x_i$$$} we can do so in $$$O(n^2)$$$

For consecutive {$$$x_i$$$} i.e. $$$x_j=x_1+(j-1)$$$, we can do so in $$$O(n)$$$ excluding binary exponentiation

Why useful?

Calculate $$$1^n+2^n+...+X^n$$$ in $$$O(n)$$$

DP optimization

With polynomials you can do cool things but not gonna cover this blog

Let's get started :)

Assume all polynomials below are $$$(n-1)$$$ degree unless specified.

disclaimer
Why (n-1) degree?

Idea:

To start let's deal with first problem:

Given $$$f(x_1)=y_1, f(x_2)=y_2, ..., f(x_n)=y_n$$$ and an integer $$$X$$$ we can calculate $$$f(X)$$$ where $$$f$$$ is the unique $$$(n-1)$$$-degree polynomial

Let's consider an easier version:

Find polynomial satisfying $$$f(x_1)=y_i, f(x_2)=0, ..., f(x_n)=0$$$

As we learn in school, $$$f(r)=0 \implies (x-r)$$$ is divisor

We can write this down:

$$$f(x)=(x-x_2)(x-x_3)\cdots (x-x_n)$$$

Now notice that $$$f(x_1)=\prod_{i=2}^n (x_1-x_i)$$$ (a constant), but we want it to be $$$y_1$$$

We can simply scale the polynomial by a correct factor i.e.

$$$f_1(x)= \prod_{i=2}^n (x-x_i) \cdot \frac{y_1}{\prod_{i=2}^n (x_1-x_i)}=y_1\prod_{i=2}^n \frac{x-x_i}{x_1-x_i}$$$

Similarly, we can find a similar polynomial such that $$$f_i(x_i)=y_i$$$ and $$$f_i(x_j)=0 \forall j \neq i$$$

$$$f_i(x) = y_i \prod_{j=1,j\neq i}^n \frac{x-x_i}{x_i-x_i}$$$

Okay now finally to answer the original question, we can notice that $$$f(x)=\sum_{i=1}^n f_i(x)$$$ satisfies the constraints. So there we have it! Lagrange interpolation:

$$$f(x)=\sum_{i=1}^n y_i \prod_{j=1,j\neq i}^n \frac{x-x_j}{x_i-x_j}$$$

$$$O(n^2)$$$ excluding inverse for division

Optimization:

Let's try to optimise it now! Assume that we're given $$$f(1), f(2), \cdots, f(n)$$$, how can we calculate $$$f(X)$$$ quickly?

Notice that here, $$$x_i=i$$$ and $$$y_i=f(i)$$$. Substitute!

$$$f(x) = \sum_{i=1}^n f(i) \prod_{j=1,j\neq i}^n \frac{x-j}{i-j} = \sum_{i=1}^n f(i) \frac{\prod_{j=1,j\neq i}^n x-j}{\prod_{j=1,j\neq i}^n i-j}$$$

Let's consider the numerator and denominator separately

Numerator

$$$\prod_{j=1,j\neq i}^n x-j = [(x-1)(x-2)\cdots (x-(i-1))] \cdot [(x-(i+1))(x-(i+2))\cdots (x-n)]$$$

We can pre-calculate prefix and suffix product of x, then we can calculate the numerator (for each $$$i$$$) in $$$O(1)$$$

Denominator

$$$\prod_{j=1,j\neq i}^n i-j = [(i-1)(i-2)(i-3)\cdots(i-(i-1))] \cdot [(i-(i+1))(i-(i+2))\cdots (i-n)]$$$
$$$=(i-1)! \cdot (-1)^{n-(i+1)+1} \cdot (n-i)!$$$
$$$=(-1)^{n-i} (n-i)! (i-1)!$$$

So we can precompute factorials (preferably their inverse directly) then $$$O(1)$$$ calculation

So now we can calculate $$$f(X)$$$ in $$$O(n)$$$

Example 0

Find_ $$$\sum_{i=1}^N i^k$$$, $$$k\leq 1e6, N\leq 1e12$$$

Solution: Notice that the required sum will be a degree $$$(k+1)$$$ polynomial, and thus we can interpolate the answer with $$$(k+2)$$$ data points (Remember, we need $$$deg(f)+1$$$ points)

To find the data points simply calculate $$$f(0)=0, f(x)=f(x-1)+x^k$$$

Code can be found below

Example 1

(BZOJ 3453, judge dead): Find_ $$$\sum_{i=0}^n \sum_{j=1}^{a+id} \sum_{k=1}^j k^x \mod 1234567891, x \leq 123, a, n, d \leq 123456789$$$

Hint

Example 2 (ADVANCED)

(Luogu P4463 submit here): Given $$$n \leq 500, k \leq 1e9$$$, we call a sequence $$${a_n}$$$ nice if $$$1\leq a_i \leq k \forall i$$$ and $$$(a_i\neq a_j) \forall i\neq j$$$. Find the sum of (products of elements in the sequence) over all nice sequences. MY code below

Hint
Hint2
Hint3
Hint4

Example 3 (ADVANCED TOO)

Problem https://www.codechef.com/problems/XETF : Find $$$\sum_{i=1,gcd(i,N)=1}^N i^k, N \leq 1e12, k \leq 256$$$. My code below

Hint

If you have any questions feel free to ask below

Practice question:

CF 622F (Example 0)

Luogu P4463 (Example 2)

Codechef XETF (Example 3)

Not in order of difficulty:

SPOJ ASUMEXTR

SPOJ ASUMHARD which is harder than EXTREME

Codechef COUNTIT

Atcoder ARC33D (Japanese)

Atcoder S8PC3G

  • Unfortunately you can't get full score, but you can get up to subtask 4

  • Still not trivial, do first 3 subtasks would give insight

  • For full solution please wait for my next blog (next year)

HDU 4059

SPOJ SOMESUMS very interesting

PE 487 Direct application

Full text and comments »

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