在小于 O(N) 的时间内找到序列的第 n 项

Question

这个问题的时间复杂度不同于被问到的类似问题。这是来自 Zauba 开发人员招聘挑战（活动于一个月前结束）的问题：

f(0) = p

f(1) = q

f(2) = r

for n > 2

f(n) = a*f(n-1) + b*f(n-2) + c*f(n-3) + g(n)

where g(n) = n*n*(n+1)

p, q, r, a, b, c, n给出。n可以大到10^18.

链接到类似的问题

在上面的链接中，没有指定时间复杂度，我已经在 中解决了这个问题O(n)，下面是伪代码（只是一种方法，所有可能的边界和边缘情况都在比赛中处理）。

if(n == 0) return p;

if(n == 1) return q;

if(n == 2) return r;

for(long i=3;i<=n;i++){

    now = a*r + b*q + c*p + i*i*(i+1);

    p = q; q = r; r = now;

}

请注意，我10^9 + 7在原始代码中的适当位置使用模数来处理溢出，在必要时处理适当的边缘情况，并且我使用了 java long 数据类型（如果它有帮助）。

但由于这仍然需要O(n)时间，我期待一个更好的解决方案来处理n ~ 10^18.

编辑

正如用户 גלעד ברקן 提到它与矩阵求幂的关系一样，我尝试这样做并停留在一个特定的点，我不确定在矩阵的第 4 行第 3 列中放置什么。请提出任何建议和更正。

| a b c  1? |   | f(n) |        | f(n+1) |

| 1 0 0  0  |   |f(n-1)|        |  f(n)  |

| 0 1 0  0  |   |f(n-2)|    =>  | f(n-1) |

| 0 0 ?! 0  |   | g(n) |        | g(n+1) |

    M               A               B

Answer 1

矩阵求幂确实是正确的方法，但还有一些工作要做。由于g(n)不是常数值，因此无法有效地（O(log n)而不是O(n)）将矩阵求幂应用于当前形式的递推关系。需要找到类似的递归关系，g(n)只有一个常数项尾随。由于g(n)是三次方，因此需要 3 个递归项：g(n) = x*g(n-1) + y*g(n-2) + z*g(n-3) + w展开它们每个的三次表达式：n³ + n² = x(n³-2n²+n) + y(n³-5n²+8n-4) + z*(n³-8n²+21n-18) + w&nbsp; &nbsp; &nbsp; &nbsp; = n³(x+y+z) + n²(-2x-5y-8z) + n(x+8y+21z) + (w-4y-18z)匹配系数以获得三个联立方程 加上x, y, z另一个来计算w：&nbsp; x +&nbsp; y +&nbsp; &nbsp;z = 1-2x - 5y -&nbsp; 8z = 1&nbsp; x + 8y + 21z = 0&nbsp; w - 4y - 18z = 0解决它们得到：x = 3&nbsp; &nbsp; y = -3&nbsp; &nbsp; z = 1&nbsp; &nbsp; w = 6方便的是，这些系数也是整数*，这意味着可以直接对递归进行模运算。* 我怀疑这是巧合——这很可能是招聘考官的意图。因此，矩阵递推方程为：|&nbsp; a&nbsp; b&nbsp; c&nbsp; 1&nbsp; 0&nbsp; 0&nbsp; 0 |&nbsp; &nbsp;| f(n-1) |&nbsp; &nbsp;|&nbsp; &nbsp;f(n) ||&nbsp; 1&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0 |&nbsp; &nbsp;| f(n-2) |&nbsp; &nbsp;| f(n-1) ||&nbsp; 0&nbsp; 1&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0 |&nbsp; &nbsp;| f(n-3) |&nbsp; &nbsp;| f(n-2) ||&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 3 -3&nbsp; 1&nbsp; 6 | x |&nbsp; &nbsp;g(n) | = | g(n+1) ||&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 1&nbsp; 0&nbsp; 0&nbsp; 0 |&nbsp; &nbsp;| g(n-1) |&nbsp; &nbsp;|&nbsp; &nbsp;g(n) ||&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 1&nbsp; 0&nbsp; 0 |&nbsp; &nbsp;| g(n-2) |&nbsp; &nbsp;| g(n-1) ||&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 1 |&nbsp; &nbsp;|&nbsp; &nbsp; &nbsp; 1 |&nbsp; &nbsp;|&nbsp; &nbsp; &nbsp; 1 |最终的矩阵求幂方程为：&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; [n-2]|&nbsp; a&nbsp; b&nbsp; c&nbsp; 1&nbsp; 0&nbsp; 0&nbsp; 0 |&nbsp; &nbsp; &nbsp; &nbsp;| f(2) |&nbsp; &nbsp;|&nbsp; &nbsp;f(n) |&nbsp; &nbsp; &nbsp; &nbsp; | f(2) |&nbsp; &nbsp;|&nbsp; r ||&nbsp; 1&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0 |&nbsp; &nbsp; &nbsp; &nbsp;| f(1) |&nbsp; &nbsp;| f(n-1) |&nbsp; &nbsp; &nbsp; &nbsp; | f(1) |&nbsp; &nbsp;|&nbsp; q ||&nbsp; 0&nbsp; 1&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0 |&nbsp; &nbsp; &nbsp; &nbsp;| f(0) |&nbsp; &nbsp;| f(n-2) |&nbsp; &nbsp; &nbsp; &nbsp; | f(0) |&nbsp; &nbsp;|&nbsp; p ||&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 3 -3&nbsp; 1&nbsp; 6 |&nbsp; &nbsp;x&nbsp; &nbsp;| g(3) | = | g(n+1) |&nbsp; &nbsp;,&nbsp; &nbsp; | g(3) | = | 36 ||&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 1&nbsp; 0&nbsp; 0&nbsp; 0 |&nbsp; &nbsp; &nbsp; &nbsp;| g(2) |&nbsp; &nbsp;|&nbsp; &nbsp;g(n) |&nbsp; &nbsp; &nbsp; &nbsp; | g(2) |&nbsp; &nbsp;| 12 ||&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 1&nbsp; 0&nbsp; 0 |&nbsp; &nbsp; &nbsp; &nbsp;| g(1) |&nbsp; &nbsp;| g(n-1) |&nbsp; &nbsp; &nbsp; &nbsp; | g(1) |&nbsp; &nbsp;|&nbsp; 2 ||&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 0&nbsp; 1 |&nbsp; &nbsp; &nbsp; &nbsp;|&nbsp; 1&nbsp; &nbsp;|&nbsp; &nbsp;|&nbsp; &nbsp; &nbsp; 1 |&nbsp; &nbsp; &nbsp; &nbsp; |&nbsp; 1&nbsp; &nbsp;|&nbsp; &nbsp;|&nbsp; 1 |（每个操作都是隐式模数10^9 + 7或提供的任何此类数字。）请注意，Java 的%运算符是remainder ，它与负数的模数不同。例子：-1 % 5 == -1&nbsp; &nbsp; &nbsp;// Java-1 = 4 (mod 5)&nbsp; &nbsp;// mathematical modulus解决方法很简单：long mod(long b, long a){&nbsp; &nbsp; // computes a mod b&nbsp; &nbsp; // assumes that b is positive&nbsp; &nbsp; return (b + (a % b)) % b;}原始迭代算法：long recurrence_original(&nbsp; &nbsp; long a, long b, long c,&nbsp; &nbsp; long p, long q, long r,&nbsp; &nbsp; long n, long m // 10^9 + 7 or whatever) {&nbsp; &nbsp; // base cases&nbsp; &nbsp; if (n == 0) return p;&nbsp; &nbsp; if (n == 1) return q;&nbsp; &nbsp; if (n == 2) return r;&nbsp; &nbsp; long f0, f1, f2;&nbsp; &nbsp; f0 = p; f1 = q; f2 = r;&nbsp; &nbsp; for (long i = 3; i <= n; i++) {&nbsp; &nbsp; &nbsp; &nbsp; long f3 = mod(m,&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; mod(m, a*f2) + mod(m, b*f1) + mod(m, c*f0) +&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; mod(m, mod(m, i) * mod(m, i)) * mod(m, i+1)&nbsp; &nbsp; &nbsp; &nbsp; );&nbsp; &nbsp; &nbsp; &nbsp; f0 = f1; f1 = f2; f2 = f3;&nbsp; &nbsp; }&nbsp; &nbsp; return f2;}模矩阵函数：long[][] matrix_create(int n){&nbsp; &nbsp; return new long[n][n];}void matrix_multiply(int n, long m, long[][] c, long[][] a, long[][] b){&nbsp; &nbsp; for (int i = 0; i < n; i++) {&nbsp; &nbsp; &nbsp; &nbsp; for (int j = 0; j < n; j++) {&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; long s = 0;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; for (int k = 0; k < n; k++)&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; s = mod(m, s + mod(m, a[i][k]*b[k][j]));&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; c[i][j] = s;&nbsp; &nbsp; &nbsp; &nbsp; }&nbsp; &nbsp; }}void matrix_pow(int n, long m, long p, long[][] y, long[][] x){&nbsp; &nbsp; // swap matrices&nbsp; &nbsp; long[][] a = matrix_create(n);&nbsp; &nbsp; long[][] b = matrix_create(n);&nbsp; &nbsp; long[][] c = matrix_create(n);&nbsp; &nbsp; // initialize accumulator to identity&nbsp; &nbsp; for (int i = 0; i < n; i++)&nbsp; &nbsp; &nbsp; &nbsp; a[i][i] = 1;&nbsp; &nbsp; // initialize base to original matrix&nbsp; &nbsp; for (int i = 0; i < n; i++)&nbsp; &nbsp; &nbsp; &nbsp; for (int j = 0; j < n; j++)&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; b[i][j] = x[i][j];&nbsp; &nbsp; // exponentiation by squaring&nbsp; &nbsp; // there are better algorithms, but this is the easiest to implement&nbsp; &nbsp; // and is still O(log n)&nbsp; &nbsp; long[][] t = null;&nbsp; &nbsp; for (long s = p; s > 0; s /= 2) {&nbsp; &nbsp; &nbsp; &nbsp; if (s % 2 == 1) {&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; matrix_multiply(n, m, c, a, b);&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; t = c; c = a; a = t;&nbsp; &nbsp; &nbsp; &nbsp; }&nbsp; &nbsp; &nbsp; &nbsp; matrix_multiply(n, m, c, b, b);&nbsp; &nbsp; &nbsp; &nbsp; t = c; c = b; b = t;&nbsp; &nbsp; }&nbsp; &nbsp; // write to output&nbsp; &nbsp; for (int i = 0; i < n; i++)&nbsp; &nbsp; &nbsp; &nbsp; for (int j = 0; j < n; j++)&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; y[i][j] = a[i][j];}最后，新算法本身：long recurrence_matrix(&nbsp; &nbsp; long a, long b, long c,&nbsp; &nbsp; long p, long q, long r,&nbsp; &nbsp; long n, long m) {&nbsp; &nbsp; if (n == 0) return p;&nbsp; &nbsp; if (n == 1) return q;&nbsp; &nbsp; if (n == 2) return r;&nbsp; &nbsp; // original recurrence matrix&nbsp; &nbsp; long[][] mat = matrix_create(7);&nbsp; &nbsp; mat[0][0] = a; mat[0][1] = b; mat[0][2] = c; mat[0][3] = 1;&nbsp; &nbsp; mat[1][0] = 1; mat[2][1] = 1;&nbsp; &nbsp; mat[3][3] = 3; mat[3][4] = -3; mat[3][5] = 1; mat[3][6] = 6;&nbsp; &nbsp; mat[4][3] = 1; mat[5][4] = 1;&nbsp; &nbsp; mat[6][6] = 1;&nbsp; &nbsp; // exponentiate&nbsp; &nbsp; long[][] res = matrix_create(7);&nbsp; &nbsp; matrix_pow(7, m, n - 2, res, mat);&nbsp; &nbsp; // multiply the first row with the initial vector&nbsp; &nbsp; return mod(m, mod(m, res[0][6])&nbsp; &nbsp; &nbsp; &nbsp; + mod(m, res[0][0]*r)&nbsp; + mod(m, res[0][1]*q)&nbsp; + mod(m, res[0][2]*p)&nbsp; &nbsp; &nbsp; &nbsp; + mod(m, res[0][3]*36) + mod(m, res[0][4]*12) + mod(m, res[0][5]*2)&nbsp; &nbsp; );}以下是上述两种算法的一些示例基准。原始迭代算法：n&nbsp; &nbsp; &nbsp; &nbsp;time (μs)-------------------10^1&nbsp; &nbsp; 9.310^2&nbsp; &nbsp; 44.910^3&nbsp; &nbsp; 401.50110^4&nbsp; &nbsp; 3882.09910^5&nbsp; &nbsp; 27940.910^6&nbsp; &nbsp; 88873.59910^7&nbsp; &nbsp; 877100.510^8&nbsp; &nbsp; 9057329.09910^9&nbsp; &nbsp; 91749994.4新矩阵算法：n&nbsp; &nbsp; &nbsp; &nbsp;time (μs)------------------10^1&nbsp; &nbsp; 69.16810^2&nbsp; &nbsp; 128.77110^3&nbsp; &nbsp; 212.69710^4&nbsp; &nbsp; 258.38510^5&nbsp; &nbsp; 318.19510^6&nbsp; &nbsp; 380.910^7&nbsp; &nbsp; 453.48710^8&nbsp; &nbsp; 560.42810^9&nbsp; &nbsp; 619.83510^10&nbsp; &nbsp;652.34410^11&nbsp; &nbsp;750.51810^12&nbsp; &nbsp;769.90110^13&nbsp; &nbsp;851.84510^14&nbsp; &nbsp;934.91510^15&nbsp; &nbsp;1016.73210^16&nbsp; &nbsp;1079.61310^17&nbsp; &nbsp;1123.41310^18&nbsp; &nbsp;1225.323旧算法用了 90 多秒来计算n = 10^9，而新算法只用了 0.6多毫秒（150,000 倍的加速）！原始算法的时间复杂度显然是线性的（正如预期的那样）；n = 10^10花了太长时间才完成，所以我没有继续。新算法的时间复杂度显然是对数的——将 的数量级加倍n导致执行时间加倍（同样，正如预期的那样，由于平方求幂）。n对于( )的“小”值，< 100矩阵分配和运算的开销掩盖了算法本身，但随着n增加而迅速变得微不足道。