euler44_v2.py

"""
Pentagonal numbers are generated by the formula, Pn=n(3n-1)/2. The first ten pentagonal numbers are:

1, 5, 12, 22, 35, 51, 70, 92, 117, 145, ...

It can be seen that P4 + P7 = 22 + 70 = 92 = P8. However, their difference, 70 - 22 = 48, is not pentagonal.

Find the pair of pentagonal numbers, Pj and Pk, for which their sum and difference are pentagonal and D = |Pk - Pj| is minimised; what is the value of D?

define pentagonal number generation:
"""

def p(n):
    return n*(3*n-1)/2

"""
p=n*(3*n-1)/2
0=1.5*n**2-.5*n-p
a=1.5
b=-.5
c=-p
n=(-b+-sqrt(b**2-4*a*c))/(2*a)
n=(.5+-sqrt(.5**2-4*1.5*-p))/(2*1.5)
n=1/6+-sqrt(.25+6*p)/3

reversed:
3n=.5+-sqrt(.25+6*p)
3n-.5=+-sqrt(.25+6p)
(3n-.5)**2=.25+6p
9n**2-3n+.25=.25+6p
9n**2-3n=6p
3n**2-n=2p
n(3n-1)/2=p

hm. so it's right, but the sqrt doesn't resolve cleanly enough to get true when it should.

so we'll round it, back in, and see if we get the same number.

now define test for pentagonality:
"""

import math

def isp(n):
    return p(round((.5+math.sqrt(.5+6*n))/3))==n
    
"""
a,b such that:
a>b
isp(p(a)+p(b))
isp(p(a)-p(b))

((.5+math.sqrt(.5+6*(a*(3*a-1)/2+b*(3*b-1)/2)))/3) %1 =0
((.5+math.sqrt(.5+6*(a*(3*a-1)/2-b*(3*b-1)/2)))/3) 
"""

def trial(i):
    for a in range(1,i+1):
        for b in range(1,a):
            pa=p(a)
            pb=p(b)
            s=pa+pb
            d=pa-pb
            sp=isp(s)
            dp=isp(d)
            success=sp and dp
            if success:
                print "{0!s}={1!s} | {2!s}={3!s} -- +{4!s} -{5!s}".format(a,pa,b,pb,s,d)