Hello,
> Mike, could you point me to your code?
I've attached the _very_ rough version I started awhile back. It just
does FormalPowerSeries and DataStream.
> I actually wonder why you would
> reprogram it in Sage and not interfacing the aldor-combinat library?
> Is there a problem with the license of the aldor compiler
> (aldor-combinat itself is GPL2).
Well, tree-like structures ( and hence species ) should be a standard
part of Sage. If we were to use aldor-combinat to do those
computations, then aldor would need to become a standard part of Sage.
This presents a number of issues:
* Aldor must build cleanly on all the platforms Sage builds on. This
is something that someone much more familiar with Aldor. The first
step would be to create an experimental spkg.
* I have no idea how (if possible) to interface an Aldor program from
C. pexpect interfaces are slow.
* licensing issues need to be sorted out since the APL2 is not
GPL-compatible. I don't know if it'd even be possible to distribute
Aldor as part of Sage.
Also, part of it was just a chance for me to learn more about species.
I definitely cannot program in Aldor, and one of the reasons why
Python has been so nice to work with is that is just sort of stays out
of the way. I'd much think about programming math instead of
programming languages.
--Mike
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~----------~----~----~----~------~----~------~--~---
from stream import Stream
from series_order import InfiniteSeriesOrder, UnknownSeriesOrder, bounded_decrement, increment
from sage.rings.all import Integer
import pdb
inf = InfiniteSeriesOrder()
unk = UnknownSeriesOrder()
"""
#Test Plus 1
sage: from sage.combinat.species.series import *
sage: from sage.combinat.species.stream import Stream
sage: gs0 = FormalPowerSeries(Stream(const=0))
sage: gs1 = FormalPowerSeries(Stream(const=1))
sage: sum1 = gs0 + gs1
sage: sum2 = gs1 + gs1
sage: sum3 = gs1 + gs0
sage: [gs0.coefficient(i) for i in range(11)]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
sage: [gs1.coefficient(i) for i in range(11)]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
sage: [sum1.coefficient(i) for i in range(11)]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
sage: [sum2.coefficient(i) for i in range(11)]
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
sage: [sum3.coefficient(i) for i in range(11)]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
#Test Plus 2
sage: l1 = [1,2,4,8,0]
sage: l2 = [-1, 0,-1,-9,22,0]
sage: gs1 = FormalPowerSeries(Stream(l1))
sage: gs2 = FormalPowerSeries(Stream(l2))
sage: sum = gs1 + gs2
sage: sum2 = gs2 + gs1
sage: [ sum.coefficient(i) for i in range(5) ]
[0, 2, 3, -1, 22]
sage: [ sum.coefficient(i) for i in range(5, 11) ]
[0, 0, 0, 0, 0, 0]
sage: [ sum2.coefficient(i) for i in range(5) ]
[0, 2, 3, -1, 22]
sage: [ sum2.coefficient(i) for i in range(5, 11) ]
[0, 0, 0, 0, 0, 0]
#Test Zero
sage: s = FormalPowerSeries(Stream([0,2,3,0]))
sage: s.zero()
False
sage: s = FormalPowerSeries(0)
sage: s.zero()
True
sage: s = FormalPowerSeries(Stream([0]))
sage: s.zero()
False
sage: s.coefficient(0)
0
sage: s.coefficient(1)
0
sage: s.zero()
True
#Test Times 1
sage: gs0 = FormalPowerSeries(Stream(const=0))
sage: gs1 = FormalPowerSeries(Stream(const=1))
sage: prod0 = gs0 * gs1
sage: prod1 = gs1 * gs0
sage: prod2 = gs1 * gs1
sage: [prod0.coefficient(i) for i in range(11)]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
sage: [prod1.coefficient(i) for i in range(11)]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
sage: [prod2.coefficient(i) for i in range(11)]
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
#Test Times 2
sage: l1 = [1,2,4,8,0]
sage: l2 = [-1, 0,-1,-9,22,0]
sage: gs1 = FormalPowerSeries(Stream(l1))
sage: gs2 = FormalPowerSeries(Stream(l2))
sage: prod1 = gs1 * gs2
sage: prod2 = gs2 * gs1
sage: [prod1.coefficient(i) for i in range(11)]
[-1, -2, -5, -19, 0, 0, 16, 176, 0, 0, 0]
sage: [prod2.coefficient(i) for i in range(11)]
[-1, -2, -5, -19, 0, 0, 16, 176, 0, 0, 0]
#Test Recursive 0
sage: s = FormalPowerSeries()
sage: s._name = 's'
sage: s.set(one+monom*s)
sage: s.aorder
0
sage: s.order
Unknown series order
sage: [s.coefficient(i) for i in range(6)]
[1, 1, 1, 1, 1, 1]
#Test Recursive 1
sage: s = FormalPowerSeries()
sage: s._name = 's'
sage: s.set(one+monom*s*s)
sage: s.aorder
0
sage: s.order
Unknown series order
sage: [s.coefficient(i) for i in range(6)]
[1, 1, 2, 5, 14, 42]
#Test Recursive 1b
sage: s = FormalPowerSeries()
sage: s._name = 's'
sage: s.set(monom + s*s)
sage: s.aorder
1
sage: s.order
Unknown series order
sage: [s.coefficient(i) for i in range(7)]
[0, 1, 1, 2, 5, 14, 42]
#Test Recursive 2
sage: s = FormalPowerSeries()
sage: s._name = 's'
sage: t = FormalPowerSeries()
sage: t._name = 't'
sage: s.set(one+monom*t*t*t)
sage: t.set(one+monom*s*s)
sage: [s.coefficient(i) for i in range(9)]
[1, 1, 3, 9, 34, 132, 546, 2327, 10191]
sage: [t.coefficient(i) for i in range(9)]
[1, 1, 2, 7, 24, 95, 386, 1641, 7150]
#Test Recursive 2b
sage: s = FormalPowerSeries()
sage: s._name = 's'
sage: t = FormalPowerSeries()
sage: t._name = 't'
sage: s.set(monom + t*t*t)
sage: t.set(monom + s*s)
sage: [s.coefficient(i) for i in range(9)]
[0, 1, 0, 1, 3, 3, 7, 30, 63]
sage: [t.coefficient(i) for i in range(9)]
[0, 1, 1, 0, 2, 6, 7, 20, 75]
#Test Recursive 3
sage: s = FormalPowerSeries()
sage: s._name = 's'
sage: s.set(one+monom*s*s*s)
sage: [s.coefficient(i) for i in range(10)]
[1, 1, 3, 12, 55, 273, 1428, 7752, 43263, 246675]
#Test PlusTimes 1
sage: s = (one + monom) * (one + monom)
sage: [s.coefficient(i) for i in range(6)]
[1, 2, 1, 0, 0, 0]
#Test Compose 1
sage: s = FormalPowerSeries(Stream(const=1))
sage: t = FormalPowerSeries(Stream([0,0,1]))
sage: u = s(t)
sage: u.coefficients(11)
[1, 0, 1, 1, 2, 3, 5, 8, 13, 21, 34]
#Test Compose 2
sage: s = FormalPowerSeries(Stream(const=1))
sage: t = FormalPowerSeries(Stream([0,0,1,0]))
sage: u = s(t)
sage: u.aorder
0
sage: u.order
Unknown series order
sage: u.coefficients(10)
[1, 0, 1, 0, 1, 0, 1, 0, 1, 0]
sage: u.aorder
0
sage: u.order
0
#Test Compose 3
# s = 1/(1-x), t = x/(1-x)
# s(t) = (1-x)/(1-2x)
sage: s = FormalPowerSeries(Stream(const=1))
sage: t = FormalPowerSeries(Stream([0,1]))
sage: u = s(t)
sage: u.coefficients(14)
[1, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
#Test Differentiate 1
sage: s = FormalPowerSeries(Stream(const=1))
sage: u = s.differentiate()
sage: u.coefficients(10)
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
#Test Differentiate 2
sage: s = FormalPowerSeries()
sage: s._name = 's'
sage: s.set(one+monom*s*s)
sage: u = s.differentiate()
sage: u.coefficients(5) #[1*1, 2*2, 3*5, 4*14, 5*42]
[1, 4, 15, 56, 210]
#Test Differentiate 3
sage: s = FormalPowerSeries(Stream(const=1))
sage: t = FormalPowerSeries(Stream([0,1]))
sage: u = s(t).differentiate()
sage: v = (s.differentiate().compose(t))*t.differentiate()
sage: u.coefficients(11)
[1, 4, 12, 32, 80, 192, 448, 1024, 2304, 5120, 11264]
sage: v.coefficients(11)
[1, 4, 12, 32, 80, 192, 448, 1024, 2304, 5120, 11264]
#Test Differentiate 4
sage: s = FormalPowerSeries(); s._name='s'
sage: t = FormalPowerSeries(); t._name='t'
sage: s.set(monom+t*t*t)
sage: t.set(monom+s*s)
sage: u = (s*t).differentiate()
sage: v = s.differentiate()*t + s*t.differentiate()
sage: u.coefficients(10)
[0, 2, 3, 4, 30, 72, 133, 552, 1791, 4260]
sage: v.coefficients(10)
[0, 2, 3, 4, 30, 72, 133, 552, 1791, 4260]
sage: u.coefficients(10) == v.coefficients(10)
True
#Test Integrate 1a
sage: s = zero
sage: t = s.integrate()
sage: t.zero()
True
#Test Integrate 1b
sage: s = zero
sage: t = s.integrate(1)
sage: t.coefficients(6)
[1, 0, 0, 0, 0, 0]
sage: t.stream.is_constant()
True
#Test Integrate 2a
sage: s = term(1, 0)
sage: t = s.integrate()
sage: t.coefficients(6)
[0, 1, 0, 0, 0, 0]
sage: t.stream.is_constant()
True
#Test Integrate 2b
sage: s = term(1,0)
sage: t = s.integrate(1)
sage: t.coefficients(6)
[1, 1, 0, 0, 0, 0]
sage: t.stream.is_constant()
True
#Test Integrate 3a
sage: s = term(1, 4)
sage: t = s.integrate()
sage: t.coefficients(10)
[0, 0, 0, 0, 0, 1/5, 0, 0, 0, 0]
#Test Integrate 3b
sage: s = term(1,4)
sage: t = s.integrate(1)
sage: t.coefficients(10)
[1, 0, 0, 0, 0, 1/5, 0, 0, 0, 0]
#Test Sum 1
sage: s = FormalPowerSeries(Stream(const=1))
sage: g = [s]*6 + [zero]
sage: t = fps_sum_generator(g)
sage: t.coefficients(10)
[1, 2, 3, 4, 5, 6, 6, 6, 6, 6]
#Test Sum 2
sage: s = FormalPowerSeries(Stream(const=1))
sage: def g():
... while True:
... yield s
sage: t = fps_sum_generator(g())
sage: t.coefficients(9)
[1, 2, 3, 4, 5, 6, 7, 8, 9]
#Test Product 1
sage: s1 = FormalPowerSeries(Stream([1,1,0]))
sage: s2 = FormalPowerSeries(Stream([1,0,1,0]))
sage: s3 = FormalPowerSeries(Stream([1,0,0,1,0]))
sage: s4 = FormalPowerSeries(Stream([1,0,0,0,1,0]))
sage: s5 = FormalPowerSeries(Stream([1,0,0,0,0,1,0]))
sage: s6 = FormalPowerSeries(Stream([1,0,0,0,0,0,1,0]))
sage: s = [s1, s2, s3, s4, s5, s6]
sage: def g():
... for a in s:
... yield a
sage: p = fps_product_generator(g())
sage: p.coefficients(26)
[1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 5, 5, 4, 4, 4, 3, 2, 2, 1, 1, 1, 0, 0, 0, 0]
#Test Product 2
sage: def m(n):
... n -= 1
... yield 1
... while True:
... for i in range(n):
... yield 0
... yield 1
...
sage: def s(n):
... q = 1/n
... n -= 1
... yield 0
... while True:
... for i in range(n):
... yield 0
... yield q
...
sage: def lhs_gen():
... n = 1
... while True:
... yield FormalPowerSeries(Stream(m(n)))
... n += 1
...
sage: def rhs_gen():
... n = 1
... while True:
... yield FormalPowerSeries(Stream(s(n)))
... n += 1
...
sage: lhs = fps_product_generator(lhs_gen())
sage: rhs = fps_sum_generator(rhs_gen()).exponentiate()
sage: lhs.coefficients(10)
[1, 1, 2, 3, 5, 7, 11, 15, 22, 30]
sage: rhs.coefficients(10)
[1, 1, 2, 3, 5, 7, 11, 15, 22, 30]
#Test exponentiate
sage: def inv_factorial():
... q = 1
... yield 0
... yield q
... n = 2
... while True:
... q = q / n
... yield q
... n += 1
sage: f = FormalPowerSeries(Stream(inv_factorial())) #e^(x)-1
sage: u = f.exponentiate()
sage: g = inv_factorial()
sage: z1 = [1,1,2,5,15,52,203,877,4140,21147,115975]
sage: l1 = [z*g.next() for z in z1]
sage: l1 = [1] + l1[1:]
sage: u.coefficients(11)
[1, 1, 1, 5/6, 5/8, 13/30, 203/720, 877/5040, 23/224, 1007/17280, 4639/145152]
sage: l1 == u.coefficients(11)
True
"""
def uninitialized():
raise RuntimeError, "we should never be here"
def FormalPowerSeries(x=None):
if x is None:
return FormalPowerSeries_class(stream=None, order=unk, aorder=unk,
aorder_changed=True, compute_aorder=uninitialized,
is_initialized=False)
if hasattr(x, 'next'):
return FormalPowerSeries(Stream(x))
if isinstance(x, Stream):
return FormalPowerSeries_class(stream=x, order=unk, aorder=0,
aorder_changed=False, compute_aorder=do_nothing,
is_initialized=True)
if x == 0:
return zero
if x == 1:
return one
if x == 'monom':
return monom
raise TypeError, "unable to create a formal power series from x (= %s)"%x
class FormalPowerSeries_class:
def __init__(self, stream=None, order=None, aorder=None, aorder_changed=True, is_initialized=False, compute_aorder=None, name=None):
self.stream = stream
if order is None:
self.order = unk
else:
self.order = order
if aorder is None:
self.aorder = unk
else:
self.aorder = aorder
if self.aorder == inf:
self.order = inf
self.aorder_changed = aorder_changed
self.is_initialized = is_initialized
self.compute_aorder = compute_aorder
self._name = name
def __repr__(self):
if self._name:
return self._name
if self.is_initialized:
l = str(self.stream.list)
else:
l = 'uninitialized'
return "Formal power series: %s"%l
def initialize(self):
if self.order != unk:
return
if self.aorder == inf:
raise ValueError, "aorder can never be infinity since order would have to be infinity"
elif self.aorder != unk and self.is_initialized:
#Try to improve the approximate order
i = self.aorder
c = self.stream
n = len(c)
if i == 0 and n > 0:
if self.stream[i] == 0:
i += 1
self.aorder += 1
#Try to recognize the zero series
if i >= n:
#For non-constant series, we cannot do anythin
if not self.stream.is_constant():
return
if self.stream[ len(self.stream) - 1] == 0:
self.aorder = inf
self.order = inf
return
while i < n and c[i] == 0:
i += 1
self.aorder = i
if i < n:
self.order = i
else:
#if self.is_initialized != False:
# raise ValueError, "the series must be uninitialized at this point"
self.compute_aorder()
#if self.is_initialized != True:
#raise ValueError, "the series must be initialized at this point"
if hasattr(self, '_references_me'):
self._references_me._copy(self)
def initialize_coefficient_stream(self,compute_coefficients):
ao = self.aorder
assert ao != unk
if ao == inf:
self.order = inf
self.stream = Stream(0)
else:
self.stream = Stream(gen=compute_coefficients(ao))
self.is_initialized = True
def coefficients(self, n):
return [self.coefficient(i) for i in range(n)]
def zero(self):
self.initialize()
return self.order == inf
def set_approximate_order(self, new_order):
self.aorder_changed = ( self.aorder != new_order )
self.aorder = new_order
return self.aorder_changed
def _copy(self, x):
self.order = x.order
self.aorder = x.aorder
self.aorder_changed = x.aorder_changed
self.compute_aorder = x.compute_aorder
self.is_initialized = x.is_initialized
self.stream = x.stream
def set(self, x):
self._copy(x)
x._references_me = self
reprx = repr(x)
if "=" in reprx:
self._name += " = (" + reprx +")"
else:
self._name += " = " + reprx
def coefficient(self, n):
#pdb.set_trace()
if self.zero():
return 0
if n < self.aorder:
return 0
assert self.is_initialized
return self.stream[n]
def get_aorder(self):
self.initialize()
return self.aorder
def get_order(self):
self.initialize()
return self.order
def get_stream(self):
self.initialize()
return self.stream
def approximate_order_two(self, x, y, new_order, compute_coefficients):
"""approximate_order_two(\n%s, \n%s, \n%s, \n%s)"""%(x,y,new_order,compute_coefficients)
if self.is_initialized:
return
ochanged = self.aorder_changed
must_initialize_coefficient_stream = ( self.aorder == unk or self.is_initialized == False)
def a_order(x,y):
#print "(%s,%s) -> %s"%(x.aorder, y.aorder, new_order(x.aorder, y.aorder))
return new_order(x.aorder, y.aorder)
ao = a_order(x,y)
#start at maximum
if ao == unk:
ao = inf
tchanged = self.set_approximate_order(ao)
if ochanged or tchanged:
x.compute_aorder()
y.compute_aorder()
ao = a_order(x,y)
tchanged = self.set_approximate_order(ao)
#assert self.is_initialized == False
if must_initialize_coefficient_stream:
self.initialize_coefficient_stream(compute_coefficients)
if hasattr(self, '_references_me'):
self._references_me._copy(self)
def approximate_order_one(self, x, new_order, compute_coefficients):
"""approximate_order_two(\n%s, \n%s, \n%s)"""%(x,new_order,compute_coefficients)
if self.is_initialized:
return
ochanged = self.aorder_changed
must_initialize_coefficient_stream = ( self.aorder == unk or self.is_initialized == False)
def a_order(z):
#print "%s -> %s"%(z.aorder, new_order(z.aorder))
return new_order(z.aorder)
ao = a_order(x)
#start at maximum
if ao == unk:
ao = inf
tchanged = self.set_approximate_order(ao)
if ochanged or tchanged:
x.compute_aorder()
ao = a_order(x)
tchanged = self.set_approximate_order(ao)
assert self.is_initialized == False
if must_initialize_coefficient_stream:
self.initialize_coefficient_stream(compute_coefficients)
if hasattr(self, '_references_me'):
self._references_me._copy(self)
def approximate_order_zero(self, new_order, compute_coefficients):
"""approximate_order_zero(\n%s, \n%s)"""%(new_order,compute_coefficients)
if self.is_initialized:
return
assert self.aorder == unk
ao = new_order()
self.set_approximate_order(ao)
self.initialize_coefficient_stream(compute_coefficients)
if hasattr(self, '_references_me'):
self._references_me._copy(self)
def new0(self,compute_coefficients, order_op, name=None):
new_fps = self.new_from_ao(name=name)
def f0():
"""new0 - f0: %s, %s"""%(order_op, compute_coefficients)
new_fps.approximate_order_zero(order_op, compute_coefficients)
new_fps.compute_aorder = f0
return new_fps
def new1(self, compute_coefficients, order_op, x, name=None):
new_fps = self.new_from_ao( name=name )
def f1():
"""new1 - f1: %s, %s, %s"""%(x, order_op, compute_coefficients)
new_fps.approximate_order_one(x, order_op, compute_coefficients)
new_fps.compute_aorder = f1
return new_fps
def new2(self,compute_coefficients, order_op, x, y, name=None):
new_fps = self.new_from_ao(name=name )
def f2():
new_fps.approximate_order_two(x, y, order_op, compute_coefficients)
f2.__doc__ = """new2 - f2: %s, %s, %s, %s"""%(x, y, order_op, compute_coefficients)
new_fps.compute_aorder = f2
return new_fps
def new_from_ao(self, name=None):
return FormalPowerSeries_class(stream=None, order=unk, aorder=self.aorder,
aorder_changed=True, compute_aorder=None,
is_initialized=False, name=name)
def __add__(self, y):
def plus(ao):
for n in range(ao):
yield 0
n = ao
while True:
yield self.coefficient(n) + y.coefficient(n)
n += 1
return self.new2(plus, min, self, y, name=repr(self)+" + " +repr(y))
def __mul__(self, y):
def times(ao):
for n in range(ao):
yield 0
n = ao
while True:
low = self.aorder
high = n - y.aorder
nth_coefficient = 0
#Handle the zero series
if low == inf or high == inf: #
yield 0 #
n += 1 #
continue #
for k in range(low, high+1):
cx = self.coefficient(k)
if cx == 0:
continue
nth_coefficient += cx * y.coefficient(n-k)
yield nth_coefficient
n += 1
reprself = repr(self)
repry = repr(y)
if "+" in reprself:
reprself = "(" + reprself + ")"
if "+" in repry:
repry = "(" + repry + ")"
return self.new2(times, lambda a,b:a+b, self, y, name=reprself+"*"+repry)
def __call__(self, y):
def compose(ao):
assert y.coefficient(0) == 0
yield self.coefficient(0)
def g():
n = 1
while True:
yield self.coefficient(n)
n += 1
rest_x = FormalPowerSeries(g())
z = rest_x.compose(y)*y
is_zero = z.zero()
n = 1
while True:
yield z.coefficient(n)
n += 1
return self.new2(compose, lambda a,b:a*b, self, y, name="(%s)o(%s)"%(repr(self), repr(y)))
def add(self, y):
return self + y
def times(self, y):
return self * y
def compose(self, y):
return self(y)
def differentiate(self):
def diff(ao):
#Fixme! (take care of approximate order stuff)
#for n in range(ao-1):
# yield 0
n = 1
while True:
yield n*self.coefficient(n)
n += 1
return self.new1(diff, bounded_decrement, self, name="D(%s)"%repr(self))
def integrate(self, integration_constant = 0):
def integrate(ao):
for n in range(ao):
yield 0
n = ao
while True:
yield (Integer(1)/Integer(n))*self.coefficient(n-1)
n += 1
if integration_constant == 0:
return self.new1(integrate, increment, self, name="I(%s)"%repr(self))
else:
return _new(0, Stream(self._integrate_nonzero(integration_constant)), name="I(%s,%s)"%(repr(self), integration_constant))
def _integrate_nonzero(self, integration_constant):
yield integration_constant
aox = self.aorder
assert aox != unk
if aox == inf:
while True:
yield 0
else:
for n in range(aox):
yield 0
n = aox + 1
while True:
yield (Integer(1)/Integer(n))*self.coefficient(n-1)
n += 1
def exponentiate(self):
s = FormalPowerSeries(); s._name = "s"
s.set( (self.differentiate()*s).integrate(1) )
s._name = "EXP(%s)"%repr(self)
return s
def do_nothing(*args, **kwargs):
## assert t.aorder != unk
## assert t.is_initialized == True
return
def new_from_ao(approximate_order):
return FormalPowerSeries_class(stream=None, order=unk, aorder=unk,
aorder_changed=True, compute_aorder=approximate_order,
is_initialized=False)
#Finite Sum
def sum_coefficients(a):
last = len(a) - 1
assert last >= 0
n = 0
while True:
r = sum( [a[i].coefficient(n) for i in range(len(a))] )
yield r
n += 1
def fps_sum(a):
return FormalPowerSeries( sum_coefficients(a) )
#Potentially infinite sum
def fps_sum_generator(g):
s = Stream(g)
def h():
n = 0
while True:
r = s[n].coefficient(n) #compute [x^n] s(n)
for i in range(len(s)-1):
r += s[i].coefficient(n)
yield r
n += 1
return FormalPowerSeries(h())
#Potentially infinite product
def fps_product_generator(g):
def h():
z = g.next()
yield z.coefficient(0)
yield z.coefficient(1)
n = 2
for x in g:
z = z * x
yield z.coefficient(n)
n += 1
while True:
yield z.coefficient(n)
n += 1
return FormalPowerSeries(h())
#################################
def _new(o, g, name=None):
return FormalPowerSeries_class(stream=g, order=o, aorder=o,
aorder_changed=False, compute_aorder=do_nothing,
is_initialized=True, name=name)
zero = _new(inf, Stream([0]), name="0")
one = _new(0, Stream([1,0]), name="1")
monom = _new(1, Stream([0,1,0]), name="monom")
def term(r, n):
if r == 0:
return zero
else:
s = [0]*n+[r,0]
return _new(int(n), Stream(s), name="term(%s,%s)"%(r,n))
from sage.rings.all import Integer
class SeriesOrderElement:
pass
class InfiniteSeriesOrder(SeriesOrderElement):
def __init__(self):
pass
def __repr__(self):
return "Infinite series order"
def __add__(self, x):
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return self
if isinstance(x, InfiniteSeriesOrder):
return self
if isinstance(x, UnknownSeriesOrder):
return x
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
__radd__ = __add__
def __mul__(self, x):
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return self
if isinstance(x, InfiniteSeriesOrder):
return self
if isinstance(x, UnknownSeriesOrder):
return x
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
def __rmul__(self, x):
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
if x == 0:
return 0
else:
return self
if isinstance(x, InfiniteSeriesOrder):
return self
if isinstance(x, UnknownSeriesOrder):
return x
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
def __lt__(self, x):
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return False
if isinstance(x, InfiniteSeriesOrder):
return False
if isinstance(x, UnknownSeriesOrder):
return False
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
def __gt__(self, x):
return True
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return False
if isinstance(x, InfiniteSeriesOrder):
return False
if isinstance(x, UnknownSeriesOrder):
return False
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
class UnknownSeriesOrder(SeriesOrderElement):
def __init__(self):
pass
def __repr__(self):
return "Unknown series order"
def __add__(self, x):
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return self
if isinstance(x, SeriesOrderElement):
return self
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
__radd__ = __add__
def __mul__(self, x):
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return self
if isinstance(x, SeriesOrderElement):
return self
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
def __rmul__(self, x):
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return self
if isinstance(x, SeriesOrderElement):
return self
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
def __lt__(self, x):
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return True
if isinstance(x, SeriesOrderElement):
return True
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
def __gt__(self, x):
return False
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return False
if isinstance(x, InfiniteSeriesOrder):
return False
if isinstance(x, UnknownSeriesOrder):
return False
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
def bounded_decrement(x):
if isinstance(x, SeriesOrderElement):
return x
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return max(0, x - 1)
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
def increment(x):
if isinstance(x, SeriesOrderElement):
return x + 1
if isinstance(x, (int, Integer)):
if x < 0:
raise ValueError, "x must be a positive integer"
return x+1
raise TypeError, "x must be a positive integer or a SeriesOrderElement"
import itertools
def _integers_from(n):
while True:
yield n
n += 1
class Stream:
"""
EXAMPLES:
sage: from itertools import izip
sage: s = Stream(0)
sage: len(s)
1
sage: [x for x in s.elements()]
[0]
sage: [x for (x,i) in izip(s, range(4))]
[0, 0, 0, 0]
sage: len(s)
1
sage: l = [4,3,5,7,4,1,9,7]
sage: s = Stream(l)
sage: s[3]
8
sage: len(s)
4
sage: s[3]
8
sage: len(s)
4
sage: s[1]
3
sage: len(s)
4
sage: s[4]
4
sage: len(s)
5
sage: [i for i in s.elements()]
[4, 3, 5, 7, 4, 1, 9, 7]
sage: len(s)
8
sage: l = [4,3,5,2,6,1]
sage: s = Stream(l)
sage: s[3]
2
sage: len(s)
4
sage: g = iter(l)
sage: s.set_gen(g)
sage: s[5]
3
sage: len(s)
6
sage: [x for x in s.elements()]
[4, 3, 5, 2, 6, 1, 4, 3, 5, 2]
sage: len(s)
10
sage: s = Stream(const=4)
sage: g = iter(s)
sage: n = s.elements()
sage: l1 = [x for (x,i) in izip(g, range(10))]
sage: l2 = [x for (x,i) in izip(n, range(10))]
sage: l = [4 for k in range(10)]
sage: l == l1
True
sage: l2
[4]
sage: t = Stream([2,3,5,7,11])
sage: g3 = t.elements(2)
sage: l3 = [x for (x,i) in izip(g3, range(10))]
sage: l3
[5, 7, 11]
sage: g4 = t.elements(10)
sage: l4 = [x for (x,i) in izip(g4, range(10))]
sage: l4
[11]
sage: l = ['Hello', ' ', 'World', '!']
sage: s = Stream(l)
sage: len(s)
0
sage: s[2]
"World"
sage: len(s)
3
sage: u = ""
sage: for i in range(len(s)): u += s[i]
sage: u
"Hello World"
sage: v = ""
sage: for i in range(10): v += s[i]
sage: v
"Hello World!!!!!!!"
sage: len(s)
4
sage: l = [2,3,5,7,11,0]
sage: s = Stream(l)
sage: s.is_constant()
False
sage: s[3]
7
sage: s.is_constant()
False
sage: s[5]
0
sage: s.is_constant()
False
sage: s[6]
0
sage: s.is_constant()
True
sage: s = Stream(const='I am constant.')
sage: s.is_constant()
True
sage: s = Stream([1,2,3,4,5,6])
sage: a = []
sage: b = []
sage: for (i,x) in enumerate(s.elements()):
a = [x] + a
b = [ s[2*i] ] + b
sage: a
[6, 5, 4, 3, 2, 1]
sage: b
[6, 6, 6, 5, 3, 1]
sage: fib = Stream()
sage: def g():
yeild 1
yield 1
n = 0
while True:
yield fib[n] + fib[n+1]
n += 1
sage: fib.set_gen(g())
sage: [x for (x,i) in izip(fib, range(11))]
[1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 69]
sage: def h(l):
if len(l) < 2:
return 1
return l[-1] + l[-2]
sage: fib = Stream(func = h)
sage: [x for (x,i) in izip(fib, range(11))]
[1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 69]
sage: s = Stream()
sage: l = Stream([1, 0])
sage: y = Stream([0,1,0])
sage: s.set_gen(iter(l + y * s * s))
sage: [x for (x,i) in izip(s, range(6))]
[1, 1, 2, 5, 14, 42]
sage: def f(l):
if len(l) == 0:
return 0
return l[-1] + 1
sage: o = Stream(func=f)
sage: [x for (x,i) in izip(o, range(10))]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
sage: double = lambda z: 2*z
sage: t = o.map(double)
sage: [x for (x,i) in izip(t, range(10))]
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18]
sage: double = lambda z: 2*z
sage: o = Stream([0,1,2,3])
sage: [x for (x,i) in izip(o, range(6))]
[0, 1, 2, 3, 3, 3]
sage: t = o.map(double)
sage: [x for (x,i) in izip(t, range(6))]
[0, 2, 4, 6, 6, 6]
sage: r = [4, 3, 5, 2, 6, 1, 1, 1, 1, 1]
sage: l = [4, 3, 5, 2, 6, 1]
sage: s = Stream(l)
sage: s[3] = -1
sage: [x for (x,i) in izip(s, r)]
[4, 3, 5, -1, 6, 1, 1, 1, 1, 1]
sage: s[5] = -2
sage: [x for (x,i) in izip(s, r)]
[4, 3, 5, -1, 6, -2, 1, 1, 1, 1]
sage: s[6] = -3
sage: [x for (x,i) in izip(s, r)]
[4, 3, 5, -1, 6, -2, -3, 1, 1, 1]
sage: s[8] = -4
sage: [x for (x,i) in izip(s, r)]
[4, 3, 5, -1, 6, -2, -3, 1, -4, 1]
sage: a = Stream(const=0)
sage: a[2] = 3
sage: [x for (x,i) in izip(a, range(4))]
[0, 0, 3, 0]
"""
def __init__(self, gen=None, const=None, list=None, constant=True, func=None):
if func != None:
if gen != None:
raise ValueError, "you cannot specify both a function and a generator"
def g():
while True:
try:
yield func(self.list)
except:
break
gen = g()
if hasattr(gen, '__iter__'):
gen = iter(gen)
#Constant stream
if const != None:
self.list = [ const ]
self.last_index = 0
self.gen = None
self.constant = const
self.end_reached = True
else:
self.list = [ ]
self.last_index = -1
self.gen = gen
self.constant = constant
self.end_reached = False
def __setitem__(self, i, t):
self[i]
if i < len(self.list):
self.list[i] = t
else:
self.list += [ self.constant ] * (i+1 - len(self.list))
self.last_index = i
self.list[i] = t
def set_gen(self, gen):
self.gen = gen
self.end_reached = False
def map(self, f):
return Stream(gen=(f(x) for x in self))
def __getitem__(self, i):
if i <= self.last_index:
return self.list[i]
elif self.end_reached:
if self.constant != "False":
return self.constant
else:
raise IndexError, "out of position"
else:
while self.last_index < i:
try:
self.list.append(self.gen.next())
self.last_index += 1
except StopIteration:
self.end_reached = True
self.constant = self.list[-1]
return self[i]
return self.list[i]
def __add__(self, ds):
if not isinstance(ds, Stream):
raise TypeError, "ds must be a Stream"
return Stream(gen=(a+b for (a,b) in itertools.izip(self,ds)))
def __mul__(self, ds):
if not isinstance(ds, Stream):
raise TypeError, "ds must be a Stream"
def times(n):
t = 0
for k in range(n+1):
sk = self[k]
if sk == 0:
continue
t += sk *ds[n-k]
return t
return Stream(gen=(times(n) for n in _integers_from(0)))
def __iter__(self):
i = 0
while True:
try:
yield self[i]
except IndexError:
break
i += 1
def __len__(self):
return len(self.list)
def number_computed(self):
return self.__len__()
def data(self):
return self.list
def is_constant(self):
return self.end_reached
def elements(self, start_index=0):
i = start_index
while not self.end_reached:
try:
self[i]
except IndexError:
break
i += 1
if start_index >= len(self.list):
return iter(self.list[-1:])
else:
return iter(self.list[start_index:])
## if start_index < 0:
## raise ValueError, "start_index must be >= 0"
## i = start_index
## while True:
## print i, self.last_index, len(self.list), self.end_reached
## z = self[i]
## print i, self.last_index, len(self.list), self.end_reached
## print "-"*20
## if self.end_reached and i >= len(self.list):
## if start_index >= len(self.list):
## yield z
## break
## yield z
## i += 1
one = Stream(const=1)
zero = Stream(const=0)
x = Stream([0,1,0])
ints = Stream(_integers_from(0))
