q =. ;
Solve[(1 - q) (1 + q)^3 == 1, q]

Print["Variant with approximate numbers up to precision prec:"]
prec = 15; steps = 30;
q = -(2/3) + (1/3) (((19 - 3 (Sqrt[33])))^(1/3)) + 
          (1/3) (((19 + 3 (Sqrt[33])))^(1/3));
a = 1; b = (1 + q) a; c = (1 + q) b; 
    d = (1 + q) c;
a = N[a, prec]; b = N[b, prec]; c = N[c, prec]; 
    d = N[d, prec];
Print[a, ",  ", b, ",  ", c, ",  ", 
        d];
For[j = 1, 
        j ≤ 
          steps, j++, 
{a, b, c, d} = {b - a, c - b, 
            d - c, d - a}; 
Print[a, ", ", b, ", ", 
          c, ", ", d]];

Print["Variant with exact numbers:"];
prec = 15; steps = 30;
q = -(2/3) + (1/3) (((19 - 3 (Sqrt[33])))^(1/3)) + 
          (1/3) (((19 + 3 (Sqrt[33])))^(1/3));
x =. ;
a = 1; b = (1 + x) a; c = (1 + x) b; 
    d = (1 + x) c;
Print[N[a /. {x -> q}, prec], ",  ", 
        N[b /. {x -> q}, prec], ",  ", 
        N[c /. {x -> q}, prec], ",  ", 
        N[d /. {x -> q}, prec]];
For[j = 1, 
        j ≤ 
          steps, j++, 
{a, b, c, 
            d} = {PolynomialMod[b - a, 2 - 2 (q^2) - (q^3)], 
            PolynomialMod[c - b, 2 - 2 (q^2) - (q^3)], 
            PolynomialMod[d - c, 2 - 2 (q^2) - (q^3)], 
            PolynomialMod[d - a, 
              2 - 2 (q^2) - (q^3)]}; 
Print[
          N[a /. {x -> q}, prec], ",  ", 
          N[b /. {x -> q}, prec], ",  ", 
          N[c /. {x -> q}, prec], ",  ", 
          N[d /. {x -> q}, prec]]];