# BIT 與 Polynomial Queries 的關聯性
## 前言
這篇文章的目的是要解釋如何使用Fenwick Tree（下稱BIT）來處理在線的Polynomial Queries（如下題），以及它的推導過程。

## 先備條件
1. 熟悉前綴和、差分、BIT的使用方法。
2. 不害怕數學抽象符號。
## 題目
來源：[zerojudge f992](https://zerojudge.tw/ShowProblem?problemid=f992)

簡易題敘：
> 給定 $n, k, q$ 與初始陣列 $a_1～a_n$ ，請執行 $q$ 次在線操作，操作皆為以下兩種之一：
> 1.給定 $l, r$ ，對於所有的 $l\leq i\leq r$ ，使 $a_i$ 增加 ${(i-l+1)}^k$ 。
> $~~~$（ $a_l$ 增加 $1^k$ 、 $a_{l+1}$ 增加 $2^k$ 、…、 $a_r$ 增加 ${(r-l+1)}^k$ ）
> 2.給定 $l, r$ ，輸出 $[l, r]$ 的區間和（對質數 ${10}^6+3$ 取模）。
> 值域：$1\leq n, q\leq{10}^5, 0 \leq k\leq 3, 1\leq a_i\leq{10}^9$
>

## 數學定義
為了方便，我們先假設一些代號。

定義 $A(0)$ 為原陣列， $A(0)_i$ 為原陣列的第 $i$ 項。

定義 $A(x)$ 為原陣列做 $x$ 次前綴和操作之後，所得到的新陣列，
而 $A(-x)$ 則為原陣列做 $x$ 次差分操作之後，所得到的新陣列， $x$ 為正整數。

對於所有的整數 $x$ 、非負整數 $i$ ，定義 $A(x)_0=0$ 、 $A(x)_{i+1}=A(x)_{i}+A(x-1)_{i+1}$ 。

舉例：

| $i$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|:-------:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:--:|
|$A(1)_i$ | 0 | 1 | 3 | 6 | 10| 15| 21| 28| 36| 36| 36 |
|$A(0)_i$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 0 | 0 |
|$A(-1)_i$| 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | -8| 0 |

## 觀察
$a_i$ 的初始值很好處理，為了方便觀察，不失一般性假設 $a_1～a_n$ 的初始值為 $0$ 。

假設 $n=10, k=3$ ，然後執行一次操作1，區間 $[l, r]$ 為 $[1, 5]$ 。

則結果為：
| $i$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|:-------:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:--:|
|$A(0)_i$ | 0 | 1 | 8 | 27| 64|125| 0 | 0 | 0 | 0 | 0 |
|$A(-1)_i$| 0 | 1 | 7 | 19| 37| 61|-125| 0| 0 | 0 | 0 |
|$A(-2)_i$| 0 | 1 | 6 | 12| 18| 24|-186|125| 0| 0 | 0 |
|$A(-3)_i$| 0 | 1 | 5 | 6 | 6 | 6 |-210|311|-125|0| 0 |
|$A(-4)_i$| 0 | 1 | 4 | 1 | 0 | 0 |-216|521|-436|125| 0|

事實上，當 $k=3$ 時，不論 $n$ 有多大，每次執行操作1最多只會改變 $A(-4)$ 陣列中的7項。

更進一步的說，對於任意的 $k$ 值，不論 $n$ 有多大，每次執行操作1最多只會改變 $A(-k-1)$ 陣列中的 $O(k)$ 項。

這個結論除了可以用來解決此題的離線版本：「所有操作2皆發生在操作1之後」，

同時也為我們找到了「如何處理更大的 $k$ 」這個問題的思考方向。

## 推導
我們先解決 $k=0$ 的情況：

$A(0)_1+A(0)_2+A(0)_3+\dots+A(0)_{x-1}+A(0)_x$

$=\displaystyle\sum_{i=1}^{1}{A(-1)_i} +\displaystyle\sum_{i=1}^{2}{A(-1)_i}+\displaystyle\sum_{i=1}^{3}{A(-1)_i}+\dots+\displaystyle\sum_{i=1}^{x-1}{A(-1)_i}+\displaystyle\sum_{i=1}^{x}{A(-1)_i}$

$=\displaystyle\sum_{i=1}^{x}{\left(A(-1)_i\times((x+1)-i)\right)}\cdots\huge①$

$=(x+1)\displaystyle\sum_{i=1}^{x}{A(-1)_i}-\displaystyle\sum_{i=1}^{x}{(A(-1)_i\times i)}$

假設 $B_i=A(-1)_i\times i$ ，那我們只要用兩個BIT分別維護 $A(-1),B$ 兩陣列的值，就可以解決 $k=0$ 情況下的問題了。

根據上個章節的結論，只要把上式的 $A(-1)$ 想辦法代換成 $A(-2)$ ，那 $k=1$ 情況下的問題就解決了。

以下為代換過程：

$A(0)_1+A(0)_2+A(0)_3+\dots+A(0)_{x-1}+A(0)_x$

$=\displaystyle\sum_{i=1}^{x}{\left(A(-1)_i\times\left(\left(x+1\right)-i\right)\right)}\cdots\huge①$

$=\displaystyle\sum_{i=1}^{x}{\left(\left(\displaystyle\sum_{j=1}^{i}{A(-2)_j}\right)\times\left(\left(x+1\right)-i\right)\right)}$

$=\displaystyle x\left(\sum_{i=1}^{1}{A(-2)_i}\right)+\left(x-1\right)\left(\displaystyle\sum_{i=1}^{2}{A(-2)_i}\right)+\left(x-2\right)\left(\displaystyle\sum_{i=1}^{3}{A(-2)_i}\right)+\dots+2\left(\displaystyle\sum_{i=1}^{x-1}{A(-2)_i}\right)+1\left(\displaystyle\sum_{i=1}^{x}{A(-2)_i}\right)$

$=\displaystyle\sum_{i=1}^{x}{\left(A(-2)_i\times\sum_{j=1}^{x-i+1}{j}\right)}$

$=\displaystyle\sum_{i=1}^{x}{\left(A(-2)_i\times\frac{\left(x-i+1\right)\left(\left(x-i+1\right)+1\right)}{2}\right)}$

$=\displaystyle\frac{1}{2}\sum_{i=1}^{x}{\left(A(-2)_i\times\left(\left(x+1\right)-i\right)\left(\left(x+2\right)-i\right)\right)}\cdots\huge②$

$=\displaystyle\frac{1}{2}\left(\left(\sum_{i=1}^{x}{A(-2)_i\times\left(x+1\right)\left(x+2\right)}\right)-\left(\sum_{i=1}^{x}{A(-2)_i\times i\times\left(\left(x+1\right)+\left(x+2\right)\right)}\right)+\left(\sum_{i=1}^{x}{A(-2)_i\times {i^2}}\right)\right)$

$=\displaystyle\frac{1}{2}\left(\left({x^2}+3x+2\right)\left(\sum_{i=1}^{x}{A(-2)_i}\right)-\left(2x+3\right)\left(\sum_{i=1}^{x}{\left(A(-2)_i\times i\right)}\right)+\left(\sum_{i=1}^{x}{\left(A(-2)_i\times {i^2}\right)}\right)\right)$

假設 $B_i=A(-2)_i\times i, C_i=A(-2)_i\times {i^2}$ ，那我們只要改成用三個BIT分別維護 $A(-2), B, C$ 三陣列的值，就可以解決 $k=1$ 情況下的問題了。

接下來我們處理 $k=2$ 的情況，也就是把上式的 $A(-2)$ 代換成 $A(-3)$ ：

$A(0)_1+A(0)_2+A(0)_3+\dots+A(0)_{x-1}+A(0)_x$

$=\displaystyle\frac{1}{2}\sum_{i=1}^{x}{\left(A(-2)_i\times\left(\left(x+1\right)-i\right)\left(\left(x+2\right)-i\right)\right)}\cdots\huge②$

$=\displaystyle\frac{1}{2}\sum_{i=1}^{x}{\left(\left(\displaystyle\sum_{j=1}^{i}{A(-3)_j}\right)\times\left(\left(x+1\right)-i\right)\left(\left(x+2\right)-i\right)\right)}$

$=\displaystyle\frac{1}{2}\left(A(-3)_1\times x\left( x+1 \right)+\left(A(-3)_1+A(-3)_2\right)\times \left( x-1 \right)\left(x\right)+\left(A(-3)_1+A(-3)_2+A(-3)_3\right)\times \left( x-2 \right)\left( x-1 \right)+\dots+\left(\displaystyle\sum_{i=1}^{n-1}{A(-3)_i}\right)\times \left(2 \times 3\right)+\left(\displaystyle\sum_{i=1}^{n}{A(-3)_i}\right)\times \left(1 \times 2\right)\right)$

$=\displaystyle\frac{1}{2}\sum_{i=1}^{x}{\left(A(-3)_i\times\sum_{j=1}^{x-i+1}{j \left( j+1 \right)}\right)}$

$=\displaystyle\frac{1}{2}\sum_{i=1}^{x}{\left(A(-3)_i\times\frac{\left( x-i+1 \right)\left( x-i+2 \right)\left( x-i+3 \right)}{3}\right)}$

$=\displaystyle\frac{1}{6}\sum_{i=1}^{x}{\left(A(-3)_i\times\left( \left( x+1 \right)-i \right)\left( \left( x+2 \right)-i \right)\left( \left( x+3 \right)-i \right)\right)}\cdots\huge③$

$=\displaystyle\frac{1}{6}\left(\left({x^3}+6{x^2}+11x+6\right)\left(\sum_{i=1}^{x}{A(-3)_i}\right)-\left(3x^2+12x+11\right)\left(\sum_{i=1}^{x}{\left(A(-3)_i\times i\right)}\right)+\left(3x+6\right)\left(\sum_{i=1}^{x}{\left(A(-3)_i\times {i^2}\right)}\right)-\left(\sum_{i=1}^{x}{\left(A(-3)_i\times {i^3}\right)}\right)\right)$

假設 $B_i=A(-3)_i\times i, C_i=A(-3)_i\times {i^2}, D_i=A(-3)_i\times {i^3}$ ，則我們能用四個BIT分別維護 $A(-3), B, C, D$ 四陣列的值來解決 $k=2$ 情況下的問題。

相信看到這裡，有人應該已經知道代換成 $A(-4)$ 的式子是什麼了，不過這個部分我們會留到下一章節，在這裡我們還是先單獨證明 $A(-4)$ 的版本：

$A(0)_1+A(0)_2+A(0)_3+\dots+A(0)_{x-1}+A(0)_x$

$=\displaystyle\frac{1}{6}\sum_{i=1}^{x}{\left(A(-3)_i\times\left( \left( x+1 \right)-i \right)\left( \left( x+2 \right)-i \right)\left( \left( x+3 \right)-i \right)\right)}\cdots\huge③$

$=\displaystyle\frac{1}{6}\sum_{i=1}^{x}{\left(\left(\displaystyle\sum_{j=1}^{i}{A(-4)_j}\right)\times\left( \left( x+1 \right)-i \right)\left( \left( x+2 \right)-i \right)\left( \left( x+3 \right)-i \right)\right)}$

$=\displaystyle\frac{1}{6}\left(A(-4)_1\times x\left( x+1 \right)\left( x+2 \right)+\left(A(-4)_1+A(-4)_2\right)\times \left( x-1 \right)\left(x\right)\left( x+1 \right)+\left(A(-4)_1+A(-4)_2+A(-4)_3\right)\times \left( x-2 \right)\left( x-1 \right)\left( x \right)+\dots+\left(\displaystyle\sum_{i=1}^{n-1}{A(-4)_i}\right)\times \left(2 \times 3\times 4\right)+\left(\displaystyle\sum_{i=1}^{n}{A(-4)_i}\right)\times \left(1 \times 2\times 3\right)\right)$

$=\displaystyle\frac{1}{6}\sum_{i=1}^{x}{\left(A(-4)_i\times\sum_{j=1}^{x-i+1}{j \left( j+1 \right)\left( j+2 \right)}\right)}$

$=\displaystyle\frac{1}{6}\sum_{i=1}^{x}{\left(A(-4)_i\times\frac{\left( x-i+1 \right)\left( x-i+2 \right)\left( x-i+3 \right)\left( x-i+4 \right)}{4}\right)}$

$=\displaystyle\frac{1}{24}\sum_{i=1}^{x}{\left(A(-4)_i\times\left( \left( x+1 \right)-i \right)\left( \left( x+2 \right)-i \right)\left( \left( x+3 \right)-i \right)\left( \left( x+4 \right)-i \right)\right)}\cdots\huge④$

$=\displaystyle\frac{1}{24}\left(\left({x^4}+10{x^3}+35{x^2}+50x+24\right)\left(\sum_{i=1}^{x}{A(-4)_i}\right)-\left(4{x^3}+30{x^2}+70{x}+50\right)\left(\sum_{i=1}^{x}{\left(A(-4)_i\times i\right)}\right)+\left(6{x^2}+30x+35\right)\left(\sum_{i=1}^{x}{\left(A(-4)_i\times {i^2}\right)}\right)-\left(4x+10\right)\left(\sum_{i=1}^{x}{\left(A(-4)_i\times {i^3}\right)}\right)+\left(\sum_{i=1}^{x}{\left(A(-4)_i\times {i^4}\right)}\right)\right)$

因此我們知道，總共用五棵BIT可以完成這題的所有限制： $0\leq k\leq 3$ 。

## 一般式

我想在這個章節之前，有人已經能猜到一般式是什麼了。

這邊公布答案：

$A(0)_1+A(0)_2+A(0)_3+\dots+A(0)_{x-1}+A(0)_x$

$=\displaystyle\frac{1}{k!}\sum_{i=1}^{x}{\left(A(-k)_i\prod_{j=1}^{k}{\left(x+j-i\right)}\right)}$

這個式子可以用數學歸納法來證明，具體來說要利用這個式子：

$\displaystyle\sum_{i=1}^{n}{\left(\prod_{j=0}^{m}{\left(i+j\right)}\right)}=\frac{1}{m+2}\prod_{j=0}^{m+1}{\left(n+j\right)}$

詳細過程請讀者自行證明，這裡不多做解釋。

## 操作的等價

我們知道，當 $k\leq 3$ 的時候，每次執行操作1都只會改變 $A(-4)$ 中的 $O(k)$ 項，但具體來說，這幾項到底會改變多少呢？

| $i$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|:-------:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:--:|
|$A(0)_i$ | 0 | 1 |$8$| 27| 64|125| 0 | 0 | 0 | 0 | 0 |
|$A(-1)_i$| 0 | 1 | 7 | 19| 37| 61|-125| 0| 0 | 0 | 0 |
|$A(-2)_i$| 0 | 1 | 6 | 12| 18| 24|-186|125| 0| 0 | 0 |
|$A(-3)_i$| 0 | 1 | 5 | 6 | 6 | 6 |-210|311|-125|0| 0 |
|$A(-4)_i$| 0 | 1 | 4 | 1 | 0 | 0 |-216|521|-436|125| 0|

如果上面這個例子的操作區間 $[l, r]$ 不是 $[1, 5]$ ，那麼 $1, 4, 1, -216, 521, -436, 125$ 這七個數字會變成什麼呢？

經過觀察（建議讀者多試幾組不同的位置與區間大小）之後可以發現，這七個數字中 $1,4,1$ 是固定的，但 $-216, 521, -436, 125$ 這四個數字的值則是跟區間大小 $(r-l+1)$ 有關。

至於有什麼關聯呢？我們可以利用代數方法來計算。

令區間大小 $m=r-l+1$ ，

假設 $n=1000,k=3$ ，操作區間$[l, r]=[1, m]$ ，且$4\leq m\leq n-4$ ：

| $i$ |$0$|$1$|$2$| $3$| $4$|$\dots$|$m-1$| $m$ |$m+1$|$m+2$|$m+3$|$m+4$|$m+5$|
|:-------:|:-:|:-:|:-:|:--:|:--:|:-----:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
|$A(0)_i$ |$0$|$1$|$8$|$27$|$64$|$\dots$|$(m-1)^3$|$m^3$|$0$|$0$|$0$|$0$|$0$|
|$A(-1)_i$|$0$|$1$|$7$|$19$|$37$|$\dots$|$3m^2-9m+7$|$3m^2-3m+1$|$-m^3$|$0$|$0$|$0$|$0$|
|$A(-2)_i$|$0$|$1$|$6$|$12$|$18$|$\dots$|$6m-12$|$6m-6$|$-m^3-3m^2+3m-1$|$m^3$|$0$|$0$|$0$|
|$A(-3)_i$|$0$|$1$|$5$| $6$| $6$|$\dots$|$6$|$6$|$-m^3-3m^2-3m+5$|$2m^3+3m^2-3m+1$|$-m^3$|$0$|$0$|
|$A(-4)_i$|$0$|$1$|$4$| $1$| $0$|$\dots$|$0$|$0$|$-m^3-3m^2-3m-1$|$3m^3+6m^2-4$|$-3m^3-3m^2+3m-1$|$m^3$|$0$|

由此可以得知：

$-216, 521, -436, 125$

分別是：

$(-m^3-3m^2-3m-1),(3m^3+6m^2-4),(-3m^3-3m^2+3m-1),(m^3)$

代入 $m=5$ 的結果。

至於 $k=0, 1, 2$ 的操作對 $A(-4)$ 的變化，就交給讀者自行證明了，在下個章節的code也能找到答案。

## 實作方法

<!-- 說明待補 -->
:::spoiler code(C++)

#include <bits/stdc++.h>
#define loop(i,a,b) for(int i=a;i<b;i++)
using namespace std;
typedef long long ll;

const int mod = (int)(1e6) + 3;
const int inv_24 = 208334;
const int mxn = (int)(1e5) + 10;
int pref[mxn];
ll i0_BIT[mxn], i1_BIT[mxn], i2_BIT[mxn], i3_BIT[mxn], i4_BIT[mxn];

inline ll neg_mod(ll x){ return ((x % mod) + mod) % mod; }
inline int lowbit(int x){ return (x) & (-x); }
inline void upd(int ind, ll val){
    for(int i = ind; i < mxn; i += lowbit(i))
        (i0_BIT[i] += val) %= mod;

    val = val * ind % mod;
    for(int i = ind; i < mxn; i += lowbit(i))
        (i1_BIT[i] += val) %= mod;

    val = val * ind % mod;
    for(int i = ind; i < mxn; i += lowbit(i))
        (i2_BIT[i] += val) %= mod;

    val = val * ind % mod;
    for(int i = ind; i < mxn; i += lowbit(i))
        (i3_BIT[i] += val) %= mod;

    val = val * ind % mod;
    for(int i = ind; i < mxn; i += lowbit(i))
        (i4_BIT[i] += val) %= mod;
}
inline ll qry(ll ind){
    int ret0 = 0, ret1 = 0, ret2 = 0, ret3 = 0, ret4 = 0;
    for(int i = ind; i > 0; i -= lowbit(i)){
        ret0 = (ret0 + i0_BIT[i] >= mod) ? (ret0 + i0_BIT[i] - mod) : (ret0 + i0_BIT[i]);
        ret1 = (ret1 + i1_BIT[i] >= mod) ? (ret1 + i1_BIT[i] - mod) : (ret1 + i1_BIT[i]);
        ret2 = (ret2 + i2_BIT[i] >= mod) ? (ret2 + i2_BIT[i] - mod) : (ret2 + i2_BIT[i]);
        ret3 = (ret3 + i3_BIT[i] >= mod) ? (ret3 + i3_BIT[i] - mod) : (ret3 + i3_BIT[i]);
        ret4 = (ret4 + i4_BIT[i] >= mod) ? (ret4 + i4_BIT[i] - mod) : (ret4 + i4_BIT[i]);
    }
    ll ind2 = ind * ind % mod;
    ll ind3 = ind2 * ind % mod;

    ll mul0 = (ind + 1) * (ind + 2) % mod * (ind + 3) * (ind + 4) % mod;
    ll mul1 = mod - ((4 * ind3 + 30 * ind2 + 70 * ind + 50) % mod);
    ll mul2 = (6 * ind2 + 30 * ind + 35) % mod;
    ll mul3 = mod - ((4 * ind + 10) % mod);
    static int mul4 = 1;
    return (ret0 * mul0 + ret1 * mul1 + ret2 * mul2 + ret3 * mul3 + ret4 * mul4) % mod;
}

int main(){
    ios::sync_with_stdio(false); cin.tie(0);
    int n, k, q;
    cin >> n >> k >> q;
    loop(i,1,n+1){
        cin >> pref[i];
        (pref[i] += pref[i-1]) %= mod;
    }
    int cmd, l, r;
    ll len, len2, len3;
    loop(queries,0,q){
        cin >> cmd >> l >> r;
        len = r - l + 1;
        if(cmd == 1){ // change value
            if(k <= 1){
                if(k == 0){
                    upd(l, 1);
                    upd(l + 1, mod - 3);
                    upd(l + 2, 3);
                    upd(l + 3, mod - 1);
                    upd(r + 1, mod - 1);
                    upd(r + 2, 3);
                    upd(r + 3, mod - 3);
                    upd(r + 4, 1);
                }
                else{
                    upd(l, 1);
                    upd(l + 1, mod - 2);
                    upd(l + 2, 1);
                    upd(r + 1, mod - len - 1);
                    upd(r + 2, len * 3 + 2);
                    upd(r + 3, mod - len * 3 - 1);
                    upd(r + 4, len);
                }
            }
            else{
                len2 = len * len % mod ;
                if(k == 2){
                    upd(l, 1);
                    upd(l + 2, mod - 1);
                    upd(r + 1, mod - (len2 + (len << 1) + 1) % mod);
                    upd(r + 2, (3 * len2 + 4 * len) % mod);
                    upd(r + 3, (3 * (mod - len2) + 2 * (mod - len) + 1) % mod);
                    upd(r + 4, len2);
                }
                else{
                    len3 = len2 * len % mod;
                    upd(l, 1);
                    upd(l + 1, 4);
                    upd(l + 2, 1);
                    upd(r + 1, ((mod - len3) + 3 * (mod - len2) + 3 * (mod - len) + (mod - 1)) % mod);
                    upd(r + 2, (3 * len3 + 6 * len2 + (mod - 4)) % mod);
                    upd(r + 3, (3 * (mod - len3) + 3 * (mod - len2) + 3 * len + (mod - 1)) % mod);
                    upd(r + 4, len3);
                }
            }
        }
        else{ // ask value
            int ans = ((qry(r) + mod - qry(l-1)) * inv_24 + (pref[r] + mod - pref[l-1])) % mod;
            cout << ans << '\n';
        }
    }
}

:::

## 更好的複雜度

從上面的code可以發現，用 BIT 單次修改的複雜度為 $O(k^2\log n)$ ，如果用線段樹來做這題，複雜度只有 $O(\log n)$ ，而且實作更簡單，只需要知道求和公式就可以實作懶標了，至於常數的部分，BIT是否比線段樹好，我就沒有試過了。