(PARI) {calc_two_diags() = /* M, n, n2 inherited; down first diag., up next */
for(r=2, n, M[r, n2r] = M[r1, n2r+1]  M[r1, n2r]);
forstep(r=n, 1, 1, M[r, n2r+1] = M[r, n2r] + M[r+1, n2r])}
/* Find pairbypair first 2k terms of A153123 and put into row 1 of M. */
/* Row n+1 of M contains 2k2n+1 terms of A153123's nth differences  */
/* for efficiency not starting until the nth term in the nth column. */
/* (Minor changes to program would calculate all nth differences.) */
/* For larger k the PARI stack size may need to be increased first. */
{k = 450; M = matrix(k+1, 2*k); M[1, 1] = 1; M[1, 2] = 2; M[2, 1] = 1;
for(n=2, k, M[n+1, n] = n; n2 = 2*n; M[1, n21] = M[1, n22] + 1;
calc_two_diags();
if(M[1, n2] <= M[1, n21], /* If first try doesn't find increasing terms, */
M[1, n21] = M[1, n21] + 1 + floor(M[2, n21]/(n1)); calc_two_diags()));
/* the above adjusts both terms just enough.*/
M[1, ] /* Print only first row of M */}
/* The above program uses these general facts: */
/* Increasing (Decreasing) one term, a(j), of a sequence by a constant c */
/* i) increases (decreases) each element of the leftmost changing */
/* diagonal of its difference array by c and */
/* ii) increases (decreases) the elements of the diagonal immediately */
/* to its right  including a(j+1)  by c, 2c, 3c, ..., dc, respectively, */
/* for the d elements above a given element on this second diagonal*/
/* which is being held fixed, where dc is the change in a(j+1). */
/* (Incidentally, a similar pattern occurs for the elements below */
/* the fixed one but the sign of the changes is reversed.) */
