[R] sign(<permutation>) in R ?

[Martin Maechler]

As alluded, I've got several off-line replies and suggestions,
notably Tobias Verbeke mentioning Matlab code {as an example of "pseudo
code" :-) } by Burkhardt on the web,
and Barry Rowlinson providing code similar to Peter Dalgaard's,
see below.
I've found myself python code on the web, unfortunately bugous
(only for *some* permutation vectors).

Of all I've got, Peter's code was clearly the best, when
n <- length(<permutation>) was not small.
For that reason, I've ``speeded'' Peter's code, and ended in
code that won hands down.

The simulation-testing of the diverse algorithms has been fun,
not the least finding that one "public" algorithm was faulty;
If I have a bit more time, I should write a small R news article
on this.
Here's the current best code
{and of course, a C version of that will be yet an order of
 magnitude faster} :

## MM: Peter's function is so compact and elegant,
##     now try to optimize for speed() :
cycLengths2 <- function(p) {
    n <- length(p)
    x <- integer(n)
    ii <- seq_len(n)
    for (i in ii) {
	z <- ii[!x][1] # index of first unmarked x[] entry
	if (is.na(z)) break
	repeat { ## mark x[] <- i  for those in i-th cycle
	    x[z] <- i
	    z <- p[z]
	    if (x[z]) break
	}
    }
    ## Now, table(x) gives the cycle lengths,
    ## where  split(seq_len(n), x)  would give the cycles list
    tabulate(x, i - 1L) ## is quite a bit faster than the equivalent
    ## table(x)
}

## MM: Here's the transformation from the cycle lengths' vector to sign(<p>):
signCycLength <- function(clen) 1L - (sum(clen %% 2 == 0) %% 2)*2L

for the *sign* of a permutation 'p'
the above has to be called as  signCycLength(cycLengths2(p))
e.g.,

   p <- sample(50)
   signCycLength(cycLengths2(p))


Thank you, again for all the feedback, and suggestions,
Martin Maechler, ETH Zurich




>>>>>> Peter Dalgaard wrote yesterday

>> Martin Maechler wrote:
>> > Ok,
>> > thanks a lot, everyone!
>> >
>> > Yes, I should rather have started thinking a bit more myself,
>> > before going the easy route to R-help....
>> >
>> > Anyway, the most obvious algorithm,
>> > just putting things into place by swapping elements,
>> > and counting how many times you have to swap, is easy and
>> > quite efficient.
>> >
>> > I'll post R code later, being busy for the next few hours.
>> > Martin
>> >
>> >   
>> It is also quite easy to generate the cycles explicitly by a marking 
>> algorithm:
>> 
>> p <- sample(1:100)
>> x <- integer(100)
>> for (i in 1:100) {
>>     z <- which(!x)[1]
>>     if (is.na(z)) break
>>     repeat {
>>         x[z] <- i
>>         z <- p[z]
>>         if (x[z]) break
>>     }
>> }
>> table(x)
>> 
>> (the which(!x)[1] bit could be optimized somewhat if this had been C. 
>> Notice that the loop is essentially O(N) because every iteration marks 
>> one element of x)
>> 
>> 
>> >   
>> >>>>>> "MM" == Martin Maechler <maechler at stat.math.ethz.ch>
>> >>>>>>     on Tue, 15 Apr 2008 18:13:43 +0200 writes:
>> >>>>>>             
>> >
>> >     MM> I am looking for an algorithm (written in R (preferably) or C,
>> >     MM> but even pseudo-code in a text book maybe fine)
>> >     MM> to determine the sign of a permutation.
>> >
>> >     MM> What is that?  Well, a permutation is either even or odd, the
>> >     MM> sign is +1 or -1, respectively, see, e.g.,
>> >     MM> http://en.wikipedia.org/wiki/Signature_of_a_permutation
>> >     MM> which also says
>> >     >>> In practice, in order to determine whether a given permutation
>> >     >>> is even or odd, one writes the permutation as a product of
>> >     >>> disjoint cycles. The permutation is odd if and only if this
>> >     >>> factorization contains an odd number of even-length cycles.
>> >
>> >     MM> but I would not know how to algorithmically
>> >     MM> "write the permutation as a product of disjoint cycles"
>> >
>> >     MM> If you start looking at R code, 
>> >     MM> let's assume the permutation {\pi(i)}_{i=1..n} is simply given
>> >     MM> as the (integer) vector (\pi(1), \pi(2), ..., \pi(n))
>> >     MM> {or equivalently, a random permutation is simply found by  'sample(n)'}
>> >
>> >     MM> Thank you in advance for further pointers,
>> >     MM> or even working R code.
>> >
>> >     MM> Best regards,
>> >     MM> Martin Maechler, ETH Zurich
>> >

>> 
>> -- 
>>    O__  ---- Peter Dalgaard             Øster Farimagsgade 5, Entr.B
>>   c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
>>  (*) \(*) -- University of Copenhagen   Denmark      Ph:  (+45) 35327918
>> ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)              FAX: (+45) 35327907