「Lagrange 插值」学习笔记

目录

• 问题引入
  • 思考
  • Lagrange 插值法
• 插值过程
  • 代码实现
• 实际应用
  • 「洛谷 P4781」「模板」拉格朗日插值
  • 「洛谷 P4463」calc
    • 题意简述
    • 数据规模
    • Solution
      • Step 1
      • Step 2
        • 证明
      • 代码
  • 「CF 995F」Cowmpany Cowmpensation
    • 题意简述
    • 数据规模
    • Solution
      • Step 1
      • Step 2
        • 证明
      • 代码
  • 「CF 662F」The Sum of the k-th Powers
    • 题意简述
    • 数据规模
    • Solution
      • 代码
  • 「BZOJ 3453」tyvj 1858 XLkxc
    • 题意简述
    • 数据规模
    • Solution
      • 结论 1
      • 结论 2
        • 证明
      • 本题
      • 代码
  • 「CF 1086F」Forest Fires
    • 题意简述
    • 数据规模
    • Solution
      • Step 0 前置-矩形面积并
      • Step 1
      • Step 2
      • 代码

问题引入

  Nephry 最近迷上了天体物理学看星星，这是她一周观察火星运行之后记录的数据（纯属虚构）：

	Mon.	Tue.	Wed.	Thu.	Fri.
太阳-火星距离	(3 imes10^4 km)	(6 imes10^4km)	刷题去了	(5 imes10^4km)	(7 imes10^4km)

  现在 Nephry 想通过其它四组数据来推测周三的数据。

思考

  显然应该抽象成函数关系。那么我们可以建立 “Dist-Date” 图如下：

  其中绿色的轨迹：(f(x)=frac{1}2x^3-frac{14}3x^2+frac{27}2x-frac{19}3) 就是一种满足已知数据的函数关系，它也是唯一过这四个点的三次函数。当然，实际物理问题中的函数关系远远不止三次这么简单，我们只是找到其中一种情况，并期望借此来估测 (f(3)) 的值。

Lagrange 插值法

  于是，为解决诸如此类问题，法国数学家 Lagrange 提出了一种构造恰好穿过平面上若干已知点的多项式函数的方法，也即是 Lagrange 插值法。所谓插值，度娘解释道：

插值是数学领域数值分析中的通过已知的离散数据求未知数据的过程或方法。

  通俗地理解成，“猜出”系列离散数据的连续函数关系，并借此求出未知数据的方法。

插值过程

  举一个最简单的例子，已知一个二次函数 (f(x)) 过三点 ((x_1,y_1),(x_2,y_2),(x_3,y_3))，求 (f(x)) 的解析式。

  代入解方程？不！上面那个函数把兔兔算傻了qwq。 Lagrange 通过构造三条二次曲线相加的方法来解决这个问题：

  第一条曲线 (f_1(x))，满足 (f_1(x_1)=1,f_1(x_2)=f_1(x_3)=0)，如图（懒得画w，搬的课件）：

  第二条曲线 (f_2(x))，满足 (f_2(x_2)=1,f_2(x_1)=f_2(x_3)=0)，如图：

  类似地，第三条曲线 (f_3(x)) 应有 (f_3(x_3)=1,f_3(x_1)=f_3(x_2)=0)，如图：

  趁你不明所以，Lagrange 写下一个等式：

[f(x)=y_1f_1(x)+y_2f_2(x)+y_3f_3(x) ]

  代入 (x_1,x_2,x_3)，发现 (f(x)) 正满足恰好过三个点的条件！（原理很简单，不赘述w。）

  推而广之，并加以严谨的数学语言描述，假设现有 (n) 个点 ((x_1,y_1),(x_2,y_2),dots,(x_n,y_n))，要找到恰好穿过这 (n) 个点的 (n-1) 次函数 (f(x))，我们只需要构造 (n) 条 (n-1) 次曲线。其中第 (i) 条有：

[f_i(x)=egin{cases}1&x=x_i\0&xin{x_1,x_2,dots,x_n}-{x_i}end{cases} ]

  如何构造出 (i) 呢？利用小学二年级的因式分解知识，我们知道 (f_i(x)) 必然含有因式 ((x-t),~tin{x_1,x_2,dots,x_n}-{x_i}) 。接下来只需要使得 (f_i(x_i)=1)，即所有以上因式都被约掉。那么就有：

[f_i(x)=prod_{jin[1,n]-{i}}frac{x-x_j}{x_i-x_j} ]

  最后，就可以根据这 (n) 条曲线构造出：

[f(x)=sum_{i=1}^ny_if_i(x) ]

代码实现

  功能：读入 (n) 个点 ((x_i,y_i))，求出过这 (n) 个点的曲线在 (x=k) 处的取值，对大素数取模。

inline int lagrange ( const vector<int> x, const vector<int> y, const int k ) {
	int n = x.size (), ret = 0;
	for ( int i = 0; i < n; ++ i ) {
		int num = y[i], den = 1;
		for ( int j = 0; j < n; ++ j ) {
			if ( i ^ j ) {
				num = 1ll * num * ( k - x[j] + p ) % p; // 直接代入 x=k 运算.
				den = 1ll * den * ( x[i] - x[j] + p ) % p;
			}
		}
		ret = ( ret + 1ll * num * qkpow ( den, p - 2 ) ) % p; // 外层 Sigma 求和.
	}
	return ret;
}

  复杂度 (O(n^2))。

实际应用

  注意下文所述的 “(f(x)) 是关于 (x) 的 (n) 次函数” 均指其最高次不超过 (n) 次。

「洛谷 P4781」「模板」拉格朗日插值

  link.

  如上文的代码实现，不赘述。

「洛谷 P4463」calc

题意简述

  link.

  给定 (k,n,p)，求：

[sum_{{a_n}}left[(forall i)left(a_iin[1,k] ight) ight]left[(forall i ot=j)left(a_i ot=a_j ight) ight]prod_{i=1}^n a_imod p ]

数据规模

  (nle500;~n+1<k<ple10^9)。

Solution

Step 1

  有一个比较 naive 的 DP：令 (f(i,j)) 表示前 (i) 个元素取值在 ([1,j]) 之间，构成的单增序列的值（即(prod_{t=1}^ia_t)）之和。易有：

[f(i,j)=jcdot f(i-1,j-1)+f(i,j-1) ]

  那么答案为：

[ans=n!f(n,k) ]

  复杂度 (O(nk))，不可过。

Step 2

  接下来，结合本文的主题，我们将会证明，当 (n) 为常数，(f(n,x)) 是关于 (x) 的 (2n+1) 次的函数。

证明

  差分 (f(i))，令 (g(i,j)=f(i,j)-f(i,j-1))。结合原题，有转移：

[g(i,j)=jsum_{t=0}^{j-1}g(i-1,t) ]

  同时，有：

[f(i,j)=sum_{t=0}^jg(i,j) ]

  由于对函数 (f(i)) 差分相当于将其次数 (-1)（最高次展开后抵消）。故只需证 (g(n,x)) 是关于 (x) 的 (2n) 次函数。

  考虑归纳：

• (1).) 当 (n=0) 时，(g(n,x)=[x=0])，成立。

• (2).) 当 (n=m) 时成立，考察 (n=m+1) 的情况：

  [g(m+1,x)=xsum_{t=0}^{x-1}g(m,t) ]

    求前缀和，次数 (+1)，再带上系数 (x)，次数总共 (+2)。由于 (n=m) 是次数为 (2m)，故 (n=m+1) 时次数为 (2m+2=2(m+1))。

    由 (1).~2).) (g(n,x)) 是关于 (x) 的 (2n) 次函数得证，故有 (f(n,x)) 是关于 (x) 的 (2n+1) 次函数，证毕。

  由此，对于 (kle2n+1) 的情况，我们直接 DP 得出答案；否则只需要 DP 出 (f(n,1..2n+1)) 的值，然后利用 Lagrange 插值法求出已被形如 (left(x,f(n,x) ight)) 的 (2n+1) 个点唯一确定的 (2n+1) 次函数，最后代入 (x=k) 求得答案即可。

  复杂度 (O(n^2))。

代码

#include <cstdio>

const int MAXN = 500;
int n, k, p, x[MAXN + 5], f[MAXN + 5][MAXN * 2 + 5];

inline int qkpow ( int a, int b ) {
	int ret = 1;
	for ( ; b; a = 1ll * a * a % p, b >>= 1 ) ret = 1ll * ret * ( b & 1 ? a : 1 ) % p;
	return ret;
}

inline int lagrange ( const int n, const int* x, const int* y, const int k ) {
	int ret = 0;
	for ( int i = 1; i <= n; ++ i ) {
		int num = y[i], den = 1;
		for ( int j = 1; j <= n; ++ j ) {
			if ( i ^ j ) {
				num = 1ll * num * ( k - x[j] + p ) % p;
				den = 1ll * den * ( x[i] - x[j] + p ) % p;
			}
		}
		ret = ( ret + 1ll * num * qkpow ( den, p - 2 ) ) % p;
	}
	return ret;
}

int main () {
	scanf ( "%d %d %d", &k, &n, &p );
	int fac = 1;
	for ( int i = 0; i <= k && i <= ( n << 1 | 1 ); ++ i ) f[0][i] = 1;
	for ( int i = 1; i <= n; fac = 1ll * fac * i ++ % p ) {
		for ( int j = 1; j <= k && j <= ( n << 1 | 1 ); ++ j ) {
			f[i][j] = ( 1ll * j * f[i - 1][j - 1] + f[i][j - 1] ) % p;
		}
	}
	if ( k <= ( n << 1 | 1 ) ) return printf ( "%d
", int ( 1ll * fac * f[n][k] % p ) ), 0;
	for ( int i = 1; i <= ( n << 1 | 1 ); ++ i ) x[i] = i;
	return printf ( "%d
", int ( 1ll * fac * lagrange ( n << 1 | 1, x, f[n], k ) % p ) ), 0;
}

「CF 995F」Cowmpany Cowmpensation

题意简述

  link.

  给定一棵 (n) 个结点的有根树，用 ([1,d]) 为每个结点染色，要求子结点颜色不大于父节点颜色，求方案数。对大素数取模。

数据规模

  (nle3000;~dle10^9)。

Solution

Step 1

  按照上一题的套路，我们还是先列出一个暴力的 DP。令 (f(u,i)) 表示以 (u) 为根的子树用 ([1,i]) 染色的方案数，转移显然有：

[f(u,i)=f(u,i-1)+prod_{vin son_u}f(v,i) ]

  直接 DP，复杂度 (O(nd))，不可过。

Step 2

  我们发现这个式子和上一题的其实有些类似，所以……

  大声说出你的猜想！

  Bingo！设 (u) 的子树大小为 (s_u)，则 (f(u,x)) 是关于 (x) 的 (s_u) 次函数。

证明

  仍考虑归纳：

• (1).) 当 (u) 为叶子结点时，(f(u,x)=1)，成立。

• (2).) 当 (vin son_u) 时成立，考虑 (u)：

  [f(u,x)-f(u,x-1)=prod_{vin son_u}f(v,x) ]

    由于 (v) 满足结论，即 (f(v,x)) 的次数为 (s_v)，则 (f(u,x)-f(u,x-1)) 的次数为 (sum_{vin son_u}s_v)，即 (s_u-1)。再还原差分，次数 (+1)，有 (f(u,x)) 是关于 (x) 的 (s_u-1+1=s_u) 次函数。

    由 (1).~2).) 证毕。

  所以，只需要处理出 (f(u,1..n))，最后利用 (left(x,f(1,x) ight)~(xin[0,n])) 求出 (f(1,x)) 的曲线方程，再代入 (d) 即可。

  连续插值部分可以优化到 (O(n))，不过这并不是唯一的瓶颈复杂度。最终复杂度仍是 (O(n^2))。

代码

#include <cstdio>

const int MAXN = 3000, MOD = 1e9 + 7;
int n, d, ecnt, head[MAXN + 5], x[MAXN + 5], f[MAXN + 5][MAXN + 5];

struct Edge { int to, nxt; } graph[MAXN + 5];

inline int qkpow ( int a, int b, const int p = MOD ) {
	int ret = 1;
	for ( ; b; a = 1ll * a * a % p, b >>= 1 ) ret = 1ll * ret * ( b & 1 ? a : 1 ) % p;
	return ret;
}

inline int lagrange ( const int n, const int* x, const int* y, const int k ) {
	int ret = 0;
	for ( int i = 1; i <= n; ++ i ) {
		int num = y[i], den = 1;
		for ( int j = 1; j <= n; ++ j ) {
			if ( i ^ j ) {
				num = 1ll * num * ( k - x[j] + MOD ) % MOD;
				den = 1ll * den * ( x[i] - x[j] + MOD ) % MOD;
			}
		}
		ret = ( ret + 1ll * num * qkpow ( den, MOD - 2 ) ) % MOD;
	}
	return ret;
}

inline void link ( const int s, const int t ) { graph[++ ecnt] = { t, head[s] }, head[s] = ecnt; }

inline void DFS ( const int u ) {
	for ( int i = 1; i <= n; ++ i ) f[u][i] = 1;
	for ( int i = head[u], v; i; i = graph[i].nxt ) {
		DFS ( v = graph[i].to );
		for ( int j = 1; j <= n; ++ j ) f[u][j] = 1ll * f[u][j] * f[v][j] % MOD;
	}
	for ( int i = 2; i <= n; ++ i ) f[u][i] = ( f[u][i] + f[u][i - 1] ) % MOD;
}

int main () {
	scanf ( "%d %d", &n, &d );
	for ( int i = 1; i <= n; ++ i ) x[i] = i;
	for ( int i = 2, f; i <= n; ++ i ) scanf ( "%d", &f ), link ( f, i );
	printf ( "%d
", ( DFS ( 1 ), d <= n ? f[1][d] : lagrange ( n + 1, x - 1, f[1] - 1, d ) ) );
	return 0;
}

「CF 662F」The Sum of the k-th Powers

题意简述

  link.

  给定 (n,k)，求 (sum_{i=1}^ni^kmod(10^9+7))。

数据规模

  (nle10^9;~kle10^6)。

Solution

  见下文结论 1。（谋篇布局的大失败qwq）暴力求出前 (O(k)) 个值，插值即可。

  复杂度 (O(klog(10^9+7)))。

代码

#include <cstdio>

#define inv( x ) qkpow ( x, MOD - 2 )

const int MAXN = 1e6, MOD = 1e9 + 7;
int n, k, y[MAXN + 5], ifac[MAXN + 5];

inline int qkpow ( int a, int b, const int p = MOD ) {
	int ret = 1;
	for ( ; b; a = 1ll * a * a % p, b >>= 1 ) ret = 1ll * ret * ( b & 1 ? a : 1 ) % p;
	return ret;
}

int main () {
	scanf ( "%d %d", &n, &k ), ifac[0] = 1;
	for ( int i = 1; i <= k + 1; ++ i ) ifac[i] = 1ll * ifac[i - 1] * inv ( i ) % MOD;
	for ( int i = 1; i <= k + 1; ++ i ) y[i] = ( y[i - 1] + qkpow ( i, k ) ) % MOD;
	if ( n <= k + 1 ) return printf ( "%d
", y[n] ), 0;
	int all = 1, ans = 0;
	for ( int i = 0; i <= k + 1; ++ i ) all = 1ll * all * ( n - i + MOD ) % MOD;
	for ( int i = 0; i <= k + 1; ++ i ) {
		ans = ( ans + ( ( k + 1 - i ) & 1 ? ( MOD - 1ll ) : 1ll )
					* y[i] % MOD * all % MOD * inv ( ( n - i + MOD ) % MOD ) % MOD
					* ifac[i] % MOD * ifac[k + 1 - i] % MOD ) % MOD;
	}
	return printf ( "%d
", ans ), 0;
}

「BZOJ 3453」tyvj 1858 XLkxc

题意简述

  给定 (k,a,n,d)，求：

[sum_{i=0}^nsum_{j=1}^{a+id}sum_{x=1}^jx^kmod1234567891 ]

数据规模

  (kle123;~a,n,dle123456789)。

  什么鬼畜的范围哟……

Solution

  别问，问就是信仰插值。

  算法本身很显，我们先来证明一个对于分析函数次数比较有用的结论吧。

结论 1

  前 (n) 个正整数的 (k) 次幂和是 (k+1) 次多项式。

  证明可参考这篇博客。

结论 2

  设 (f(x),g(x)) 是关于 (x) 的多项式函数，(f(x)) 的次数为 (p_f)，(g(x)) 的次数为 (p_g)，对于函数：

[h(x)=sum_{i=0}^xf(g(i)) ]

  其中 (xinmathbb N)，则 (h(x)) 的次数 (p_h=p_f+p_g+1)。

证明

  既然只考虑最高次，不妨 (f(x)=x^{p_f},g(x)=x^{p_g})，代入有 (h(x)=sum_{i=0}^xi^{p_f+p_g})，利用结论 1 即可证明。

本题

  记 (f(x)=sum_{x=1}^jx^k)，(g(x)=sum_{j=1}^{a+id}sum_{x=1}^jx^k=sum_{j=1}^{a+id}f(j))，(h(x)=sum_{i=0}^nsum_{j=1}^{a+id}sum_{x=1}^jx^k=sum_{i=0}^ng(a+id))。由结论 2，(f(x),g(x),h(x)) 的次数分别为 (k+1,k+2,k+3)，直接插值即可。

  复杂度 (O(Tk^2))。

代码

#include <cstdio>

#define inv( x ) qkpow ( ( x ) % MOD, MOD - 2 )

typedef long long LL;

const LL MAXK = 123, MOD = 1234567891;
LL k, a, n, d, ifac[MAXK + 5], s[MAXK + 5], t[MAXK + 5], x[MAXK + 5], y[MAXK + 5];

inline LL qkpow ( LL a, LL b, const LL p = MOD ) {
	LL ret = 1;
	for ( ; b; a = a * a % p, b >>= 1 ) ret = ret * ( b & 1 ? a : 1 ) % p;
	return ret;
}

inline LL calc ( const LL k, const LL* y, const LL n ) {
	if ( n <= k ) return y[n];
	LL all = 1, ret = 0;
	for ( int i = 0; i <= k; ++ i ) all = all * ( n - i + MOD ) % MOD;
	for ( int i = 0; i <= k; ++ i ) {
		ret = ( ret + ( ( k - i ) & 1 ? ( MOD - 1 ) : 1 )
					* y[i] % MOD * all % MOD * inv ( ( n - i + MOD ) % MOD ) % MOD
					* ifac[i] % MOD * ifac[k - i] % MOD ) % MOD;
	}
	return ret;
}

int main () {
	ifac[0] = 1;
	for ( int i = 1; i <= MAXK + 3; ++ i ) ifac[i] = inv ( i ) * ifac[i - 1] % MOD;
	int T;
	for ( scanf ( "%d", &T ); T --; ) {
		scanf ( "%lld %lld %lld %lld", &k, &a, &n, &d );
		for ( int i = 1; i <= k + 2; ++ i ) s[i] = ( s[i - 1] + qkpow ( i, k ) ) % MOD;
		for ( int i = 1; i <= k + 2; ++ i ) t[i] = ( t[i - 1] + s[i] ) % MOD;
		y[0] = calc ( k + 2, t, a );
		for ( int i = 1; i <= k + 3; ++ i ) y[i] = ( y[i - 1] + calc ( k + 2, t, ( a + i * d % MOD ) % MOD ) ) % MOD;
		printf ( "%lld
", calc ( k + 3, y, n ) );
	}
	return 0;
}

「CF 1086F」Forest Fires

题意简述

  link.

  一个无限大的网格图上在第 (0) 时刻有 (n) 个格子着火，每一时刻，每个着火的格子会让周围八连通的格子着火。一个格子的权值定义为此格子最早的着火时刻，没有着火的格子权值为 (0)。求第 (t) 时刻所有格子的权值之和。

数据规模

  (nle50;~|x_i|,|y_i|,tle10^8)。

Solution

Step 0 前置-矩形面积并

  求网格图上多个可交矩形的覆盖面积和，用扫描线+线段树 (O(nlog n)) 求解。这里不赘述。

Step 1

  令 (f(i)) 表示 (t) 时刻时权值为 (i) 的格子的个数。那么答案就是 (sum_{i=0}^ticdot f(i)) 如果强行求每一个 (f(i)) 是很复杂的，因为权值为 (i) 的格子的分布是许多个环。

  换种思路——对 (f(i)) 求前缀和，令 (g(i)=sum_{j=0}^if(i))。可以发现 (g(i)) 实际上就是 (i) 时刻网格图上被覆盖的面积之和。并且 (g(t)-g(i)) 实质上就是权值大于 (i) 的格子数量。所以答案又可以表示为 (sum_{i=0}^{t-1}g(t)-g(i)=tg(t)-sum_{i=0}^{t-1}g(i))。

  如果暴力求每个 (g(i))，复杂度 (O(tnlog n))。原地爆炸。

Step 2

  再次结合本文主题，我们来研究一下关于 (x) 的函数 (g(x)) 的次数性质。

  首先，假设矩形都不交。单独考虑每个矩形，它的面积按时刻递增呈现匀加速增长，这是典型的二次函数。所以，常数（指与 (x) 无关）个二次函数之和仍然是二次函数。

  接着，考虑两矩形相交后的情形。我们可以把这两个矩形整体看做一个不规则图形。单独考虑上下左右每一面，其随时刻增加，仍是匀加速增涨，面积亦为二次函数。

  综合上面两个结论，可以得到：当没有新的矩形相交时，(g(x)) 是关于 (x) 的二次函数。也即是，(g(x)) 是一个分段函数，每一段都是一个二次函数。

  显然段数为 (O(n^2))，对于每一段，先暴力算出三个 (g(x)) 的值，以此解出 (g(x)=ax^2+bx+c)（可以插值，也可以手玩解方程），再求一个 (g(x)) 的前缀和即可。

  复杂度 (O(n^3log n))。

代码

#include <cstdio>
#include <vector>
#include <algorithm>

inline int rint () {
	int x = 0, f = 1; char s = getchar ();
	for ( ; s < '0' || '9' < s; s = getchar () ) f = s == '-' ? -f : f;
	for ( ; '0' <= s && s <= '9'; s = getchar () ) x = x * 10 + ( s ^ '0' );
	return x * f;
}

template<typename Tp>
inline void wint ( Tp x ) {
	if ( x < 0 ) putchar ( '-' ), x = ~ x + 1;
	if ( 9 < x ) wint ( x / 10 );
	putchar ( x % 10 ^ '0' );
}

inline int abs_ ( const int x ) { return x < 0 ? -x : x; }

inline int max_ ( const int a, const int b ) { return a < b ? b : a; }

const int MAXN = 50, MOD = 998244353, INV2 = 499122177, INV6 = 166374059;
int n, t, x[MAXN + 5], y[MAXN + 5];

namespace MatrixInter {

int tmp[MAXN * 4 + 5];

class SegmentTree {
private:
	int cover[MAXN << 4], cnt[MAXN << 4];

public:
	inline void clear ( const int rt, const int l, const int r ) {
		cover[rt] = cnt[rt] = 0;
		if ( l == r ) return ;
		int mid = l + r >> 1;
		clear ( rt << 1, l, mid ), clear ( rt << 1 | 1, mid + 1, r );
	}
	inline void pushup ( const int rt, const int l, const int r ) {
		cover[rt] = cnt[rt] ? tmp[r] - tmp[l - 1] : cover[rt << 1] + cover[rt << 1 | 1];
	}
	inline void add ( const int rt, const int al, const int ar, const int l, const int r, const int v ) {
		if ( al <= l && r <= ar ) {
			if ( ! ( cnt[rt] += v ) ) cover[rt] = l == r ? 0 : cover[rt << 1] + cover[rt << 1 | 1];
			else cover[rt] = tmp[r] - tmp[l - 1];
			return ;
		}
		int mid = l + r >> 1;
		if ( al <= mid ) add ( rt << 1, al, ar, l, mid, v );
		if ( mid < ar ) add ( rt << 1 | 1, al, ar, mid + 1, r, v );
		pushup ( rt, l, r );
	}
	inline int getTop () { return cover[1]; }
} st;

struct Event { int l, r, opt; };
std :: vector<Event> evt[MAXN * 4 + 5];

inline int calcS ( const int t ) {
	int cnt = 0;
	for ( int i = 1; i <= n; ++ i ) {
		tmp[++ cnt] = x[i] - t - 1, tmp[++ cnt] = x[i] + t;
		tmp[++ cnt] = y[i] - t - 1, tmp[++ cnt] = y[i] + t;
	}
	std :: sort ( tmp + 1, tmp + cnt + 1 );
	cnt = std :: unique ( tmp + 1, tmp + cnt + 1 ) - tmp - 1;
	st.clear ( 1, 1, cnt );
	for ( int i = 1; i <= cnt; ++ i ) evt[i].clear ();
	for ( int i = 1, up, dn, le, ri; i <= n; ++ i ) {
		up = std :: lower_bound ( tmp + 1, tmp + cnt + 1, y[i] - t - 1 ) - tmp + 1;
		dn = std :: lower_bound ( tmp + 1, tmp + cnt + 1, y[i] + t ) - tmp;
		le = std :: lower_bound ( tmp + 1, tmp + cnt + 1, x[i] - t - 1 ) - tmp;
		ri = std :: lower_bound ( tmp + 1, tmp + cnt + 1, x[i] + t ) - tmp;
		evt[le].push_back ( { up, dn, true } ), evt[ri].push_back ( { up, dn, false } );
	}
	int ret = 0;
	for ( int i = 1; i <= cnt; ++ i ) {
		ret = ( ret + 1ll * st.getTop () * ( tmp[i] - tmp[i - 1] ) % MOD ) % MOD;
		for ( Event e: evt[i] ) st.add ( 1, e.l, e.r, 1, cnt, e.opt ? 1 : -1 );
	}
	return ret;
}

} // namespace MatrixInter.

int itime[MAXN * MAXN + 5];

inline int preS ( const int n ) { return n * ( n + 1ll ) % MOD * INV2 % MOD; }

inline int sqrS ( const int n ) { return n * ( n + 1ll ) % MOD * ( n * 2ll + 1 ) % MOD * INV6 % MOD; }

int main () {
#define MI MatrixInter
	n = rint (), t = rint ();
	for ( int i = 1; i <= n; ++ i ) x[i] = rint (), y[i] = rint ();
	int icnt = 0;
	for ( int i = 1; i < n; ++ i ) {
		for ( int j = i + 1, tmpt; j <= n; ++ j ) {
			tmpt = max_ ( abs_ ( x[i] - x[j] ), abs ( y[i] - y[j] ) ) >> 1;
			if ( tmpt <= t ) itime[++ icnt] = tmpt;
			if ( ++ tmpt <= t ) itime[++ icnt] = tmpt;
		}
	}
	itime[++ icnt] = 0, itime[++ icnt] = t;
	std :: sort ( itime + 1, itime + icnt + 1 );
	icnt = std :: unique ( itime + 1, itime + icnt + 1 ) - itime - 1;
	int ans = 1ll * t * MI :: calcS ( t ) % MOD;
	for ( int i = 1, tl, tr, a, b, c, ta, tb; i < icnt; ++ i ) {
		if ( ( tl = itime[i] ) + 3 > ( tr = itime[i + 1] ) ) {
			for ( ; tl < tr; ++ tl ) ans = ( ans - MI :: calcS ( tl ) + MOD ) % MOD;
			continue;
		}
		c = MI :: calcS ( tl ), ta = MI :: calcS ( tl + 1 ), tb = MI :: calcS ( tl + 2 );
		a = ( ( tb - ta * 2ll % MOD + c ) % MOD + MOD ) % MOD * INV2 % MOD;
		b = ( ( ta - c - a ) % MOD + MOD ) % MOD;
		ans = ( ans - 1ll * a * sqrS ( tr - tl - 1 ) % MOD + MOD ) % MOD;
		ans = ( ans - 1ll * b * preS ( tr - tl - 1 ) % MOD + MOD ) % MOD;
		ans = ( ans - 1ll * c * ( tr - tl ) % MOD + MOD ) % MOD;
	}
	return wint ( ans ), putchar ( '
' ), 0;
}