#!/usr/bin/env python
###############################################################
####### machine proof of 6th order symplectic factorization ###
###############################################################

#number of operators
nh = 3
#check the terms up to #th order
norder = 6

def op_mult(x,y):
    """multiply two polynomials (of operators) x and y"""
    z = {}
    for key1,value1 in x.iteritems():
        for key2,value2 in y.iteritems():
            key = key1+key2
            if(len(key) <= norder):
                if z.has_key(key):
                    z[key]=z[key]+value1*value2
                else:
                    z[key]=value1*value2
    for key in z.keys():
        if abs(z[key])<1.e-6:
            del z[key]
    return z

def op_add(x,y):
    """add two polynomials x and y"""
    z = x
    for key,value in y.iteritems():
        if z.has_key(key):
            z[key] = z[key] + value
            if abs(z[key])<1.e-6:
                del z[key]
        else:
            z[key] = value
    return z

def op_scale(x, s):
    """multiply a polynomial (of operators) x by a c-number s"""  
    z = {}
    if s==0:
        return {}
    for key,value in x.iteritems():
        z[key]=x[key]*s
    return z

def op_minus(x,y):
    """x-y"""
    z = x
    for key,value in y.iteritems():
        if z.has_key(key):
            z[key] = z[key] - value
            if abs(z[key])<1.e-6:
                del z[key]
        else:
            z[key] = -value
    return z


def op_exp(x):
    """polynomial expansion of e^x up to norder-th order"""
    y = {'':1}
    z = y
    for i in range(1,norder+1):
        y = op_scale(op_mult(y,x),1./i)
        z = op_add(z, y)
    return z

def op_expscale(x,s):
    """return e^{sx} up to norder-th order"""
    z = op_exp(op_scale(x,s))
    return z

si_k = 4
si_p = 1./(si_k-1.)
c1 = 1./(2. - pow(2.,si_p))
c0 = 1.- 2.*c1

#the operator A+B+C
x = {'A':1.,'B':1.,'C':1}

#the operator A
x1 = {'A':1.}

#the operator B
x2 = {'B':1.}

#the operator C
x3 = {'C':1.}


w1 = -1.17767998417887
w2 = 0.235573213359357
w3 = 0.784513610477560
w0 = 1.-2.*(w1+w2+w3)

z = op_expscale(x1,w3/2.)
z = op_mult(z, op_expscale(x2,w3/2.))
z = op_mult(z, op_expscale(x3,w3))
z = op_mult(z, op_expscale(x2,w3/2.))
z = op_mult(z, op_expscale(x1,(w3+w2)/2.))
z = op_mult(z, op_expscale(x2,w2/2.))
z = op_mult(z, op_expscale(x3,w2))
z = op_mult(z, op_expscale(x2,w2/2.))
z = op_mult(z, op_expscale(x1,(w2+w1)/2.))
z = op_mult(z, op_expscale(x2,w1/2.))
z = op_mult(z, op_expscale(x3,w1))
z = op_mult(z, op_expscale(x2,w1/2.))
z = op_mult(z, op_expscale(x1,(w1+w0)/2.))
z = op_mult(z, op_expscale(x2,w0/2.))
z = op_mult(z, op_expscale(x3,w0))
z = op_mult(z, op_expscale(x2,w0/2.))
z = op_mult(z, op_expscale(x1,(w1+w0)/2.))
z = op_mult(z, op_expscale(x2,w1/2.))
z = op_mult(z, op_expscale(x3,w1))
z = op_mult(z, op_expscale(x2,w1/2.))
z = op_mult(z, op_expscale(x1,(w2+w1)/2.))
z = op_mult(z, op_expscale(x2,w2/2.))
z = op_mult(z, op_expscale(x3,w2))
z = op_mult(z, op_expscale(x2,w2/2.))
z = op_mult(z, op_expscale(x1,(w3+w2)/2.))
z = op_mult(z, op_expscale(x2,w3/2.))
z = op_mult(z, op_expscale(x3,w3))
z = op_mult(z, op_expscale(x2,w3/2.))
z = op_mult(z, op_expscale(x1,w3/2.))
if len(op_minus(z,op_exp(x))) == 0:
    print "All terms vanish."
else:
    print "Remain "+str(len(op_minus(z,op_exp(x))))+" terms"