# Maple code for Gilbreath transform and its underlying array
N. J. A. Sloane, May 09 2023

# Computes first M terms of Gilbreath transform of sequence [f(1), f(2), f(3), ...]

# get next row
nexr:=proc(a) local b,i; b:=[seq(abs(a[i+1]-a[i]),i=1..nops(a)-1)]; end;

# the main program
GA := proc(f,M) local U,A,s,b,i,j; global nexr;
A:=Array(1..M, 1..M):
s:=[seq(f(i),i=1..M)];
for j from 1 to M do A[1,j]:=f(j); od:
for i from 2 to M do 
s := nexr(s);  
for j from 1 to M+1-i do A[i,j]:=s[j]; od:
od:
lprint("Gilbreath array:");
for i from 1 to M do
lprint(seq(A[i,j],j=1..M+1-i),"...");
od;
lprint("...");
lprint("antidiagonals of Gilbreath array:"); 
for i from 1 to M do
lprint(seq(A[i+1-j,j],j=1..i));
od;
lprint("...");
lprint("Gilbreath array formatted for OEIS:");
U:=[]; 
for i from 1 to M do
U:=[op(U), seq(A[i+1-j,j],j=1..i)];
od;
lprint(U);
lprint("Gilbreath transform:");
lprint([seq(A[i,1],i=1..M)]);
end;

# Example: the Gilbreath transform of the sum-of-divisors function sigma (A000203)

with(numtheory);
GA(sigma,8);
This produces the following output:

"Gilbreath array:"
1, 3, 4, 7, 6, 12, 8, 15, "..."
2, 1, 3, 1, 6, 4, 7, "..."
1, 2, 2, 5, 2, 3, "..."
1, 0, 3, 3, 1, "..."
1, 3, 0, 2, "..."
2, 3, 2, "..."
1, 1, "..."
0, "..."
"..."
"Gilbreath antidiagonals:"
1
2, 3
1, 1, 4
1, 2, 3, 7
1, 0, 2, 1, 6
2, 3, 3, 5, 6, 12
1, 3, 0, 3, 2, 4, 8
0, 1, 2, 2, 1, 3, 7, 15
"..."
"Gilbreath array for OEIS:"
[1, 2, 3, 1, 1, 4, 1, 2, 3, 7, 1, 0, 2, 1, 6, 2, 3, 3, 5, 6, 12, 1, 3, 0, 3, 2, 4, 8, 0, 1, 2, 2, 1, 3, 7, 15]   # This is A362464
"Gilbreath transform:"
[1, 2, 1, 1, 1, 2, 1, 0]  # This is A362451