 # PARI sourcecode for demo in A055881

Dec 18th, 2012
621
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. perm(n)=
13. /* Compute the first n! terms of A055881 by counting
14.   in rising factorial base.
15.   Generate all permutations of n along the way.
16. */
17. {
18.     my( fc=vector(n) ); /* mixed radix numbers (rising factorial base) */
19.     my( ct=0, j, a=0 );
20.     my( p=vector(n,k,k-1) ); /* permutation */
21.     my( k, t );  /* for permutations */
22.     while ( 1,
23.         print1(ct, ":  ", p); /* permutation */
24.         print1("  ", a );  /* A055881(ct): position of rightmost change in fc[] */
25.         print1("    ", fc);  /* factorial number */
26.         print();
27.         ct += 1;
28.
29.         /* increment factorial number fc[]: */
30.         j = 1;
31.         while ( fc[j] == j,  fc[j]=0;  j+=1; );  /* scan over nines */
32.         if ( j==n,  return() );  /* current is last */
33.         fc[j] += 1;
34.         /* update permutation p[], reverse prefix of length j+1: */
35.         a = j;  /* next term of A055881 */
36.         j += 1;  k = 1;
37.         while ( k < j,
38.             t=p[j]; p[j]=p[k]; p[k]=t;
39.             k+=1; j-=1;
40.         );
41.     );
42. }
43.
44. print("Demonstration: Reproduces the example given in the comments section for A055881");
45. perm(4);
46. print(" ");
47. print("You are free to call perm() and experiment with it. When finished, type \"quit\" to exit.");
48. print(" ");
RAW Paste Data