stochastic_matmul.py

import random
import matplotlib.pyplot as plt

# Our matrix multiplication
A = [[8,10],[4,6]]
B = [6,10]

c1_vals = []
c2_vals = []

num_trials  = 2**12
num_samples = 2**10

for j in range(num_trials):
    # Results
    c = [0,0]

    # Multiply intermediates
    mul  = [[0,0],[0,0]]

    # Add intermediates
    atff = [0,0]

    for i in range(num_samples):
        # Make a random value
        rand = random.randint(0, 15)
        randb = random.randint(0, 15)

        # Sample all values according to LFSR
        a = [[_ > rand for _ in row] for row in A]
        b = [_ > randb for _ in B]

        # Perform multiplication
        for n in range(2):
            for m in range(2):
                # First, output based on current state.
                mul[n][m] = 1 if (a[n][m] == 1 and b[m] == 1) else 0

        # Peform additon
        for n in range(2):
            # First, sum the current values
            xor = 1 if (mul[n][0] != mul[n][1]) else 0
            add = atff[n] if (xor == 1) else mul[n][1]
            c[n] += add

            # The, update flip-flops
            atff[n] = (1 - atff[n]) if (xor == 1) else atff[n]

    c1_vals.append(c[0]*16/num_samples)
    c2_vals.append(c[1]*16/num_samples)

c1_expected = (A[0][0]*B[0]/16 + A[0][1]*B[1]/16)/2
c2_expected = (A[1][0]*B[0]/16 + A[1][1]*B[1]/16)/2

print("C1 Expected: " + str(c1_expected) + " | C1 Actual: " + str(sum(c1_vals)/len(c1_vals)))
print("C2 Expected: " + str(c2_expected) + " | C2 Actual: " + str(sum(c2_vals)/len(c2_vals)))

plt.hist(c1_vals, bins=25)
plt.hist(c2_vals, bins=25)
plt.show()