def lcmfirstints(n):
    return lcm([1..n])

def LPF(n):
    if n==1:
        return 1
    return factor(n)[-1][0]

def LPP(n):
    if n==1:
        return 0
    f=factor(n)
    return max([p[0]^p[1] for p in f])

def CountPrimSets(V): #Returns the number of primitive subsets of the list of integers V
    if len(V)==0:
        return 1
    components=ListConnectedComponents(V)
    count=1
    for c in components:
        if len(c)==1:
            count*=2
        elif len(c)==2:
            count*=3
        elif len(c)==3:
            if not c[1]%c[0]==0 or not c[2]%c[1]==0:
                count*=5
            else:
                count*=4
        else:
            count*=sum(reduceCountPrimSubsets(c))
    return count

def ListConnectedComponents(V): #Break the list of integers V into connected components of the divisibility poset
    a=V[0]
    n=V[-1]
    L=dict()
    Components=[]
    for i in V:
        if i in L: #i has been visited, so its multiples
            #have already been taken into account
            continue
        L[i]=[i]
        Components.append(L[i])
        if i>n/2:
            continue
        for j in range(2,int(floor(n/i))+1):
            h=i*j
            if not h in V:
                continue
            if not h in L:
                L[i].append(i*j)
                L[h]=L[i]
                continue
            elif L[i]==L[h]:
                continue
            L[h].extend(L[i])
            Components.remove(L[i])
            for k in reversed(L[i]): #Reversed so that i is processed last...
                L[k]=L[h]
    if len(Components)==1:
        return [V[:]]
    for c in Components:
        c.sort()
    return Components



def reduceCountPrimSubsets(L): #Returns a pair, the count of primitive sets containing the first element of L, and the number not containing the first element of L
    V=L[:]
    if len(V)==0:
        return (0,1)
    a=V[0]
    n=V[-1]
    del V[0]
    countnothasa=CountPrimSets(V)
    if a==1:
	counthasa=1
    else:
	for i in range(2,int(floor(n/a))+1):
            if (a*i) in V:
                V.remove(a*i)
        counthasa=CountPrimSets(V)
    return (counthasa,countnothasa)


def seqgcd(L):
    g=L[0]
    for i in L:
        if not i%g==0:
            g=gcd(Integer(g),Integer(i))
            if g==1:
                return 1
    return g

def is_pow(p,n):
    if (n^(1/p)).is_integer():
        return true
    return false

def largestPwrDiv(p,n):
    prod=1
    for f in factor(n):
        prod*=f[0]**floor(f[1]/p)
    return prod


def validDs(n):
    P=[]
    for p in Primes():
        if p<n:
            P.append(p)
        else:
            break
    d=1
    L=[1..n]
    grow=true
    while grow:
        grow=false
        for p in P:
            for l in L:
                if not p*l in L and p*l<d*(n+1):
                    L.append(p*l)
                    grow=true
                    break
        #print L
        for p in P:
            for l in L:
                if not p.divides(l) and l/p > d:
                    d*=p
                    L=[k*p for k in L]
                    L.append(l)
                    grow=true
                    break
        #print L
        #print d,
    return d

def validOddDs(n):   #Returns an integer, supported on the primes p<=n which is the largest at which new structure is added to the graph.
                     #WARNING: this is only finite for small values of n.  When n is sufficiently large this will run forever...
    P=[]
    for p in Primes():
        if p<n:
            if not p==2:
                P.append(p)
        else:
            break
    d=1
    L=[1..n]
    grow=true
    while grow:
        grow=false
        for p in P:
            for l in L:
                if not p*l in L and p*l<d*(n+1):
                    L.append(p*l)
                    grow=true
                    break
        #print L
        for p in P:
            for l in L:
                if not p.divides(l) and l/p > d:
                    d*=p
                    L=[k*p for k in L]
                    L.append(l)
                    grow=true
                    break
        #print L
        #print d,
    return d

[validOddDs(n) for n in [1..10]]

newprodlog=0
accountedfor=0   # Total size of the subset of [0,1] that has been taken into account
R=RealField(300)
for n in range(1,11): #process the interval (1/(n+1),1/n]
    D=validOddDs(n)   # For small values of n the potential values of d is a finite set, this returns the largest value of d that needs to be considered
    dlist=[d for d in divisors(D)]
    for d in dlist:
        valstoprocess=[[d,i] for i in range((d*n),(d*(n+1))) if i==(d*n) or LPF(i)<=n]
        vals=[i for i in range((d*n)+1,(d*(n+1)+1)) if i==(d*(n+1)) or LPF(i)<=n]
        prev=d*n
        L=[]
        for p in valstoprocess:
            v=p[1]
            if len(L)==0:
                I=[i for i in range(d,d*n+1) if not 2.divides(i)]
                L=sorted(ListConnectedComponents(I))[0]
                (counthasd,countnothasd)=reduceCountPrimSubsets(L[:])
            else:
                vDivs=[j for j in divisors(v) if j>d and not j in L]
                if len(vDivs)>1:
                    CurList=sorted(ListConnectedComponents([i for i in [d..v] if i%2==1]))[0]
                    NewElements=[i for i in CurList if not i in L and not i==v]
                    CountPrimSubsNewElements=sum(reduceCountPrimSubsets(NewElements))
                    counthasd*=CountPrimSubsNewElements
                    countnothasd*=CountPrimSubsNewElements
                else:
                    CurList=L+vDivs
                    CurList.sort()
                (x,y)=reduceCountPrimSubsets([i for i in CurList if not Integer(i).divides(v)]) #Count sets which contain v
                counthasd+=x
                countnothasd+=y
                L=CurList
            nextv=vals.pop(0)
            r=1+counthasd/countnothasd
            print '\t',d,prev,r,N(r),prev,nextv,
            for pr in r.support():
                newprodlog+=(d*euler_phi(D/d)*(nextv-prev)/(D*(prev)*(nextv)))*log(R(pr))*r.valuation(pr)
            print RR(newprodlog)
            accountedfor+=(d*euler_phi(D/d)*(nextv-prev)/(D*(prev)*(nextv)))
            prev=nextv
        print (n+1),d,d*n,RR(exp(newprodlog/2)), RR(exp(newprodlog/2)*sqrt(2)^(1.0-accountedfor))
        print " "