Source code in the scilab programming language

clear
xdel(winsid())

stacksize('max')
sz=stacksize();

n=7; //For the expansion up to power of n
g=50; //For test of primes up to g

function X = pascal(g) //Pascal´s triangle
    X(1,1)=1; //Zeroth power
    X(2,1)=1; //First power
    X(2,2)=1; 
    for q=3:1:g+1 //From second power use this loop
       X(q,1)=1;
       X(q,q)=1;
        for p=2:1:q-1
            X(q,p)=X(q-1,p-1)+X(q-1,p);
        end
    end
endfunction

Z=pascal(g); //Generate Pascal's triangle up to g

Q(0+1)="(x-1)^0 = 1"; //For nicer display
Q(1+1)="(x-1)^1 = x^1-1"; //For nicer display

disp(Q(1))
disp(Q(2))

function cf=coef(Z,q,p) //Return coeffiecents for nicer display of expansion without "ones"
    if Z(q,p)==1 then
        cf="";
    else
        cf=string(Z(q,p));
    end
endfunction

for q=3:n+1  //Generate and display the expansions
    Q(q)=strcat(["(x-1)^",string(q-1)," = "]);
    sing=""; //Sign of coeff.
        for p=1:q-1 //Number of coefficients equals power minus 1
            Q(q)=strcat([Q(q),sing,coef(Z,q,p),"x^",string(q-p)]);
            if sing=="-" then sing="+"; else sing="-"; end
        end
        Q(q)=strcat([Q(q),sing,string(1)]);
    disp(Q(q))
    clear Q
end

function prime=prime(Z,g)
    prime="true";
    for p=2:g
        if abs(floor(Z(g+1,p)/g)-Z(g+1,p)/g)>0 then
            prime="false";
            break;
        end
    end
endfunction

R="2"; //For nicer display
for r=3:g
    if prime(Z,r)=="true" then 
        R=strcat([R, ", ",string(r)]);
    end
end
disp(R)