response to http://math.stackexchange.com/questions/1159073/ about GCD function

gcd.py

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# http://stackoverflow.com/questions/11175131/code-for-greatest-common-divisor-in-python
def gcd(x, y):
    while y != 0:
        (x, y) = (y, x % y)
    return x

# a hack
def gcdsteps(x, y):
    c = 0
    while y != 0:
        (x, y) = (y, x % y)
        c += 1
    return c

# my work

X = [gcdsteps(int(1000000*np.random.random()),int(1000000*np.random.random())) for s in range(100) for t in range(100)]

plt.hist( X, bins = np.arange(30))
plt.show()

smooth.py

bits = np.random.random(len(P))

P = [2,3,5,7,11,13,17,19]
P = np.array(P)
P1, P2 = P[bits > 0.5], P[bits < 0.5]

def smooth(P):

    X = [1]

    for p in P:
        B = [p**k for k in np.arange(6*np.log(10)/np.log(p)).astype(int)]
        X = [x*b for x in X for b in B if(x*b) < 10**6]

    return X

plt.hist([gcdsteps(s,t) for s in smooth(P1) for t in smooth(P2)], bins = np.arange(30))
plt.show()

# or print smooth relatively prime pairs
#
# (540225, 131072, 12),
# (385875, 524288, 13),
# (275625, 524288, 11),
# (196875, 262144, 10),

A = [a for a in smooth(P1)[:] if a > 100000]
B = [b for b in smooth(P1)[:] if b > 100000]

relativelyPrime = [(s,t, gcdsteps(s,t)) for s in A for t in B if gcd(s,t) == 1 ]