Z<-SHUFFLE_CARDS STR;CNT;N;JX;A;QIN;QV
|Shuffle pack of N cards
|Compute the cycle of shuffle
|See Ian Stewart. Scientific American , Nov 1998
QIN<-QV<-0 & STR<-"52" #ifndef "STR"
STR<-"-cycle QV<-1|-in QIN<-1" GETOPTS STR
N<-#fi STR
Z<-y(2,N%2)riN    | out shuffle
IF QIN;Z<-mZ      | in Shuffle
Z<-,Z & IF ~QV; ->0
A<-iN & CNT<-1
WHILE QV;BREAKIF (iN) #equiv A<-A[Z];CNT++;WEND
Z<-CNT

Z<-GPERM S;K;N
|Permutations of the set S
N<-rS & K<-0 & Z<-i0
WHILE K<N; Z<-Z,:S[K], GPERM (K #ne iN)/S & K++ & WEND

Z<-FC X2SPLIT X;C;J;K;N;S;T
|Components of the system i<-X[i]
|X range of function
|Use this to express permutation as a product of cycles
|if FC then give length only
 N<-(rX)c 1+c/X
 Z<-NrC<- 0 & T<-"" & S<-i0
WHILE 0<rJ<-(Z=0)/iN
C<-C+1 & J<-1tJ
BB:Z[J]<-C
IF 0=K<-Z[J<-X[J]]; ->BB
J<-(Z=C)/iN & Z[J]<- K |Adjust component number
T<-T,"(",(zJ),")"
S<-S,rJ
WEND
IF 0 < FC <-0 #ifndef "FC"; Z<-S; ->0
Z<-T

Z<-P_MULT STR;IDX;NR;TAB
|Mapping composition of permutations
Z<-#deb STR<-"" #ifndef "STR"
NR<-1tr:TAB<-#sstomat " ", STR
Z<-$NL[#deb TAB[IDX<-NR-1 :]]
WHILE 0 < IDX--; Z<-($NL[#deb TAB[IDX:]])[Z]; WEND

ZA<-JMULT STR
|Mapping composition of permutations
|Each permutation is entered as a single letter.
IF 0=r ZA<- STR<-"" #ifndef "STR"; ->0
ZA<-$NL[1tSTR] & STR<-1dSTR
WHILE 0 < r STR; ZA<-($NL[1tSTR])[ZA]; STR<-1dSTR; WEND

Z<-G_ORDER X;T;N
|Order Z of permutation X. X^Z=1
Z<-1; N<-irT<-X
WHILE ~ T #equiv N; Z++; T<-T[X]; WEND

Z<-M12 STR
|Two generators for M12
Z<-STR<-"AB" #ifndef "STR"
$NL[STR[0]]<-mi12
$NL[STR[1]]<- ,y 2 6 r (i6),6+mi6

Z<-N GEN2 STR;K
|Two generators : reverse + shuffle
N<-12 #ifndef "N"
Z<-STR<-"AB" #ifndef "STR"
K<-N%2
$NL[STR[0]]<-mi2*K
$NL[STR[1]]<- ,y (2,K) r (iK),K+miK

Z<-X CV Y;I;J;N;T
|Convolution or polynomial multiplication
Z<-iI<-0 & N<-1+c/J<-,(irX) #outer+ irY & T<-,X #outer* Y
AA: Z<-Z,+/(I=J)/T & I++ & ->(I<N)/AA

Z<-CV_XSQ Y;T
|First N terms of Y*Y
Z<-i0;K<-1
WHILE (rZ)<rY; Z<-Z,+/T*mT<-KtY; K++; WEND

Z<-CV_PRODUCT IX;N
|Product (1+x^i[0])(1+x^i[1])...
|Double precision for USA Electoral College.
|This is the real time saving part of the Algorithm
Z<-1.0 & ->(0=rIX)/0
->(1<rIX)/AA
|One term
Z<-1,((IX-1)r0),1.0 & ->0
AA: |Split
N<-(rIX)%2
Z<-(CV_PRODUCT NtIX) CV CV_PRODUCT NdIX

Z<-MAX CV_ALT_PRODUCT SIX;N
|Product (1+s[0]*x^i[0])(1+s[1]*x^i[1])...
|s[k] + or - signs
|i[k] indices
MAX<-50 #ifndef "MAX"   | Truncate series
Z<-1; IF 0 = r SIX; ->0
IF 2 = r SIX; Z<-(1, (SIX[1]-1), 1)/1 0, SIX[0]; ->0
N<-(rSIX)%2; N<-2*N%2
Z<-(MAX CV_ALT_PRODUCT N t SIX) CV MAX CV_ALT_PRODUCT N d SIX
IF MAX < rZ; Z<-MAX t Z | Truncate

Z<-N P_INV A
|Inverse of series a[0]+a[1]x+a[2]x*x....
|assume a[0]=1
Z<-1 & N<-20 #ifndef"N"
WHILE N>rZ;Z<-Z,-+/(mZ)*(rZ)t1dA;WEND

Z<-N LI_INV A;K
|Inverse of series a[0]+a[1]x+a[2]x*x....
|assume a[0]=1
|Give Z as a string of Long Int.
Z<-"1" & K<-1; N<-20 #ifndef"N"
WHILE K<N
Z<-Z," ",(- m K t1dA) R_SUM Z
K++
WEND

Z<-K R_SUM STR;X
|Sum k[i]*STR[i] where STR is Long Int vector
Z<-"0"
WHILE 0<rK; Z<-Z + (z1tK) * SHIFT; K<-1dK; WEND

Z<-PARTITION_SERIES N;NP
|Compute p[k] 0<k<=N where p[k] is number of partitions of k.
|Use Euler's generating function.
|sum (p[k]x^k) = 1/product(1-x^k)
NP<-N+8 | Take a few extra terms
Z<-NP CV_ALT_PRODUCT ,_1,y1+iNP
IF N<110; Z<-(N+1) t NP P_INV Z; ->0
Z<-NP LI_INV Z
Z<-(#sstomat " ",Z)[iN+1]

Z<-TAU N
|Ramanujan numbers
|Expand (Product (1-x^2n))^24
Z<-N CV_ALT_PRODUCT ,_1, y1+iN
Z<-N t Z CV Z CV Z      |P^3
Z<-CV_XSQ 1.0*Z
Z<-CV_XSQ Z
Z<-CV_XSQ Z

Z<-S_NEXT X;K
|Next partition of N = +/X in lexicographic order
IF 1 >= rZ<-X; ->0
IF X[_2+rX] = K<-X[0]; Z<-K, ((+/X)-K<-1+X[0])/1; ->0
Z<-K, S_NEXT 1dX;

Z<-LIST_PARTS N;X;Y
|List partitions of N as a P(N) x N matrix
IF N>25; Z<-"sanity check"; ->0
Z<-X<-Nr1
WHILE ~ X #equiv Y<- S_NEXT X
Z<-Z,:N t X<-Y
WEND

Z<-DISP_PARTS N;H;X;Y
|List partitions of N as a P(N) histograms
IF N>20; Z<-"sanity check"; ->0
Z<-X<-Nr1
WHILE ~ X #equiv Y<- S_NEXT X
H<-".O" [(N,N) t #theta(iN) #outer < X]
H #wput 0 0,rH
0 r #sleep 0 800000
X<-Y
WEND

"shuffle.txt" << EOF
    Marcus Du Sautoy and moonshine.
    Groups of order 3^3.
    Search for a function f:N->N which enumerates finite groups
of order n. In the case of Abelian groups a partial solution is
possible. The number of groups of order p^n is P(n) the partition
function. (Sloane, encyclopedia of integer sequences.)

For abelian group enumeration try dirichlet series
    Product (Zeta(ns)) n=1 to infinity.

There are also functions to evaluate coefficients of powers of
A(x)=Product(1+x^n)    B(x)=Product(1-x^n)
C(x)=Product(1+x^2n-1) D(x)=Product(1-x^2n-1)
Try (C^8-D^8)%16 against x A(x^2)

EOF

Z<-SEQ_8 STR;E;N;OFILE;T
|Products (1+x^n)^8
OFILE<-""
N<-50
STR<-"-o OFILE<-SHIFT" GETOPTS STR
IF 0<rSTR; N<-#fi STR
T<-N CV_ALT_PRODUCT ,1,y1+iN
E<-N CV_XSQ 1.0*T
E<-N CV_XSQ E
E<-N CV_XSQ E
Z<-(zy1+iN)," ",16 0 zyE
T<-N CV_ALT_PRODUCT ,_1,y1+iN
E<-N CV_XSQ 1.0*T
E<-N CV_XSQ E
E<-N CV_XSQ E
Z<-Z,16 0 zyE
T<-N CV_ALT_PRODUCT ,1,y1+2*iN
E<-N CV_XSQ 1.0*T
E<-N CV_XSQ E
E<-N CV_XSQ E
Z<-Z,16 0 zyE
T<-N CV_ALT_PRODUCT ,_1,y1+2*iN
E<-N CV_XSQ 1.0*T
E<-N CV_XSQ E
E<-N CV_XSQ E
Z<-Z,16 0 zyE


Z<-PROD_1 STR;A;B;C;D;N;OFILE
|Products A=1+x^n B=1-x^n C=1+x^2n-1 D=1+x^2n-1
OFILE<-""
N<-50
STR<-"-o OFILE<-SHIFT" GETOPTS STR
IF 0<rSTR; N<-#fi STR
A<-N CV_ALT_PRODUCT ,1,y1+iN
B<-N CV_ALT_PRODUCT ,_1,y1+iN
C<-N CV_ALT_PRODUCT ,1,y1+2*iN
D<-N CV_ALT_PRODUCT ,_1,y1+2*iN
Z<-yiN
Z<-Z,(yA),(yB),(yC),yD
Z<-zZ

Z<-M SD_PRODUCT N;K;SIG
|Use sigma function to evaluate product series
|evaluate (1-x^n)^m
M<-1 #ifndef "M"
SIG<-1dSD_SIGMA N+5
Z<-1.0
|Product 1-x^n with use logarithmic derivative
|Sign change gives 1/product 1-x^n
WHILE (rZ)<N
Z<-Z,-(M*+/Z*mKtSIG)%K<-rZ
WEND

Z<-SD_PARTITION N;K;SIG
|Use sigma function to evaluate product series
SIG<-1dSD_SIGMA N+5
Z<-1
|Product 1/1-x^n with use logarithmic derivative
WHILE (rZ)<N
Z<-Z,(+/Z*mKtSIG)%K<-rZ
WEND

Z<-SD_ODD N;T
|Sum of odd divisors of numbers up to N
T<-i1+N
T<-T*(N+1)r0 1
Z<-SD_CV T,:0,N/1

Z<-M SD_ODDP N;K;SIG
|Product (1-x^(2n-1))^m
|Product 1/1-x^(2n-1) enumerate partitions into odd numbers  (M=-1)
|Same as Product 1+x^n
M<-1 #ifndef "M"
SIG<-1 d SD_ODD N+5
Z<-1.0
WHILE (rZ)<N
Z<-Z,- (M* +/Z*mKtSIG)%K<-rZ
WEND

Z<-SD_TAU N;K;SIG
|Use sigma function to evaluate product series
SIG<-1dSD_SIGMA N+5
Z<-1.0
|Product (1-x^n)^24 with use logarithmic derivative
|Sign change gives 1/product 1-x^n
WHILE (rZ)<N
Z<-Z,-(24* +/Z*mKtSIG)%K<-rZ
WEND

Z<-NJAS STR;CMD
CMD<-"firefox http://www.research.att.com/~njas/sequences/index.html?q="
Z<-#cmd CMD, (URL_ENCODE STR), " &"

Z<-LI_TAU N;K;MA;SIG;T;V
|Use sigma function to evaluate product (1-x^n)^24
SIG<-1dSD_SIGMA N+5
K<-1
MA<-(N,30)r" "
MA[0]<-30t"1"
|Product (1-x^n)^24 with use logarithmic derivative
WHILE K<N
T<-YR_SUM m KtSIG
V<-("-24"*T)%zK
("TAU[####]= " zK) , V
MA[K]<-30tV
K++
WEND
Z<-MA

Z<-YR_SUM K;I
|Import matrix MA
I<-0; Z<-"0"
WHILE I<rK
Z<-Z+(zK[I])* #deb MA[I:]
I++
WEND

Z<-EP_TABLE STR;FW;LP;K;M;N;T;LC;X
|Tabulate coefficients of Product (1-x^n)^k for values of K
N<-50
STR<-"-n N<-#fi SHIFT" GETOPTS STR
K<-#fi STR
Z<-(zy0,iN)
WHILE 0<rK; M<-1 t K; K<-1dK
X<-M SD_PRODUCT N
X<-"/-0/0/*" #rxsubs #deb 16 0 z X
X<-(#sstomat " ",X)," "
FW<-1dr: X<-(- "-" #ne , X[:0]) m X
X<-(FWt("P^", zM)),:X
Z<-Z," ",X
WEND

Z<-SD_TABLE STR;A;B;C;D;N;OFILE
|Products A=1+x^n B=1-x^n C=1+x^2n-1 D=1+x^2n-1
OFILE<-""
N<-50
STR<-"-o OFILE<-SHIFT" GETOPTS STR
IF 0<rSTR; N<-#fi STR
B<-B_SERIES N
D<-D_SERIES N
A<-A_SERIES N
C<- #abs D
Z<-(5 z yiN)," ", 14 0 z (yA),(yB),(yC),yD

Z<-M A_SERIES N;ASUM;K
|Direct computation of (Product (1+x^n))^m
|Use logarithmic derivative and alternating sum of divisor functions
M<-1 #ifndef "M"
ASUM<-1 d SD_CV (0,Nr1 _1),: i 1+N
Z<-1.0
WHILE (rZ)<N
Z<-Z, (M* +/Z*mKtASUM)%K<-rZ
WEND

Z<-M B_SERIES N;K;SIG
|Use sigma function to evaluate series Product ((1-x^n))^m
M<-1 #ifndef "M"
SIG<-1dSD_SIGMA N
Z<-1.0
|use logarithmic derivative
|M = -1 gives 1/product 1-x^n
WHILE (rZ)<N
Z<-Z,-(M*+/Z*mKtSIG)%K<-rZ
WEND

Z<-M D_SERIES N;DSUM;K
|Direct computation of (Product (1-x^2n-1))^m
|Use logarithmic derivative and sum of odd divisor function
M<-1 #ifndef "M"
DSUM<-1 d SD_CV (0,Nr1),: ((N+1)r0 1)*i 1+N
Z<-1.0
WHILE (rZ)<N
Z<-Z, (- M* +/Z*mKtDSUM)%K<-rZ
WEND

Z<-M C_SERIES N
|Compute (Product (1+x^2n-1))^m
M<-1 #ifndef"M"
Z<-#abs M D_SERIES N

Z<-K S_OPEN X;J
|Open a series f(x)=1+a[1]x+a[2]x^2... via x->x^k
K<-2 #ifndef "K"
Z<-K/X
Z<-(rX) r Z* (rZ) r Kt1

Z<-FXY_SERIES STR;CNT;EXPR;M;T;X;Y
|Count points for which f(x,y)=0 (mod M)
Z<-1 1;M<-2;CNT<-50
EXPR<-"-cnt CNT<-#fi SHIFT" GETOPTS STR
WHILE CNT >= rZ
T<-iM,M
Y<-T[:0]; X<-T[:1]
0 r #do "T<-",EXPR
Z<-Z,+/0=T #mod M
M++
WEND