1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
|
#!/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()
|
