' Some basic combinatorial subroutines for VBA (Excel).
' Version: 2.7.2007, 9.2.2008.
' Author: Thomas Wieder
' http://www.thomas-wieder.privat.t-online.de/default.html
' http://homepages.tu-darmstadt.de/%7Ewieder/Welcome.html
' A totally free library.
'
' Some subroutines are taken from the free Fortran library "combo".
' "combo" was written by John Burkardt and can be found at: 
' http://orion.math.iastate.edu/burkardt/
' Some subroutines are taken from the free Fortran library "subset".
' "subset" was written by John Burkardt and can be found at: 
' http://orion.math.iastate.edu/burkardt/
' I converted these Fortran subroutines to VBA by hand.
'
' All subroutines listed here are provided with no guarantee at all.
' Use these subroutines at your own risk.
'
' Reason for this library:
' VBA comes with no functions for many basic combinatorial quantities
' like the Stirling numbers and the Bell numbers.
' Unfortunately, only very few of them are compiled here.
' VBA = Visual Basic for Applications 

Option Explicit



'1111111111111111111111111111111111111111111111111111111111111111111111111111111


Sub TestBell()
Dim i As Long, n As Long, Result As Long
n = 9
For i = 1 To n
Result = Bell(i)
Debug.Print "TestBell:", n, Result
'Cells(5, 3 + i) = Result
Next i
End Sub

Function Bell(n As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Berechne die Bellschen Zahlen.

Dim imsgbox As Integer
Dim i As Long

If n > 2147483648# Then
imsgbox = MsgBox("n needs to be smaller than 2147483648 !", vbOKOnly, "Bell")
End
End If
If (n < 0) Then
imsgbox = MsgBox("n needs to be greater than 0 !", vbOKOnly, "Bell")
End
End If

Bell = 0

For i = 1 To n
Bell = Bell + Stirling2(n, i)
Next i

End Function


'2222222222222222222222222222222222222222222222222222222222222222222222222222222


Sub Test_Stirling2()
Dim n As Long, k As Long, dummy As Long, dummy2 As Long
k = 3
For n = 1 To 10
dummy = Stirling2(n, k)
Debug.Print n, k, dummy
Next n
End Sub

Function Stirling2(n As Long, k As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Berechne die Stirlingschen Zahlen der zweiten Art.

    Dim i As Long, Summe As Long
    Summe = 0
    For i = 0 To k
    Summe = Summe + ((-1) ^ i) * binomial(k, i) * (k - i) ^ n
    Next i
Stirling2 = (1 / faculty(k)) * Summe
End Function


' ----------------------------------------------------------------------------


Sub test_stirling1()
Dim Resultat As Long
Resultat = Stirling1(5, 2)
Debug.Print Resultat
End Sub

Function Stirling1(m As Long, n As Long)
'
'*******************************************************************************
'
'! STIRLING_NUMBERS1 computes Stirling numbers of the first kind.
'
'
'  Reference:
'
'    Algorithm 3.17,
'    Donald Kreher and Douglas Simpson,
'    Combinatorial Algorithms,
'    CRC Press, 1998, page 89.
'
'  Modified:
'
'    20 January 1999
'
'  Author:
'
'    John Burkardt
'
'  Parameters:
'
'    Input, integer M, the maximum row to compute.
'    M must be nonnegative.
'
'    Input, integer N, the maximum column to compute.
'    N must be nonnegative.
'
'    Output, integer S(0:M,0:N), the first M+1 rows and N+1 columns
'    of the table of Stirling numbers of the first kind.
'
'  implicit none
'
  Dim i As Long, j As Long, s(0 To 100, 0 To 100) As Long
'
  s(0, 0) = 1
  For i = 1 To m
  s(i, 0) = 0
  Next i
'
'  This loop may be extraneous.
'
  For i = 0 To longMinimum(m, n - 1)
    s(i, i + 1) = 0
  Next i

  For i = 1 To m
    For j = 1 To n
      If (j <= i) Then
        s(i, j) = s(i - 1, j - 1) - (i - 1) * s(i - 1, j)
      Else
        s(i, j) = 0
      End If
    Next j
  Next i

' Wieder
Stirling1 = s(m, n)

End Function


'3333333333333333333333333333333333333333333333333333333333333333333333333333333


Sub Test_binomial()
Dim n As Long, dummy As Long
n = 6
dummy = binomial(n, 4)
Debug.Print n, dummy
End Sub

Function binomial(m As Long, n As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Berechne die Binomialzahlen.

Dim Result As Long
Result = faculty(m) / (faculty(m - n) * faculty(n))
binomial = Result
End Function

Public Function Binomialkoeffizient(n As Long, k As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Berechne die Binomialzahlen.

Dim imsgbox As Long

If n < 0 Then
imsgbox = MsgBox("n kleiner 0! Stop!", vbOKOnly, "Binomialkoeffizient")
Stop
ElseIf k < 0 Then
imsgbox = MsgBox("k kleiner 0! Stop!", vbOKOnly, "Binomialkoeffizient")
Stop
ElseIf n < k Then
Stop
imsgbox = MsgBox("n kleiner k! Stop!", vbOKOnly, "Binomialkoeffizient")
End If

Binomialkoeffizient = Fakultät(n) / (Fakultät(n - k) * Fakultät(k))

End Function


'4444444444444444444444444444444444444444444444444444444444444444444444444444444


Sub Test_faculty()
Dim dummy As Long
dummy = faculty(6)
Debug.Print dummy
End Sub

Function faculty(n As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Berechne die Fakultät (ohne Rekursion).

Dim Result As Long, i As Long
Result = 1
For i = 2 To n
Result = Result * i
Next i
faculty = Result
End Function

Function Fakultät(n As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Berechne die Fakultät rekursiv.

    If n <= 1 Then    ' Rekursionsende ist erreicht.
        Fakultät = 1    ' (N = 0) Keine weiteren Aufrufe.
    Else ' Fakultät wird erneut aufgerufen
          ' wenn N > 0.
        Fakultät = Fakultät(n - 1) * n
    End If
End Function


'5555555555555555555555555555555555555555555555555555555555555555555555555555555


Sub TestZahlPartitonenInTeile()
Dim n As Long, k As Long, Resultat As Long
n = 8
k = 4
Resultat = ZahlPartitionen(n, k)
Debug.Print "TestZahlPartitonen: n,k,Resultat:", n, k, Resultat
End Sub

Public Function ZahlPartitionenInTeile(n As Long, k As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Berechne rekursiv die Zahl der ganzahligen Partitionen von n in k Teile.

Dim imsgbox As Integer

If n > 2147483648# Or k > 2147483648# Then
imsgbox = MsgBox("n and k need to be smaller than 2147483648 !", vbOKOnly, "ZahlPartitionen")
End
End If
If (n < 0 Or k < 0) Then
imsgbox = MsgBox("n and k need to be greater than 0 !", vbOKOnly, "ZahlPartitionen")
End
End If
'If k > n Then
'imsgbox = MsgBox("k needs to be <= n !", vbOKOnly, "ZahlPartitionen")
'End
'End If

If k = 1 Then
ZahlPartitionen = 1
Exit Function
ElseIf k = n Then
ZahlPartitionen = 1
Exit Function
ElseIf k > n Then
ZahlPartitionen = 0
Exit Function
End If

ZahlPartitionen = ZahlPartitionen(n - 1, k - 1) + ZahlPartitionen(n - k, k)

End Function


'6666666666666666666666666666666666666666666666666666666666666666666666666666666


Sub test_Produkt()
Dim Faktoren(10) As Long, Ergebnis As Long, ierror As Long
Faktoren(1) = 1
Faktoren(2) = 2
Faktoren(3) = 3
Faktoren(4) = 4
Faktoren(5) = 5
Ergebnis = longProdukt(1, 5, Faktoren, 2, ierror)
Debug.Print "Ergebnis=", Ergebnis, ierror
End Sub

Function longProdukt(UntereGrenze As Long, ObereGrenze As Long, Faktoren() As Long, _
imodus As Long, ierror As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Berechne ein Produkt.

Dim i As Long, Ergebnis As Long

ierror = 0
If ObereGrenze < UntereGrenze Then
ierror = 1
Debug.Print "Function Product: UntereGrenze, ObereGrenze", UntereGrenze, ObereGrenze
Exit Function
End If

Ergebnis = 1
For i = UntereGrenze To ObereGrenze
If imodus = 1 Then
Ergebnis = Ergebnis * Faktoren(i)
ElseIf imodus = 2 Then
If Faktoren(i) <> 0 Then Ergebnis = Ergebnis * Faktoren(i)
End If
Next i

longProdukt = Ergebnis

End Function


'7777777777777777777777777777777777777777777777777777777777777777777777777777777


Function longSum(UntereGrenze As Long, ObereGrenze As Long, Summanden() As Long, ierror As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Berechne eine Summe.

Dim i As Long, Ergebnis As Long
Dim iverbose As Integer
Static zeile

iverbose = 0
ierror = 0
zeile = zeile + 1

If ObereGrenze < UntereGrenze Then
ierror = 1
Debug.Print "Function longSum: UntereGrenze, ObereGrenze", UntereGrenze, ObereGrenze
Exit Function
End If

Ergebnis = 0
For i = UntereGrenze To ObereGrenze
If iverbose = 1 Then
Debug.Print "longSum: i, summand(i)", i, Summanden(i)
zeile = zeile + 1
Cells(zeile, 12 + i) = Summanden(i)
End If
Ergebnis = Ergebnis + Summanden(i)
Next i

longSum = Ergebnis

End Function


'8888888888888888888888888888888888888888888888888888888888888888888888888888888


Function kronecker_delta(arg1, arg2)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Berechne das Kronecker-delta.

If arg1 = arg2 Then
kronecker_delta = 1
Else
kronecker_delta = 0
End If

End Function


'9999999999999999999999999999999999999999999999999999999999999999999999999999999


Sub test_ListeAllePartitionenVon_n_Auf()

Dim index_partitionen As Long, index_part As Long, ZahlPartitionen As Long
Dim ZahlTeileDerPartitionen(1000) As Long, ListeDerTeileDerPartitionen(100, 1000) As Long

Call ListeAllePartitionenVon_n_Auf(5, ZahlPartitionen, ZahlTeileDerPartitionen, ListeDerTeileDerPartitionen)

For index_partitionen = 1 To ZahlPartitionen
Debug.Print "Nächste Partition:"
For index_part = 1 To ZahlTeileDerPartitionen(index_partitionen)
Debug.Print "test_ListeAllePartitionenVon_n_Auf: i-te Partition, Teile", index_partitionen, ListeDerTeileDerPartitionen(index_partitionen, index_part)
Next index_part
Next index_partitionen

End Sub

Sub ListeAllePartitionenVon_n_Auf( _
n As Long, ZahlPartitionen As Long, ZahlTeileDerPartitionen() As Long, ListeDerTeileDerPartitionen() As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Liste alle ganzzahligen Zerlegungen auf.

Dim iverbose As Integer
Dim done As Boolean
Dim A(100) As Long, mult(100) As Long, npart As Long, ndiffpart As Long
Dim index As Long, index_hilf As Long, index_partition As Long, index_part As Long

index_partition = 0

While done = False

If index_partition = 0 Then done = True

Call i_partition_next(n, npart, A, ndiffpart, mult, done)

If (done) Then
ZahlPartitionen = index_partition
Exit Sub
End If

index_partition = index_partition + 1

' Call i_partition_print(n, npart, a, mult)
Debug.Print ' '
index_part = 0
For index = 1 To npart
If iverbose = 1 Then
Debug.Print "ListeAllePartitionenVon_n_Auf: index_partition, A(index), mult(index)", index_partition, A(index), mult(index)
End If
For index_hilf = 1 To mult(index)
index_part = index_part + 1
' Liste der Teile enthält die Partition in ausführlicher Form.
ListeDerTeileDerPartitionen(index_partition, index_part) = A(index)
ZahlTeileDerPartitionen(index_partition) = index_part
Next index_hilf
Next index

Wend

Exit Sub

End

End Sub

Sub i_partition_next(n As Long, npart As Long, A() As Long, ndiffpart As Long, mult() As Long, done As Boolean)

'!*******************************************************************************
'!
'!! I_PARTITION_NEXT generates the partitions of an intiger, one at a time.
'!
'!
'!  Modified:
'!
'!    15 April 1999
'!
'!  Author:
'!
'!    Albert Nijenhuis, Herbert Wilf
'!
'!    FORTRAN90 version by John Burkardt
'!
'!  Reference:
'!
'!    Albert Nijenhuis, Herbert Wilf,
'!    Combinatorial Algorithms,
'!    Academic Press, 1978, second edition,
'!    ISBN 0-12-519260-6.
'!
'!  Parameters:
'!
'!    Input, long N, the integer to be partitioned.
'!
'!    Input/output, integer NPART, the number of "parts" in the partition.
'!
'!    Input/output, integer A(N).  A contains the parts of
'!    the partition.  The value of N is represented by
'!
'!      N = sum ( 1 <= I <= NPART ) MULT(I) * A(I).
'!
'!    Input/output, integer MULT(N).  MULT counts the multiplicity of
'!    the parts of the partition.  MULT(I) is the multiplicity
'!    of the part A(I), for I = 1 to NPART.
'!
'!    Input/output, logical DONE.
'!    On first call, the user should set DONE to TRUE to signal
'!    that the program should initialize data.
'!    On each return, the programs sets DONE to FALSE if it
'!    has another partition to return.  If the program returns
'!    with DONE TRUE, then there are no more partitions.
'!
'  implicit none

  Dim iss As Long
  Dim iu As Long
  Dim iv As Long
  Dim iw As Long
  Dim k As Long
  Dim k1 As Long
  Dim index As Long

  If (n <= 0) Then
  Debug.Print ' '
  Debug.Print 'I_PARTITION_NEXT - Fatal error!'
  Debug.Print '  N must be positive.'
  Debug.Print '  The input value of N was ', n
  Stop
  End If

  If (done = True) Then

    A(1) = n
    mult(1) = 1
    For index = 2 To n
    A(index) = 0
    mult(index) = 0
    Next index
    npart = 1
    done = False

  Else

    If (1 < A(npart) Or 1 < npart) Then

      done = False

      If (A(npart) = 1) Then
        iss = A(npart - 1) + mult(npart)
        k = npart - 1
      Else
        iss = A(npart)
        k = npart
      End If

      iw = A(k) - 1
      iu = iss \ iw
      iv = modulo(iss, iw)
      mult(k) = mult(k) - 1

      If (mult(k) = 0) Then
        k1 = k
      Else
        k1 = k + 1
      End If

      mult(k1) = iu
      A(k1) = iw
      
      'Wieder
      ndiffpart = k1

      If (iv = 0) Then
        npart = k1
      Else
        mult(k1 + 1) = 1
        A(k1 + 1) = iv
        npart = k1 + 1
      End If

    Else
      done = True
    End If

  End If

  Exit Sub
  
End Sub


'10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 


Sub test_npart_rsf_lex_successor()

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Aufruf des Programmes *** npart_rsf_lex_successor *** von John Burkardt.

Dim n As Long, npart As Long, rank As Long, rankold As Long, ihilf As Long
Dim A(100) As Long
n = 12
npart = 3
rank = -1
Do
    Debug.Print "Next Partition:"
    rankold = rank
    Call npart_rsf_lex_successor(n, npart, A, rank)
    If (rank <= rankold) Then
      Exit Do
    End If
    For ihilf = 1 To npart
    Debug.Print rank, A(ihilf)
    Next ihilf
Loop
End Sub

Sub npart_rsf_lex_successor(n As Long, npart As Long, A() As Long, rank As Long)

'!*******************************************************************************
'!
'!! NPART_RSF_LEX_SUCCESSOR computes the RSF lex successor for NPART partitions.
'!
'!  Modified:
'!
'!    28 February 2001
'!
'!  Author:
'!
'!    John Burkardt
'!
'!  Reference:
'!
'!    Donald Kreher and Douglas Stinson,
'!    Algorithm 3.7,
'!    Combinatorial Algorithms,
'!    CRC Press, 1998, page 76.
'!
'!  Parameters:
'!
'!    Input, integer N, the integer to be partitioned.
'!    N must be at least 1.
'!
'!    Input, integer NPART, the number of parts of the partition.
'!    1 <= NPART <= N.
'!
'!    Input/output, integer A(NPART), contains the partition.
'!    A(1) through A(NPART) contain the nonzero longs which
'!    sum to N.
'!
'!    Input/output, integer RANK, the rank.
'!
'!    If RANK = -1 on input, then the routine understands that this is
'!    the first call, and that the user wishes the routine to supply
'!    the first element in the ordering, which has RANK = 0.
'!
'!    In general, the input value of RANK is increased by 1 for output,
'!    unless the very last element of the ordering was input, in which
'!    case the output value of RANK is 0.
'!
'  implicit none

  Dim d As Long, i As Long, ierror As Long, j As Long, ihilf As Long
  Dim imsgbox As Long, mldg As String

On Error GoTo Fehler

'!
'!  Return the first element.
'!
  If (rank = -1) Then

    If (npart < 1) Then
      Debug.Print " "
      Debug.Print "NPART_RSF_LEX_SUCCESSOR - Fatal error!"
      Debug.Print "  NPART < 1."
      Stop
    End If
    
    For ihilf = 1 To npart - 1
    A(ihilf) = 1
    Next ihilf
    A(npart) = n - (npart - 1)

    rank = 0
    Exit Sub    'Return

  End If
'!
'!  Check.
'!
  Call part_rsf_check(n, npart, A, ierror)

  If (ierror <> 0) Then
    Debug.Print " "
    Debug.Print "NPART_RSF_LEX_SUCCESSOR - Fatal error!"
    Debug.Print "  The input array is illegal."
    Debug.Print "  IERROR = ", ierror
    Stop
  End If
'!
'!  Find the first index I for which A(NPART+1-I) + 1 < A(NPART).
'!
  i = 2

  Do
  
    If (npart < i) Then
      Exit Do
    End If

    If (A(npart + 1 - i) + 1 < A(npart)) Then
      Exit Do
    End If

    i = i + 1

  Loop
'!
'!  If no such index, we've reached the end of the line.
'!
  If (i = npart + 1) Then
    
    For ihilf = 1 To npart - 1
    A(ihilf) = 1
    Next ihilf
    A(npart) = n - (npart - 1)

    rank = 0
    Exit Sub 'Return
'!
'!  Otherwise, increment A(NPART+1-I), and adjust other entries.
'!
  Else

    A(npart + 1 - i) = A(npart + 1 - i) + 1
    d = -1

    For j = i - 1 To 2 Step -1
      d = d + A(npart + 1 - j) - A(npart + 1 - i)
      A(npart + 1 - j) = A(npart + 1 - i)
    Next j

    A(npart) = A(npart) + d

  End If

  rank = rank + 1

  Exit Sub
  
Fehler:
imsgbox = MsgBox("Ein Fehler ist aufgetreten!", vbOKOnly, "npart_ref_lex_succesor")
Debug.Print ("Fehler in SumProdBinomial")
If Err.Number <> 0 Then
    mldg = "Fehler # " & Str(Err.Number) & " wurde ausgelöst von " _
            & Err.Source & Chr(13) & Err.Description
    MsgBox mldg, , "Fehler", Err.HelpFile, Err.HelpContext
End If

End Sub

Sub part_rsf_check(n As Long, npart As Long, A() As Long, ierror As Long)

'!*******************************************************************************
'!
'!! PART_RSF_CHECK checks a reverse standard form partition of an long.
'!
'!  Modified:
'!
'!    24 January 1999
'!
'!  Author:
'!
'!    John Burkardt
'!
'!  Reference:
'!
'!    Donald Kreher and Douglas Stinson,
'!    Combinatorial Algorithms,
'!    CRC Press, 1998.
'!
'!  Parameters:
'!
'!    Input, integer N, the long to be partitioned.
'!    N must be positive.
'!
'!    Input, integer NPART, the number of parts of the partition.
'!    1 <= NPART <= N.
'!
'!    Input, integer A(NPART), contains the partition.
'!    A(1) through A(NPART) contain the nonzero integers which
'!    sum to N.  The entries must be in ASCENDING order.
'!
'!    Output, long IERROR, error flag.
'!    0, no error.
'!    -1, N is illegal.
'!    -2, NPART is illegal.
'!    -3, the entries do not add up to N.
'!    I, the I-th entry of A is illegal.
'!
'  implicit none

  Dim i As Long, Summe As Long, imsgbox As Long, mldg As String

On Error GoTo Fehler

  ierror = 0

  If (n < 1) Then
    ierror = -1
  Exit Sub ' Return
  End If

  If (npart < 1 Or n < npart) Then
    ierror = -2
    Exit Sub
  End If
'!
'!  Every entry must lie between 1 and N.
'!
  For i = 1 To npart
    If (A(i) < 1 Or n < A(i)) Then
      ierror = i
    Exit Sub
    End If
  Next i
'!
'!  The entries must be in ascending order.
'!
  For i = 2 To npart
    If (A(i) < A(i - 1)) Then
      ierror = i
      Exit Sub
    End If
  Next i
'!
'!  The entries must add up to N.
'!
  Summe = 0
  For i = 1 To npart
  Summe = Summe + A(i)
  Next i
  If (Summe <> n) Then
    ierror = -3
  End If

  Exit Sub
  
Fehler:
imsgbox = MsgBox("Ein Fehler ist aufgetreten!", vbOKOnly, "part_rsf_checl")
Debug.Print ("Fehler in SumProdBinomial")
If Err.Number <> 0 Then
    mldg = "Fehler # " & Str(Err.Number) & " wurde ausgelöst von " _
            & Err.Source & Chr(13) & Err.Description
    MsgBox mldg, , "Fehler", Err.HelpFile, Err.HelpContext
End If
  
End Sub


'11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 


Sub test_perm_next()

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Aufruf des Programmes *** perm_next *** von John Burkardt.

Dim n As Long, P(1000) As Long, more As Boolean, even As Boolean, i As Long, j As Long

n = 3
more = False
i = 0
Do
Call perm_next(n, P(), more, even)
Debug.Print "Next Permutation:"
i = i + 1
For j = 1 To n
Debug.Print P(j)
Cells(i, j) = P(j)
Next j
Loop Until more = False

End Sub

Sub perm_next(n As Long, P() As Long, more As Boolean, even As Boolean)

'!*******************************************************************************
'!
'!! PERM_NEXT computes all of the permutations of N objects, one at a time.
'!
'!  Discussion:
'!
'!    The routine is initialized by calling with MORE = TRUE, in which case
'!    it returns the identity permutation.
'!
'!    If the routine is called with MORE = FALSE, then the successor of the
'!    input permutation is computed.
'!
'!  Modified:
'!
'!    30 March 2001
'!
'!  Author:
'!
'!    Albert Nijenhuis, Herbert Wilf
'!
'!    FORTRAN90 version by John Burkardt
'!
'!  Reference:
'!
'!    Albert Nijenhuis, Herbert Wilf,
'!    Combinatorial Algorithms,
'!    Academic Press, 1978, second edition,
'!    ISBN 0-12-519260-6.
'!
'!  Parameters:
'!
'!    Input, integer N, the number of objects being permuted.
'!
'!    Input/output, integer P(N), the permutation, in standard index form.
'!    On the first call, the input value is unimportant.
'!    On subsequent calls, the input value should be the same
'!    as the output value from the previous call.  In other words, the
'!    user should just leave P alone.
'!    On output, contains the "next" permutation.
'!
'!    Input/output, logical MORE.
'!    Set MORE = FALSE before the first call.
'!    MORE will be reset to TRUE and a permutation will be returned.
'!    Each new call produces a new permutation until
'!    MORE is returned FALSE.
'!
'!    Input/output, logical EVEN.
'!    The input value of EVEN should simply be its output value from the
'!    previous call; (the input value on the first call doesn't matter.)
'!    On output, EVEN is TRUE if the output permutation is even, that is,
'!    involves an even number of transpositions.
'!
'  implicit none

  Dim i As Long, i1 As Long, ia As Long, id As Long, ist As Long, j As Long, l As Long, m  As Long

  If Not more Then

    Call ivec_indicator(n, P)
    more = True
    even = True

    If (n = 1) Then
      more = False
'      Return
    Exit Sub
    End If

    If (P(n) <> 1 Or P(1) <> 2 + modulo(n, 2)) Then
'      Return
     Exit Sub
    End If

    For i = 1 To n - 3
      If (P(i + 1) <> P(i) + 1) Then
        Return
      End If
    Next i

    more = False

  Else

    If (n = 1) Then
      P(1) = 0
      more = False
'      Return
    Exit Sub
    End If

    If (even) Then

      ia = P(1)
      P(1) = P(2)
      P(2) = ia
      even = False

      If (P(n) <> 1 Or P(1) <> 2 + modulo(n, 2)) Then
'        Return
        Exit Sub
      End If

      For i = 1 To n - 3
        If (P(i + 1) <> P(i) + 1) Then
'          Return
        Exit Sub
        End If
      Next i

      more = False
'      Return
      Exit Sub

    Else

      more = False

      ist = 0

      For i1 = 2 To n

        ia = P(i1)
        i = i1 - 1
        id = 0

        For j = 1 To i
          If (ia < P(j)) Then
            id = id + 1
          End If
        Next j

        ist = id + ist
        If (id <> i * modulo(ist, 2)) Then
          more = True
          Exit For
        End If

      Next i1

      If Not more Then
        P(1) = 0
        Return
      End If

    End If

    m = modulo(ist + 1, 2) * (n + 1)

    For j = 1 To i

'      if ( sign ( 1, p(j)-ia ) /= sign ( 1, p(j)-m ) ) then
       If (Sgn(P(j) - ia) <> Sgn(P(j) - m)) Then
        m = P(j)
        l = j
      End If

    Next j

    P(l) = ia
    P(i1) = m
    even = True

  End If

'  Return
End Sub

Sub ivec_indicator(n As Long, A() As Long)

'!*******************************************************************************
'!
'!! IVEC_INDICATOR sets an integer vector to the indicator vector.
'!
'!  Modified:
'!
'!    09 November 2000
'!
'!  Author:
'!
'!    John Burkardt
'!
'!  Parameters:
'!
'!    Input, integer N, the number of elements of A.
'!
'!    Output, integer A(N), the array to be initialized.
'!
'  implicit none
  
  Dim i As Long

  For i = 1 To n
    A(i) = i
  Next i

'  Return
End Sub


'12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 


Sub test_perm_index_to_cycle()

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Aufruf des Programmes *** perm_index_to_cycle *** von John Burkardt.

Dim j As Long, n As Long, p1(10) As Long, p2(10) As Long

n = 9
p1(1) = 2
p1(2) = 3
p1(3) = 9
p1(4) = 6
p1(5) = 7
p1(6) = 8
p1(7) = 5
p1(8) = 4
p1(9) = 1

Call perm_index_to_cycle(n, p1, p2)

For j = 1 To n
Cells(1, j) = p2(j)
Next j

End Sub

Sub perm_index_to_cycle(n As Long, p1() As Long, p2() As Long)

'!*******************************************************************************
'!
'!! PERM_INDEX_TO_CYCLE converts a permutation from standard index to cycle form.
'!
'!  Example:
'!
'!    Input:
'!
'!n = 9
'!      P1 = 2, 3, 9, 6, 7, 8, 5, 4, 1
'!
'!    Output:
'!
'!      P2 = -1, 2, 3, 9, -4, 6, 8, -5, 7
'!
'!  Modified:
'!
'!    31 July 2000
'!
'!  Author:
'!
'!    John Burkardt
'!
'!  Reference:
'!
'!    Albert Nijenhuis, Herbert Wilf,
'!    Combinatorial Algorithms,
'!    Academic Press, 1978, second edition,
'!    ISBN 0-12-519260-6.
'!
'!  Parameters:
'!
'!    Input, integer N, the number of objects being permuted.
'!
'!    Input, integer P1(N), the permutation, in standard index form.
'!
'!    Output, integer P2(N), the permutation, in cycle form.
'!
'  implicit none

  Dim i As Long, j As Long, k As Long

  i = 0
  j = 1

  Do While (j <= n)

    If (p1(j) < 0) Then

      j = j + 1

    Else

      k = j

      i = i + 1
      p2(i) = -k

      Do While (p1(k) <> j)
        i = i + 1
        p2(i) = p1(k)
        p1(k) = -p1(k)
        k = Abs(p1(k))
      Loop

      p1(k) = -p1(k)

    End If

  Loop

  For i = 1 To n
  p1(i) = Abs(p1(i))
  Next i

  Exit Sub
End Sub


'13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 


Sub test_perm_cycle()

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Aufruf des Programmes *** perm_cycle *** von John Burkardt.

Dim n As Long, iopt As Long, P(1000) As Long, isgn As Long, ncycle As Long, i As Long
Dim LengthOfCycle(1000) As Long, MultiplicitiesOfCycles(1000) As Long
iopt = 0
n = 9

P(1) = 2
P(2) = 3
P(3) = 9
P(4) = 6
P(5) = 7
P(6) = 8
P(7) = 5
P(8) = 4
P(9) = 1

Call perm_cycle(n, iopt, P(), isgn, ncycle, LengthOfCycle(), MultiplicitiesOfCycles())
Debug.Print "Zahl Zyklen =", ncycle

For i = 1 To n
Debug.Print i, LengthOfCycle(i), MultiplicitiesOfCycles(i)
Next i

P(1) = 1
P(2) = 2
P(3) = 3
P(4) = 4
P(5) = 5
P(6) = 6
P(7) = 7
P(8) = 8
P(9) = 9

Call perm_cycle(n, iopt, P(), isgn, ncycle, LengthOfCycle(), MultiplicitiesOfCycles())
Debug.Print "Zahl Zyklen =", ncycle

For i = 1 To n
Debug.Print i, LengthOfCycle(i), MultiplicitiesOfCycles(i)
Next i

P(1) = 9
P(2) = 8
P(3) = 7
P(4) = 6
P(5) = 5
P(6) = 4
P(7) = 3
P(8) = 2
P(9) = 1

Call perm_cycle(n, iopt, P(), isgn, ncycle, LengthOfCycle(), MultiplicitiesOfCycles())
Debug.Print "Zahl Zyklen =", ncycle

For i = 1 To n
Debug.Print i, LengthOfCycle(i), MultiplicitiesOfCycles(i)
Next i

End Sub

Sub perm_cycle(n As Long, iopt As Long, P() As Long, isgn As Long, ncycle As Long, _
LengthOfCycle() As Long, MultiplicitiesOfCycles() As Long)

'!*******************************************************************************
'!
'!! PERM_CYCLE analyzes a permutation.
'!
'!  Discussion:
'!
'!    The routine will count cycles, find the sign of a permutation,
'!    and tag a permutation.
'!
'!  Example:
'!
'!    Input:
'!
'!n = 9
'!iopt = 1
'!      P = 2, 3, 9, 6, 7, 8, 5, 4, 1
'!
'!    Output:
'!
'!ncycle = 3
'!isgn = 1
'!      P = -2, 3, 9, -6, -7, 8, 5, 4, 1
'!
'!  Modified:
'!
'!    15 April 1999
'!
'!  Author:
'!
'!    Albert Nijenhuis, Herbert Wilf
'!
'!    FORTRAN90 version by John Burkardt
'!
'!  Reference:
'!
'!    Albert Nijenhuis, Herbert Wilf,
'!    Combinatorial Algorithms,
'!    Academic Press, 1978, second edition,
'!    ISBN 0-12-519260-6.
'!
'!  Parameters:
'!
'!    Input, integer N, the number of objects being permuted.
'!
'!    Input, integer IOPT, requests tagging.
'!    0, the permutation will not be tagged.
'!    1, the permutation will be tagged.
'!
'!    Input/output, integer P(N).  On input, P describes a
'!    permutation, in the sense that entry I is to be moved to P(I).
'!    If IOPT = 0, then P will not be changed by this routine.
'!    If IOPT = 1, then on output, P will be "tagged".  That is,
'!    one element of every cycle in P will be negated.  In this way,
'!    a user can traverse a cycle by starting at any entry I1 of P
'!    which is negative, moving to I2 = ABS(P(I1)), then to
'!    P(I2), and so on, until returning to I1.
'!
'!    Output, integer ISGN, the "sign" of the permutation, which is
'!    +1 if the permutation is even, -1 if odd.  Every permutation
'!    may be produced by a certain number of pairwise switches.
'!    If the number of switches is even, the permutation itself is
'!    called even.
'!
'!    Output, integer NCYCLE, the number of cycles in the permutation.
'!
'  implicit none

Dim i As Long, i1 As Long, i2 As Long, ierror As Long, ist As Long
Dim icycle As Long, ilength As Long, ncycle_old As Long

For ilength = 1 To n
LengthOfCycle(ilength) = 1
MultiplicitiesOfCycles(ilength) = 0
Next ilength

Call perm_check(n, P, ierror)

If (ierror <> 0) Then
  Debug.Print "PERM_CYCLE - Fatal error!"
  Debug.Print "The input array does not represent a proper permutation. In particular, the array is missing the value ", ierror
  Stop
End If

ist = 1
ncycle = n
icycle = 1
ncycle_old = n

For i = 1 To n

''  Test
'   If p(i) < 0 Then
'   Debug.Print i, p(i)
'   End If

  i1 = P(i)

  Do While (i < i1)
    ncycle = ncycle - 1
    i2 = P(i1)
    P(i1) = -i2
    i1 = i2
    
  LengthOfCycle(icycle) = LengthOfCycle(icycle) + 1
  
  Loop
    
    For ilength = 1 To n
        If LengthOfCycle(icycle) = ilength Then
        MultiplicitiesOfCycles(ilength) = MultiplicitiesOfCycles(ilength) + 1
        End If
    Next ilength
    
    If ncycle <> ncycle_old Then
    MultiplicitiesOfCycles(1) = MultiplicitiesOfCycles(1) - 1
    icycle = icycle + 1
    ncycle_old = ncycle
    End If

  If (iopt <> 0) Then
    ist = -Sgn(P(i))
  End If

'    p(i) = sign(p(i), ist)
   P(i) = -Sgn(ist) * P(i)

Next i

isgn = 1 - 2 * modulo(n - ncycle, 2)

  MultiplicitiesOfCycles(1) = ncycle - icycle + 1
  For ilength = ncycle + 1 To n
  LengthOfCycle(ilength) = 0
  Next ilength

Exit Sub
End Sub

Sub perm_check(n As Long, P() As Long, ierror As Long)

'!*******************************************************************************
'!
'!! PERM_CHECK checks that a vector represents a permutation.
'!
'!  Discussion:
'!
'!    The routine verifies that each of the longs from 1
'!    to N occurs among the N entries of the permutation.
'!
'!  Modified:
'!
'!    06 August 2000
'!
'!  Author:
'!
'!    John Burkardt
'!
'!  Parameters:
'!
'!    Input, integer N, the number of entries.
'!
'!    Input, integer P(N), the permutation, in standard index form.
'!
'!    Output, integer IERROR, error flag.
'!    0, the array does represent a permutation.
'!    nonzero, the array does not represent a permutation.  The smallest
'!    missing value is equal to IERROR.
'!
'  implicit none

  Dim ifind  As Long, iseek  As Long
  
  ierror = 0

  For iseek = 1 To n

    ierror = iseek

    For ifind = 1 To n
      If (P(ifind) = iseek) Then
        ierror = 0
        Exit For
      End If
    Next ifind

    If (ierror <> 0) Then
      Exit Sub
    End If

  Next iseek

  Exit Sub
End Sub


' 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 


Sub TestMultiplizitäten()

Dim i As Long, j As Long
Dim Feld(12) As Long, Multiplizität(200, 200, 2) As Long, LängeMultiplizität As Long

Feld(1) = 1
Feld(2) = 1
Feld(3) = 2
Feld(4) = 5
Feld(5) = 5
Feld(6) = 5
Feld(7) = 2
Feld(8) = 1
Feld(9) = 4
Feld(10) = 8
Feld(11) = 2
Feld(12) = 10

i = 1
Call Multiplizitäten(i, Feld(), 12, Multiplizität(), LängeMultiplizität)

For j = 1 To 12
Debug.Print i, "TestMultiplizitäten", Multiplizität(i, j, 1), Multiplizität(i, j, 2), LängeMultiplizität
Next j

End Sub

Sub Multiplizitäten(iprttn As Long, Feld() As Long, LängeFeld As Long, Multiplizität() As Long, LängeMultiplizität As Long)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' Bestimme die Multiplizitäten eines Feldes von ganzen Zahlen.

Dim h As Long, i As Long, j As Long
Dim Hilfe(200) As Long
Dim MinWert As Long, MaxWert As Long

LängeMultiplizität = 200
MinWert = 0

MaxWert = LängeFeld
For i = 1 To MaxWert
Hilfe(i) = 0
Next i

For j = 1 To MaxWert
For i = 1 To LängeFeld
If Feld(i) = j Then
Hilfe(j) = Hilfe(j) + 1
' Ende Abfrage *** Feld(i) = i ***.
End If
Next i
Next j

j = 0
For i = 1 To MaxWert
If Hilfe(i) <> 0 Then
j = j + 1
LängeMultiplizität = j
' Speichere den Wert
Multiplizität(iprttn, j, 1) = i
' Speichere die Multiplizität des Werts.
Multiplizität(iprttn, j, 2) = Hilfe(i)
End If
Next i

End Sub


' 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 


Public Function modulo(Zahl1 As Variant, Zahl2 As Variant)

' Begonnen am: 1.7.2007
' Letzte Änderung am: 1.7.2007
' thomas.wieder@t-online.de
' zahl1 modulo zahl2 = Rest von zahl1 - INT(zahl1/zahl2)*zahl2

If Zahl2 = 0 Then
modulo = 0
Exit Function
End If

modulo = Zahl1 - Int(Zahl1 / Zahl2) * Zahl2

End Function


' 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 


Sub TEST028a()

'*******************************************************************************
'
'! TEST028a tests EQUIV_NEXT2.
'
'  implicit none

  Dim i As Integer, rank As Integer
  Dim n As Integer, a(100) As Integer
  Dim done As Boolean
  
  n = 4

  Debug.Print " "
  Debug.Print "TEST028a"
  Debug.Print "  EQUIV_NEXT2 generates all partitions of a set."
  Debug.Print "  Here, N = ", n
  rank = 0
  done = True
  Debug.Print " "
  'Debug.Print "rank / element: ", ( i, i = 1, n )
  Debug.Print " "

  Do

    Call equiv_next2(n, a(), done)

    If (done) Then
      Exit Sub
    End If

    rank = rank + 1
    'write ( *, '(2x,i4,10x,6i4)' ) rank, a(1:n)
    Debug.Print rank, a(1), a(2), a(3), a(4)

  Loop

  Return
End Sub

Sub equiv_next2(n As Integer, a() As Integer, done As Boolean)

'*******************************************************************************
'
'! EQUIV_NEXT2 computes, one at a time, the partitions of a set.
'
'  Discussion:
'
'    A partition of a set assigns each element to exactly one subset.
'
'    The number of partitions of a set of size N is the Bell number B(N).
'
'    The entries of A are the partition subset to which each
'    element of the original set belongs.  If there are NPART distinct
'    parts of the partition, then each entry of A will be a
'    number between 1 and NPART.  Every number from 1 to NPART will
'    occur somewhere in the list.  If the entries of A are
'    examined in order, then each time a new partition subset occurs,
'    it will be the next unused integer.
'
'    For instance, for N = 4, the program will describe the set
'    where each element is in a separate subset as 1, 2, 3, 4,
'    even though such a partition might also be described as
'    4, 3, 2, 1 or even 1, 5, 8, 19.
'
'  Modified:
'
'    15 April 1999
'
'  Parameters:
'
'    Input, integer N, the number of elements in the set.
'
'    Input/output, integer A(N), contains the information
'    defining the current partition.  The user should not alter
'    A between calls.  Except for the very first
'    call, the routine uses the previous output value of A to compute
'    the next value.
'
'    Input/output, logical DONE.  Before the very first call, the
'    user should set DONE to TRUE, which prompts the program
'    to initialize its data, and return the first partition.
'    Thereafter, the user should call again, for the next
'    partition, and so on, until the routine returns with DONE
'    equal to TRUE, at which point there are no more partitions
'    to compute.
'
'  implicit none

  Dim i As Integer
  Dim imax As Integer, j As Integer, jmax As Integer

  If (done) Then
    done = False
    
  For i = 1 To n
    a(i) = 1
  Next i
    
  Else
'
'  Find the last element J that can be increased by 1.
'  This is the element that is not equal to its maximum possible value,
'  which is the maximum value of all preceding elements +1.
'
    jmax = a(1)
    imax = 1

    For j = 2 To n

      If (jmax < a(j)) Then
        jmax = a(j)
      Else
        imax = j
      End If

    Next j
'
'  If no element can be increased by 1, we are done.
'
    If (imax = 1) Then
      done = True
    Exit Sub
    End If
'
'  Increase the value of the IMAX-th element by 1, set its successors to 1.
'
    done = False
    a(imax) = a(imax) + 1
    For i = imax + 1 To n
    a(i) = 1
    Next i
    
  End If

End Sub