# PARI sourcecode for demo in A055881

Dec 18th, 2012
1,173
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. /*
2. This demonstrative program should be used for either documentation or instructive purposes only,
3. and goes without any warranty whatsoever.
4.
5. This is part of: The On-Line Encyclopedia Of Integer Sequences OEIS(TM)
6.
7. Visit: http://www.oeis.org/A055881
8.
9. Author: Joerg Arndt (www.jjj.de), Dec 12 2012
10. Pastebin by: R. J. Cano (remy@ula.ve), Dec 17 2012
11.
12. */
13.
14. perm(n)=
15. /* Compute the first n! terms of A055881 by counting
16.   in rising factorial base.
17.   Generate all permutations of n along the way.
18. */
19. {
20.     my( fc=vector(n) ); /* mixed radix numbers (rising factorial base) */
21.     my( ct=0, j, a=0 );
22.     my( p=vector(n,k,k-1) ); /* permutation */
23.     my( k, t );  /* for permutations */
24.     while ( 1,
25.         print1(ct, ":  ", p); /* permutation */
26.         print1("  ", a );  /* A055881(ct): position of rightmost change in fc[] */
27.         print1("    ", fc);  /* factorial number */
28.         print();
29.         ct += 1;
30.
31.         /* increment factorial number fc[]: */
32.         j = 1;
33.         while ( fc[j] == j,  fc[j]=0;  j+=1; );  /* scan over nines */
34.         if ( j==n,  return() );  /* current is last */
35.         fc[j] += 1;
36.         /* update permutation p[], reverse prefix of length j+1: */
37.         a = j;  /* next term of A055881 */
38.         j += 1;  k = 1;
39.         while ( k < j,
40.             t=p[j]; p[j]=p[k]; p[k]=t;
41.             k+=1; j-=1;
42.         );
43.     );
44. }
45.
46. print("Demonstration: Reproduces the example given in the comments section for A055881");
47. perm(4);
48. print(" ");
49. print("You are free to call perm() and experiment with it. When finished, type \"quit\" to exit.");
50. print(" ");