# Cauchy Product

The Cauchy Product is the discrete equivalent to the convolution operation. It is important in combinatorics (e.g. generating functions) and will cross your path because convolution gives the probability of a sum of independent random variables.


```python
import numpy as np
from time import time

def timer(func):
    def wrapperfunc(*args):
        t0 = time()
        res = func(*args)
        print("this took %1.3fs" % (time()-t0))
        return res
    return wrapperfunc


def memoize(func):
    """
    memoizing function wrapper
    """
    S = {}
    def wrapperfunction(*args):
        if args not in S:
            S[args] = func(*args)
        return S[args]
    return wrapperfunction

    

def constant(const = 1, min_n = 0, max_n=np.inf):
    """
    return function that gives series coefficent for the ith term

    c = const if min_n< = n<=max_n
    
    c = 0 else
    """
    def coeff(i):
        if i >= min_n and i <= max_n:
            return const
        else:
            return 0
    
    return coeff



def cauchy_2(a,b):
    """
    returns function that gives coefficient of the cauchy product of 
    two series with coefficients a(n) and b(n)
    
    (it's just the discrete convolution)
    """
    
    @memoize
    def coeff(i):
        return np.sum([a(h)*b(i-h) for h in range(i+1)])
    
    return coeff


@memoize
def cauchy(*args):
    """
    returns function that gives coefficient of the cauchy product of 
    many series with coefficients a(n) and b(n)
    
    (it's just the discrete convolution)
    
    Maximum recursion depth will be reached for products of too many series.
    
    Evaluation will also become slow.
    """
    if len(args) == 1:
        return args[0]
    
    else:
        return cauchy_2(args[0],cauchy(*args[1:]))
```


```python
a = constant()
b = constant()

c = cauchy(a,b)

d = cauchy(*tuple(500*[a]))

for i in range(10):
    print(i,a(i),b(i),c(i),d(i))
```

    0 1 1 1 1
    1 1 1 2 500
    2 1 1 3 125250
    3 1 1 4 20958500
    4 1 1 5 2635531375
    5 1 1 6 265661562600
    6 1 1 7 22359848185500
    7 1 1 8 1616297597409000
    8 1 1 9 102432860235795375
    9 1 1 10 5781765888864894500


# Dice 

Using the Cauchy Product to calculate the probability distribution of the sum of $n$ 6-sided dice.


```python
import matplotlib.pyplot as plt 



d6 = constant(const = 1/6, min_n=1,max_n=6)


N = 150

t0=time()
plt.figure(figsize=(12,10))
ss = list(range(-10,int(3.5*N)))
for n in range(1,N):

    d = cauchy(*tuple(n*[d6]))

    p = [d(s) for s in ss]
    
    plt.plot(ss,p)
    plt.text(x=3.5*n,y=max(p),s=n)
    
    print(n,end='\r')
    
print(time()-t0)
```

# Integer Partitons

## NEEDS REVISION -- NOT CORRECT 

The famed partition funcions P(k,n), which give the number of ways in which an integer k can divided across n subsets.

Previously, the coefficients of the product of the generating functions was:

\begin{equation}
c_k = \sum_{l=0}^{k} a_l b_{k-l}
\end{equation}

This counts the number of multisets with $k$ elements. I.e., it counts $[a_k,b_0]$ as well as $[a_0,b_k]$. In order to only count the partitions, the unordered pair $[a_k,b_0]$ should only show up once. This is accomplished by modifying the limits of the generating function:

\begin{equation}
\sum_{i=0}^k a_i x^i \sum_{j=0}^k b_j x^j \rightarrow \sum_{i=0}^k a_i x^i \sum_{j=0}^i b_j x^j
\end{equation}


```python
"""
The problem here is that the partition_n operation can't be built up from partition_2 operations the way I am trying here.
"""

def partition_2(a,b):
    """
    returns function that gives coefficient of the cauchy product of 
    two series with coefficients a(n) and b(n) modified so that the second coefficient is always
    zero whenever its subscript is larger than the subscript of the first coefficient
    """
    
    f = lambda i,h: 1 if i-h <= h else 0
    
    @memoize
    def coeff(i):
        return np.sum([a(h)*b(i-h)*f(i,h) for h in range(i+1)])
    
    return coeff




@memoize
def partition_n(*args):
    """
    returns function that gives coefficient of the cauchy product of 
    many series with coefficients a(n) and b(n)
    
    (it's just the discrete convolution)
    
    Maximum recursion depth will be reached for products of too many series.
    
    Evaluation will also become slow.
    """
    if len(args) == 1:
        return args[0]
    
    else:
        return partition_2(args[0],partition_n(*args[1:]))


        
        
@memoize
def partition(k,n):
    if n == 0:
        return 0
    else:
        d = partition_n(*tuple(n*[constant()]))
        return d(k)
```


```python
N = 16
print(''.join(['\t']+['k=%i\t' % n for n in range(1,N)]))
for n in range(1,N):
    print('\t'.join(['n=%i' % n]+['%i' % partition(k,n) for k in range(1,N)]))
```

    	k=1	k=2	k=3	k=4	k=5	k=6	k=7	k=8	k=9	k=10	k=11	k=12	k=13	k=14	k=15	
    n=1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1
    n=2	1	2	2	3	3	4	4	5	5	6	6	7	7	8	8
    n=3	1	2	2	4	4	6	6	9	9	12	12	16	16	20	20
    n=4	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=5	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=6	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=7	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=8	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=9	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=10	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=11	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=12	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=13	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=14	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26
    n=15	1	2	2	4	4	6	6	10	10	14	14	20	20	26	26



```python
N = 50
P = np.zeros([int(N/2),N])

for i in range(int(N/2)):
    n = i
    for j in range(N):
        k = j
        
        P[i,j] = partition(k,n)
```


```python
plt.figure(figsize=(18,9))
plt.pcolor(P,cmap='gist_rainbow')
plt.xlabel('k')
plt.ylabel('n')

for i in range(int(N/2)):
    for j in range(N):
        plt.text(x=j+0.125,y=i+0.25,s=str(int(P[i,j])))
```


```python
#@memoize
def partitions(k,n,i_max=None):
    """
    number partitions 
    """
    
    if not i_max:
        i_max = k
        
    if k == 0:
        return 1
    elif n == 1:
        return 1
    elif n == 0:
        return 0
    else:
        i_min = int(np.ceil(k/n))
        return int(np.sum([partitions(k-i,n-1,i_max=i) for i in range(i_min,i_max+1)]))
```


```python
for i in range(10):
    print(partitions(5,i))
```

    0
    1
    3
    9
    25
    69
    172
    414
    951
    2115


1175800c9ed8738f7a5c1679d72c502dcf5d1cda
ipynb
Jupyter Notebook
Math - Cauchy Product.ipynb
jpbm/probabilism
a2f5c1595aed616236b2b889195604f365175899
MIT
Math - Cauchy Product.ipynb
jpbm/probabilism
a2f5c1595aed616236b2b889195604f365175899
MIT
Math - Cauchy Product.ipynb
jpbm/probabilism
a2f5c1595aed616236b2b889195604f365175899
MIT
Qwen/Qwen-72B
1. YES
2. YES