如何在之字形排序中找到第 i 项？
How to find ith item in zigzag ordering?
上周的一个问题定义了 n x m 矩阵上的之字形排序并提出了 how to list the elements in that order。
我的问题是如何快速找到之字形排序中的第i项？也就是说，不遍历矩阵（对于大 n 和 m 太慢了）。
例如图片中的 n=m=8 和 (x, y) 描述（行、列）
f(0) = (0, 0)
f(1) = (0, 1)
f(2) = (1, 0)
f(3) = (2, 0)
f(4) = (1, 1)
...
f(63) = (7, 7)
具体问题：百万乘百万矩阵之字形排序中的第十亿（1e10）项是什么？
假设所需元素位于矩阵的上半部分。对角线的长度是1, 2, 3 ..., n。
让我们找到所需的对角线。它满足以下属性：
sum(1, 2 ..., k) >= pos 但 sum(1, 2, ..., k - 1) < pos。 1, 2, ..., k 的总和是 k * (k + 1) / 2。所以我们只需要找到最小的整数 k 使得 k * (k + 1) / 2 >= pos。我们可以使用二分搜索或明确地解决这个二次不等式。
当我们知道了k，我们只需要找到这条对角线的pos - (k - 1) * k / 2元素即可。我们知道它从哪里开始以及我们应该移动到哪里（向上或向下，取决于 k 的奇偶校验），因此我们可以使用一个简单的公式找到所需的单元格。
此解的时间复杂度为O(1)或O(log n)（这取决于我们是使用二分查找还是在步骤2中显式求解不等式）。
如果想要的元素位于矩阵的下半部分，我们可以将这个问题求解一个pos' = n * n - pos + 1，然后利用对称性得到原问题的解。
我在这个解决方案中使用了基于 1 的索引，使用基于 0 的索引可能需要在某处添加 +1 或 -1，但解决方案的思路是相同的。
如果矩阵是矩形而不是正方形，我们需要考虑这样一个事实：当我们搜索 k，因此如果 k <= m 则总和变为 k * (k + 1) / 2，否则 k * (k + 1) / 2 + m * (k - m)。
import math, random
def naive(n, m, ord, swap = False):
dx = 1
dy = -1
if swap:
dx, dy = dy, dx
cur = [0, 0]
for i in range(ord):
cur[0] += dy
cur[1] += dx
if cur[0] < 0 or cur[1] < 0 or cur[0] >= n or cur[1] >= m:
dx, dy = dy, dx
if cur[0] >= n:
cur[0] = n - 1
cur[1] += 2
if cur[1] >= m:
cur[1] = m - 1
cur[0] += 2
if cur[0] < 0: cur[0] = 0
if cur[1] < 0: cur[1] = 0
return cur
def fast(n, m, ord, swap = False):
if n < m:
x, y = fast(m, n, ord, not swap)
return [y, x]
alt = n * m - ord - 1
if alt < ord:
x, y = fast(n, m, alt, swap if (n + m) % 2 == 0 else not swap)
return [n - x - 1, m - y - 1]
if ord < (m * (m + 1) / 2):
diag = int((-1 + math.sqrt(1 + 8 * ord)) / 2)
parity = (diag + (0 if swap else 1)) % 2
within = ord - (diag * (diag + 1) / 2)
if parity: return [diag - within, within]
else: return [within, diag - within]
else:
ord -= (m * (m + 1) / 2)
diag = int(ord / m)
within = ord - diag * m
diag += m
parity = (diag + (0 if swap else 1)) % 2
if not parity:
within = m - within - 1
return [diag - within, within]
if __name__ == "__main__":
for i in range(1000):
n = random.randint(3, 100)
m = random.randint(3, 100)
ord = random.randint(0, n * m - 1)
swap = random.randint(0, 99) < 50
na = naive(n, m, ord, swap)
fa = fast(n, m, ord, swap)
assert na == fa, "(%d, %d, %d, %s) ==> (%s), (%s)" % (n, m, ord, swap, na, fa)
print fast(1000000, 1000000, 9999999999, False)
print fast(1000000, 1000000, 10000000000, False)
所以第100亿个元素（序号为9999999999的）和第100亿个元素（序号为10^10的）是：
[20331, 121089]
[20330, 121090]
解析解
在一般情况下，您的矩阵将分为 3 个区域：
- 初始三角形t1
- 偏斜部分mid 其中对角线长度恒定
- 最后一个三角形t2
我们将 p 称为对角线的索引 运行。
我们想定义两个函数 x(p) 和 y(p) 来给出 p 的列和行第个单元格。
初始三角形
让我们看一下初始三角形部分 t1，其中每条新对角线都比前一条长一个单位。
现在让我们调用 d 包含单元格的对角线的索引，并且
Sp = sum(di) for i in [0..p -1]
我们有 p = Sp + k，其中 0 <=k <= d 和
Sp = d(d+1)/2
如果我们求解 d，它带来
d²+d-2p = 0，一个只保留正根的二次方程：
d = (-1+sqrt(1+8*p))/2
现在我们想要最接近d的最高整数值，即floor(d).
最后，我们有
p = d + k 和 d = floor((-1+sqrt(1+8*p))/2) 和 k = p - d(d+1)/2
我们打电话
o(d) 如果 d 等于 1 的函数是 odd 和 0 否则，和
e(d) 如果 d 等于 1 的函数是 even 否则为 0。
我们可以这样计算 x(p) 和 y(p)：
d = floor((-1+sqrt(1+8*p))/2)
k = p - d(d+1)/2
o = d % 2
e = 1 - o
x = e*d + (o-e)*k
y = o*d + (e-o)*k
even and odd functions are used to try to salvage some clarity, but you can replace
e(p) with 1 - o(p) and have slightly more efficient but less symetric formulaes for x and y.
中间部分
让我们考虑最小的矩阵维度s，即s = min (m,n).
前面的公式一直有效，直到 x 或 y（以先到者为准）达到值 s。
p 的上限，例如 x(i) <= s 和 y(i) <= s 对于所有 i 在 [0..p]
（即由 p 索引的单元格在初始三角形 t1 内）由
给出
pt1 = s(s+1)/2.
对于 p >= pt1，对角线长度保持等于 s 直到我们到达第二个三角形 t2.
当在 mid 内时，我们有：
p = s(s+1)/2 + ds + k with k in [0..s[。
产生：
d = floor ((p - s(s+1)/2)/s)
k = p - ds
然后我们可以使用相同的 even/odd 技巧来计算 x(p) 和 y(p):
p -= s(s+1)/2
d = floor (p / s)
k = p - d*s
o = (d+s) % 2
e = 1 - o
x = o*s + (e-o)*k
y = e*s + (o-e)*k
if (n > m)
x += d+e
y -= e
else
y += d+o
x -= o
最终三角形
利用对称性，我们可以计算出pt2 = m*n - s(s+1)/2
我们现在面临与 t1 几乎相同的问题，只是对角线可能 运行 与 t1[= 的方向相同163=] 或反方向（若n+m为奇数）.
使用对称技巧，我们可以像这样计算 x(p) 和 y(p)：
p = n*m -1 - p
d = floor((-1+sqrt(1+8*p))/2)
k = p - d*(d+1)/2
o = (d+m+n) % 2
e = 1 - $o;
x = n-1 - (o*d + (e-o)*k)
y = m-1 - (e*d + (o-e)*k)
综合起来
这是一个示例 C++ 实现。
- 出于纯粹的懒惰，我使用了 64 位整数。大多数可以替换为 32 位值。
- 可以通过预先计算更多的系数来提高计算效率。
- 大部分代码都可以分解，但我怀疑这样做是否值得。
由于这只是一个快速而肮脏的概念证明，我没有对其进行优化。
#include <cstdio> // printf
#include <algorithm> // min
using namespace std;
typedef long long tCoord;
void panic(const char * msg)
{
printf("PANIC: %s\n", msg);
exit(-1);
}
struct tPoint {
tCoord x, y;
tPoint(tCoord x = 0, tCoord y = 0) : x(x), y(y) {}
tPoint operator+(const tPoint & p) const { return{ x + p.x, y + p.y }; }
bool operator!=(const tPoint & p) const { return x != p.x || y != p.y; }
};
class tMatrix {
tCoord n, m; // dimensions
tCoord s; // smallest dimension
tCoord pt1, pt2; // t1 / mid / t2 limits for p
public:
tMatrix(tCoord n, tCoord m) : n(n), m(m)
{
s = min(n, m);
pt1 = (s*(s + 1)) / 2;
pt2 = n*m - pt1;
}
tPoint diagonal_cell(tCoord p)
{
tCoord x, y;
if (p < pt1) // inside t1
{
tCoord d = (tCoord)floor((-1 + sqrt(1 + 8 * p)) / 2);
tCoord k = p - (d*(d + 1)) / 2;
tCoord o = d % 2;
tCoord e = 1 - o;
x = o*d + (e - o)*k;
y = e*d + (o - e)*k;
}
else if (p < pt2) // inside mid
{
p -= pt1;
tCoord d = (tCoord)floor(p / s);
tCoord k = p - d*s;
tCoord o = (d + s) % 2;
tCoord e = 1 - o;
x = o*s + (e - o)*k;
y = e*s + (o - e)*k;
if (m > n) // vertical matrix
{
x -= o;
y += d + o;
}
else // horizontal matrix
{
x += d + e;
y -= e;
}
}
else // inside t2
{
p = n * m - 1 - p;
tCoord d = (tCoord)floor((-1 + sqrt(1 + 8 * p)) / 2);
tCoord k = p - (d*(d + 1)) / 2;
tCoord o = (d + m + n) % 2;
tCoord e = 1 - o;
x = n - 1 - (o*d + (e - o)*k);
y = m - 1 - (e*d + (o - e)*k);
}
return{ x, y };
}
void check(void)
{
tPoint move[4] = { { 1, 0 }, { -1, 1 }, { 1, -1 }, { 0, 1 } };
tPoint pos;
tCoord dir = 0;
for (tCoord p = 0; p != n * m ; p++)
{
tPoint dc = diagonal_cell(p);
if (pos != dc) panic("zot!");
pos = pos + move[dir];
if (dir == 0)
{
if (pos.y == m - 1) dir = 2;
else dir = 1;
}
else if (dir == 3)
{
if (pos.x == n - 1) dir = 1;
else dir = 2;
}
else if (dir == 1)
{
if (pos.y == m - 1) dir = 0;
else if (pos.x == 0) dir = 3;
}
else
{
if (pos.x == n - 1) dir = 3;
else if (pos.y == 0) dir = 0;
}
}
}
};
void main(void)
{
const tPoint dim[] = { { 10, 10 }, { 11, 11 }, { 10, 30 }, { 30, 10 }, { 10, 31 }, { 31, 10 }, { 11, 31 }, { 31, 11 } };
for (tPoint d : dim)
{
printf("Checking a %lldx%lld matrix...", d.x, d.y);
tMatrix(d.x, d.y).check();
printf("done\n");
}
tCoord p = 10000000000;
tMatrix matrix(1000000, 1000000);
tPoint cell = matrix.diagonal_cell(p);
printf("Coordinates of %lldth cell: (%lld,%lld)\n", p, cell.x, cell.y);
}
根据 "manual" 扫描矩阵检查结果。
这种 "manual" 扫描是一种丑陋的技巧，不适用于单行或单列矩阵，尽管 diagonal_cell() 对任何矩阵都有效（ "diagonal" 在这种情况下扫描变为线性）。
为 1.000.000x1.000.000 矩阵的第 10.000.000.000 个单元格找到的坐标似乎是一致的，因为单元格所在的对角线 d 大约为 sqrt(2 *1e10), 约。 141421，单元格坐标之和约等于d（121090+20330 = 141420）。此外，这也是另外两位海报的报道。
我想说这块混淆代码很有可能实际上为您的问题生成一个复杂度为 O(1) 的解决方案。
上周的一个问题定义了 n x m 矩阵上的之字形排序并提出了 how to list the elements in that order。
我的问题是如何快速找到之字形排序中的第i项？也就是说，不遍历矩阵（对于大 n 和 m 太慢了）。
例如图片中的 n=m=8 和 (x, y) 描述（行、列）
f(0) = (0, 0)
f(1) = (0, 1)
f(2) = (1, 0)
f(3) = (2, 0)
f(4) = (1, 1)
...
f(63) = (7, 7)
具体问题：百万乘百万矩阵之字形排序中的第十亿（1e10）项是什么？
假设所需元素位于矩阵的上半部分。对角线的长度是
1, 2, 3 ..., n。让我们找到所需的对角线。它满足以下属性：
sum(1, 2 ..., k) >= pos但sum(1, 2, ..., k - 1) < pos。1, 2, ..., k的总和是k * (k + 1) / 2。所以我们只需要找到最小的整数k使得k * (k + 1) / 2 >= pos。我们可以使用二分搜索或明确地解决这个二次不等式。当我们知道了
k，我们只需要找到这条对角线的pos - (k - 1) * k / 2元素即可。我们知道它从哪里开始以及我们应该移动到哪里（向上或向下，取决于k的奇偶校验），因此我们可以使用一个简单的公式找到所需的单元格。
此解的时间复杂度为O(1)或O(log n)（这取决于我们是使用二分查找还是在步骤2中显式求解不等式）。
如果想要的元素位于矩阵的下半部分，我们可以将这个问题求解一个pos' = n * n - pos + 1，然后利用对称性得到原问题的解。
我在这个解决方案中使用了基于 1 的索引，使用基于 0 的索引可能需要在某处添加 +1 或 -1，但解决方案的思路是相同的。
如果矩阵是矩形而不是正方形，我们需要考虑这样一个事实：当我们搜索 k，因此如果 k <= m 则总和变为 k * (k + 1) / 2，否则 k * (k + 1) / 2 + m * (k - m)。
import math, random
def naive(n, m, ord, swap = False):
dx = 1
dy = -1
if swap:
dx, dy = dy, dx
cur = [0, 0]
for i in range(ord):
cur[0] += dy
cur[1] += dx
if cur[0] < 0 or cur[1] < 0 or cur[0] >= n or cur[1] >= m:
dx, dy = dy, dx
if cur[0] >= n:
cur[0] = n - 1
cur[1] += 2
if cur[1] >= m:
cur[1] = m - 1
cur[0] += 2
if cur[0] < 0: cur[0] = 0
if cur[1] < 0: cur[1] = 0
return cur
def fast(n, m, ord, swap = False):
if n < m:
x, y = fast(m, n, ord, not swap)
return [y, x]
alt = n * m - ord - 1
if alt < ord:
x, y = fast(n, m, alt, swap if (n + m) % 2 == 0 else not swap)
return [n - x - 1, m - y - 1]
if ord < (m * (m + 1) / 2):
diag = int((-1 + math.sqrt(1 + 8 * ord)) / 2)
parity = (diag + (0 if swap else 1)) % 2
within = ord - (diag * (diag + 1) / 2)
if parity: return [diag - within, within]
else: return [within, diag - within]
else:
ord -= (m * (m + 1) / 2)
diag = int(ord / m)
within = ord - diag * m
diag += m
parity = (diag + (0 if swap else 1)) % 2
if not parity:
within = m - within - 1
return [diag - within, within]
if __name__ == "__main__":
for i in range(1000):
n = random.randint(3, 100)
m = random.randint(3, 100)
ord = random.randint(0, n * m - 1)
swap = random.randint(0, 99) < 50
na = naive(n, m, ord, swap)
fa = fast(n, m, ord, swap)
assert na == fa, "(%d, %d, %d, %s) ==> (%s), (%s)" % (n, m, ord, swap, na, fa)
print fast(1000000, 1000000, 9999999999, False)
print fast(1000000, 1000000, 10000000000, False)
所以第100亿个元素（序号为9999999999的）和第100亿个元素（序号为10^10的）是：
[20331, 121089]
[20330, 121090]
解析解
在一般情况下，您的矩阵将分为 3 个区域：
- 初始三角形t1
- 偏斜部分mid 其中对角线长度恒定
- 最后一个三角形t2
我们将 p 称为对角线的索引 运行。
我们想定义两个函数 x(p) 和 y(p) 来给出 p 的列和行第个单元格。
初始三角形
让我们看一下初始三角形部分 t1，其中每条新对角线都比前一条长一个单位。
现在让我们调用 d 包含单元格的对角线的索引，并且
Sp = sum(di) for i in [0..p -1]
我们有 p = Sp + k，其中 0 <=k <= d 和
Sp = d(d+1)/2
如果我们求解 d，它带来
d²+d-2p = 0，一个只保留正根的二次方程：
d = (-1+sqrt(1+8*p))/2
现在我们想要最接近d的最高整数值，即floor(d).
最后，我们有 p = d + k 和 d = floor((-1+sqrt(1+8*p))/2) 和 k = p - d(d+1)/2
我们打电话
o(d) 如果 d 等于 1 的函数是 odd 和 0 否则，和
e(d) 如果 d 等于 1 的函数是 even 否则为 0。
我们可以这样计算 x(p) 和 y(p)：
d = floor((-1+sqrt(1+8*p))/2)
k = p - d(d+1)/2
o = d % 2
e = 1 - o
x = e*d + (o-e)*k
y = o*d + (e-o)*k
even and odd functions are used to try to salvage some clarity, but you can replace
e(p) with 1 - o(p) and have slightly more efficient but less symetric formulaes for x and y.
中间部分
让我们考虑最小的矩阵维度s，即s = min (m,n).
前面的公式一直有效，直到 x 或 y（以先到者为准）达到值 s。
p 的上限，例如 x(i) <= s 和 y(i) <= s 对于所有 i 在 [0..p]
（即由 p 索引的单元格在初始三角形 t1 内）由
给出
pt1 = s(s+1)/2.
对于 p >= pt1，对角线长度保持等于 s 直到我们到达第二个三角形 t2.
当在 mid 内时，我们有：
p = s(s+1)/2 + ds + k with k in [0..s[。
产生：
d = floor ((p - s(s+1)/2)/s)
k = p - ds
然后我们可以使用相同的 even/odd 技巧来计算 x(p) 和 y(p):
p -= s(s+1)/2
d = floor (p / s)
k = p - d*s
o = (d+s) % 2
e = 1 - o
x = o*s + (e-o)*k
y = e*s + (o-e)*k
if (n > m)
x += d+e
y -= e
else
y += d+o
x -= o
最终三角形
利用对称性，我们可以计算出pt2 = m*n - s(s+1)/2
我们现在面临与 t1 几乎相同的问题，只是对角线可能 运行 与 t1[= 的方向相同163=] 或反方向（若n+m为奇数）.
使用对称技巧，我们可以像这样计算 x(p) 和 y(p)：
p = n*m -1 - p
d = floor((-1+sqrt(1+8*p))/2)
k = p - d*(d+1)/2
o = (d+m+n) % 2
e = 1 - $o;
x = n-1 - (o*d + (e-o)*k)
y = m-1 - (e*d + (o-e)*k)
综合起来
这是一个示例 C++ 实现。
- 出于纯粹的懒惰，我使用了 64 位整数。大多数可以替换为 32 位值。
- 可以通过预先计算更多的系数来提高计算效率。
- 大部分代码都可以分解，但我怀疑这样做是否值得。
由于这只是一个快速而肮脏的概念证明，我没有对其进行优化。
#include <cstdio> // printf
#include <algorithm> // min
using namespace std;
typedef long long tCoord;
void panic(const char * msg)
{
printf("PANIC: %s\n", msg);
exit(-1);
}
struct tPoint {
tCoord x, y;
tPoint(tCoord x = 0, tCoord y = 0) : x(x), y(y) {}
tPoint operator+(const tPoint & p) const { return{ x + p.x, y + p.y }; }
bool operator!=(const tPoint & p) const { return x != p.x || y != p.y; }
};
class tMatrix {
tCoord n, m; // dimensions
tCoord s; // smallest dimension
tCoord pt1, pt2; // t1 / mid / t2 limits for p
public:
tMatrix(tCoord n, tCoord m) : n(n), m(m)
{
s = min(n, m);
pt1 = (s*(s + 1)) / 2;
pt2 = n*m - pt1;
}
tPoint diagonal_cell(tCoord p)
{
tCoord x, y;
if (p < pt1) // inside t1
{
tCoord d = (tCoord)floor((-1 + sqrt(1 + 8 * p)) / 2);
tCoord k = p - (d*(d + 1)) / 2;
tCoord o = d % 2;
tCoord e = 1 - o;
x = o*d + (e - o)*k;
y = e*d + (o - e)*k;
}
else if (p < pt2) // inside mid
{
p -= pt1;
tCoord d = (tCoord)floor(p / s);
tCoord k = p - d*s;
tCoord o = (d + s) % 2;
tCoord e = 1 - o;
x = o*s + (e - o)*k;
y = e*s + (o - e)*k;
if (m > n) // vertical matrix
{
x -= o;
y += d + o;
}
else // horizontal matrix
{
x += d + e;
y -= e;
}
}
else // inside t2
{
p = n * m - 1 - p;
tCoord d = (tCoord)floor((-1 + sqrt(1 + 8 * p)) / 2);
tCoord k = p - (d*(d + 1)) / 2;
tCoord o = (d + m + n) % 2;
tCoord e = 1 - o;
x = n - 1 - (o*d + (e - o)*k);
y = m - 1 - (e*d + (o - e)*k);
}
return{ x, y };
}
void check(void)
{
tPoint move[4] = { { 1, 0 }, { -1, 1 }, { 1, -1 }, { 0, 1 } };
tPoint pos;
tCoord dir = 0;
for (tCoord p = 0; p != n * m ; p++)
{
tPoint dc = diagonal_cell(p);
if (pos != dc) panic("zot!");
pos = pos + move[dir];
if (dir == 0)
{
if (pos.y == m - 1) dir = 2;
else dir = 1;
}
else if (dir == 3)
{
if (pos.x == n - 1) dir = 1;
else dir = 2;
}
else if (dir == 1)
{
if (pos.y == m - 1) dir = 0;
else if (pos.x == 0) dir = 3;
}
else
{
if (pos.x == n - 1) dir = 3;
else if (pos.y == 0) dir = 0;
}
}
}
};
void main(void)
{
const tPoint dim[] = { { 10, 10 }, { 11, 11 }, { 10, 30 }, { 30, 10 }, { 10, 31 }, { 31, 10 }, { 11, 31 }, { 31, 11 } };
for (tPoint d : dim)
{
printf("Checking a %lldx%lld matrix...", d.x, d.y);
tMatrix(d.x, d.y).check();
printf("done\n");
}
tCoord p = 10000000000;
tMatrix matrix(1000000, 1000000);
tPoint cell = matrix.diagonal_cell(p);
printf("Coordinates of %lldth cell: (%lld,%lld)\n", p, cell.x, cell.y);
}
根据 "manual" 扫描矩阵检查结果。
这种 "manual" 扫描是一种丑陋的技巧，不适用于单行或单列矩阵，尽管 diagonal_cell() 对任何矩阵都有效（ "diagonal" 在这种情况下扫描变为线性）。
为 1.000.000x1.000.000 矩阵的第 10.000.000.000 个单元格找到的坐标似乎是一致的，因为单元格所在的对角线 d 大约为 sqrt(2 *1e10), 约。 141421，单元格坐标之和约等于d（121090+20330 = 141420）。此外，这也是另外两位海报的报道。
我想说这块混淆代码很有可能实际上为您的问题生成一个复杂度为 O(1) 的解决方案。