#!/usr/bin/python
#
# I really like this problem. Let K be any field containing the rationals
# and two elements a and c, and let A be the K-algebra generated by
# the X[i], where i runs over the integers, with the relations
#
#   X[i]^2 = X[i+1] - c + a.X[i-1]  for all i in Z.
#
# Clearly, any element of A can be written as a polynomial in X[i] so that
# all variables have exponents less than two -- otherwise, you could use
# the relation above to reduce the global degree. It's less easy to see that
# the reduced representation is canonic; in particular, it doesn't depend on
# the order of the reductions.
#
# For instance, you have:
#
# X[0]^2.X[1] = (X[1]-c+a.X[-1]).X[1]
#             = X[1]^2 - c.X[1] + a.X[-1].X[1]
#             = X[2] - c + a.X[0] - c.X[1] + a.X[-1].X[1]
#
# This first phase is called "melting". Now, it turns out that there is
# a linear form on A (let's call it the integral) which is particularly
# easy to describe on the reduced form. Each term is (up to a factor)
# a product of X[i] where i runs over some finite subset S of the integers.
# Then we define integral(prod X[i]) = (3.a/4)^(card(S)/2) if all the
# connected components of S have even cardinality, and 0 otherwise.
# In the above case the integral is (-c).
#
# If you want to know where all of this comes from, get the preprint at
# http://topo.math.u-psud.fr/~bousch/preprints/ called "Algebres de Henon"
# (in french). Also, have a look at the file called "fission.form" if
# you're interested in the contorsions necessary to do this in FORM.
# And if you have FORM, don't forget to do some performance comparisons.

from saml1 import *
import sys, string, regex

def extract_literals(poly, rx):
  """Extract all literals matching a certain regular expression"""
  from string import splitfields, joinfields
  prog = regex.compile(rx)
  str = repr(poly)
  found = {}
  while 1:
    position = prog.search(str)
    if position < 0:
      break
    matched = prog.match(str, position)
    lit = str[position:position+matched]
    found[lit] = 1
    # Now replace all the occurences of this literal with zeroes
    str = joinfields(splitfields(str,lit),'0')
  return found.keys()

def melt(poly):
  sys.stderr.write('Literals: ');
  candidates = extract_literals(poly, "X\\[-?[0-9]+\\]")
  sys.stderr.write(string.join(candidates,',')+'\n')
  sys.stderr.write('melting... ')
  for x in candidates:
    p = string.atoi(x[2:-1])
    X = poly.parse('X[%d]^2' % p)
    Y = poly.parse('a.X[%d]-c+X[%d]' % (p-1,p+1))
    poly = poly.subs(X,Y)
  sys.stderr.write('done\n')
  return poly

def combine(poly):
  literals = extract_literals(poly, "X\\[-?[0-9]+\\]")
  indices = map(lambda s:string.atoi(s[2:-1]), literals)
  for i in range(min(indices),max(indices)-min(indices)):
    X = poly.parse('X[%d].X[%d]' % (i,i+1))
    Y = poly.parse('3.a/4')
    poly = poly.subs(X,Y)
  for l in literals:
    poly = poly.subs(poly.parse(l), poly.zero())
  return poly

def test(expression):
  model = Mathnode(ST_RATIONAL,'0').cast(ST_APOLY)
  poly = model.parse(expression)
  while 1:
    new = melt(poly)
    if new==poly: break
    poly = new
  print 'After melting:'
  print poly
  print 'Integral:'
  poly = combine(poly)
  print poly
  print 'And when a=0 it reduces to:'
  print poly.subs(poly.parse('a'), poly.zero())

if __name__ == '__main__':
  if len(sys.argv) == 1:
    expr = '(X[-1]+X[0]+X[1])^12'
    print 'Initial expression:', expr
  else:
    expr = sys.argv[1]
  test(expr)
  print MnAllocStats()