[Search for users] [Overall Top Noters] [List of all Conferences] [Download this site]

Conference rusure::math

Title:Mathematics at DEC
Moderator:RUSURE::EDP
Created:Mon Feb 03 1986
Last Modified:Fri Jun 06 1997
Last Successful Update:Fri Jun 06 1997
Number of topics:2083
Total number of notes:14613

842.0. "Stirling numbers" by ELWD2::CHINNASWAMY () Sat Mar 12 1988 11:33

I happened to look into this conference only a few days ago. I missed a
golden opportunity to participate in it for the last several years in DEC.

A few months ago, I was working on the behavior of a cross bar switch. 
I got the following expression which I found very difficult to simlify to a
closed form. I wrote a program to get some numerical results for different
values of N, M and B. ( for instance N = 40, M = 24 and B = 16). It took for
over a week on a 8650 to finish the computation. I then simplified the
expression somewhat. With this simplification, I was able to get the results
I believe, in a day. I was not happy with my simplification. I got busy with
other things. I haven't touched it since. But if any of you are more familiar
with stirling numbers and their properties, may be this is an easy problem. 
I would like to get a simple closed form.

                   N                    Min(n,M)
                  ==                    ==
                  \ (N,n)p^n(1-p)^(N-n) \ (k-B)*k!S(n,k)(M,k)/M^n 
                  /                     / 
                  ==                    ==
                 n = B+1                k=B+1

where S(N,k) is the Stirling number of the second kind and (M,k) denotes
the combination of N things taken k at a time and p^n denotes p raised to the
power n etc.


T.RTitleUserPersonal
Name
DateLines
842.1PLDVS2::YOUNGBack from the Shadows Again,Sat Mar 12 1988 22:212
    What is "p" in your equation?  Also some of the closure seems ambiguous
    to me, could you clarify?
842.2ELWD2::CHINNASWAMYSun Mar 13 1988 21:2222
>    What is "p" in your equation?  Also some of the closure seems ambiguous
>    to me, could you clarify?

    I thought I explained everything. Here is some more explanation:


           N                     Min(n,M)
          ==                      ==
          \ (N,n)*p^n*(1-p)^(N-n) \ (k-B)*k!*S(n,k)*(M,k)/M^n 
          /                       / 
          ==                      ==
         n = B+1                 k=B+1

    where S(N,k) is the Stirling number of the second kind and (M,k) denotes
    the combination of N things taken k at a time and p^n denotes p raised to
    the power n , p is probability ( Assume any value, if a numerical value is
    desired. But an algebraic expression is needed), k! denotes the factorial
    of k, etc.

    I have added a few asterisks to denote multiplication.

842.3Faster Program:PLDVS2::YOUNGBack from the Shadows Again,Mon Mar 14 1988 02:06161
    Well, if all you want or need is to get an answer faster, the following
    program seems to complete in less than 1/10 th of a second on an
    VAX 8250.
    
    I cannot check it for correctness, as I do not know the correct
    answers.  If this has bugs or is not what you need let me know.
    
    --  Barry
    
-----------------CROSS_BAR.FOR-------------------------------------------
        OPTIONS /G_FLOAT
        Program CROSS_BAR

        Implicit none
        Include 'COMMONDEF.FOR'
        External    Chins_problem
        Real*8      Chins_problem, Temp, p
        Integer     N, M, B, I
        
        
        !--

        Type *, 'Enter p (probabilty): '
        Accept *, p

        Do while (1.eq.1)
        
            Type *, 'Enter N,M,B : '
100         Format (I,I,I)
            Accept 100, N, M, B

            Type *, N, M, B

            Temp = Chins_problem ( N, M, B, p )

            Type *, temp

        End do


        End


!----

        OPTIONS /G_FLOAT /extend_source
        Real*8 function CHINS_PROBLEM ( N, M, B, p )

!++++
!
!   This subroutine calculates the solution to the solution to the formula
! posted int MATH::842 by CHINNASWAMMY, which is:
!
!           N                    Min(i,M)
!          ==                    ==
!          \ (N,i)p^i(1-p)^(N-i) \ (k-B)*k!S2(i,k)(M,k)/M^i 
!          /                     / 
!          ==                    ==
!         i = B+1                k=B+1
!
!   Where  (M,k)  is the number of combinations of M things taken k at a 
! time, and  S2(i,k)  is a stirling number of the second kind.
!
!----

        Implicit none
        Include 'COMMONDEF.FOR'
        Real*8      Fact, S2, Comb, Sum1, Sum2
        Integer     First_time  /1/
        Integer     I, J, K, L
        Integer     M, N, B
        Real*8      p
        
        Fact(l)  =  Factl(l+1)
        S2(l,j)  =  stirl(j,l)
        Comb(l,j)=  fact(l) / ( fact(j)*fact(l-j) )
        
        !--

        If ( first_time ) then
            Call Init_factorials
            Call Init_stirling2
            First_time  =   0
        End if

        Sum1    =   0

        Do i = B+1, N

            Sum2    =   0

            Do k = B+1, Min(M,i)
                Sum2    =   Sum2  +  (k-B) * Fact(k) * S2(i,k) * Comb(M,k)
            End do

            Sum1    =   Sum1 + ( Comb(N,i) * p**i * ((1-p)**(N-i)) * Sum2 ) 
     +              /   ( Dble(M)**i )

        End do

        Chins_problem   =   Sum1


        End

!----

        OPTIONS /G_FLOAT
        Subroutine INIT_FACTORIALS

        Implicit none
        Include 'COMMONDEF.FOR'

        Integer     I

        !--

        Factl ( 1 ) =   1

        Do i = 2, 101

            Factl (i)   =   Factl (i-1) * i-1       ! Offset 1 to account
                                                    ! for (0!).
        End do

        End

!----


        OPTIONS /G_FLOAT
        Subroutine INIT_STIRLING2

        Implicit none
        Include 'COMMONDEF.FOR'

        Integer     I, J

        !--

        Do i = 1, 100

            Stirl ( 1, i )  =   1
            Stirl ( i, i )  =   1

            Do j = 2, i-1

                Stirl ( j, i )  =   j * stirl(j,i-1) + stirl(j-1,i-1)

            End do

        End do

        End

-------------------------COMMONDEF.FOR-------------------------------------
        Real*8      Factl(101), Stirl(100,100)

        Common      /Funct_vals/
     +              Factl,
     +              Stirl
842.4Typo in Factorial routineSSDEVO::LARYMon Mar 14 1988 12:434
		Factl(i) = factl(i-1) * i-1

should be:	Factl(i) = factl(i-1) * (i-1)

842.5Gaaaahhhhh!!!!PLDVS2::YOUNGBack from the Shadows Again,Mon Mar 14 1988 13:12158
    Thanks for the sanity-check.
    
    Corrected version follows:
    
-----------------CROSS_BAR.FOR-------------------------------------------
        OPTIONS /G_FLOAT
        Program CROSS_BAR

        Implicit none
        Include 'COMMONDEF.FOR'
        External    Chins_problem
        Real*8      Chins_problem, Temp, p
        Integer     N, M, B, I
        
        
        !--

        Type *, 'Enter p (probabilty): '
        Accept *, p

        Do while (1.eq.1)
        
            Type *, 'Enter N,M,B : '
100         Format (I,I,I)
            Accept 100, N, M, B

            Type *, N, M, B

            Temp = Chins_problem ( N, M, B, p )

            Type *, temp

        End do


        End


!----

        OPTIONS /G_FLOAT /extend_source
        Real*8 function CHINS_PROBLEM ( N, M, B, p )

!++++
!
!   This subroutine calculates the solution to the solution to the formula
! posted int MATH::842 by CHINNASWAMMY, which is:
!
!           N                    Min(i,M)
!          ==                    ==
!          \ (N,i)p^i(1-p)^(N-i) \ (k-B)*k!S2(i,k)(M,k)/M^i 
!          /                     / 
!          ==                    ==
!         i = B+1                k=B+1
!
!   Where  (M,k)  is the number of combinations of M things taken k at a 
! time, and  S2(i,k)  is a stirling number of the second kind.
!
!----

        Implicit none
        Include 'COMMONDEF.FOR'
        Real*8      Fact, S2, Comb, Sum1, Sum2
        Integer     First_time  /1/
        Integer     I, J, K, L
        Integer     M, N, B
        Real*8      p
        
        Fact(l)  =  Factl(l+1)
        S2(l,j)  =  stirl(j,l)
        Comb(l,j)=  fact(l) / ( fact(j)*fact(l-j) )
        
        !--

        If ( first_time ) then
            Call Init_factorials
            Call Init_stirling2
            First_time  =   0
        End if

        Sum1    =   0

        Do i = B+1, N

            Sum2    =   0

            Do k = B+1, Min(M,i)
                Sum2    =   Sum2  +  (k-B) * Fact(k) * S2(i,k) * Comb(M,k)
            End do

            Sum1    =   Sum1 + ( Comb(N,i) * p**i * ((1-p)**(N-i)) * Sum2 ) 
     +              /   ( Dble(M)**i )

        End do

        Chins_problem   =   Sum1


        End

!----

        OPTIONS /G_FLOAT
        Subroutine INIT_FACTORIALS

        Implicit none
        Include 'COMMONDEF.FOR'

        Integer     I

        !--

        Factl ( 1 ) =   1

        Do i = 2, 101

            Factl (i)   =   Factl (i-1) * (i-1)     ! Offset 1 to account
                                                    ! for (0!).
        End do

        End

!----


        OPTIONS /G_FLOAT
        Subroutine INIT_STIRLING2

        Implicit none
        Include 'COMMONDEF.FOR'

        Integer     I, J

        !--

        Do i = 1, 100

            Stirl ( 1, i )  =   1
            Stirl ( i, i )  =   1

            Do j = 2, i-1

                Stirl ( j, i )  =   j * stirl(j,i-1) + stirl(j-1,i-1)

            End do

        End do

        End

-------------------------COMMONDEF.FOR-------------------------------------
        Real*8      Factl(101), Stirl(100,100)

        Common      /Funct_vals/
     +              Factl,
     +              Stirl


842.6it's not "times (i-1)"ZFC::DERAMOperforming without a netMon Mar 14 1988 22:5013
    Re .4:
    
>>                       -< Typo in Factorial routine >-
>>
>>		Factl(i) = factl(i-1) * i-1
>>
>>   should be:	Factl(i) = factl(i-1) * (i-1)
    
    factorial(i) = 1 * 2 * ... * (i - 1) * i
    
                 = factorial(i - 1) * i
    
    Dan
842.7Some hand-waving here...CHOVAX::YOUNGBack from the Shadows Again,Tue Mar 15 1988 03:2713
    Actually Dan, I think that this is explained in the comments which
    are in the program.
    
    Since FORTRAN has 1's based arrays, the array Factl() could not
    return the value of 0! if the mapping of array index and function
    argument were exact.  Therefore I bulit in a 1 place shift so that
    Factl(1) = 0!, and in general Factl(i) = (i-1)! .  For clarity I
    define a FORTRAN statement function elsewhere called Fact() that
    makes the 1 step translation so that  Fact(i) = Factl(i+1) = i!
    
    So this should be correct.
    
    --  Barry
842.8ELWD2::CHINNASWAMYTue Mar 15 1988 10:3816
	Barry,

	Thanks for your solution. The only difference, if I recall correctly
	between my program and yours is that I wrote it in Pascal, there was
	one more term in the calculation, I calculated for several values of
	p and submitted it as a batch job. But these things should not make
	that much of a difference in compute time. I will look at the
	difference more carefully this weekend.

	It would be nice if we can find an algebraic expression.

	Again, thanks for your help.

	Swamy

842.9ZFC::DERAMOperforming without a netTue Mar 15 1988 14:257
    Re .7, factorials etc.
    
    I suspected that it might have been "out of context", but I
    had skipped over the program. I was going to suggest that you
    rename the function to Gamma. (-:
    
    Dan
842.10PLDVS1::YOUNGBack from the Shadows Again,Tue Mar 15 1988 15:057
    Re .8:
    
    Did your Pascal program also precalculate the Factorial and Strirling
    values?  This change in my program made it an order O(n^2) program
    instead of the order O(n^4) program that it was originally.
    
    --  Barry
842.11Hey, what the heck is a Cross-bar switch?CHOVAX::YOUNGBack from the Shadows Again,Sat Mar 19 1988 10:071
    
842.12Stirling numbers and cross barELWD2::CHINNASWAMYMon Mar 28 1988 00:1417
    The following three references will give an idea of what a cross bar
    switch is and where the formulas came from:

    1. Yu-cheng Liu and Chi-jiunn Jou, "Effective Memory Bandwidth and
       Processor Blocking Probability in Multiple-Bus Systems", IEEE Trans.
       Comput., Vol. C-36, No. 6, pp 761 - 764, June 1987

    2. D. Y. Chang, D. J. Kuck,, and D. H. Lawrie, "On the Effective Bandwidth
       of Parallel Memories", IEEE Trans. Comput., vol. c-26, pp. 480-490, May
       1977.

    3. T. Lang, M. Valero, and I. Alegre, "Bandwidth of crossbar and multiple-
       bus connections for multiprocessors", IEEE Trans. Comput., vol. C-31,
       pp. 1227-1234, CDec. 1982.

    swamy