#!/usr/bin/env python
"""
combino.py -- combinatory stuff related to squozen_code.html and combos.html .
The table() function creates a table used on combos.html .
"""

from factorial import *
from math import log, floor


def right0( s ):
    """
    Estimate of "linear" part of nbitsC(s) function based on combinations
    nbitsC(s) = log( (1000000+2**32-1)/1000000, 2 )          # 12.0687673
              + log(   (999999 +2**32-1  )! * (2**32-1-s)!
                    / ((999999 +2**32-1-s)! * (2**32-1  )!), 2 )
    and Stirling's formula.  Using Stirlings formula adds wrinkles that aren't
    really there.
    """
    y = log2_factorial(999999+2**32-1) + log2_factorial(2**32-1-s)
    y -= log2_factorial(999999+2**32-1-s) + log2_factorial(2**32-1)
    return y


def right1(s):
    """
    estimate of "linear" part of nbits(s) function based on combinations
    nbitsC(s) = log( (1000000+2**32-1)/1000000, 2 )          # 12.0687673
              + log(   (999999 +2**32-1  )! * (2**32-1-s)!
                    / ((999999 +2**32-1-s)! * (2**32-1  )!), 2 )
    Exact but slow.  Doesn't like s to be a float.
    """
    result = 0
    for i in xrange(1,s+1):
        result += log( (999999.0+2**32-i) / (2**32-i), 2 )
    return result


def right2( s, npieces=10000 ):
    """
    Estimate "linear" part of nbits(s) function based on combinations.
    estimate of "linear" part of nbits(s) function based on combinations
    nbitsC(s) = log( (1000000+2**32-1)/1000000, 2 )          # 12.0687673
              + log(   (999999 +2**32-1  )! * (2**32-1-s)!
                    / ((999999 +2**32-1-s)! * (2**32-1  )!), 2 )
    Medium-accurate.  Medium-fast.  Doesn't like s to be a float.
    npieces is ballpark how many pieces to break the sum into.
    """
    result = 0
    # i will range from 2**32-1 down to 2**32-s
    # Because i is the denominator, average over smaller and smaller
    # groups of terms as i gets smaller.
    i = 2**32 - 1
    bottom = 2**32 - s
    while i > bottom + npieces:
        j = max( bottom + npieces, int( i - i / npieces - 1 ) )
        n = i - j + 1
        j = ( i + j ) * .5
        result += log( (999999.0 + j) / j, 2 ) * n
        i -= n
    while i >= bottom:
        result += log( (999999.0 + i) / i, 2 )
        i -= 1
    return result


def r2t():
    """
    This compares the result of right2() to that of right1() for accuracy.
    """
    for i in [ 1, 2, 10, 100, 200, 500, 1000, 4000, 9000, 100000, 1000000 ]:
        r2 = right2(i)/i
        if i <= 1000000:
            r1 = right1(i)/i
            r1s = "%16.14f" % r1
            d = (r2-r1)/r1
        else:
            r1s = "%16s" % " "
            d = ""
        print "%10d" % i, r1s, "%16.14f" % r2, d


def table():
    print """
    <style type="text/css">
       .smallsville {
          font-size: 90%;
       }
    </style>
    """

    print '<table border="0" cellpadding="3" class="smallsville">'
    print '  <tr>'
    for heading in [ '<i>s</i>', 'nbitsC(<i>s</i>)', 
                     '(nbitsC(<i>s</i>) -&nbsp;nbitsC(0)) / <i>s</i>',
                     'nbitsC(M) -&nbsp;nbitsC(s)' ]:
        print '    <td align="center" valign="bottom">' + heading + '<br><hr></td>'
    print '  </tr>'
    nbitsC232_1 = right2(2**32-1) + 12.0687673
    for i in ( [ 1, 10, 100, 1000, 4295 ] +
               [ 10**i for i in range(4,10) ] + 
               [ 2*10**9, 3*10**9 ] +
               [ 2**32-1 - 10**i for i in range(9,0,-1) ] +
               [ 2**32-2, 2**32-1 ] ):
        r2 = right2(i)
        nbitsCs = 12.0687673+r2
        print '  <tr>',
	if i == 4295:
	    astr = "<sup>*</sup>"
  	else:
	    astr = "&nbsp;"
	if i == 4295  or  i == 2**32-1:
	    b1, b0 = "<b>", "</b>"
  	else:
	    b1, b0 = "", ""
        if 2**32-1 - 10**9 <= i <= 2**32-1:
            istr = "2<sup>32</sup>-%1d" % (2**32-i)
	else:
	    istr = "%10d" % i

        # s
        print '    <td align=right><code>' + istr + astr + '</code></td>'

        # nbitsC(s)
        print '    <td align=right><code>'+b1, "%12.3f  " % (nbitsCs), 
	print                      b0+'</code></td>'

        # ( nbitsC(s)-nbitsC(0) ) / s
        print '    <td align=right><code>', "%16.14f" % ( r2 / i ), '</code></td>',

        # nbitsC( 2*32-1 ) - nbitsC( s )
        print '    <td align=right><code>', "%12.3f  " % ( nbitsC232_1 - nbitsCs ), 
	print                      '</code></td>'

        print '  </tr>'
    print '</table>'
    print '<span class="smallsville">'
    print '    <sup>*</sup>the average gap width'
    print '</span>'


table()