(* Mathematica program for A000620-A000625, 
   Jean-Francois Alcover, Jul 18 2022, 
   translated and adapted from N. J. A. Sloane's Maple code *) 

Ps = {0, 0, 0}; Pn = {1, 1, 1}; Ss = {0, 0, 0}; 
Sn = {0, 0, 1}; Ts = {0, 0, 0}; Tn = {0, 0, 0};
As = {0, 0, 0}; An = {1, 1, 2}; T = {1, 1, 2}; 
Do[AppendTo[Ps, As[[n - 1]]]; AppendTo[Pn, An[[n - 1]]]; 
  t1 = Sum[2*T[[n - 1 - j]]*T[[j]], {j, 1, Floor[(n - 2)/2]}]; 
  If[Mod[n, 2] == 1, i = (n - 1)/2; t1 = t1 + T[[i]]^2 - An[[i]]]; 
  AppendTo[Ss, t1]; t2 = 0; If[Mod[n, 2] == 1, i = (n - 1)/2; 
   t2 = An[[i]]]; AppendTo[Sn, t2]; t3 = 0; 
  Do[Do[k = n - 1 - i - j; If[j < k, t3 = t3 + 2*T[[i]]*T[[j]]*T[[k]]],
       {j, i + 1, (n - 2)/2}], {i, 1, (n - 1)/3}]; 
  t4 = 0; t5 = 0; 
  Do[j = n - 1 - 2*i; If[j > 0 && i != j, 
       t4 = t4 + (T[[i]]^2 - An[[i]])*T[[j]] + An[[i]]*As[[j]]; 
       t5 = t5 + An[[i]]*An[[j]]],
      {i, 1, (n - 2)/2}];
  t6 = 0; t7 = 0; If[Mod[n, 3] == 1, i = (n - 1)/3; 
   t6 = (2*T[[i]] + T[[i]]^3)/3 - An[[i]]^2; 
   t7 = An[[i]]^2]; 
  AppendTo[Ts, t3 + t4 + t6]; 
  AppendTo[Tn, t5 + t7]; 
  AppendTo[As, Ps[[n]] + Ss[[n]] + Ts[[n]]]; 
  AppendTo[An, Pn[[n]] + Sn[[n]] + Tn[[n]]]; 
  AppendTo[T, As[[n]] + An[[n]]], {n, 4, 60}]; 

A000620 = Ps; A000621 = Pn; A000622 = Ss; 
A000623 = Ts; A000624 = Tn; A000625 = T;