1942.sagews

a=1

a+2

3

eqn(x)=x^2+x+1

eqn(2)

7

plot(eqn)

min = 0.01
n = 1000000
sum = [0.0] * n
for i in xrange(n):
    a = 1.0
    while a > min:
        a *= random()
        sum[i] += a
sum.sort()

hist_w = 0.1
hist_n = int(sum[-1] / hist_w) + 1
hist = [0] * hist_n
j = 0
for i in xrange(hist_n):
    while j < n and sum[j] < (i+1) * hist_w:
        hist[i] += 1
        j += 1
        
line2d(zip([hist_w * i for i in xrange(hist_n)], hist))

# coolest environment ever!

# goodbye hairy mathematica syntax!

# (oh except when I need to do some heavier stuff, maybe)

line2d(zip([hist_w * i for i in xrange(hist_n)], [y / n / hist_w for y in hist]))