Author: Zafeirakis Zafeirakopoulos
Views : 95
Description: Polyhedral Omega Demo
# Import Polyhedral Omega
from lds import *
from util import *
from geometry import *

# Import time for timing
from time import *

# Define a function to help choose data for the demo

def runPO(A, b, E, decomposition=True, parallelepipeds=True, rationalfunction=True, output=False):
lds = LinearDiophantineSystem(A, b, E)
print "A:"
for row in lds._A:
print row

if decomposition == True:
print "Computing cones"
t1=time()
cones = lds.symbolic_cones()
t2=time()
print "Cones computation took: ", t2-t1, " secs for ", len(cones), " cones"

if parallelepipeds == True:
print "Computing parallelepipeds"
t3=time()
pis = lds.fundamental_parallelepipeds()
t4=time()
print "Parallelepipeds computation took: ", t4-t3, " secs"

if rationalfunction == True:
print "Computing rational functions"
t5=time()
r = lds.rational_function_string()
t6=time()
print "Rational function computation took: ", t6-t5, " secs"

if output == True:
print "Symbolic Cones:"
for c in cones:
print c

return cones, pis,r

A = ( (1,-1) )
b = ( 0, )
E = [ 0 ]
runPO(A,b,E)

A: 1 -1 Computing cones Cones computation took: 0.000853061676025 secs for 2 cones Computing parallelepipeds Parallelepipeds computation took: 0.0496678352356 secs Computing rational functions Rational function computation took: 0.000174045562744 secs (1 * [ Cone: generators = ((0, 1), (1, 0)), apex = (0, 0), openness = (0, 0) ] + -1 * [ Cone: generators = ((0, 1), (1, 1)), apex = (0, 0), openness = (1, 0) ] , {[ Cone: generators = ((0, 1), (1, 0)), apex = (0, 0), openness = (0, 0) ]: [(0, 0)], [ Cone: generators = ((0, 1), (1, 1)), apex = (0, 0), openness = (1, 0) ]: [(0, 1)]}, '1*(z0**0*z1**0)/((1-z0**0*z1**1)*(1-z0**1*z1**0))+-1*(z0**0*z1**1)/((1-z0**0*z1**1)*(1-z0**1*z1**1))')
A = ( (1,1,1),(1,1,1),(1,1,1) )
b = ( 10,10,10)
E = [ 1,1,1]
cones, pis,r = runPO(A,b,E)
P = PolynomialRing(QQ,["z0","z1","z2"])
print P
F = P.fraction_field()
print F
print r
rr = F(r)
sr = SR(rr)
sr.parent()
var("z0,z1,z2")
print sr
print sr(1,1,1)
#sr.series(z0==0,2)


A: (1, 1, 1) (1, 1, 1) (1, 1, 1) Computing cones Cones computation took: 0.00151515007019 secs for 3 cones Computing parallelepipeds Parallelepipeds computation took: 0.00966191291809 secs Computing rational functions Rational function computation took: 0.000119209289551 secs Multivariate Polynomial Ring in z0, z1, z2 over Rational Field Fraction Field of Multivariate Polynomial Ring in z0, z1, z2 over Rational Field 1*(z0**12*z1**-1*z2**-1)/((1-z0**1*z1**-1*z2**0)*(1-z0**1*z1**0*z2**-1))+1*(z0**0*z1**0*z2**10)/((1-z0**0*z1**1*z2**-1)*(1-z0**1*z1**0*z2**-1))+-1*(z0**0*z1**11*z2**-1)/((1-z0**0*z1**1*z2**-1)*(1-z0**1*z1**-1*z2**0)) Symbolic Ring (z0, z1, z2) z0^10 + z0^9*z1 + z0^8*z1^2 + z0^7*z1^3 + z0^6*z1^4 + z0^5*z1^5 + z0^4*z1^6 + z0^3*z1^7 + z0^2*z1^8 + z0*z1^9 + z1^10 + z0^9*z2 + z0^8*z1*z2 + z0^7*z1^2*z2 + z0^6*z1^3*z2 + z0^5*z1^4*z2 + z0^4*z1^5*z2 + z0^3*z1^6*z2 + z0^2*z1^7*z2 + z0*z1^8*z2 + z1^9*z2 + z0^8*z2^2 + z0^7*z1*z2^2 + z0^6*z1^2*z2^2 + z0^5*z1^3*z2^2 + z0^4*z1^4*z2^2 + z0^3*z1^5*z2^2 + z0^2*z1^6*z2^2 + z0*z1^7*z2^2 + z1^8*z2^2 + z0^7*z2^3 + z0^6*z1*z2^3 + z0^5*z1^2*z2^3 + z0^4*z1^3*z2^3 + z0^3*z1^4*z2^3 + z0^2*z1^5*z2^3 + z0*z1^6*z2^3 + z1^7*z2^3 + z0^6*z2^4 + z0^5*z1*z2^4 + z0^4*z1^2*z2^4 + z0^3*z1^3*z2^4 + z0^2*z1^4*z2^4 + z0*z1^5*z2^4 + z1^6*z2^4 + z0^5*z2^5 + z0^4*z1*z2^5 + z0^3*z1^2*z2^5 + z0^2*z1^3*z2^5 + z0*z1^4*z2^5 + z1^5*z2^5 + z0^4*z2^6 + z0^3*z1*z2^6 + z0^2*z1^2*z2^6 + z0*z1^3*z2^6 + z1^4*z2^6 + z0^3*z2^7 + z0^2*z1*z2^7 + z0*z1^2*z2^7 + z1^3*z2^7 + z0^2*z2^8 + z0*z1*z2^8 + z1^2*z2^8 + z0*z2^9 + z1*z2^9 + z2^10 66


binomial(12,2)

66


Symbolic Ring




(z0, z1, z2) -1/(z0*z1*z2 - z0*z1 - z0*z2 - z1*z2 + z0 + z1 + z2 - 1)
Error in lines 3-3 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute flags=compile_flags), namespace, locals) File "", line 1, in <module> File "sage/symbolic/expression.pyx", line 5467, in sage.symbolic.expression.Expression.__call__ (build/cythonized/sage/symbolic/expression.cpp:31337) return self._parent._call_element_(self, *args, **kwds) File "sage/symbolic/ring.pyx", line 987, in sage.symbolic.ring.SymbolicRing._call_element_ (build/cythonized/sage/symbolic/ring.cpp:11144) return _the_element.subs(d, **kwds) File "sage/symbolic/expression.pyx", line 5321, in sage.symbolic.expression.Expression.substitute (build/cythonized/sage/symbolic/expression.cpp:30459) res = self._gobj.subs_map(smap, 0) ValueError: power::eval(): division by zero
sr

-1/(z0*z1*z2 - z0*z1 - z0*z2 - z1*z2 + z0 + z1 + z2 - 1)
sr.series?

File: /ext/sage/sage-8.8_1804/local/lib/python2.7/site-packages/sage/symbolic/expression.pyx
Signature : sr.series(self, symbol, order=None)
Docstring :
Return the power series expansion of self in terms of the given
variable to the given order.

INPUT:

* "symbol" - a symbolic variable or symbolic equality such as "x
== 5"; if an equality is given, the expansion is around the value
on the right hand side of the equality

* "order" - an integer; if nothing given, it is set to the global
default ("20"), which can be changed using
"set_series_precision()"

OUTPUT:

A power series.

To truncate the power series and obtain a normal expression, use
the "truncate()" command.

EXAMPLES:

We expand a polynomial in x about 0, about 1, and also truncate it
back to a polynomial:

sage: var('x,y')
(x, y)
sage: f = (x^3 - sin(y)*x^2 - 5*x + 3); f
x^3 - x^2*sin(y) - 5*x + 3
sage: g = f.series(x, 4); g
3 + (-5)*x + (-sin(y))*x^2 + 1*x^3 + Order(x^4)
sage: g.truncate()
x^3 - x^2*sin(y) - 5*x + 3
sage: g = f.series(x==1, 4); g
(-sin(y) - 1) + (-2*sin(y) - 2)*(x - 1) + (-sin(y) + 3)*(x - 1)^2 + 1*(x - 1)^3 + Order((x - 1)^4)
sage: h = g.truncate(); h
(x - 1)^3 - (x - 1)^2*(sin(y) - 3) - 2*(x - 1)*(sin(y) + 1) - sin(y) - 1
sage: h.expand()
x^3 - x^2*sin(y) - 5*x + 3

We computer another series expansion of an analytic function:

sage: f = sin(x)/x^2
sage: f.series(x,7)
1*x^(-1) + (-1/6)*x + 1/120*x^3 + (-1/5040)*x^5 + Order(x^7)
sage: f.series(x)
1*x^(-1) + (-1/6)*x + ... + Order(x^20)
sage: f.series(x==1,3)
(sin(1)) + (cos(1) - 2*sin(1))*(x - 1) + (-2*cos(1) + 5/2*sin(1))*(x - 1)^2 + Order((x - 1)^3)
sage: f.series(x==1,3).truncate().expand()
-2*x^2*cos(1) + 5/2*x^2*sin(1) + 5*x*cos(1) - 7*x*sin(1) - 3*cos(1) + 11/2*sin(1)

Expressions formed by combining series can be expanded by applying
series again:

sage: (1/(1-x)).series(x, 3)+(1/(1+x)).series(x,3)
(1 + 1*x + 1*x^2 + Order(x^3)) + (1 + (-1)*x + 1*x^2 + Order(x^3))
sage: _.series(x,3)
2 + 2*x^2 + Order(x^3)
sage: (1/(1-x)).series(x, 3)*(1/(1+x)).series(x,3)
(1 + 1*x + 1*x^2 + Order(x^3))*(1 + (-1)*x + 1*x^2 + Order(x^3))
sage: _.series(x,3)
1 + 1*x^2 + Order(x^3)

Following the GiNaC tutorial, we use John Machin's amazing formula
pi = 16 tan^{-1}(1/5) - 4 tan^{-1}(1/239) to compute digits of
pi. We expand the arc tangent around 0 and insert the fractions
1/5 and 1/239.

sage: x = var('x')
sage: f = atan(x).series(x, 10); f
1*x + (-1/3)*x^3 + 1/5*x^5 + (-1/7)*x^7 + 1/9*x^9 + Order(x^10)
sage: float(16*f.subs(x==1/5) - 4*f.subs(x==1/239))
3.1415926824043994