\\ Pari functions for Permutations
\\ Created by Michael Hortmann, Michael.Hortmann@math.uni-Bremen.de

\\ Permutations are represented by vectors of n elements.
\\ Generally, inputs are not type checked.
\\ Can be directly imported into Pari, e.g. by copy/paste.



\\ multiplication of two permutations of same length
perm_mult(p,q)={
    local(n=length(p),r=vector(n));
    if(n!=length(q),return("p und q are not of same length"));

    for(i=1,n,r[i]=p[q[i]]);

    return(r)
}

\\ this function returns the identity permutation
perm_id(n)=return(vector(n,i,i))


\\ Inverse permutation
perm_inverse(p)={
    local(n=length(p),q=vector(n));
    for(i=1,n,q[p[i]]=i);
    return(q)
}

\\ square of a permutation
perm_square(p)=return(perm_mult(p,p))


\\ n-th power of a permuation, inefficient for large n
perm_smallpow(p,n)={
    if(n==0,return(vector(length(p),i,i)));

    m=divrem(n,2);
    if(n>0,return(perm_mult(p,perm_pow(p,n-1))));
    if(n<0,return(perm_pow(perm_inverse(p),-n)))
}

\\ create a random permutation of length n.
perm_random(n)={
    local(p=vector(n,i,i),tmp,j);

    for(i=1,n,tmp=p[i];p[i]=p[j=random(n)+1];p[j]=tmp);
    return(p)
}

\\ real neat algorithm for creating a random permutation of length n
perm_random_recursive(n)={
    local(i=random(n)+1,p);

    if(n==1,return([1]));
    
    p=concat(perm_random_recursive(n-1),[0]);
    p[n]=p[i];p[i]=n;

    return(p)
}

\\ auxiliary function: returns the binary representation of n as a vector
bits(n)={
    local(x=[],m);
    
    while(n,n=(m=divrem(n,2))[1];x=concat(m[2],x));
    return(x)
}

\\ efficient computation of exponentiation by repeated squaring
perm_pow(p,n)={
    local(m,b,l,v=p);

    if(n==0,return(vector(length(p),i,i)));
    if(n==1,return(p));


    if(n>1,
        b=bits(n);
        l=length(b);
        
        for(i=2,l,if(b[i],v=perm_mult(perm_square(v),p),v=perm_square(v)))
    );
        
    if(n<0,v=perm_pow(perm_inverse(p),-n));

    return(v)
}

\\ auxiliary function: reverses factor: factor_inverse(factor(n))==n
factor_inverse(v)=return(prod(i=1,matsize(v)[1],v[i,1]^v[i,2]))


\\ inefficient computation of the order of a permutation
\\ by repeated division of primefactors out of the factorial or length(p)
perm_order(p)={
    local(n=length(p),m=n!,v=factor(m),k=length(v[,2]),t);

    for(i=1,k,while(perm_pow(p,t=m/v[i,1])==perm_id(n),m=t;v[i,2]-=1));
    return(factor_inverse(v))
}

\\ computes the least common multiple of the components of a vector
lcmv(v)=if(length(v)==2,return(lcm(v[1],v[2])),return(lcm(v[1],lcmv(cut(v)))))

\\ computes the cycle structure of a permutation
\\ e.g. computes [[1,4],[5,2,3]] from [4,3,5,1,2]
perm_cycles(p)={
    local(v,u=[],w=[],i,count,n=length(p));
    
    while(count<n,    
        v=vector(1); i=1;count++;
        v[1]=p[eval(setminus(Set(p),Set(u)))[1]];
        while(p[v[i]]!=v[1],v=concat(v,p[v[i]]);i++;count++);

        w=concat(w,[v]);u=concat(u,v)
    );

    return(w)
}

\\ the cycle structure of a permutation makes computation of order effective:
\\ the order is simply the least common multiple of the length of the cycles.
perm_order_cycles(p)={
    local(q=perm_cycles(p));
    return(lcmv(vector(length(q),i,length(q[i]))));
}

\\ recomputes the original permuation from the output of perm_cycles
perm_cycles_reverse(w)={
    local(l,p,count);

    l=sum(i=1,length(w),length(w[i]));
    p=vector(l);

    for(i=1,length(w),for(j=1,length(w[i])-1, p[w[i][j]]=w[i][j+1]); p[w[i][length(w[i])]]=w[i][1]);

    return(p)
}

\\ unused auxiliary functions
cut(v)=vecextract(v,concat("2..",Str(length(v))))
cut2(v)=vector(length(v)-1,i,v[i+1])
cut3(v,i)=vector(length(v)-i,j,v[i+j])