# T is the 3x+1 function discussed in equation (2.1)
restart:  p := 50: n := 2^p -1: # Change input p for T




 print(`n = `.``(2)^p-1 =  n);

 T := proc(n:posint)
 if type(n,odd) then (3*n+1)/2 else n/2 fi
 end:

 s := proc(n:posint) # calculates the expansion factor.
 local m,c ;

 c := n; m := n;

 while c > 1 do
 c := T(c); m := max(m,c);
 od;

 m/n;
 end:


 print(`s(`.n.`)` =  s(n));

 s(27);