| 1 | """ |
| 2 | Cached formal powerseries and formal Laurant series in one variable. |
| 3 | |
| 4 | Author: Henryk Trappmann |
| 5 | """ |
| 6 | |
| 7 | from sage.structure.sage_object import SageObject |
| 8 | from sage.rings.arith import factorial |
| 9 | from sage.rings.arith import binomial |
| 10 | from sage.rings.integer import Integer |
| 11 | from sage.calculus.calculus import var |
| 12 | from sage.symbolic.expression import Expression |
| 13 | from sage.calculus.functional import diff |
| 14 | from sage.rings.ring import Ring |
| 15 | from sage.rings.ring_element import RingElement |
| 16 | from sage.rings.rational_field import QQ, RationalField |
| 17 | from sage.rings.rational import Rational |
| 18 | from sage.rings.real_mpfr import RR, RealField |
| 19 | from sage.rings.complex_field import ComplexField_class |
| 20 | from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing |
| 21 | from sage.rings.polynomial.polynomial_ring import PolynomialRing_field |
| 22 | from sage.misc.misc_c import prod |
| 23 | from sage.misc.functional import n |
| 24 | from sage.rings.infinity import Infinity |
| 25 | from sage.rings.power_series_ring_element import PowerSeries |
| 26 | from sage.rings.polynomial.polynomial_element import Polynomial |
| 27 | from sage.matrix.constructor import matrix |
| 28 | |
| 29 | def decidable0(K): |
| 30 | """ |
| 31 | Returns true if K has a decidable 0. |
| 32 | For example powerseries have no decidable 0. |
| 33 | |
| 34 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing,decidable0 |
| 35 | sage: decidable0(QQ) |
| 36 | True |
| 37 | sage: decidable0(FormalPowerSeriesRing(QQ)) |
| 38 | False |
| 39 | """ |
| 40 | if K == Integer or K == int: |
| 41 | return True |
| 42 | if K == Rational: |
| 43 | return True |
| 44 | if isinstance(K,RationalField): |
| 45 | return True |
| 46 | if isinstance(K,RealField): |
| 47 | return True |
| 48 | if isinstance(K,PolynomialRing_field): |
| 49 | return True |
| 50 | # if isinstance(K,RealField): |
| 51 | # return true |
| 52 | # if isinstance(K,ComplexField_class): |
| 53 | # return true |
| 54 | return False |
| 55 | |
| 56 | def _isnat(n): |
| 57 | """ |
| 58 | Returns True if n is a non-negative int or Integer. |
| 59 | |
| 60 | sage: from sage.rings.formal_powerseries import _isnat |
| 61 | sage: _isnat(5r) |
| 62 | True |
| 63 | sage: _isnat(5) |
| 64 | True |
| 65 | sage: _isnat(0) |
| 66 | True |
| 67 | sage: _isnat(-1r) |
| 68 | False |
| 69 | """ |
| 70 | if isinstance(n,int) or isinstance(n,Integer): |
| 71 | if n >= 0: |
| 72 | return True |
| 73 | return False |
| 74 | |
| 75 | def _assert_nat(n): |
| 76 | """ |
| 77 | Throws an exception if not _isnat(n). |
| 78 | sage: from sage.rings.formal_powerseries import _assert_nat |
| 79 | sage: #try: |
| 80 | sage: # _assert_nat(-1) |
| 81 | sage: #except AssertionError: |
| 82 | sage: # print 'bang' |
| 83 | """ |
| 84 | assert _isnat(n), repr(n)+ " must be natural number." |
| 85 | |
| 86 | class FormalPowerSeriesRing(Ring): |
| 87 | def by_lambda(self,f,min_index=0): |
| 88 | """ |
| 89 | Returns the powerseries with coefficients f(n). |
| 90 | |
| 91 | Alternative expression: |
| 92 | self.by_lambda(f) == self(f) |
| 93 | self.by_lambda(f,min_index) == self(f,min_index) |
| 94 | |
| 95 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 96 | sage: P = FormalPowerSeriesRing(QQ) |
| 97 | sage: P.by_lambda(lambda n: 1/factorial(n)) |
| 98 | [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...] |
| 99 | |
| 100 | Initialization as formal Laurant series: |
| 101 | sage: P.by_lambda(lambda n: n,-3) |
| 102 | [-3, -2, -1; 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, ...] |
| 103 | |
| 104 | Corruptly setting min_index=3 |
| 105 | sage: P.by_lambda(lambda n: n,3) |
| 106 | [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...] |
| 107 | |
| 108 | Note that functions can not be serialized/pickled. |
| 109 | If you want to have a serializable/picklable object, you can derive |
| 110 | from FormalPowerSeries and define the method coeffs. |
| 111 | |
| 112 | sage: from sage.rings.formal_powerseries import FormalPowerSeries |
| 113 | sage: #class F(FormalPowerSeries):^J def coeffs(self,n): return n |
| 114 | sage: #F(P) |
| 115 | """ |
| 116 | return FormalPowerSeries(self,f,min_index) |
| 117 | |
| 118 | |
| 119 | def by_iterator(self,g,min_index=0): |
| 120 | """ |
| 121 | Returns a powerseries from the generator g. |
| 122 | The powerseries coefficients start at min_index |
| 123 | which is allowed be negative obtaining a formal Laurant series. |
| 124 | |
| 125 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 126 | sage: P = FormalPowerSeriesRing(QQ) |
| 127 | sage: P.by_iterator(iter(ZZ),-2) |
| 128 | [0, 1; -1, 2, -2, 3, -3, 4, -4, 5, -5, 6, -6, 7, -7, 8, -8, 9, -9, 10, -10, ...] |
| 129 | """ |
| 130 | |
| 131 | return Iterated(self,g,min_index) |
| 132 | |
| 133 | def by_undefined(self,min_index=0): |
| 134 | """ |
| 135 | Returns an undefined powerseries. This is useful to use method `define'. |
| 136 | |
| 137 | Alternative expressions: |
| 138 | self.by_undefined() == self() |
| 139 | self.by_undefined(m) == self(None,m) |
| 140 | |
| 141 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 142 | sage: P = FormalPowerSeriesRing(QQ) |
| 143 | sage: a = P.by_undefined() |
| 144 | sage: a |
| 145 | Undefined |
| 146 | sage: P.by_undefined(2).min_index |
| 147 | 2 |
| 148 | """ |
| 149 | return FormalPowerSeries(self,min_index=min_index) |
| 150 | |
| 151 | def by_list(self,list,start=0): |
| 152 | """ |
| 153 | Returns the powerseries with coefficients p[n] where |
| 154 | p[n]==0 for 0<=n<start, p[m+start]==list[m] for all list indices m, |
| 155 | and p[n]==0 for all later indices n. |
| 156 | |
| 157 | Alternative expression: |
| 158 | self.by_list(list) == self(list) |
| 159 | self.by_list(list,start) == self(list,start) |
| 160 | |
| 161 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 162 | sage: P = FormalPowerSeriesRing(QQ) |
| 163 | sage: f = P.by_list([1,2,3,4,5]) |
| 164 | sage: f |
| 165 | [1, 2, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 166 | sage: P.by_list([1,2,3],5) |
| 167 | [0, 0, 0, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 168 | """ |
| 169 | return List(self,list,start).reclass() |
| 170 | |
| 171 | def by_polynomial(self,p): |
| 172 | """ |
| 173 | Returns the FormalPowerSeries from the given Polynomial. |
| 174 | Alternative expression: self.by_polynomial(p) == self(p) |
| 175 | |
| 176 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 177 | sage: PP = PolynomialRing(QQ,x) |
| 178 | sage: P = FormalPowerSeriesRing(QQ) |
| 179 | sage: pp = PP([1,2,3]) |
| 180 | sage: P.by_polynomial(pp) |
| 181 | [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 182 | sage: x = PP(x) |
| 183 | sage: P(2*x+x**2) |
| 184 | [0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 185 | """ |
| 186 | return self.by_list(p.padded_list()) |
| 187 | |
| 188 | def by_powerseries(self,p): |
| 189 | """ |
| 190 | Returns the FormalPowerSeries from the given PowerSeries. |
| 191 | Alternative expression: self.py_powerseries(ps)==self(ps) |
| 192 | |
| 193 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 194 | sage: PS = PowerSeriesRing(QQ,x) |
| 195 | sage: FPS = FormalPowerSeriesRing(QQ) |
| 196 | sage: FPS.by_powerseries(PS([0,1,2,3])) |
| 197 | [0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 198 | """ |
| 199 | return self.by_polynomial(p.polynomial()) |
| 200 | |
| 201 | def by_taylor(self,expr,v,at=0,**kwargs): |
| 202 | """ |
| 203 | Returns the taylor series of `expr' with respect to `v' at `at'. |
| 204 | |
| 205 | Alternative expressions: |
| 206 | self.by_taylor(expr,v) == self(expr,v) |
| 207 | self.by_taylor(expr,v,at) == self(expr,v,at) |
| 208 | |
| 209 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 210 | sage: P = FormalPowerSeriesRing(QQ) |
| 211 | sage: f = P.by_taylor(ln(x),x,1) |
| 212 | sage: f |
| 213 | [0, 1, -1/2, 1/3, -1/4, 1/5, -1/6, 1/7, -1/8, 1/9, -1/10, 1/11, -1/12, ...] |
| 214 | """ |
| 215 | |
| 216 | #print "after",min_index |
| 217 | return Taylor(self,expr,v,at).reclass() |
| 218 | |
| 219 | def by_constant(self,c,**kwargs): |
| 220 | """ |
| 221 | Returns the powerseries with coefficients [c,0,0,...]. |
| 222 | |
| 223 | Alternative expression: self.by_constant(c) == self(c) |
| 224 | |
| 225 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 226 | sage: P = FormalPowerSeriesRing(QQ) |
| 227 | sage: f = P.by_constant(42) |
| 228 | sage: f |
| 229 | [42, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 230 | """ |
| 231 | |
| 232 | return Constant(self,c) |
| 233 | |
| 234 | |
| 235 | def __call__(self,p1=None,p2=None,p3=None,**kwargs): |
| 236 | """ |
| 237 | Initialization by finite sequence of coefficients: |
| 238 | Examples: |
| 239 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 240 | sage: PQ = FormalPowerSeriesRing(QQ) |
| 241 | sage: PQ([1,2,3]) |
| 242 | [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 243 | sage: PQ([]) |
| 244 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 245 | |
| 246 | Initialization by coefficient function: |
| 247 | Example: |
| 248 | sage: PQ(lambda n: 1/factorial(n)) |
| 249 | [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...] |
| 250 | |
| 251 | Initialization by expresion: |
| 252 | Examples: |
| 253 | sage: PQ(1+2*x+3*x^2,x) |
| 254 | [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 255 | sage: PQ(exp(x),x) |
| 256 | [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...] |
| 257 | sage: PQ(ln(x),x,1) |
| 258 | [0, 1, -1/2, 1/3, -1/4, 1/5, -1/6, 1/7, -1/8, 1/9, -1/10, 1/11, -1/12, ...] |
| 259 | |
| 260 | Note: This is much slower than directly providing the coefficient function. |
| 261 | |
| 262 | See also methods: by_const, by_undef, by_list, by_taylor, by_lambda |
| 263 | """ |
| 264 | |
| 265 | if isinstance(p1,Integer) or isinstance(p1,int) or isinstance(p1,Rational): |
| 266 | return self.by_constant(p1,**kwargs) |
| 267 | |
| 268 | if isinstance(p1,list): |
| 269 | if p2 == None: |
| 270 | return self.by_list(p1,**kwargs) |
| 271 | return self.by_list(p1,p2,**kwargs) |
| 272 | |
| 273 | if isinstance(p1,Expression): |
| 274 | if p3 == None: |
| 275 | return self.by_taylor(p1,p2,**kwargs) |
| 276 | return self.by_taylor(p1,p2,p3,**kwargs) |
| 277 | |
| 278 | if isinstance(p1,Polynomial): |
| 279 | return self.by_polynomial(p1) |
| 280 | |
| 281 | if isinstance(p1,PowerSeries): |
| 282 | return self.by_powerseries(p1) |
| 283 | |
| 284 | #TODO generator if isinstance(p1, |
| 285 | |
| 286 | if type(p1) is type(lambda n: 0): |
| 287 | if p2 == None: |
| 288 | return self.by_lambda(p1,**kwargs) |
| 289 | return self.by_lambda(p1,p2,**kwargs) |
| 290 | |
| 291 | if p1 == None: |
| 292 | return self.by_undefined(p2) |
| 293 | |
| 294 | raise TypeError, "unrecognized initialization input" + repr(type(p1)) |
| 295 | |
| 296 | def is_field(self): |
| 297 | """ |
| 298 | Returns True if self is a field, i.e. if it can be used as |
| 299 | formal laurant series. |
| 300 | |
| 301 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 302 | sage: FormalPowerSeriesRing(IntegerRing()).is_field() |
| 303 | False |
| 304 | sage: FormalPowerSeriesRing(QQ).is_field() |
| 305 | True |
| 306 | """ |
| 307 | return self.K.is_field() |
| 308 | |
| 309 | def is_commutative(self): |
| 310 | """ |
| 311 | The powerseries is commutative if the base ring is. |
| 312 | |
| 313 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 314 | sage: FormalPowerSeriesRing(IntegerRing()).is_commutative() |
| 315 | True |
| 316 | """ |
| 317 | return self.K.is_commutative() |
| 318 | |
| 319 | def is_exact(self): |
| 320 | """ |
| 321 | The powerseries is exact if the base ring is. |
| 322 | |
| 323 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 324 | sage: FormalPowerSeriesRing(RR).is_exact() |
| 325 | False |
| 326 | sage: FormalPowerSeriesRing(QQ).is_exact() |
| 327 | True |
| 328 | """ |
| 329 | return self.K.is_exact() |
| 330 | |
| 331 | def is_finite(self): |
| 332 | """ |
| 333 | Powerseries ring is never finite (except the base ring has only 1 element). |
| 334 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 335 | sage: FormalPowerSeriesRing(QQ).is_finite() |
| 336 | False |
| 337 | """ |
| 338 | return False |
| 339 | |
| 340 | def zero_element(self): |
| 341 | """ |
| 342 | Returns the zero element of this power series ring. |
| 343 | |
| 344 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 345 | sage: P = FormalPowerSeriesRing(QQ) |
| 346 | sage: P.zero_element() |
| 347 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 348 | """ |
| 349 | return self.Zero |
| 350 | |
| 351 | def one_element(self): |
| 352 | """ |
| 353 | Returns the one element of this power series ring. |
| 354 | |
| 355 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 356 | sage: P = FormalPowerSeriesRing(QQ) |
| 357 | sage: P.one_element() |
| 358 | [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 359 | """ |
| 360 | return self.One |
| 361 | |
| 362 | def __init__(self,base_ring): |
| 363 | """ |
| 364 | Returns the powerseries ring over base_ring. |
| 365 | |
| 366 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 367 | sage: FormalPowerSeriesRing(QQ) |
| 368 | FormalPowerSeriesRing over Rational Field |
| 369 | """ |
| 370 | if base_ring == None: |
| 371 | return |
| 372 | |
| 373 | self.K = base_ring |
| 374 | |
| 375 | def PSS(seq): |
| 376 | """ sage: None # indirect doctest """ |
| 377 | return self.by_list(seq) |
| 378 | |
| 379 | if self.K == int: |
| 380 | self.K = Integer |
| 381 | |
| 382 | K = self.K |
| 383 | K1 = K.one_element() |
| 384 | self.K1 = K1 |
| 385 | |
| 386 | self.Zero = Zero(self) |
| 387 | self.One = One(self) |
| 388 | self.Id = Id(self,min_index=1) |
| 389 | self.Inc = Inc(self) |
| 390 | self.Dec = Dec(self) |
| 391 | self.Exp = Exp(self) |
| 392 | self.Dec_exp = Dec_exp(self,min_index=1) |
| 393 | self.Log_inc = Log_inc(self,min_index=1) |
| 394 | self.Sin = Sin(self,min_index=1) |
| 395 | self.Cos = Cos(self) |
| 396 | self.Arcsin = Arcsin(self,min_index=1) |
| 397 | self.Arctan = Arctan(self,min_index=1) |
| 398 | self.Sinh = Sinh(self,min_index=1) |
| 399 | self.Cosh = Cosh(self) |
| 400 | self.Arcsinh = Arcsinh(self,min_index=1) |
| 401 | self.Arctanh = Arctanh(self,min_index=1) |
| 402 | self.Bernoulli = (self.Id / self.Exp.dec()).mul_fact() |
| 403 | self.Bernoulli.__doc__ = """ |
| 404 | The n-th Bernoulli number is equal to |
| 405 | the n-th derivative of 1/(exp(x)-1) at 0. |
| 406 | """ |
| 407 | self.Tan = Tan(self,min_index=1) |
| 408 | self.Tanh = Tanh(self,min_index=1) |
| 409 | self.Xexp = Xexp(self,min_index=1) |
| 410 | self.Lambert_w = Lambert_w(self,min_index=1) |
| 411 | self.Sqrt_inc = Sqrt_inc(self) |
| 412 | |
| 413 | #dont go into a recursion defining stirling1 |
| 414 | if isinstance(K,FormalPowerSeriesRing): |
| 415 | return |
| 416 | |
| 417 | self.Stirling1 = Stirling1(self) |
| 418 | |
| 419 | # def lehmer_comtet(n,k): #A008296 |
| 420 | # """ sage: None # indirect doctest """ |
| 421 | # r = 0 |
| 422 | # for l in range(k,n+1): |
| 423 | # r += binomial(l, k)*k**(l-k)*self.Stirling1[n][l] |
| 424 | # return K(r) |
| 425 | |
| 426 | self.Lehmer_comtet = Lehmer_comtet(self) |
| 427 | self.A000248 = self.Lehmer_comtet |
| 428 | |
| 429 | #self.selfpower_inc = PSF(lambda n: K(sum([ lehmer_comtet(n,k) for k in range(0,n+1))/factorial(n),K0)) |
| 430 | self.Selfpower_inc = Selfpower_inc(self) |
| 431 | |
| 432 | self.Superroot_inc = Superroot_inc(self) |
| 433 | |
| 434 | self.A003725 = A003725(self) |
| 435 | |
| 436 | self.Selfroot_inc = Selfroot_inc(self) |
| 437 | |
| 438 | self.Inv_selfroot_inc = Inv_selfroot_inc(self) |
| 439 | |
| 440 | def _repr_(self): |
| 441 | """ |
| 442 | Description of this FormalPowerSeriesRing. |
| 443 | |
| 444 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 445 | sage: FormalPowerSeriesRing(QQ)._repr_() |
| 446 | 'FormalPowerSeriesRing over Rational Field' |
| 447 | """ |
| 448 | return "FormalPowerSeriesRing over " + repr(self.K) |
| 449 | |
| 450 | def _coerce_map_from_(self,T): |
| 451 | """ |
| 452 | Returns true if type T can be coerced to a FormalPowerSeries |
| 453 | with self as parent. This can be done always when |
| 454 | T can be coerced to self.base_ring(). |
| 455 | This is used for operations like lmul and rmul. |
| 456 | |
| 457 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 458 | sage: P = FormalPowerSeriesRing(RR) |
| 459 | sage: P._coerce_map_from_(QQ) |
| 460 | True |
| 461 | sage: P._coerce_map_from_(int) |
| 462 | True |
| 463 | sage: P = FormalPowerSeriesRing(QQ) |
| 464 | sage: P._coerce_map_from_(int) |
| 465 | True |
| 466 | sage: P._coerce_map_from_(RR) |
| 467 | False |
| 468 | """ |
| 469 | #print self.K, T,not self.K.coerce_map_from(T) == None |
| 470 | if not self.K.coerce_map_from(T) == None: |
| 471 | return True |
| 472 | if isinstance(T,FormalPowerSeriesRing): |
| 473 | return not self.K.coerce_map_from(T.K) == None |
| 474 | return False |
| 475 | |
| 476 | |
| 477 | def base_ring(self): |
| 478 | """ |
| 479 | Returns the base ring of the powerseries ring. |
| 480 | |
| 481 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 482 | sage: FormalPowerSeriesRing(QQ).base_ring() == QQ |
| 483 | True |
| 484 | """ |
| 485 | return self.K |
| 486 | |
| 487 | def _test(self): |
| 488 | """ |
| 489 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 490 | sage: P = FormalPowerSeriesRing(QQ) |
| 491 | |
| 492 | sage: P.Bernoulli |
| 493 | [1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66, 0, -691/2730, 0, 7/6, ...] |
| 494 | |
| 495 | # #takes too long: |
| 496 | # sage: P(exp(x),x)-P.Exp |
| 497 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 498 | # sage: P(log(x+1),x) - P.Log_inc |
| 499 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 500 | # sage: P(cos(x),x)-P.Cos |
| 501 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 502 | # sage: P(sin(x),x)-P.Sin |
| 503 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 504 | # sage: P(arcsin(x),x) - P.Arcsin |
| 505 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 506 | # sage: P(arctan(x),x) - P.Arctan |
| 507 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 508 | # sage: P(sinh(x),x)-P.Sinh |
| 509 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 510 | # sage: P(cosh(x),x)-P.Cosh |
| 511 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 512 | # sage: P(arcsinh(x),x)-P.Arcsinh |
| 513 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 514 | # sage: P(arctanh(x),x)-P.Arctanh |
| 515 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 516 | # sage: P(sqrt(x+1),x)-P.Sqrt_inc |
| 517 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 518 | # sage: P(x*exp(x),x)-P.Xexp |
| 519 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 520 | # sage: P(exp(x)-1,x)-P.Exp.dec() |
| 521 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 522 | |
| 523 | sage: P.Exp.dec().reclass() | P.Log_inc |
| 524 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 525 | sage: P.Sin | P.Arcsin |
| 526 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 527 | sage: P.Tan | P.Arctan |
| 528 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 529 | sage: P.Tanh | P.Arctanh |
| 530 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 531 | sage: P.Xexp | P.Lambert_w |
| 532 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 533 | sage: P.Superroot_inc.dec().reclass() | P.Selfpower_inc |
| 534 | [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 535 | sage: P.Inv_selfroot_inc.dec().reclass() | P.Selfroot_inc |
| 536 | [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 537 | |
| 538 | sage: P.Superroot_inc ** P.Superroot_inc |
| 539 | [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 540 | sage: P.Tan - P.Sin / P.Cos |
| 541 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 542 | sage: P.Sin*P.Sin + P.Cos*P.Cos |
| 543 | [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 544 | sage: p = P([3,2,1]) |
| 545 | sage: p.rcp()*p |
| 546 | [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 547 | |
| 548 | |
| 549 | sage: P._test() |
| 550 | """ |
| 551 | pass |
| 552 | |
| 553 | # class UndefinedFormalPowerSeries(RingElement): |
| 554 | # """ |
| 555 | # Undefined powerseries. |
| 556 | # """ |
| 557 | # coeffs = None |
| 558 | # def _repr_(a): |
| 559 | # return "Undefined" |
| 560 | |
| 561 | class FormalPowerSeries(RingElement): |
| 562 | """ |
| 563 | Formal power series: |
| 564 | |
| 565 | A powerseries p is basically seen as an infinite sequence of coefficients |
| 566 | The n-th coefficient is retrieved by p[n]. |
| 567 | |
| 568 | EXAMPLES: |
| 569 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 570 | sage: PQ = FormalPowerSeriesRing(QQ) |
| 571 | sage: #Predefined PowerSeries |
| 572 | sage: expps = PQ.Exp |
| 573 | sage: expps.polynomial(10,x) |
| 574 | 1/362880*x^9 + 1/40320*x^8 + 1/5040*x^7 + 1/720*x^6 + 1/120*x^5 + 1/24*x^4 + 1/6*x^3 + 1/2*x^2 + x + 1 |
| 575 | sage: expps |
| 576 | [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...] |
| 577 | sage: PQ.Zero |
| 578 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 579 | sage: PQ.One |
| 580 | [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 581 | sage: PQ.Id |
| 582 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 583 | sage: #finite powerseries |
| 584 | sage: p = PQ([1,2,3,4]) |
| 585 | sage: p.polynomial(10,x) |
| 586 | 4*x^3 + 3*x^2 + 2*x + 1 |
| 587 | sage: p |
| 588 | [1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 589 | sage: one = PQ([1]) |
| 590 | sage: id = PQ([0,1]) |
| 591 | |
| 592 | sage: #power series are just functions that map the index to the coefficient |
| 593 | sage: expps[30] |
| 594 | 1/265252859812191058636308480000000 |
| 595 | |
| 596 | Formal Laurant Series can have negative minimum index |
| 597 | sage: PQ(lambda n: n,-5) |
| 598 | [-5, -4, -3, -2, -1; 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, ...] |
| 599 | |
| 600 | Power series operations |
| 601 | sage: p+p |
| 602 | [2, 4, 6, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 603 | sage: p-p |
| 604 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 605 | sage: p*p |
| 606 | [1, 4, 10, 20, 25, 24, 16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 607 | sage: one/p |
| 608 | [1, -2, 1, 0, 5, -14, 13, -4, 25, -90, 121, -72, 141, -550, 965, -844, 993, ...] |
| 609 | sage: p.rcp()/p*p*p |
| 610 | [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 611 | sage: p**2 |
| 612 | [1, 4, 10, 20, 25, 24, 16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 613 | sage: #composition only works for coefficient 0 being 0 in the second operand |
| 614 | sage: dexp = (expps - one).reclass() |
| 615 | sage: expps(dexp) |
| 616 | [1, 1, 1, 5/6, 5/8, 13/30, 203/720, 877/5040, 23/224, 1007/17280, ...] |
| 617 | |
| 618 | sage: #we come into interesting regions ... |
| 619 | sage: dexp(dexp) |
| 620 | [0, 1, 1, 5/6, 5/8, 13/30, 203/720, 877/5040, 23/224, 1007/17280, ...] |
| 621 | sage: dexp.nit(2) |
| 622 | [0, 1, 1, 5/6, 5/8, 13/30, 203/720, 877/5040, 23/224, 1007/17280, ...] |
| 623 | sage: dexp.it(-1) |
| 624 | [0, 1, -1/2, 1/3, -1/4, 1/5, -1/6, 1/7, -1/8, 1/9, -1/10, 1/11, -1/12, ...] |
| 625 | sage: dexp.it(-1)(dexp) |
| 626 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 627 | |
| 628 | sage: hdexp = dexp.it(1/2) |
| 629 | sage: hdexp |
| 630 | [0, 1, 1/4, 1/48, 0, 1/3840, -7/92160, 1/645120, 53/3440640, -281/30965760, ...] |
| 631 | sage: #verifying that shoudl be Zero |
| 632 | sage: hdexp.it(2) - dexp |
| 633 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 634 | |
| 635 | sage: #symbolic (parabolic) iteration |
| 636 | sage: dexp.it(x)[5].coefficients() |
| 637 | [[-1/180, 1], [1/24, 2], [-13/144, 3], [1/16, 4]] |
| 638 | sage: q = dexp.it(1/x).it(x) |
| 639 | |
| 640 | sage: expand(q[3]) |
| 641 | 1/6 |
| 642 | sage: dexp[3] |
| 643 | 1/6 |
| 644 | |
| 645 | sage: #you can initialize power series with arbitrary functions on natural numbers |
| 646 | sage: #for example the power series of sqrt(2)^x can be given as |
| 647 | sage: bsrt = FormalPowerSeriesRing(SR)(sqrt(2)^x,x) |
| 648 | |
| 649 | sage: #making the 0-th coefficient 0 to get the decremented exponential |
| 650 | sage: dbsrt = bsrt.set_item(0,0) |
| 651 | |
| 652 | sage: #and now starting hyperbolic iteration |
| 653 | sage: dbsrt2 = dbsrt.it(x).it(1/x) |
| 654 | sage: #Sage is not able to simplify |
| 655 | sage: #simplify(dbsrt2[3]) |
| 656 | sage: #but numerically we can verify equality |
| 657 | sage: RR(dbsrt2[3](x=0.73)-dbsrt[3]) < 10**(-18) |
| 658 | True |
| 659 | """ |
| 660 | |
| 661 | def __init__(self,parent,f=None,min_index=None,base_ring=None): |
| 662 | """ |
| 663 | Returns the formal powerseries. |
| 664 | |
| 665 | This method should only be called from FormalPowerSeriesRing. |
| 666 | See the description there how to obtain a FormalPowerSeries. |
| 667 | |
| 668 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing,FormalPowerSeries |
| 669 | sage: FormalPowerSeries(FormalPowerSeriesRing(QQ)) |
| 670 | Undefined |
| 671 | """ |
| 672 | if parent == None: |
| 673 | if base_ring==None: |
| 674 | parent = FormalPowerSeriesRing(QQ) |
| 675 | else: |
| 676 | parent = FormalPowerSeriesRing(base_ring) |
| 677 | |
| 678 | RingElement.__init__(self, parent) |
| 679 | |
| 680 | if not self.__class__.__dict__.has_key('coeffs'): |
| 681 | self.coeffs = f |
| 682 | |
| 683 | if min_index == None: |
| 684 | min_index = 0 |
| 685 | self.min_index = min_index #the minimal non-zero index |
| 686 | self._memo = {} |
| 687 | self._powMemo = {} |
| 688 | self._itMemo = {} |
| 689 | |
| 690 | self.K = self.parent().K |
| 691 | self.K1 = self.parent().K1 |
| 692 | |
| 693 | self.min_index = min_index |
| 694 | |
| 695 | if self.min_index > 0: |
| 696 | self._subclass(FormalPowerSeries0) |
| 697 | #if not f==None: |
| 698 | # self.reclass() |
| 699 | |
| 700 | # def new(self,f=None,min_index=0,**kwargs): |
| 701 | # return type(self)(self.parent(),f,**kwargs) |
| 702 | |
| 703 | def _subclass(self,T): |
| 704 | """ |
| 705 | Makes the methods in T available to self. |
| 706 | |
| 707 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 708 | sage: P = FormalPowerSeriesRing(QQ) |
| 709 | sage: #class B: |
| 710 | sage: # def f(self): return 'B' |
| 711 | sage: #b = P.Exp._subclass(B) |
| 712 | sage: #b.f() |
| 713 | """ |
| 714 | if isinstance(self,T): |
| 715 | return self |
| 716 | |
| 717 | # if issubclass(T,self.__class__): |
| 718 | # self.__class__ = T |
| 719 | # return self |
| 720 | # |
| 721 | # bs = () |
| 722 | # for C in self.__class__.__bases__: |
| 723 | # if issubclass(T,C): |
| 724 | # assert not T == C |
| 725 | # bs += (T,) |
| 726 | # else: |
| 727 | # bs += (C,) |
| 728 | # self.__class__.__bases__ = bs |
| 729 | |
| 730 | class Sub(T,self.__class__): pass |
| 731 | |
| 732 | self.__class__ = Sub |
| 733 | |
| 734 | return self |
| 735 | |
| 736 | def new(self,f=None,min_index=None,**kwargs): |
| 737 | """ |
| 738 | Returns a FormalPowerSeries from the same parent and class. |
| 739 | |
| 740 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 741 | sage: p = FormalPowerSeriesRing(QQ)([1,2]) |
| 742 | sage: p.new(lambda n: n).parent() == p.parent() |
| 743 | True |
| 744 | sage: p = FormalPowerSeriesRing(QQ)([0,2]) |
| 745 | sage: type(p.new(lambda n: n)) == type(p) |
| 746 | True |
| 747 | sage: p = FormalPowerSeriesRing(QQ)([0,1]) |
| 748 | sage: type(p.new()) == type(p) |
| 749 | True |
| 750 | """ |
| 751 | res = FormalPowerSeries(self.parent(),f,min_index,**kwargs) |
| 752 | if min_index == None: |
| 753 | res.__class__ = self.__class__ |
| 754 | if isinstance(self,FormalPowerSeries0): |
| 755 | res.min_index = 1 |
| 756 | return res |
| 757 | |
| 758 | # def __nonzero__(self): |
| 759 | # """ |
| 760 | # Returns always None, |
| 761 | # as it is not decidable whether a powerseries is 0. |
| 762 | |
| 763 | # sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 764 | # sage: FormalPowerSeriesRing(QQ)(0).is_zero() == None |
| 765 | # """ |
| 766 | # return None |
| 767 | |
| 768 | # def is_one(self): |
| 769 | # """ |
| 770 | # Returns always None, |
| 771 | # as it is not decidable whether a powerseries is 1. |
| 772 | # """ |
| 773 | # return None |
| 774 | |
| 775 | def define(self,p): |
| 776 | """ |
| 777 | Defines the power series self by another powerseries p |
| 778 | which is allowed to refer to self. |
| 779 | |
| 780 | For example one can define exp by taking the integral of its |
| 781 | derivative which is again exp: |
| 782 | |
| 783 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 784 | sage: P = FormalPowerSeriesRing(QQ) |
| 785 | sage: f = P.by_undefined() |
| 786 | sage: f.define(integral(f,1)) |
| 787 | sage: f |
| 788 | [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...] |
| 789 | |
| 790 | Or generally one can define g = exp o f by taking the integral of |
| 791 | g * f'. For example for f(x)=x**2, [0,0,1]: |
| 792 | sage: g = P() |
| 793 | sage: f = P([0,0,1]) |
| 794 | sage: fd = f.diff() |
| 795 | sage: g.define(integral(g*fd,1)) |
| 796 | sage: g - (f | P.Exp) |
| 797 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 798 | """ |
| 799 | self.coeffs = p.coeffs |
| 800 | |
| 801 | def __getitem__(self,n): |
| 802 | """ |
| 803 | Returns the n-th coefficient of this powerseries: self[n]. |
| 804 | This is the coefficient of x^n. |
| 805 | |
| 806 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 807 | sage: P = FormalPowerSeriesRing(QQ) |
| 808 | sage: P([1,2,3])[1] == 2 |
| 809 | True |
| 810 | """ |
| 811 | |
| 812 | if not self._memo.has_key(n): |
| 813 | #self._memo[n] = simplify(expand(self.coeffs(n))) |
| 814 | self._memo[n] = self.coeffs(n) |
| 815 | return self._memo[n] |
| 816 | |
| 817 | def __getslice__(self,i,j): # [i:j] |
| 818 | """ |
| 819 | Returns the list of coefficients with indices in range(i,j). |
| 820 | |
| 821 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 822 | sage: P = FormalPowerSeriesRing(QQ) |
| 823 | sage: P(lambda n: n)[1:3] |
| 824 | [1, 2] |
| 825 | """ |
| 826 | return [self[k] for k in range(i,j)] |
| 827 | |
| 828 | def set_item(a, index, value): |
| 829 | """ |
| 830 | Returns the powerseries that has a[index] replaced by value. |
| 831 | |
| 832 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 833 | sage: P = FormalPowerSeriesRing(QQ) |
| 834 | sage: P([1,2,3]).set_item(1,42) |
| 835 | [1, 42, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 836 | sage: P([1,2]).set_item(0,0).min_index |
| 837 | 1 |
| 838 | """ |
| 839 | min_index = a.min_index |
| 840 | if min_index == index and value == 0: |
| 841 | min_index += 1 |
| 842 | return a.new(lambda n: value if n == index else a[n],min_index) |
| 843 | |
| 844 | def __setitem__(a,index,value): |
| 845 | """ |
| 846 | Replaces a[index] by value, returns None. |
| 847 | |
| 848 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 849 | sage: P = FormalPowerSeriesRing(QQ) |
| 850 | sage: p = P([1,2,3]) |
| 851 | sage: p[1] = 42 |
| 852 | sage: p |
| 853 | [1, 42, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 854 | sage: p = P([1,2]) |
| 855 | sage: p[0] = 0 |
| 856 | sage: p.min_index |
| 857 | 1 |
| 858 | """ |
| 859 | min_index = a.min_index |
| 860 | if min_index == index and value == 0: |
| 861 | min_index += 1 |
| 862 | |
| 863 | a.min_index = min_index |
| 864 | f = a.coeffs |
| 865 | a.coeffs = lambda n: value if n == index else f(n) |
| 866 | |
| 867 | def set_slice(a,i,j,seq): |
| 868 | """ |
| 869 | Returns the powerseries that has a[i:j] replaced by seq, returns None. |
| 870 | |
| 871 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 872 | sage: P = FormalPowerSeriesRing(QQ) |
| 873 | sage: P(lambda n: n).set_slice(5,10,[42,43,44,45,46]) |
| 874 | [0, 1, 2, 3, 4, 42, 43, 44, 45, 46, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, ...] |
| 875 | sage: P([1,2]).set_slice(0,1,[0]).min_index |
| 876 | 1 |
| 877 | """ |
| 878 | |
| 879 | min_index = a.min_index |
| 880 | min_s=j-i |
| 881 | for k in range(0,j-i): |
| 882 | if not seq[k] == 0: |
| 883 | min_s = k |
| 884 | break |
| 885 | |
| 886 | if i <= min_index and min_index <= i+min_s: |
| 887 | min_index = i+min_s |
| 888 | else: |
| 889 | min_index = min(min_index,i+min_s) |
| 890 | return a.new(lambda n: seq[n-i] if i<=n and n<j else a[n],min_index) |
| 891 | |
| 892 | def __setslice__(a, i, j, seq): |
| 893 | """ |
| 894 | Replaces a[i:j] by seq. |
| 895 | |
| 896 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 897 | sage: P = FormalPowerSeriesRing(QQ) |
| 898 | sage: p = P(lambda n: n) |
| 899 | sage: p[5:10] = [42,43,44,45,46] |
| 900 | sage: p |
| 901 | [0, 1, 2, 3, 4, 42, 43, 44, 45, 46, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, ...] |
| 902 | sage: p = P([1,2]) |
| 903 | sage: p[0:1] = [0] |
| 904 | sage: p.min_index |
| 905 | 1 |
| 906 | """ |
| 907 | |
| 908 | min_index = a.min_index |
| 909 | min_s=j-i |
| 910 | for k in range(0,j-i): |
| 911 | if not seq[k] == 0: |
| 912 | min_s = k |
| 913 | break |
| 914 | |
| 915 | if i <= min_index and min_index <= i+min_s: |
| 916 | min_index = i+min_s |
| 917 | else: |
| 918 | min_index = min(min_index,i+min_s) |
| 919 | |
| 920 | a.min_index = min_index |
| 921 | f = a.coeffs |
| 922 | a.coeffs = lambda n: seq[n-i] if i<=n and n<j else f(n) |
| 923 | |
| 924 | def extinct_before(a,min_index): |
| 925 | """ |
| 926 | Returns the powerseries with elements before min_index replaced by 0. |
| 927 | |
| 928 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 929 | sage: P = FormalPowerSeriesRing(QQ) |
| 930 | sage: P(lambda n: n).extinct_before(5) |
| 931 | [0, 0, 0, 0, 0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...] |
| 932 | """ |
| 933 | return ExtinctBefore(a,min_index) |
| 934 | |
| 935 | def reclass(p): |
| 936 | """ |
| 937 | Recalculates the class of this object which possibly changes |
| 938 | to a subclass having more (and possibly different) operations |
| 939 | available. Returns p. |
| 940 | |
| 941 | Reclass queries p[0] and p[1], so for the sake of a lazy `define' |
| 942 | this is not automatically done on creation. |
| 943 | |
| 944 | Due to implementation effort a reclassed object can not be |
| 945 | pickled/serialized with dumps. |
| 946 | |
| 947 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing, FormalPowerSeries, FormalPowerSeries0, FormalPowerSeries01 |
| 948 | sage: P = FormalPowerSeriesRing(QQ) |
| 949 | sage: isinstance(P([0,1]).reclass(),FormalPowerSeries01) |
| 950 | True |
| 951 | sage: isinstance(P([0,2]).reclass(),FormalPowerSeries0) |
| 952 | True |
| 953 | sage: isinstance(P([1,1]).reclass(),FormalPowerSeries) |
| 954 | True |
| 955 | """ |
| 956 | |
| 957 | # if not decidable0(p.K): |
| 958 | # if p.min_index > 0 and not isinstance(p,FormalPowerSeries0): |
| 959 | # p._subclass(FormalPowerSeries0) |
| 960 | # return p |
| 961 | |
| 962 | min_index = max(2,p.min_index) |
| 963 | for n in range(p.min_index,2): |
| 964 | if not p[n] == 0: |
| 965 | min_index = n |
| 966 | break |
| 967 | |
| 968 | p.min_index = min_index |
| 969 | |
| 970 | if p.min_index < 0: |
| 971 | return p |
| 972 | |
| 973 | if min_index > 0: |
| 974 | if p[1] == 1: |
| 975 | return p._subclass(FormalPowerSeries01) |
| 976 | return p._subclass(FormalPowerSeries0) |
| 977 | return p |
| 978 | |
| 979 | def mul_fact(a): |
| 980 | """ |
| 981 | The sequence a[n]*n! |
| 982 | |
| 983 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 984 | sage: P = FormalPowerSeriesRing(QQ) |
| 985 | sage: P.Exp.mul_fact() |
| 986 | [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...] |
| 987 | """ |
| 988 | return MulFact(a) |
| 989 | |
| 990 | def div_fact(a): |
| 991 | """ |
| 992 | Returns the sequence a[n]/n!. |
| 993 | |
| 994 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 995 | sage: P = FormalPowerSeriesRing(QQ) |
| 996 | sage: P(lambda n: 1).div_fact() - P.Exp |
| 997 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 998 | """ |
| 999 | return DivFact(a) |
| 1000 | |
| 1001 | def inc(a): |
| 1002 | """ |
| 1003 | Increment: a + One. |
| 1004 | |
| 1005 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1006 | sage: P = FormalPowerSeriesRing(QQ) |
| 1007 | sage: P.Zero.inc() |
| 1008 | [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1009 | sage: P.Zero + P.One |
| 1010 | [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1011 | """ |
| 1012 | return IncMethod(a) |
| 1013 | |
| 1014 | def dec(a): |
| 1015 | """ |
| 1016 | Decrement: a - One. |
| 1017 | |
| 1018 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1019 | sage: P = FormalPowerSeriesRing(QQ) |
| 1020 | |
| 1021 | sage: P.Zero.dec() |
| 1022 | [-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1023 | """ |
| 1024 | return DecMethod(a) |
| 1025 | |
| 1026 | def rmul(a,s): |
| 1027 | """ |
| 1028 | Scalar multiplication from right with scalar s |
| 1029 | |
| 1030 | Alternative expression: a.rmul(s) == a * s |
| 1031 | |
| 1032 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1033 | sage: P = FormalPowerSeriesRing(QQ) |
| 1034 | sage: p = P([1,2,3]) |
| 1035 | sage: p.rmul(2) |
| 1036 | [2, 4, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1037 | sage: p * (2/3) |
| 1038 | [2/3, 4/3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1039 | """ |
| 1040 | return RMul(a,s) |
| 1041 | |
| 1042 | _rmul_ = rmul |
| 1043 | |
| 1044 | def lmul(a,s): |
| 1045 | """ |
| 1046 | Scalar multiplication from left with scalar s |
| 1047 | |
| 1048 | Alternative expression: a.lmul(s) == s * a |
| 1049 | |
| 1050 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1051 | sage: P = FormalPowerSeriesRing(QQ) |
| 1052 | sage: p = P([1,2,3]) |
| 1053 | sage: p.lmul(2) |
| 1054 | [2, 4, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1055 | sage: (2/3) * p |
| 1056 | [2/3, 4/3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1057 | """ |
| 1058 | return LMul(a,s) |
| 1059 | |
| 1060 | _lmul_ = lmul |
| 1061 | |
| 1062 | def add(a,b): |
| 1063 | """ |
| 1064 | Addition of two powerseries. |
| 1065 | |
| 1066 | Alternative expression: a.add(b) == a+b |
| 1067 | |
| 1068 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1069 | sage: P = FormalPowerSeriesRing(QQ) |
| 1070 | sage: P([1,2,3]) + P([4,5,6]) |
| 1071 | [5, 7, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1072 | """ |
| 1073 | |
| 1074 | return Add(a,b) |
| 1075 | |
| 1076 | _add_ = add |
| 1077 | |
| 1078 | def sub(a,b): |
| 1079 | """ |
| 1080 | Subtraction of two powerseries: a-b. |
| 1081 | |
| 1082 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1083 | sage: P = FormalPowerSeriesRing(QQ) |
| 1084 | sage: P(lambda n: n) - P(lambda n: n) |
| 1085 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1086 | sage: P([0,1]).sub(P([1,0])) |
| 1087 | [-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1088 | """ |
| 1089 | return Sub(a,b) |
| 1090 | |
| 1091 | _sub_ = sub |
| 1092 | |
| 1093 | def neg(a): |
| 1094 | """ |
| 1095 | Negation of the powerseries a. |
| 1096 | |
| 1097 | Alternative expression a.neg() == -a. |
| 1098 | |
| 1099 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1100 | sage: P = FormalPowerSeriesRing(QQ) |
| 1101 | sage: -P(lambda n: 1) |
| 1102 | [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...] |
| 1103 | """ |
| 1104 | return Neg(a) |
| 1105 | |
| 1106 | _neg_ = neg |
| 1107 | |
| 1108 | def mul(a,b): |
| 1109 | """ |
| 1110 | Multiplication of two powerseries. |
| 1111 | This a lazy algorithm: for initial 0 in a[k] and initial 0 in b[n-k] |
| 1112 | the corresponding b[n-k] or a[k] is not evaluated. |
| 1113 | |
| 1114 | Alternative expression: a.mul(b) == a*b |
| 1115 | |
| 1116 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1117 | sage: P = FormalPowerSeriesRing(QQ) |
| 1118 | sage: P(lambda n: 1) * P(lambda n:1) |
| 1119 | [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, ...] |
| 1120 | """ |
| 1121 | #multiplication of two powerseries |
| 1122 | |
| 1123 | return Mul(a,b) |
| 1124 | |
| 1125 | _mul_ = mul |
| 1126 | |
| 1127 | def div(c,b): |
| 1128 | """ |
| 1129 | Division. |
| 1130 | Alternative expression: c.div(b) == c/b |
| 1131 | |
| 1132 | Precondition: b != Zero |
| 1133 | It satisfies: a/b*b=a, a*b/b=a |
| 1134 | |
| 1135 | If c[0]==0 it returns a formal Laurant series, i.e. min_index < 0. |
| 1136 | |
| 1137 | This operation is not lazy in b, it retrieves values starting from |
| 1138 | i=min_index until it finds b[i]!=0 |
| 1139 | |
| 1140 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1141 | sage: P = FormalPowerSeriesRing(QQ) |
| 1142 | sage: P([1,2,3])/P([2,3,4]) |
| 1143 | [1/2, 1/4, 1/8, -11/16, 25/32, 13/64, -239/128, 613/256, 73/512, ...] |
| 1144 | sage: P.One/P(lambda n: (-1)**n) |
| 1145 | [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1146 | """ |
| 1147 | b.min_index = b.val() |
| 1148 | |
| 1149 | return Div(c,b) |
| 1150 | |
| 1151 | _div_ = div |
| 1152 | |
| 1153 | def rcp(a): |
| 1154 | """ |
| 1155 | Returns the reciprocal power series, i.e. One/a. |
| 1156 | a must not be Zero. |
| 1157 | If a[0] == 0 it returns a fromal Laurant series (min_index < 0). |
| 1158 | |
| 1159 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1160 | sage: P = FormalPowerSeriesRing(QQ) |
| 1161 | sage: P(lambda n:1).rcp() |
| 1162 | [1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1163 | """ |
| 1164 | return a.parent().One/a |
| 1165 | |
| 1166 | def npow_mult(a,n): |
| 1167 | """ |
| 1168 | Power with natural exponent n computed in the most simple way by |
| 1169 | multiplying the powerseries n times. |
| 1170 | This function is cached, it remembers the power for each n. |
| 1171 | |
| 1172 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1173 | sage: P = FormalPowerSeriesRing(QQ) |
| 1174 | sage: P([1,2,3]).npow(2)/P([1,2,3]) |
| 1175 | [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1176 | """ |
| 1177 | |
| 1178 | _assert_nat(n) |
| 1179 | |
| 1180 | if n==0: |
| 1181 | return a.parent().One |
| 1182 | if n==1: |
| 1183 | return a |
| 1184 | if not a._powMemo.has_key(n): |
| 1185 | n1 = int(n) / 2 |
| 1186 | n2 = int(n) / 2 + n % 2 |
| 1187 | a._powMemo[n] = a.npow_mult(n1) * a.npow_mult(n2) |
| 1188 | return a._powMemo[n] |
| 1189 | |
| 1190 | def _s(f,k,m,n): |
| 1191 | """ |
| 1192 | Returns the sum over |
| 1193 | m_1+2*m_2+...+k*m_k = n, m_0+m_1+...+m_k = m |
| 1194 | of |
| 1195 | (m over m_1,m_2,...,m_k)* f[0]^{m_0} ... f[k]^{m_k} |
| 1196 | |
| 1197 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1198 | sage: P = FormalPowerSeriesRing(QQ) |
| 1199 | sage: P.Exp._s(3,5,2) |
| 1200 | 25/2 |
| 1201 | """ |
| 1202 | # print k,m,n |
| 1203 | if f._powMemo.has_key((k,m,n)): |
| 1204 | return f._powMemo[(k,m,n)] |
| 1205 | |
| 1206 | if m == 0: |
| 1207 | if n == 0: |
| 1208 | return 1 |
| 1209 | return 0 |
| 1210 | if k < f.min_index: |
| 1211 | return 0 |
| 1212 | # if n < f.min_index*m: #case only for speed optimizing |
| 1213 | # return 0 |
| 1214 | # if k == f.min_index: #case only for speed optimizing |
| 1215 | # if n == k*m: |
| 1216 | # return f[k]**m |
| 1217 | # return 0 |
| 1218 | # print 'p' |
| 1219 | res = 0 |
| 1220 | mk = 0 |
| 1221 | while min(f.min_index*m,0)+k*mk <= n and mk <= m: |
| 1222 | v = f._s(k-1,m-mk,n-k*mk) |
| 1223 | if not v == 0: #be lazy in evaluating f[k] |
| 1224 | if not mk == 0: #be lazy in evaluating f[k] |
| 1225 | v *= f[k]**mk |
| 1226 | res += v * binomial(m,mk) |
| 1227 | mk+=1 |
| 1228 | |
| 1229 | f._powMemo[(k,m,n)] = res |
| 1230 | return res |
| 1231 | |
| 1232 | def npow_combin(f,m): |
| 1233 | """ |
| 1234 | Power with natural exponent m. |
| 1235 | Computed via the cominatorial way similar to Faa di Bruno's formula. |
| 1236 | |
| 1237 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1238 | sage: P = FormalPowerSeriesRing(QQ) |
| 1239 | sage: P([1,2,3]).npow_combin(2)/P([1,2,3]) |
| 1240 | [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1241 | """ |
| 1242 | |
| 1243 | _assert_nat(m) |
| 1244 | return Npow(f,m) |
| 1245 | |
| 1246 | npow = npow_mult |
| 1247 | |
| 1248 | def nipow(a,t): |
| 1249 | """ |
| 1250 | Non-integer power of a (works also for integers but less efficiently). |
| 1251 | a[0] must be nonzero and a[0]**t must be defined |
| 1252 | as well multiplication of t with all coefficients. |
| 1253 | |
| 1254 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1255 | sage: P = FormalPowerSeriesRing(QQ) |
| 1256 | sage: P.Exp.nipow(1/2) |
| 1257 | [1, 1/2, 1/8, 1/48, 1/384, 1/3840, 1/46080, 1/645120, 1/10321920, ...] |
| 1258 | |
| 1259 | sage: PR = FormalPowerSeriesRing(RR) |
| 1260 | sage: PR.Exp.nipow(0.5).n(20) |
| 1261 | [1.0000, 0.50000, 0.12500, 0.020833, 0.0026042, 0.00026042, 0.000021701, ...] |
| 1262 | """ |
| 1263 | return Nipow(a,t) |
| 1264 | |
| 1265 | def pow(a,t): |
| 1266 | """ |
| 1267 | The t-th (possibly non-integer) power of a. |
| 1268 | a[0]**t has to be defined. |
| 1269 | |
| 1270 | Alternative expression: a.pow(t) == a ** t |
| 1271 | |
| 1272 | It satisfies: |
| 1273 | a.nipow(0) = One |
| 1274 | a.nipow(1) = a |
| 1275 | a.nipow(s+t) == a.nipow(s)*a.nipow(t) |
| 1276 | |
| 1277 | See also: npow, nipow |
| 1278 | |
| 1279 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1280 | sage: P = FormalPowerSeriesRing(QQ) |
| 1281 | sage: P(lambda n: 1).pow(-1) |
| 1282 | [1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1283 | sage: P = FormalPowerSeriesRing(RealField(16)) |
| 1284 | sage: P([1,2,3])**(-0.37) |
| 1285 | [1.000, -0.7400, -0.09619, 1.440, -2.228, 0.6642, 4.092, -9.079, 6.390, ...] |
| 1286 | """ |
| 1287 | |
| 1288 | if isinstance(t,FormalPowerSeries): |
| 1289 | P = a.parent() |
| 1290 | return ( a.log() * t ) | P.Exp |
| 1291 | if isinstance(t,Integer) or isinstance(t,int): |
| 1292 | if t >= 0: |
| 1293 | return a.npow(t) |
| 1294 | return a.rcp().npow(-t) |
| 1295 | return a.nipow(t) |
| 1296 | |
| 1297 | __pow__ = pow |
| 1298 | |
| 1299 | def sqrt(a): |
| 1300 | """ |
| 1301 | Square root of a. a.sqrt() * a.sqrt() == a |
| 1302 | |
| 1303 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1304 | sage: P = FormalPowerSeriesRing(QQ) |
| 1305 | sage: P([1,2,1]).sqrt() |
| 1306 | [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1307 | sage: (P.One - P.Sin ** 2).sqrt() - P.Cos |
| 1308 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1309 | """ |
| 1310 | return a.rt(2) |
| 1311 | |
| 1312 | def rt(a,n): |
| 1313 | """ |
| 1314 | n-th root of a. a.rt(n)**n == a |
| 1315 | |
| 1316 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1317 | sage: P = FormalPowerSeriesRing(QQ) |
| 1318 | |
| 1319 | sage: P([1,3,3,1]).rt(3) |
| 1320 | [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1321 | """ |
| 1322 | return a.nipow(1/Integer(n)) |
| 1323 | |
| 1324 | def compose(b,a): |
| 1325 | """ |
| 1326 | Composition (left after right), in mathematical notation b o a. |
| 1327 | Alternative expressions: b.compose(a) == b(a) == a | b |
| 1328 | |
| 1329 | Precondition: |
| 1330 | a[0] == 0. |
| 1331 | Satisfies: |
| 1332 | a(b).polynomial(m*n,x) == a.polynomial(m,b.polynomial(n,x)) |
| 1333 | |
| 1334 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1335 | sage: P = FormalPowerSeriesRing(QQ) |
| 1336 | sage: P([1,2,3]).compose(P([0,1,2])) |
| 1337 | [1, 2, 7, 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1338 | sage: P([1,2,3])(P([0,1,2])) |
| 1339 | [1, 2, 7, 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1340 | """ |
| 1341 | a._assertp0() |
| 1342 | |
| 1343 | return Compose(b,a) |
| 1344 | |
| 1345 | __call__ = compose |
| 1346 | |
| 1347 | def log(a): |
| 1348 | """ |
| 1349 | Logarithm of powerseries a with a[0]==1. |
| 1350 | |
| 1351 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1352 | sage: P = FormalPowerSeriesRing(QQ) |
| 1353 | sage: P.Exp.log() |
| 1354 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1355 | """ |
| 1356 | |
| 1357 | P = a.parent() |
| 1358 | |
| 1359 | dec_a = a.set_item(0,0) |
| 1360 | |
| 1361 | # if decidable0(a.K): |
| 1362 | # assert a[0] == 1 |
| 1363 | return dec_a | P.Log_inc |
| 1364 | |
| 1365 | # def __xor__(a,t): # ^ |
| 1366 | # #Not recognized as it seems to be mapped to ** in sage |
| 1367 | # return NotImplemented |
| 1368 | |
| 1369 | def bell_polynomials(a,k): |
| 1370 | """ |
| 1371 | The sequence of Bell polynomials (partial/of the second kind) |
| 1372 | [B_{0,k}(a[1],a[2],...),B_{1,k}(a[1],a[2],...),...] |
| 1373 | |
| 1374 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1375 | sage: P = PolynomialRing(QQ,['c1','c2','c3','c4','c5']) |
| 1376 | sage: c1 = P('c1'); c2= P('c2'); c3 = P('c3'); c4 = P('c4'); c5 = P('c5') |
| 1377 | sage: c = FormalPowerSeriesRing(QQ)([0,c1,c2,c3,c4,c5]) |
| 1378 | sage: c.bell_polynomials(2)[6] == 6*c5*c1 + 15*c4*c2 + 10*c3**2 |
| 1379 | True |
| 1380 | """ |
| 1381 | return (a.div_fact()**k).mul_fact().rmul(Integer(1)/factorial(k)) |
| 1382 | |
| 1383 | def bell_polynomial(a,n,k): |
| 1384 | """ |
| 1385 | The Bell polynomial (partial/of the second kind). |
| 1386 | |
| 1387 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1388 | sage: P = PolynomialRing(QQ,['c1','c2','c3','c4','c5']) |
| 1389 | sage: c1 = P('c1'); c2= P('c2'); c3 = P('c3'); c4 = P('c4'); c5 = P('c5') |
| 1390 | sage: c = FormalPowerSeriesRing(QQ)([0,c1,c2,c3,c4,c5]) |
| 1391 | sage: c.bell_polynomial(6,3) == 15*c4*c1**2 + 60*c3*c2*c1 + 15*c2**3 |
| 1392 | True |
| 1393 | sage: c.bell_polynomial(6,2) == 6*c5*c1 + 15*c4*c2 + 10*c3**2 |
| 1394 | True |
| 1395 | """ |
| 1396 | return a.bell_polynomials(k)[n] |
| 1397 | |
| 1398 | def bell_complete(a,n): |
| 1399 | """ |
| 1400 | The complete Bell polynomial. |
| 1401 | |
| 1402 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1403 | sage: P = PolynomialRing(QQ,['c1','c2','c3','c4','c5']) |
| 1404 | sage: c1 = P('c1');c2 = P('c2');c3 = P('c3');c4 = P('c4');c5 = P('c5') |
| 1405 | sage: PS = FormalPowerSeriesRing(P) |
| 1406 | sage: c = PS([0,c1,c2,c3,c4,c5]) |
| 1407 | sage: (PS.Exp(c.div_fact()) - c.bell_complete(5).div_fact())[1:6] |
| 1408 | [0, 0, 0, 0, 0] |
| 1409 | """ |
| 1410 | if n <= 0: |
| 1411 | return a.parent().Zero |
| 1412 | |
| 1413 | res = a.bell_polynomials(1) |
| 1414 | for k in range(2,n+1): |
| 1415 | res += a.bell_polynomials(k) |
| 1416 | return res |
| 1417 | |
| 1418 | |
| 1419 | def genfunc(a,n): |
| 1420 | """ |
| 1421 | Returns the generating function of this powerseries up to term |
| 1422 | a[n-1]*x**(n-1). |
| 1423 | |
| 1424 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1425 | sage: P = FormalPowerSeriesRing(QQ) |
| 1426 | sage: f = P.by_list([1,2,3],-2).polynomial(5) |
| 1427 | sage: g = P.by_list([1,2,3],-2).genfunc(5) |
| 1428 | sage: f(3.7)==g(3.7) |
| 1429 | True |
| 1430 | """ |
| 1431 | m = a.min_index |
| 1432 | return lambda x: sum([a[k]*x**k for k in range(m,n)],a.K(0)) |
| 1433 | |
| 1434 | def polynomial(a,n,x=var('x')): |
| 1435 | """ |
| 1436 | Returns the associated polynomial for the first n coefficients. |
| 1437 | f_0 + f_1*x + f_2*x^2 + ... + f_{n-1}*x^{n-1} |
| 1438 | |
| 1439 | In case of a Laurant series with e.g. min_index=-2: |
| 1440 | f_{-2} x^(-2) + f_{-1} x^(-1) + ... + f_{n-1}*x^{n-1} |
| 1441 | |
| 1442 | You can adjust the variable by setting x. |
| 1443 | |
| 1444 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1445 | sage: P = FormalPowerSeriesRing(QQ) |
| 1446 | sage: P.by_list([0,1,2]).polynomial(5).padded_list() |
| 1447 | [0, 1, 2] |
| 1448 | """ |
| 1449 | |
| 1450 | # return PolynomialRing(a.K,x)(sum([a[k]*x**k for k in range(n)],a.K(0))) |
| 1451 | P = PolynomialRing(a.K,x) |
| 1452 | m = a.min_index |
| 1453 | if m >= 0: |
| 1454 | return P(a[:n]) |
| 1455 | return P(a[m:n])/P(x**(-m)) |
| 1456 | |
| 1457 | # def subs(a,*args,**kwargs): |
| 1458 | # def f(n): |
| 1459 | # if n < a.min_index: |
| 1460 | # return 0 |
| 1461 | # if a[n] == 0: |
| 1462 | # return 0 |
| 1463 | # if a[n] == 1: |
| 1464 | # return 1 |
| 1465 | # return a[n].subs(*args,**kwargs) |
| 1466 | # return FormalPowerSeries(f,a.min_index) |
| 1467 | |
| 1468 | def _assertp0(a): |
| 1469 | """ |
| 1470 | Asserts a[0]==0. |
| 1471 | |
| 1472 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1473 | sage: P = FormalPowerSeriesRing(QQ) |
| 1474 | sage: P([0,2])._assertp0() |
| 1475 | """ |
| 1476 | assert a.min_index > 0, "min_index must be > 0, but is " + repr(a.min_index) + ". Use reclass() if necessary." |
| 1477 | |
| 1478 | def _repr_(a): |
| 1479 | """ |
| 1480 | Returns a string representation of a, as list of the first |
| 1481 | coefficients. |
| 1482 | If it is a formal Laurant series |
| 1483 | the coefficients before index 0 are seperated by ";" |
| 1484 | |
| 1485 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1486 | sage: FormalPowerSeriesRing(QQ).by_list([1,2,3],-2) |
| 1487 | [1, 2; 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1488 | sage: FormalPowerSeriesRing(QQ).by_list([1,2,3],2) |
| 1489 | [0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1490 | sage: FormalPowerSeriesRing(QQ).by_list([1,2,3],-2)._repr_() |
| 1491 | '[1, 2; 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]' |
| 1492 | """ |
| 1493 | # res = "" |
| 1494 | # for n in range(80): |
| 1495 | # coeff = a(n) |
| 1496 | # s = "" |
| 1497 | # if coeff != 0: |
| 1498 | # if coeff != 1: |
| 1499 | # s += repr(a(n)) + "*" |
| 1500 | # if n != 0: |
| 1501 | # s += "x" |
| 1502 | # if n != 1: |
| 1503 | # s += "^" + repr(n) |
| 1504 | # else: |
| 1505 | # s += repr(a(n)) |
| 1506 | # s += " + " |
| 1507 | # if len(res)+len(s) > 75: break |
| 1508 | # else: res += s |
| 1509 | |
| 1510 | # res += "O(x^" + repr(n) + ")" |
| 1511 | |
| 1512 | if a.coeffs == None: |
| 1513 | return "Undefined" |
| 1514 | |
| 1515 | res = "[" |
| 1516 | if a.min_index < 0: |
| 1517 | for k in range(a.min_index,0): |
| 1518 | res += repr(a[k]) |
| 1519 | if k==-1: |
| 1520 | res += "; " |
| 1521 | else: |
| 1522 | res += ", " |
| 1523 | for n in range(80): |
| 1524 | coeff = a[n] |
| 1525 | s = repr(a[n]) + ", " |
| 1526 | if len(res)+len(s) > 76: break |
| 1527 | else: res += s |
| 1528 | |
| 1529 | res += "...]"; |
| 1530 | #res = repr([ a(m) for m in range(10)]) |
| 1531 | return res |
| 1532 | |
| 1533 | # def complex_contour(a,N,fp=0): |
| 1534 | # """ |
| 1535 | # Returns a contour plot of this powerseries. |
| 1536 | # Experimental yet. |
| 1537 | # """ |
| 1538 | # r = abs(a[N])**(-1/Integer(N)) |
| 1539 | # l = r/sqrt(2.0) |
| 1540 | # f = a.polynomial(N) |
| 1541 | # x0=real(fp) |
| 1542 | # y0=imag(fp) |
| 1543 | # return contour_plot(lambda x,y: real(f(CC(x+i*y-fp))),(x0-l,x0+l),(y0-l,y0+l),fill=false) + contour_plot(lambda x,y: imag(f(CC(x+i*y-fp))),(x0-l,x0+l),(y0-l,y0+l),fill=false) |
| 1544 | |
| 1545 | def n(a,*args,**kwargs): |
| 1546 | """ |
| 1547 | Applies n to the coefficients |
| 1548 | |
| 1549 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1550 | sage: P = FormalPowerSeriesRing(QQ) |
| 1551 | sage: P.Exp.n(digits=3) |
| 1552 | [1.00, 1.00, 0.500, 0.167, 0.0417, 0.00833, 0.00139, 0.000198, 0.0000248, ...] |
| 1553 | """ |
| 1554 | return N(a,*args,**kwargs) |
| 1555 | |
| 1556 | def val(a): |
| 1557 | """ |
| 1558 | Returns the first index i such that a[i] != 0 |
| 1559 | Does not terminate if a == 0 |
| 1560 | |
| 1561 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1562 | sage: P = FormalPowerSeriesRing(QQ) |
| 1563 | sage: P([0,0,42]).val() |
| 1564 | 2 |
| 1565 | sage: P.by_list([1],-42).val() |
| 1566 | -42 |
| 1567 | """ |
| 1568 | n = a.min_index |
| 1569 | while a[n] == 0: |
| 1570 | n += 1 |
| 1571 | return n |
| 1572 | |
| 1573 | def __lshift__(a,m=1): |
| 1574 | """ |
| 1575 | Returns the powerseries with coefficients shiftet to the left by m. |
| 1576 | The coefficients a[0],...,a[m-1] get lost. |
| 1577 | |
| 1578 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1579 | sage: P = FormalPowerSeriesRing(QQ) |
| 1580 | sage: (P.Exp.mul_fact() << 1).div_fact() - P.Exp |
| 1581 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1582 | """ |
| 1583 | return Lshift(a,m) |
| 1584 | |
| 1585 | def __rshift__(a,m=1): |
| 1586 | """ |
| 1587 | Returns the powerseries with coefficients shifted to the right by m. |
| 1588 | The resulting coefficients b[0],...,b[m-1] are zero. |
| 1589 | |
| 1590 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1591 | sage: P = FormalPowerSeriesRing(QQ) |
| 1592 | sage: (P.Exp.mul_fact() >> 1).div_fact() - P.Exp |
| 1593 | [-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1594 | """ |
| 1595 | return Rshift(a,m) |
| 1596 | |
| 1597 | def diff(a,m=1): |
| 1598 | """ |
| 1599 | Differentiates the powerseries m times. |
| 1600 | |
| 1601 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1602 | sage: P = FormalPowerSeriesRing(QQ) |
| 1603 | sage: P.Exp.diff(3) - P.Exp |
| 1604 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1605 | sage: P.by_list([1],5).diff(5)[0] == factorial(5) |
| 1606 | True |
| 1607 | """ |
| 1608 | return Diff(a,m) |
| 1609 | |
| 1610 | # def integral(a,c=0,m=1): |
| 1611 | # """ |
| 1612 | # Integrates the powerseries m times with integration constant being 0 |
| 1613 | # """ |
| 1614 | # for i in range(-m,0): |
| 1615 | # if a[i] != 0: |
| 1616 | # if m == 1: |
| 1617 | # raise ValueError, "Coefficient -1 must be 0, but it is: " + repr(a[-1]) |
| 1618 | # else: |
| 1619 | # raise ValueError, "Coefficients from -"+repr(m)+" upto -1 must be 0, however a[" +repr(i)+"]=" + repr(a[i]) |
| 1620 | |
| 1621 | # def f(n): |
| 1622 | # if 0 <= n and n < m: |
| 1623 | # return c |
| 1624 | # return a[n-m]/prod(Integer(k) for k in range(n-m+1,n+1)) |
| 1625 | # return a.new(f,a.min_index+m) |
| 1626 | |
| 1627 | def integral(a,c=0): |
| 1628 | """ |
| 1629 | Integrates the powerseries with integration constant c |
| 1630 | |
| 1631 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1632 | sage: P = FormalPowerSeriesRing(QQ) |
| 1633 | sage: P.Exp.integral(1) - P.Exp |
| 1634 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1635 | sage: P.Exp.integral(0) + P.One - P.Exp |
| 1636 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1637 | """ |
| 1638 | |
| 1639 | return Integral(a,c) |
| 1640 | |
| 1641 | # def indefinite_sum(f,c=0): |
| 1642 | # def ids(n,m): |
| 1643 | # N = m+1 |
| 1644 | # if n > N: |
| 1645 | # return 0 |
| 1646 | # if n < 0: |
| 1647 | # return 0 |
| 1648 | # if n == 0: |
| 1649 | # return c |
| 1650 | # if n == N: |
| 1651 | # return 1/QQ(n) |
| 1652 | # return - sum([ f[k]*binomial(k,n) for k in range(n+1,N+1)])/QQ(n) |
| 1653 | # print ids(1,2), ids(2,2), ids(3,2) |
| 1654 | |
| 1655 | ### finite approximate operations |
| 1656 | |
| 1657 | def carleman_matrix(p,N,M=None): |
| 1658 | """ |
| 1659 | The carleman_matrix has as nth row the coefficients of p^n. |
| 1660 | It has N rows and N columns, except M specifies a different number of |
| 1661 | rows. |
| 1662 | |
| 1663 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1664 | sage: P = FormalPowerSeriesRing(QQ) |
| 1665 | sage: P([1,1]).carleman_matrix(4) == matrix([[1,0,0,0],[1,1,0,0],[1,2,1,0],[1,3,3,1]]) |
| 1666 | True |
| 1667 | """ |
| 1668 | if M == None: |
| 1669 | M=N |
| 1670 | return matrix([[p.npow(m)[n] for n in range(N) ] for m in range(M)]) |
| 1671 | |
| 1672 | def bell_matrix(a,N,M=None): |
| 1673 | """ |
| 1674 | The Bell matrix with N (or M) rows and N columns of this power series. |
| 1675 | The n-th column of the Bell matrix is the sequence of coefficients |
| 1676 | of the n-th power of the power series. |
| 1677 | The Bell matrix is the transpose of the Carleman matrix. |
| 1678 | |
| 1679 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1680 | sage: P = FormalPowerSeriesRing(QQ) |
| 1681 | sage: P([1,1]).bell_matrix(4) == matrix([[1,1,1,1],[0,1,2,3],[0,0,1,3],[0,0,0,1]]) |
| 1682 | True |
| 1683 | """ |
| 1684 | if M == None: |
| 1685 | M=N |
| 1686 | return matrix([[a.npow(n)[m] for n in range(N)] for m in range(M)]) |
| 1687 | |
| 1688 | class FormalPowerSeries0(FormalPowerSeries): |
| 1689 | """ |
| 1690 | The formal powerseries f with f[0]==0. |
| 1691 | """ |
| 1692 | # def __init__(self,f=None,min_index=1,**kwargs): |
| 1693 | # """ |
| 1694 | # Initializes this FormalPowerSeries0. |
| 1695 | # Should be called only from FormalPowerSeriesRing. |
| 1696 | |
| 1697 | # sage: None |
| 1698 | # """ |
| 1699 | # assert min_index >= 1 |
| 1700 | # super(FormalPowerSeries0,self).__init__(f,min_index,**kwargs) |
| 1701 | |
| 1702 | def __or__(a,b): |
| 1703 | """ |
| 1704 | Composition (right after left): a|b. |
| 1705 | For documentation see: `compose'. |
| 1706 | |
| 1707 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1708 | sage: P = FormalPowerSeriesRing(QQ) |
| 1709 | sage: P([0,1,2]) | P([1,2,3]) |
| 1710 | [1, 2, 7, 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1711 | """ |
| 1712 | return b(a) |
| 1713 | |
| 1714 | def nit(a,n): |
| 1715 | """ |
| 1716 | n times iteration (n times composition), n is a natural number. |
| 1717 | |
| 1718 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1719 | sage: P = FormalPowerSeriesRing(QQ) |
| 1720 | sage: P(1/(x+1)-1,x).nit(2) |
| 1721 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1722 | """ |
| 1723 | _assert_nat(n) |
| 1724 | |
| 1725 | if n == 0: |
| 1726 | return a.parent().Id |
| 1727 | if n == 1: |
| 1728 | return a |
| 1729 | if not a._itMemo.has_key(n): |
| 1730 | n1 = int(n) / 2 |
| 1731 | n2 = int(n) / 2 + n % 2 |
| 1732 | a._itMemo[n] = a.nit(n1).compose(a.nit(n2)) |
| 1733 | return a._itMemo[n] |
| 1734 | |
| 1735 | def it(a,t): |
| 1736 | """ |
| 1737 | The t-th (possibly non-integer) iteration. |
| 1738 | |
| 1739 | It satisfies: |
| 1740 | a.it(0) = Id |
| 1741 | a.it(1) = a |
| 1742 | a.it(s+t) == a.it(s)(a.it(t)) |
| 1743 | |
| 1744 | Alternative expression: a.it(t) == a & t |
| 1745 | |
| 1746 | See also: nit, regit (for restrictions depending t's kind) |
| 1747 | |
| 1748 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1749 | sage: P = FormalPowerSeriesRing(QQ) |
| 1750 | sage: P([0,2]).it(-2) |
| 1751 | [0, 1/4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1752 | sage: p = P([0,2]).regit_b(1/2) |
| 1753 | sage: (p[0],p[1],p[2]) |
| 1754 | (0, sqrt(2), 0) |
| 1755 | """ |
| 1756 | a._assertp0() |
| 1757 | |
| 1758 | if isinstance(t,Integer) or isinstance(t,int): |
| 1759 | if t >= 0: |
| 1760 | return a.nit(t) |
| 1761 | return a.inv().nit(-t) |
| 1762 | return a.regit(t) |
| 1763 | |
| 1764 | def regit(a,t): |
| 1765 | """ |
| 1766 | Regular (non-integer) iteration (works also for integer t). |
| 1767 | Preconditions: |
| 1768 | a[1]**n != a[1] for all n. |
| 1769 | Particularly a[1]!=0 and a[1]!=1. |
| 1770 | a[1]**t must be defined. |
| 1771 | |
| 1772 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1773 | sage: P = FormalPowerSeriesRing(QQ) |
| 1774 | sage: P([0,2]).regit(1/2)[:3] |
| 1775 | [0, sqrt(2), 0] |
| 1776 | """ |
| 1777 | |
| 1778 | a._assertp0() |
| 1779 | |
| 1780 | return Regit(a,t) |
| 1781 | |
| 1782 | def regit_b(a,t): |
| 1783 | """ |
| 1784 | Regular iteration via the schroeder function. |
| 1785 | |
| 1786 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1787 | sage: P = FormalPowerSeriesRing(QQ) |
| 1788 | sage: p = P([0,2]) |
| 1789 | sage: p.regit_b(1/2)[0:3] |
| 1790 | [0, sqrt(2), 0] |
| 1791 | """ |
| 1792 | s = a.schroeder() |
| 1793 | return s.inv()(s.rmul(a[1]**t)) |
| 1794 | |
| 1795 | |
| 1796 | |
| 1797 | __and__ = it |
| 1798 | |
| 1799 | def inv(a): |
| 1800 | """ |
| 1801 | The inverse. |
| 1802 | |
| 1803 | Precondition: a[1] != 0 |
| 1804 | Satisfies: a.inv()(a)=Id |
| 1805 | Alternative expression: a.inv() == ~a |
| 1806 | |
| 1807 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1808 | sage: P = FormalPowerSeriesRing(QQ) |
| 1809 | sage: P.Xexp.inv() - P.Lambert_w |
| 1810 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1811 | sage: ~P.Inv_selfroot_inc.dec().reclass() - P.Selfroot_inc.dec() |
| 1812 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1813 | """ |
| 1814 | a._assertp0() |
| 1815 | |
| 1816 | return Inv(a) |
| 1817 | |
| 1818 | __invert__ = inv |
| 1819 | |
| 1820 | # def julia_b(a): |
| 1821 | # """ |
| 1822 | # diff(it(a,t),t)(t=0) == ln(a[1])*julia(a) |
| 1823 | # """ |
| 1824 | # a._assertp0() |
| 1825 | |
| 1826 | # Poly=PolynomialRing(a.K,'x') |
| 1827 | # b = FormalPowerSeriesRing(Poly)() |
| 1828 | # b.min_index = 1 |
| 1829 | |
| 1830 | # def f(n): |
| 1831 | # """ sage: None # indirect doctest """ |
| 1832 | # if decidable0(a.K): |
| 1833 | # assert a[1] != 0 |
| 1834 | |
| 1835 | # if n == 0: |
| 1836 | # return Poly([0]) |
| 1837 | # if n == 1: |
| 1838 | # return Poly([0,1]) |
| 1839 | # res = a[n]*(b[1]**n)-b[1]*a[n] |
| 1840 | # res += sum([a[m]*b.npow(m)[n] - b[m]*a.npow(m)[n] for m in range(2,n)],a.K(0)) |
| 1841 | # res /= a[1]**n - a[1] |
| 1842 | # return res |
| 1843 | # b.coeffs = f |
| 1844 | |
| 1845 | # def h(p): |
| 1846 | # """ sage: None # indirect doctest """ |
| 1847 | # return sum([p.coeffs()[n]*n for n in range(p.degree()+1)],a.K(0)) |
| 1848 | |
| 1849 | # return a.new(lambda n: h(b[n])) |
| 1850 | |
| 1851 | def julia(a): |
| 1852 | """ |
| 1853 | Iterative logarithm or Julia function. |
| 1854 | Has different equivalent definitions: |
| 1855 | 1. Solution j of: j o a = a' * j, j[1]=1 |
| 1856 | 2. j = diff(f.it(t),t)(t=0) |
| 1857 | |
| 1858 | Precondition: a[1]**n!=a[1] for all n>1 |
| 1859 | |
| 1860 | It has similar properties like the logarithm: |
| 1861 | itlog(f^t) == t*itlog(f) |
| 1862 | |
| 1863 | It can be used to define the regular Abel function abel(f) by |
| 1864 | abel(f)' = 1/itlog(f) |
| 1865 | |
| 1866 | Refs: |
| 1867 | Eri Jabotinsky, Analytic iteration (1963), p 464 |
| 1868 | Jean Ecalle, Theorie des Invariants Holomorphes (1974), p 19 |
| 1869 | |
| 1870 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1871 | sage: P = FormalPowerSeriesRing(QQ) |
| 1872 | sage: a = P([0,2,1]) |
| 1873 | sage: j = a.julia() |
| 1874 | sage: j(a) - a.diff()*j |
| 1875 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1876 | """ |
| 1877 | |
| 1878 | return Julia(a) |
| 1879 | |
| 1880 | |
| 1881 | # def itlog(a): |
| 1882 | # """ |
| 1883 | # """ |
| 1884 | |
| 1885 | # #TODO this should be possible directly |
| 1886 | # _t = var('_t') |
| 1887 | # g = a.it(_t) |
| 1888 | # def f(n): |
| 1889 | # """ sage: None # indirect doctest """ |
| 1890 | # return diff(g[n],_t)(_t=0) |
| 1891 | # res = a.new(f) |
| 1892 | # res.min_index = res.val() |
| 1893 | # return res |
| 1894 | |
| 1895 | def schroeder(a): |
| 1896 | """ |
| 1897 | Returns the Schroeder powerseries s with s[1]=1 |
| 1898 | for a powerseries a with a[0]=0 and a[1]^n!=a[1] for all n |
| 1899 | |
| 1900 | The Schroeder powerseries s is determined up to a multiplicative |
| 1901 | constant by the functional equation: |
| 1902 | s(a(z))=a[1]*s(z) |
| 1903 | |
| 1904 | Let f(x) = 2*x + x**2, let s(x)=log(1+x), then s(f(x))=2*s(x): |
| 1905 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1906 | sage: P = FormalPowerSeriesRing(QQ) |
| 1907 | sage: P([0,2,1]).schroeder() - P.Log_inc |
| 1908 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1909 | """ |
| 1910 | a._assertp0() |
| 1911 | |
| 1912 | return Schroeder(a) |
| 1913 | |
| 1914 | def inv_schroeder(a): |
| 1915 | """ |
| 1916 | Returns the inverse Schroeder powerseries. |
| 1917 | |
| 1918 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1919 | sage: P = FormalPowerSeriesRing(QQ) |
| 1920 | sage: p = P([0,2,3]) |
| 1921 | sage: p.inv_schroeder() | p.schroeder() |
| 1922 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1923 | """ |
| 1924 | return InvSchroeder(a) |
| 1925 | |
| 1926 | def abel(f): |
| 1927 | """ |
| 1928 | The regular Abel function of a powerseries f (f[1]**n != f[1]) |
| 1929 | has the form a(x)=(ln(x)+ps(x))/ln(q) |
| 1930 | where q=f[1]!=0,1 and ps is a powerseries |
| 1931 | |
| 1932 | This method returns ps. |
| 1933 | |
| 1934 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing,FormalPowerSeries0 |
| 1935 | sage: a = var('a') |
| 1936 | sage: p = FormalPowerSeriesRing(PolynomialRing(QQ,a))(exp(a*x)-1,x) |
| 1937 | sage: pah = p.abel() |
| 1938 | sage: pac = p.abel2() |
| 1939 | sage: pah |
| 1940 | [0, 1/2*a/(-a + 1), (5/24*a^3 + 1/24*a^2)/(a^3 - a^2 - a + 1), ...] |
| 1941 | sage: [pac[k] - pah[k]==0 for k in range(0,5)] |
| 1942 | [True, True, True, True, True] |
| 1943 | """ |
| 1944 | f._assertp0() |
| 1945 | |
| 1946 | P = f.parent() |
| 1947 | return (f.schroeder()<<1).dec().reclass() | P.Log_inc |
| 1948 | |
| 1949 | def abel2(a): |
| 1950 | """ |
| 1951 | A different implementation of the regular Abel function via |
| 1952 | integration of 1/julia(a). |
| 1953 | |
| 1954 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing, FormalPowerSeries0 |
| 1955 | sage: a = var('a') |
| 1956 | sage: p = FormalPowerSeriesRing(PolynomialRing(QQ,a))(exp(a*x)-1,x) |
| 1957 | sage: p.abel2() |
| 1958 | [0, -1/2*a/(a - 1), (5/12*a^3 + 1/12*a^2)/(2*a^3 - 2*a^2 - 2*a + 2), ...] |
| 1959 | |
| 1960 | """ |
| 1961 | |
| 1962 | return a.julia().rcp().extinct_before(0).integral() |
| 1963 | |
| 1964 | |
| 1965 | class FormalPowerSeries01(FormalPowerSeries0): |
| 1966 | """ |
| 1967 | The FormalPowerSeriess p with p[0]==0 and p[1]==1. |
| 1968 | """ |
| 1969 | |
| 1970 | def valit(a): |
| 1971 | """ |
| 1972 | Returns the first index i such that a[i] != Id[i] |
| 1973 | |
| 1974 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 1975 | sage: P = FormalPowerSeriesRing(QQ) |
| 1976 | sage: P([0,1,0,0,1]).valit() |
| 1977 | 4 |
| 1978 | """ |
| 1979 | if not a[0] == 0: |
| 1980 | return 0 |
| 1981 | if not a[1] == 1: |
| 1982 | return 1 |
| 1983 | n = 2 |
| 1984 | while a[n] == 0: |
| 1985 | n+=1 |
| 1986 | return n |
| 1987 | |
| 1988 | # def it_b(p,t): |
| 1989 | # """ |
| 1990 | # A different direct implementation for `it'. |
| 1991 | |
| 1992 | # sage: p = P.Dec_exp.it_b(1/2) |
| 1993 | # sage: (p | p) - P.Dec_exp |
| 1994 | # [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 1995 | # """ |
| 1996 | # N=p.valit() |
| 1997 | # P = p.parent() |
| 1998 | # q = P() |
| 1999 | # def f(n): |
| 2000 | # """ sage: None # indirect doctest """ |
| 2001 | # if n < N: |
| 2002 | # return P.Id[n] |
| 2003 | # if n == N: |
| 2004 | # return t * p[N] |
| 2005 | # if n > N: |
| 2006 | # r=p[n] |
| 2007 | # r+=sum([p[m]*(q**m)[n] - q[m]*(p**m)[n] for m in range(N,n)]) |
| 2008 | # return r |
| 2009 | |
| 2010 | # q.coeffs = f |
| 2011 | # return q |
| 2012 | |
| 2013 | def regit(a,t): |
| 2014 | """ |
| 2015 | Regular iteration for powerseries with a[0]==0 and a[1]==1. |
| 2016 | The iteration index t needs not to be an integer. |
| 2017 | t should be a number in the base ring of this powerseries, |
| 2018 | but at least have a defined multiplication with the elements |
| 2019 | of the base ring. |
| 2020 | |
| 2021 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2022 | sage: P = FormalPowerSeriesRing(QQ) |
| 2023 | sage: p = P.Dec_exp.regit(1/2) |
| 2024 | sage: (p | p) - P.Dec_exp |
| 2025 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 2026 | """ |
| 2027 | |
| 2028 | return Regit01(a,t) |
| 2029 | |
| 2030 | def julia(a): |
| 2031 | """ |
| 2032 | Iterative logarithm or Julia function. |
| 2033 | Has different equivalent definitions: |
| 2034 | 1. Solution j of: j o a = a' * j, j[1]=1 |
| 2035 | 2. j = diff(f.it(t),t)(t=0) |
| 2036 | |
| 2037 | It has similar properties like the logarithm: |
| 2038 | itlog(f^t) == t*itlog(f) |
| 2039 | |
| 2040 | It can be used to define the regular Abel function abel(f) by |
| 2041 | abel(f)' = 1/itlog(f) |
| 2042 | |
| 2043 | Refs: |
| 2044 | Eri Jabotinsky, Analytic iteration (1963), p 464 |
| 2045 | Jean Ecalle, Theorie des Invariants Holomorphes (1974), p 19 |
| 2046 | |
| 2047 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2048 | sage: P = FormalPowerSeriesRing(QQ) |
| 2049 | sage: j = P.Dec_exp.julia() |
| 2050 | sage: j(P.Dec_exp) - P.Dec_exp.diff() * j |
| 2051 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 2052 | """ |
| 2053 | #diff(,t)(t=0) is the first coefficient of binomial(t,m) |
| 2054 | #Stirling1(m)[k] is the kth coefficient of m!*binomial(t,m) |
| 2055 | |
| 2056 | res = Julia01(a) |
| 2057 | #res.min_index = res.val() |
| 2058 | return res |
| 2059 | |
| 2060 | |
| 2061 | def abel_coeffs(a): |
| 2062 | """ |
| 2063 | The Abel function a of a power series f with f[0]=0 and f[1]=1 |
| 2064 | generally has the form |
| 2065 | F(z) + p(z), where |
| 2066 | F(z)=ji[-m]*z^{-m+1}/(-m+1)+...+ji[-2]*z^{-1}/(-1)+ji[-1]*log(z) |
| 2067 | and p is a powerseries (with positive indices only) |
| 2068 | |
| 2069 | ji[-1] is called the iterative residue (Ecalle) |
| 2070 | ji is the reciprocal of the Julia function |
| 2071 | (also called iterative logarithm) -which is meromorphic- of f |
| 2072 | |
| 2073 | The method returns the sequence [ji[-1],[F[-m+1],...,F[-1],p[0],p[1],p[2],...] |
| 2074 | |
| 2075 | These coefficients can be useful to compute the Abel function, |
| 2076 | for which the powerseries is a non-converging asymptotic development. |
| 2077 | |
| 2078 | The Abel function can then be gained by |
| 2079 | lim_{n->oo} F(f&n(z))-n |
| 2080 | |
| 2081 | |
| 2082 | sage: from sage.rings.formal_powerseries import * |
| 2083 | sage: P = FormalPowerSeriesRing(QQ) |
| 2084 | sage: P([0,1,0,2]).abel_coeffs() |
| 2085 | [3/2, [-1/4, 0; 0, 0, -5/4, 0, 21/8, 0, -35/4, 0, 2717/80, 0, -13429/100, 0, ...]] |
| 2086 | sage: p = P([0,1,0,0,1,2,3]) |
| 2087 | sage: a = p.abel_coeffs() |
| 2088 | sage: a |
| 2089 | [6, [-1/3, 1, -1; 0, -10, 11/2, 17/9, -169/12, 349/30, 13/18, -544/21, 1727/24, ...]] |
| 2090 | sage: ((p << 1).log().rmul(a[0]) + (p | a[1]) - a[1]).reclass() |
| 2091 | [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 2092 | """ |
| 2093 | |
| 2094 | juli = a.julia().rcp() |
| 2095 | # m = jul.val() |
| 2096 | # juli = (jul << m).rcp() |
| 2097 | # return [[ juli[m+i]/(i+1) for i in range(-m,-1) ],juli[m-1], (juli<<m).integral()] |
| 2098 | resit = juli[-1] |
| 2099 | #juli[-1]=0 |
| 2100 | return [resit,juli.set_item(-1,0).integral()] |
| 2101 | |
| 2102 | ### Constructors ### |
| 2103 | |
| 2104 | class Constant(FormalPowerSeries): |
| 2105 | def __init__(self,parent,c): |
| 2106 | """ |
| 2107 | Description and tests at FormalPowerSeriesRing.by_constant |
| 2108 | sage: None # indirect doctest |
| 2109 | """ |
| 2110 | FormalPowerSeries.__init__(self,parent) |
| 2111 | self.c = parent.K(c) |
| 2112 | |
| 2113 | def coeffs(self,n): |
| 2114 | """ sage: None # indirect doctest """ |
| 2115 | if n == 0: |
| 2116 | return self.c |
| 2117 | return 0 |
| 2118 | |
| 2119 | |
| 2120 | class Iterated(FormalPowerSeries): |
| 2121 | def __init__(self,parent,g,min_index): |
| 2122 | """ |
| 2123 | Description and tests at FormalPowerSeriesRing.by_generator |
| 2124 | sage: None # indirect doctest |
| 2125 | """ |
| 2126 | FormalPowerSeries.__init__(self,parent,min_index=min_index) |
| 2127 | self.g = g |
| 2128 | |
| 2129 | def coeffs(res,n): |
| 2130 | """ sage: None # indirect doctest """ |
| 2131 | g = res.g |
| 2132 | |
| 2133 | if n<res.min_index: |
| 2134 | return 0 |
| 2135 | if n==res.min_index: |
| 2136 | return g.next() |
| 2137 | x = res[n-1] #dummy to compute the prev value |
| 2138 | return g.next() |
| 2139 | |
| 2140 | class List(FormalPowerSeries): |
| 2141 | def __init__(self,parent,list,start): |
| 2142 | """ |
| 2143 | Description and tests at FormalPowerSeriesRing.by_list |
| 2144 | sage: None # indirect doctest |
| 2145 | """ |
| 2146 | FormalPowerSeries.__init__(self,parent) |
| 2147 | |
| 2148 | l = len(list) |
| 2149 | M=0 |
| 2150 | for M in range(l): |
| 2151 | if list[M] != 0: |
| 2152 | break |
| 2153 | N=-1 |
| 2154 | for i in range(l): |
| 2155 | N = l-i-1 |
| 2156 | if list[N] != 0: |
| 2157 | break |
| 2158 | |
| 2159 | self.min_index = start + M |
| 2160 | self.max_index = start + N |
| 2161 | self.start = start |
| 2162 | self.list = list |
| 2163 | |
| 2164 | |
| 2165 | def coeffs(self,k): |
| 2166 | """ sage: None # indirect doctest """ |
| 2167 | if k<self.min_index or k>self.max_index: |
| 2168 | return 0 |
| 2169 | return self.list[k-self.start] |
| 2170 | |
| 2171 | class Taylor(FormalPowerSeries): |
| 2172 | def __init__(self,parent,expr,v,at): |
| 2173 | """ |
| 2174 | Description and tests at FormalPowerSeriesRing.by_taylor |
| 2175 | sage: None # indirect doctest |
| 2176 | """ |
| 2177 | assert not v == None |
| 2178 | assert isinstance(v,Expression) #this should be Variable |
| 2179 | |
| 2180 | si = FormalPowerSeries.__init__ |
| 2181 | #coeffs always returns non-empty list, at least [0,0] is contained |
| 2182 | min_index = expr.taylor(v,at,2).substitute({v:v+at}).coeffs(v)[0][1] |
| 2183 | si(self,parent,min_index=min_index) |
| 2184 | |
| 2185 | self.expr = expr |
| 2186 | self.v = v |
| 2187 | self.at = at |
| 2188 | |
| 2189 | def coeffs(self,n): |
| 2190 | """ sage: None # indirect doctest """ |
| 2191 | #too slow |
| 2192 | #if not at == 0 and n==0: |
| 2193 | # return expr({v:at})-at |
| 2194 | #return simplify(diff(expr,v,n).substitute({v:at})/factorial(n)) |
| 2195 | expr = self.expr |
| 2196 | v = self.v |
| 2197 | at = self.at |
| 2198 | return self.K(expr.taylor(v,at,n).substitute({v:v+at}).coeff(v,n)) |
| 2199 | |
| 2200 | ### Constants ### |
| 2201 | |
| 2202 | class Zero(FormalPowerSeries): |
| 2203 | """ |
| 2204 | The zero element power series. |
| 2205 | |
| 2206 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2207 | sage: P = FormalPowerSeriesRing(QQ) |
| 2208 | sage: loads(dumps(P.Zero)) |
| 2209 | [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 2210 | """ |
| 2211 | def coeffs(self,n): |
| 2212 | """ sage: None # indirect doctest """ |
| 2213 | return 0 |
| 2214 | |
| 2215 | class One(FormalPowerSeries): |
| 2216 | """ |
| 2217 | The one element power series. |
| 2218 | |
| 2219 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2220 | sage: P = FormalPowerSeriesRing(QQ) |
| 2221 | sage: loads(dumps(P.One)) |
| 2222 | [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 2223 | """ |
| 2224 | def coeffs(self,n): |
| 2225 | """ sage: None # indirect doctest """ |
| 2226 | if n == 0: |
| 2227 | return self.K1 |
| 2228 | return 0 |
| 2229 | |
| 2230 | class Id(FormalPowerSeries01): |
| 2231 | """ |
| 2232 | The powerseries of the identity function id(x)=x. |
| 2233 | |
| 2234 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2235 | sage: P = FormalPowerSeriesRing(QQ) |
| 2236 | sage: loads(dumps(P.Id)) |
| 2237 | [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 2238 | """ |
| 2239 | def coeffs(self,n): |
| 2240 | """ sage: None # indirect doctest """ |
| 2241 | if n == 1: |
| 2242 | return self.K1 |
| 2243 | return 0 |
| 2244 | |
| 2245 | class Inc(FormalPowerSeries): |
| 2246 | """ |
| 2247 | The powerseries of the increment function inc(x)=x+1. |
| 2248 | |
| 2249 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2250 | sage: P = FormalPowerSeriesRing(QQ) |
| 2251 | sage: loads(dumps(P.Inc)) |
| 2252 | [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 2253 | """ |
| 2254 | def coeffs(self,n): |
| 2255 | """ sage: None # indirect doctest """ |
| 2256 | if n == 0 or n == 1: |
| 2257 | return self.K1 |
| 2258 | return 0 |
| 2259 | |
| 2260 | class Dec(FormalPowerSeries): |
| 2261 | """ |
| 2262 | The powerseries of the decrement function dec(x)=x-1. |
| 2263 | |
| 2264 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2265 | sage: P = FormalPowerSeriesRing(QQ) |
| 2266 | sage: loads(dumps(P.Dec)) |
| 2267 | [-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...] |
| 2268 | """ |
| 2269 | def coeffs(self,n): |
| 2270 | """ sage: None # indirect doctest """ |
| 2271 | if n == 0: |
| 2272 | return -self.K1 |
| 2273 | if n == 1: |
| 2274 | return self.K1 |
| 2275 | return 0 |
| 2276 | |
| 2277 | class Exp(FormalPowerSeries): |
| 2278 | """ |
| 2279 | The powerseries of the exponential function. |
| 2280 | |
| 2281 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2282 | sage: P = FormalPowerSeriesRing(QQ) |
| 2283 | sage: loads(dumps(P.Exp)) |
| 2284 | [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...] |
| 2285 | """ |
| 2286 | def coeffs(self,n): |
| 2287 | """ sage: None # indirect doctest """ |
| 2288 | return self.K1/factorial(n) |
| 2289 | |
| 2290 | class Dec_exp(FormalPowerSeries01): |
| 2291 | """ |
| 2292 | The powerseries of the decremented exponential dec_exp(x)=exp(x)-1. |
| 2293 | |
| 2294 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2295 | sage: P = FormalPowerSeriesRing(QQ) |
| 2296 | sage: loads(dumps(P.Dec_exp)) |
| 2297 | [0, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...] |
| 2298 | """ |
| 2299 | def coeffs(self,n): |
| 2300 | """ sage: None # indirect doctest """ |
| 2301 | if n == 0: return 0 |
| 2302 | return self.K1/factorial(n) |
| 2303 | |
| 2304 | class Log_inc(FormalPowerSeries01): |
| 2305 | """ |
| 2306 | The powerseries of the logarithm at 1 |
| 2307 | or the powerseries of f(x)=log(x+1) at 0. |
| 2308 | |
| 2309 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2310 | sage: P = FormalPowerSeriesRing(QQ) |
| 2311 | sage: loads(dumps(P.Log_inc)) |
| 2312 | [0, 1, -1/2, 1/3, -1/4, 1/5, -1/6, 1/7, -1/8, 1/9, -1/10, 1/11, -1/12, ...] |
| 2313 | """ |
| 2314 | def coeffs(self,n): |
| 2315 | """ sage: None # indirect doctest """ |
| 2316 | if n == 0: return 0 |
| 2317 | return self.K((-1)**(n+1)/Integer(n)) |
| 2318 | |
| 2319 | class Sin(FormalPowerSeries01): |
| 2320 | """ |
| 2321 | The powerseries of the sine. |
| 2322 | |
| 2323 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2324 | sage: P = FormalPowerSeriesRing(QQ) |
| 2325 | sage: loads(dumps(P.Sin)) |
| 2326 | [0, 1, 0, -1/6, 0, 1/120, 0, -1/5040, 0, 1/362880, 0, -1/39916800, 0, ...] |
| 2327 | """ |
| 2328 | def coeffs(self,n): |
| 2329 | """ sage: None # indirect doctest """ |
| 2330 | if n % 2 == 0: return 0 |
| 2331 | return self.K((-1)**((n-1)/2)/factorial(n)) |
| 2332 | |
| 2333 | class Cos(FormalPowerSeries): |
| 2334 | """ |
| 2335 | The powerseries of the cosine. |
| 2336 | |
| 2337 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2338 | sage: P = FormalPowerSeriesRing(QQ) |
| 2339 | sage: loads(dumps(P.Cos)) |
| 2340 | [1, 0, -1/2, 0, 1/24, 0, -1/720, 0, 1/40320, 0, -1/3628800, 0, 1/479001600, ...] |
| 2341 | """ |
| 2342 | def coeffs(self,n): |
| 2343 | """ sage: None # indirect doctest """ |
| 2344 | if n % 2 == 1: return 0 |
| 2345 | return self.K((-1)**(n/2)/factorial(n)) |
| 2346 | |
| 2347 | class Arcsin(FormalPowerSeries01): |
| 2348 | """ |
| 2349 | The powerseries of arcsin. |
| 2350 | |
| 2351 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2352 | sage: P = FormalPowerSeriesRing(QQ) |
| 2353 | sage: loads(dumps(P.Arcsin)) |
| 2354 | [0, 1, 0, 1/6, 0, 3/40, 0, 5/112, 0, 35/1152, 0, 63/2816, 0, 231/13312, 0, ...] |
| 2355 | """ |
| 2356 | def coeffs(self,n): |
| 2357 | """ sage: None # indirect doctest """ |
| 2358 | |
| 2359 | if n % 2 == 0: |
| 2360 | return 0 |
| 2361 | evenprod = Integer(1) |
| 2362 | oddprod = Integer(1) |
| 2363 | for k in range(2,n): |
| 2364 | if k % 2 == 0: |
| 2365 | evenprod *= k |
| 2366 | else: |
| 2367 | oddprod *=k |
| 2368 | return self.K(oddprod/evenprod/n) |
| 2369 | |
| 2370 | class Arctan(FormalPowerSeries01): |
| 2371 | """ |
| 2372 | The powerseries of arctan. |
| 2373 | |
| 2374 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2375 | sage: P = FormalPowerSeriesRing(QQ) |
| 2376 | sage: loads(dumps(P.Arctan)) |
| 2377 | [0, 1, 0, -1/3, 0, 1/5, 0, -1/7, 0, 1/9, 0, -1/11, 0, 1/13, 0, -1/15, 0, ...] |
| 2378 | """ |
| 2379 | def coeffs(self,n): |
| 2380 | """ sage: None # indirect doctest """ |
| 2381 | if n % 2 == 0: return 0 |
| 2382 | return self.K((-1)**(n/2)/Integer(n)) |
| 2383 | |
| 2384 | class Sinh(FormalPowerSeries01): |
| 2385 | """ |
| 2386 | The powerseries of sinh. |
| 2387 | |
| 2388 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2389 | sage: P = FormalPowerSeriesRing(QQ) |
| 2390 | sage: loads(dumps(P.Sinh)) |
| 2391 | [0, 1, 0, 1/6, 0, 1/120, 0, 1/5040, 0, 1/362880, 0, 1/39916800, 0, ...] |
| 2392 | """ |
| 2393 | def coeffs(self,n): |
| 2394 | """ sage: None # indirect doctest """ |
| 2395 | if n % 2 == 0: return 0 |
| 2396 | return self.K(1/factorial(n)) |
| 2397 | |
| 2398 | class Cosh(FormalPowerSeries): |
| 2399 | """ |
| 2400 | The powerseries of cosh. |
| 2401 | |
| 2402 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2403 | sage: P = FormalPowerSeriesRing(QQ) |
| 2404 | sage: loads(dumps(P.Cosh)) |
| 2405 | [1, 0, 1/2, 0, 1/24, 0, 1/720, 0, 1/40320, 0, 1/3628800, 0, 1/479001600, 0, ...] |
| 2406 | """ |
| 2407 | def coeffs(self,n): |
| 2408 | """ sage: None # indirect doctest """ |
| 2409 | if n % 2 == 1: return 0 |
| 2410 | return self.K(1/factorial(n)) |
| 2411 | |
| 2412 | class Arcsinh(FormalPowerSeries01): |
| 2413 | """ |
| 2414 | The powerseries of arcsinh. |
| 2415 | |
| 2416 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2417 | sage: P = FormalPowerSeriesRing(QQ) |
| 2418 | sage: loads(dumps(P.Arcsinh)) |
| 2419 | [0, 1, 0, -1/6, 0, 3/40, 0, -5/112, 0, 35/1152, 0, -63/2816, 0, 231/13312, ...] |
| 2420 | """ |
| 2421 | def coeffs(self,n): |
| 2422 | """ sage: None # indirect doctest """ |
| 2423 | if n % 2 == 0: |
| 2424 | return 0 |
| 2425 | evenprod = Integer(1) |
| 2426 | oddprod = Integer(1) |
| 2427 | for k in range(2,n): |
| 2428 | if k % 2 == 0: |
| 2429 | evenprod *= k |
| 2430 | else: |
| 2431 | oddprod *= k |
| 2432 | return self.K((-1)**(n/2)*oddprod/evenprod/n) |
| 2433 | |
| 2434 | class Arctanh(FormalPowerSeries01): |
| 2435 | """ |
| 2436 | The powerseries of arctanh. |
| 2437 | |
| 2438 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2439 | sage: P = FormalPowerSeriesRing(QQ) |
| 2440 | sage: loads(dumps(P.Arctanh)) |
| 2441 | [0, 1, 0, 1/3, 0, 1/5, 0, 1/7, 0, 1/9, 0, 1/11, 0, 1/13, 0, 1/15, 0, 1/17, ...] |
| 2442 | """ |
| 2443 | def coeffs(self,n): |
| 2444 | """ sage: None # indirect doctest """ |
| 2445 | if n % 2 == 0: return 0 |
| 2446 | return self.K(1/Integer(n)) |
| 2447 | |
| 2448 | class Tan(FormalPowerSeries01): |
| 2449 | """ |
| 2450 | The powerseries of tan. |
| 2451 | |
| 2452 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2453 | sage: P = FormalPowerSeriesRing(QQ) |
| 2454 | sage: loads(dumps(P.Tan)) |
| 2455 | [0, 1, 0, 1/3, 0, 2/15, 0, 17/315, 0, 62/2835, 0, 1382/155925, 0, ...] |
| 2456 | """ |
| 2457 | def coeffs(self,N): |
| 2458 | """ sage: None # indirect doctest """ |
| 2459 | if N % 2 == 0: |
| 2460 | return 0 |
| 2461 | n = (N + 1) / 2 |
| 2462 | P = self.parent() |
| 2463 | return P.K(P.Bernoulli[2*n] * (-4)**n * (1-4**n) / factorial(2*n)) |
| 2464 | |
| 2465 | class Tanh(FormalPowerSeries01): |
| 2466 | """ |
| 2467 | The powerseries of tanh. |
| 2468 | |
| 2469 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2470 | sage: P = FormalPowerSeriesRing(QQ) |
| 2471 | sage: loads(dumps(P.Tanh)) |
| 2472 | [0, 1, 0, -1/3, 0, 2/15, 0, -17/315, 0, 62/2835, 0, -1382/155925, 0, ...] |
| 2473 | """ |
| 2474 | def coeffs(self,N): |
| 2475 | """ sage: None # indirect doctest """ |
| 2476 | if N % 2 == 0: |
| 2477 | return 0 |
| 2478 | n = (N+1)/2 |
| 2479 | P = self.parent() |
| 2480 | return P.K(P.Bernoulli[2*n] * (-1)**(2*n) * 4**n * (4**n-1) / factorial(2*n)) |
| 2481 | |
| 2482 | class Xexp(FormalPowerSeries01): |
| 2483 | """ |
| 2484 | The powerseries of x*exp(x) |
| 2485 | |
| 2486 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2487 | sage: P = FormalPowerSeriesRing(QQ) |
| 2488 | sage: loads(dumps(P.Xexp)) |
| 2489 | [0, 1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, ...] |
| 2490 | """ |
| 2491 | def coeffs(self,n): |
| 2492 | """ sage: None # indirect doctest """ |
| 2493 | if n==0: return 0 |
| 2494 | return self.K(1/factorial(n-1)) |
| 2495 | |
| 2496 | class Lambert_w(FormalPowerSeries01): |
| 2497 | """ |
| 2498 | The Lambert W function is the inverse of f(x)=x*e^x |
| 2499 | |
| 2500 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2501 | sage: P = FormalPowerSeriesRing(QQ) |
| 2502 | sage: loads(dumps(P.Lambert_w)) |
| 2503 | [0, 1, -1, 3/2, -8/3, 125/24, -54/5, 16807/720, -16384/315, 531441/4480, ...] |
| 2504 | """ |
| 2505 | def coeffs(self,n): |
| 2506 | """ sage: None # indirect doctest """ |
| 2507 | if n==0: return 0 |
| 2508 | return self.K((-n)**(n-1)/factorial(n)) |
| 2509 | |
| 2510 | class Sqrt_inc(FormalPowerSeries): |
| 2511 | """ |
| 2512 | The powerseries of sqrt at 1, or of sqrt(x+1) at 0. |
| 2513 | |
| 2514 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2515 | sage: P = FormalPowerSeriesRing(QQ) |
| 2516 | sage: loads(dumps(P.Sqrt_inc)) |
| 2517 | [1, 1/2, -1/8, 1/16, -5/128, 7/256, -21/1024, 33/2048, -429/32768, ...] |
| 2518 | """ |
| 2519 | def coeffs(self,n): |
| 2520 | """ sage: None # indirect doctest """ |
| 2521 | evenprod=Integer(1) |
| 2522 | oddprod=Integer(1) |
| 2523 | for k in range(2,2*n+1): |
| 2524 | if k%2 == 0: |
| 2525 | evenprod *= k |
| 2526 | else: |
| 2527 | oddprod *= k |
| 2528 | return self.K((-1)**n *oddprod/evenprod/(1-2*n)) |
| 2529 | |
| 2530 | # class BernoulliVar(FormalPowerSeries): |
| 2531 | # """ |
| 2532 | # The variant of the Bernoulli numbers where B[1]=1/2 |
| 2533 | # instead of B[1]=-1/2 |
| 2534 | # """ |
| 2535 | # def coeffs(self,n): |
| 2536 | # if n == 1: |
| 2537 | # return self.K(1)/self.K(2) |
| 2538 | # return self.Bernoulli[n] |
| 2539 | |
| 2540 | # class Faulhaber(FormalPowerSeries): |
| 2541 | # """ |
| 2542 | # Faulhaber's formula. |
| 2543 | # """ |
| 2544 | # def coeffs(self,n): |
| 2545 | |
| 2546 | |
| 2547 | class Stirling1(FormalPowerSeries): |
| 2548 | """ |
| 2549 | Returns the sequence of Stirling numbers of the first kind. |
| 2550 | These are the coefficients of the polynomial x(x-1)(x-2)...(x-n+1). |
| 2551 | stirling1[n][k] is the coefficient of x^k in the above polynomial. |
| 2552 | """ |
| 2553 | def coeffs(self,n): |
| 2554 | """ sage: None # indirect doctest """ |
| 2555 | P = self.parent() |
| 2556 | if n==0: |
| 2557 | return P.One |
| 2558 | |
| 2559 | res = Stirling1_Succ(P,self[n-1],n) |
| 2560 | return res |
| 2561 | |
| 2562 | class Stirling1_Succ(FormalPowerSeries): |
| 2563 | """ |
| 2564 | Constructs the sequence of Stirling numbers stirling1[n] |
| 2565 | from g = stirling1[n-1] |
| 2566 | """ |
| 2567 | def __init__(self,parent,g,n): |
| 2568 | """ sage: None # indirect doctest """ |
| 2569 | FormalPowerSeries.__init__(self,parent,min_index=1) |
| 2570 | self.g = g |
| 2571 | self.n = n |
| 2572 | |
| 2573 | def coeffs(self,k): |
| 2574 | """ sage: None # indirect doctest """ |
| 2575 | g = self.g |
| 2576 | n = self.n |
| 2577 | return g[k-1]-(n-1)*g[k] |
| 2578 | |
| 2579 | |
| 2580 | class Lehmer_comtet(FormalPowerSeries): |
| 2581 | """ |
| 2582 | The n-th Lehmer-Comtet number is the n-th derivative of x^x at 1. |
| 2583 | See Sloane A008296. |
| 2584 | |
| 2585 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2586 | sage: P = FormalPowerSeriesRing(QQ) |
| 2587 | sage: loads(dumps(P.Lehmer_comtet)) |
| 2588 | [1, 1, 3, 10, 41, 196, 1057, 6322, 41393, 293608, 2237921, 18210094, ...] |
| 2589 | """ |
| 2590 | def coeffs(self,n): |
| 2591 | """ sage: None # indirect doctest """ |
| 2592 | return self.K(sum([k**(n-k)*binomial(n,k) for k in range(n+1)])) |
| 2593 | |
| 2594 | class Selfpower_inc(FormalPowerSeries): |
| 2595 | """ |
| 2596 | The powerseries of x^x at 1. |
| 2597 | |
| 2598 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2599 | sage: P = FormalPowerSeriesRing(QQ) |
| 2600 | sage: loads(dumps(P.Selfpower_inc)) |
| 2601 | [1, 1, 1, 1/2, 1/3, 1/12, 3/40, -1/120, 59/2520, -71/5040, 131/10080, ...] |
| 2602 | """ |
| 2603 | def coeffs(self,n): |
| 2604 | """ sage: None # indirect doctest """ |
| 2605 | P = self.parent() |
| 2606 | return P.K(sum([P.Stirling1[n][k]*P.A000248[k] for k in range(n+1)]))/factorial(n) |
| 2607 | |
| 2608 | class Superroot_inc(FormalPowerSeries): |
| 2609 | """ |
| 2610 | The powerseries of the inverse of x^x developed at 1. |
| 2611 | |
| 2612 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2613 | sage: P = FormalPowerSeriesRing(QQ) |
| 2614 | sage: loads(dumps(P.Superroot_inc)) |
| 2615 | [1, 1, -1, 3/2, -17/6, 37/6, -1759/120, 13279/360, -97283/1008, ...] |
| 2616 | """ |
| 2617 | def coeffs(self,n): |
| 2618 | """ sage: None # indirect doctest """ |
| 2619 | P = self.parent() |
| 2620 | return P.K(sum([P.Stirling1[n][k]*Integer(1-k)**(k-1) for k in range(n+1)]))/factorial(n) |
| 2621 | |
| 2622 | class A003725(FormalPowerSeries): |
| 2623 | """ |
| 2624 | The derivatives of exp(x*e^(-x)) at 0. |
| 2625 | |
| 2626 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2627 | sage: P = FormalPowerSeriesRing(QQ) |
| 2628 | sage: loads(dumps(P.A003725)) |
| 2629 | [1, 1, -1, -2, 9, -4, -95, 414, 49, -10088, 55521, -13870, -2024759, ...] |
| 2630 | """ |
| 2631 | def coeffs(self,n): |
| 2632 | """ sage: None # indirect doctest """ |
| 2633 | return self.K(sum([ (-k)**(n-k)*binomial(n, k) for k in range(n+1)])) |
| 2634 | |
| 2635 | class Selfroot_inc(FormalPowerSeries): |
| 2636 | """ |
| 2637 | The powerseries of x^(1/x) at 1. |
| 2638 | |
| 2639 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2640 | sage: P = FormalPowerSeriesRing(QQ) |
| 2641 | sage: loads(dumps(P.Selfroot_inc)) |
| 2642 | [1, 1, -1, 1/2, 1/6, -3/4, 131/120, -9/8, 1087/1260, -271/720, -2291/10080, ...] |
| 2643 | """ |
| 2644 | def coeffs(self,n): |
| 2645 | """ sage: None # indirect doctest """ |
| 2646 | P = self.parent() |
| 2647 | return P.K(sum([P.Stirling1[n][k]*P.A003725[k] for k in range(n+1)]))/factorial(n) |
| 2648 | |
| 2649 | class Inv_selfroot_inc(FormalPowerSeries): |
| 2650 | """ |
| 2651 | The inverse of the self root x^(1/x) at 1. |
| 2652 | The inverse of the self root for x=b is the fixed point of f(y)=b^y. |
| 2653 | |
| 2654 | sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing |
| 2655 | sage: P = FormalPowerSeriesRing(QQ) |
| 2656 | sage: loads(dumps(P.Inv_selfroot_inc)) |
| 2657 | [1, 1, 1, 3/2, 7/3, 4, 283/40, 4681/360, 123101/5040, 118001/2520, ...] |
| 2658 | """ |
| 2659 | def coeffs(self,n): |
| 2660 | """ sage: None # indirect doctest """ |
| 2661 | P = self.parent() |
| 2662 | if n<0: |
| 2663 | return 0 |
| 2664 | if n==0: |
| 2665 | return P.K1 |
| 2666 | r = 0 |
| 2667 | for k in range(1,n+1): |
| 2668 | r += P.Stirling1[n][k]*P.K((k+1))**(k-1) |
| 2669 | |
| 2670 | return P.K(r)/factorial(n) |
| 2671 | |
| 2672 | ### Methods ### |
| 2673 | |
| 2674 | class ExtinctBefore(FormalPowerSeries): |
| 2675 | def __init__(self,a,min_index): |
| 2676 | """ |
| 2677 | Description and tests at FormalPowerSeries.extinct_before |
| 2678 | sage: None # indirect doctest |
| 2679 | """ |
| 2680 | FormalPowerSeries.extinct_before.__doc__ |
| 2681 | |
| 2682 | si = FormalPowerSeries.__init__ |
| 2683 | si(self,a.parent()) |
| 2684 | self.a=a |
| 2685 | self.min_index=min_index |
| 2686 | |
| 2687 | def coeffs(self,n): |
| 2688 | """ sage: None # indirect doctest """ |
| 2689 | if n < self.min_index: |
| 2690 | return 0 |
| 2691 | return self.a[n] |
| 2692 | |
| 2693 | class MulFact(FormalPowerSeries): |
| 2694 | def __init__(self,a): |
| 2695 | """ |
| 2696 | Description and tests at FormalPowerSeries.mul_fact |
| 2697 | sage: None # indirect doctest |
| 2698 | """ |
| 2699 | FormalPowerSeries.__init__(self,a.parent(),min_index=a.min_index) |
| 2700 | self.a = a |
| 2701 | def coeffs(self,n): |
| 2702 | """ sage: None # indirect doctest """ |
| 2703 | return self.a[n]*self.K(factorial(n)) |
| 2704 | |
| 2705 | class DivFact(FormalPowerSeries): |
| 2706 | def __init__(self,a): |
| 2707 | """ |
| 2708 | Description and tests at FormalPowerSeries.div_fact |
| 2709 | sage: None # indirect doctest |
| 2710 | """ |
| 2711 | FormalPowerSeries.__init__(self,a.parent(),min_index=a.min_index) |
| 2712 | self.a = a |
| 2713 | def coeffs(self,n): |
| 2714 | """ sage: None # indirect doctest """ |
| 2715 | return self.a[n]/self.K(factorial(n)) |
| 2716 | |
| 2717 | class IncMethod(FormalPowerSeries): |
| 2718 | def __init__(self,a): |
| 2719 | """ |
| 2720 | Description and tests at FormalPowerSeries.inc |
| 2721 | sage: None # indirect doctest |
| 2722 | """ |
| 2723 | FormalPowerSeries.__init__(self,a.parent()) |
| 2724 | self.a = a |
| 2725 | def coeffs(self,n): |
| 2726 | """ sage: None # indirect doctest """ |
| 2727 | if n == 0: |
| 2728 | return self.a[0] + self.K1 |
| 2729 | return self.a[n] |
| 2730 | |
| 2731 | class DecMethod(FormalPowerSeries): |
| 2732 | def __init__(self,a): |
| 2733 | """ |
| 2734 | Description and tests at FormalPowerSeries.dec |
| 2735 | sage: None # indirect doctest |
| 2736 | """ |
| 2737 | FormalPowerSeries.__init__(self,a.parent()) |
| 2738 | self.a = a |
| 2739 | def coeffs(self,n): |
| 2740 | """ sage: None # indirect doctest """ |
| 2741 | if n == 0: |
| 2742 | return self.a[0] - self.K1 |
| 2743 | return self.a[n] |
| 2744 | |
| 2745 | class RMul(FormalPowerSeries): |
| 2746 | def __init__(self,a,s): |
| 2747 | """ |
| 2748 | Description and tests at FormalPowerSeries.rmul |
| 2749 | sage: None # indirect doctest |
| 2750 | """ |
| 2751 | FormalPowerSeries.__init__(self,a.parent(),min_index=a.min_index) |
| 2752 | self.a = a |
| 2753 | self.s = s |
| 2754 | |
| 2755 | def coeffs(self,n): |
| 2756 | """ sage: None # indirect doctest """ |
| 2757 | return self.a[n] * self.s |
| 2758 | |
| 2759 | class LMul(FormalPowerSeries): |
| 2760 | def __init__(self,a,s): |
| 2761 | """ |
| 2762 | Description and tests at FormalPowerSeries.lmul |
| 2763 | sage: None # indirect doctest |
| 2764 | """ |
| 2765 | FormalPowerSeries.__init__(self,a.parent(),min_index=a.min_index) |
| 2766 | self.a = a |
| 2767 | self.s = s |
| 2768 | |
| 2769 | def coeffs(self,n): |
| 2770 | """ sage: None # indirect doctest """ |
| 2771 | return self.s * self.a[n] |
| 2772 | |
| 2773 | class Add(FormalPowerSeries): |
| 2774 | def __init__(self,a,b): |
| 2775 | """ |
| 2776 | Description and tests at FormalPowerSeries.add |
| 2777 | sage: None # indirect doctest |
| 2778 | """ |
| 2779 | si = FormalPowerSeries.__init__ |
| 2780 | si(self,a.parent(),min_index=min(a.min_index,b.min_index)) |
| 2781 | self.a = a; self.b = b |
| 2782 | |
| 2783 | def coeffs(self,n): |
| 2784 | """ sage: None # indirect doctest """ |
| 2785 | a = self.a; b = self.b |
| 2786 | if n < a.min_index: |
| 2787 | if n < b.min_index: |
| 2788 | return 0 |
| 2789 | return b[n] |
| 2790 | if n < b.min_index: |
| 2791 | return a[n] |
| 2792 | return a[n]+b[n] |
| 2793 | |
| 2794 | class Sub(FormalPowerSeries): |
| 2795 | def __init__(self,a,b): |
| 2796 | """ |
| 2797 | Description and tests at FormalPowerSeries.sub |
| 2798 | sage: None # indirect doctest |
| 2799 | """ |
| 2800 | si = FormalPowerSeries.__init__ |
| 2801 | si(self,a.parent(),min_index=min(a.min_index,b.min_index)) |
| 2802 | self.a = a; self.b = b |
| 2803 | |
| 2804 | def coeffs(self,n): |
| 2805 | """ sage: None # indirect doctest """ |
| 2806 | a = self.a; b = self.b |
| 2807 | if n < a.min_index: |
| 2808 | if n < b.min_index: |
| 2809 | return 0 |
| 2810 | #a[0]==0 |
| 2811 | return -b[n] |
| 2812 | if n < b.min_index: |
| 2813 | #b[0]==0 |
| 2814 | return a[n] |
| 2815 | return a[n]-b[n] |
| 2816 | |
| 2817 | class Neg(FormalPowerSeries): |
| 2818 | def __init__(self,a): |
| 2819 | """ |
| 2820 | Description and tests at FormalPowerSeries.neg |
| 2821 | sage: None # indirect doctest |
| 2822 | """ |
| 2823 | si = FormalPowerSeries.__init__ |
| 2824 | si(self,a.parent(),min_index=a.min_index) |
| 2825 | self.a = a |
| 2826 | |
| 2827 | def coeffs(self,n): |
| 2828 | """ sage: None # indirect doctest """ |
| 2829 | a = self.a |
| 2830 | if n < a.min_index: |
| 2831 | return 0 |
| 2832 | return -a[n] |
| 2833 | |
| 2834 | class Mul(FormalPowerSeries): |
| 2835 | def __init__(self,a,b): |
| 2836 | """ |
| 2837 | Description and tests at FormalPowerSeries.mul |
| 2838 | sage: None # indirect doctest |
| 2839 | """ |
| 2840 | si = FormalPowerSeries.__init__ |
| 2841 | si(self,a.parent(),min_index=a.min_index+b.min_index) |
| 2842 | self.a,self.b = a,b |
| 2843 | |
| 2844 | def ab(self,m,n): |
| 2845 | """ |
| 2846 | Lazy product of a[m] and b[n] |
| 2847 | sage: None # indirect doctest |
| 2848 | """ |
| 2849 | a = self.a |
| 2850 | b = self.b |
| 2851 | if a[m] == 0: |
| 2852 | return 0 |
| 2853 | return a[m]*b[n] |
| 2854 | |
| 2855 | def coeffs(self,n): |
| 2856 | """ sage: None # indirect doctest """ |
| 2857 | a,b = self.a,self.b |
| 2858 | return sum([self.ab(k,n-k) for k in range(a.min_index,n+1-b.min_index)]) |
| 2859 | |
| 2860 | class Div(FormalPowerSeries): |
| 2861 | def __init__(self,c,b): |
| 2862 | """ |
| 2863 | Description and tests at FormalPowerSeries.div |
| 2864 | sage: None # indirect doctest |
| 2865 | """ |
| 2866 | si = FormalPowerSeries.__init__ |
| 2867 | si(self,c.parent(),min_index=c.min_index - b.min_index) |
| 2868 | self.c=c |
| 2869 | self.b=b |
| 2870 | |
| 2871 | def _ab(a,m,n): |
| 2872 | """ |
| 2873 | Lazy product of b[n] and a[m]. |
| 2874 | sage: None # indirect doctest |
| 2875 | """ |
| 2876 | b = a.b |
| 2877 | if b[n] == b.K(0): |
| 2878 | return b.K(0) |
| 2879 | return a[m]*b[n] |
| 2880 | |
| 2881 | def coeffs(a,n): |
| 2882 | """ sage: None # indirect doctest """ |
| 2883 | c = a.c |
| 2884 | b = a.b |
| 2885 | if n<a.min_index: |
| 2886 | return 0 |
| 2887 | r = c[n+b.min_index] |
| 2888 | for k in range(a.min_index,n): |
| 2889 | r -= a._ab(k,n+b.min_index-k) |
| 2890 | return r/b[b.min_index] |
| 2891 | |
| 2892 | class Npow(FormalPowerSeries): |
| 2893 | def __init__(self,f,m): |
| 2894 | """ |
| 2895 | Description and tests at FormalPowerSeries.npow |
| 2896 | sage: None # indirect doctest |
| 2897 | """ |
| 2898 | si = FormalPowerSeries.__init__ |
| 2899 | si(self,f.parent(),min_index=f.min_index*m) |
| 2900 | self.f = f |
| 2901 | self.m = m |
| 2902 | |
| 2903 | def coeffs(self,n): |
| 2904 | """ sage: None # indirect doctest """ |
| 2905 | m = self.m |
| 2906 | if n<self.min_index: |
| 2907 | return 0 |
| 2908 | return self.f._s(n-self.f.min_index*(m-1),m,n) |
| 2909 | |
| 2910 | class Nipow(FormalPowerSeries): |
| 2911 | def __init__(self,a,t): |
| 2912 | """ |
| 2913 | Description and tests at FormalPowerSeries.nipow |
| 2914 | sage: None # indirect doctest |
| 2915 | """ |
| 2916 | si = FormalPowerSeries.__init__ |
| 2917 | si(self,a.parent()) |
| 2918 | self.a = a |
| 2919 | self.t = t |
| 2920 | |
| 2921 | def coeffs(s,n): |
| 2922 | """ sage: None # indirect doctest """ |
| 2923 | a = s.a |
| 2924 | t = s.t |
| 2925 | if decidable0(a.K): |
| 2926 | assert a[0] != 0, "0th coefficient is " + repr(a[0]) + ", but must be non-zero for non-integer powers" |
| 2927 | |
| 2928 | da = a.set_item(0,0) |
| 2929 | |
| 2930 | if n>=0 and a[0] == 1: |
| 2931 | #dont use symbolic arithmetic for ratonal powers of 1 |
| 2932 | return sum([binomial(t,k) * da.npow(k)[n] for k in range(n+1)],a.K(0)) |
| 2933 | return sum([binomial(t,k) * a[0]**t/a[0]**k * da.npow(k)[n] for k in range(n+1)],a.K(0)) |
| 2934 | |
| 2935 | class Compose(FormalPowerSeries): |
| 2936 | def __init__(self,b,a): |
| 2937 | """ |
| 2938 | Description and tests at FormalPowerSeries.compose |
| 2939 | sage: None # indirect doctest |
| 2940 | """ |
| 2941 | si = FormalPowerSeries.__init__ |
| 2942 | si(self,b.parent(),min_index=b.min_index*a.min_index) |
| 2943 | self.b = b |
| 2944 | self.a = a |
| 2945 | |
| 2946 | def coeffs(self,n): |
| 2947 | """ sage: None # indirect doctest """ |
| 2948 | b = self.b |
| 2949 | a = self.a |
| 2950 | res = sum([b[k]*(a.npow(k)[n]) for k in range(n+1)]) |
| 2951 | if b.min_index < 0: |
| 2952 | bi = a.rcp() |
| 2953 | res += sum([b[k]*(bi.npow(-k)[n]) for k in range(b.min_index,0)],b.K(0)) |
| 2954 | return res |
| 2955 | class Lshift(FormalPowerSeries): |
| 2956 | def __init__(self,a,m=1): |
| 2957 | """ |
| 2958 | Description and tests at FormalPowerSeries.__lshift__ |
| 2959 | sage: None # indirect doctest |
| 2960 | """ |
| 2961 | si = FormalPowerSeries.__init__ |
| 2962 | si(self,a.parent()) |
| 2963 | self.a = a |
| 2964 | self.m = m |
| 2965 | |
| 2966 | def coeffs(self,n): |
| 2967 | """ sage: None # indirect doctest """ |
| 2968 | return self.a[n+self.m] |
| 2969 | |
| 2970 | class Rshift(FormalPowerSeries): |
| 2971 | def __init__(self,a,m=1): |
| 2972 | """ |
| 2973 | Description and tests at FormalPowerSeries.__rshift__ |
| 2974 | sage: None # indirect doctest |
| 2975 | """ |
| 2976 | si = FormalPowerSeries.__init__ |
| 2977 | si(self,a.parent()) |
| 2978 | self.a = a |
| 2979 | self.m = m |
| 2980 | |
| 2981 | def coeffs(self,n): |
| 2982 | """ sage: None # indirect doctest """ |
| 2983 | a = self.a |
| 2984 | m = self.m |
| 2985 | if n < m: |
| 2986 | return 0 |
| 2987 | return a[n-m] |
| 2988 | |
| 2989 | class Diff(FormalPowerSeries): |
| 2990 | def __init__(self,a,m=1): |
| 2991 | """ |
| 2992 | Description and tests at FormalPowerSeries.diff |
| 2993 | sage: None # indirect doctest |
| 2994 | """ |
| 2995 | si = FormalPowerSeries.__init__ |
| 2996 | si(self,a.parent(),min_index=self.deg(a.min_index,m)) |
| 2997 | self.a = a |
| 2998 | self.m = m |
| 2999 | |
| 3000 | def deg(self,v,m): |
| 3001 | """ sage: None # indirect doctest """ |
| 3002 | if v >= 0: |
| 3003 | return max(v-m,0) |
| 3004 | return v-m |
| 3005 | |
| 3006 | def coeffs(self,n): |
| 3007 | """ sage: None # indirect doctest """ |
| 3008 | a = self.a |
| 3009 | m = self.m |
| 3010 | if -m <= n and n < 0: |
| 3011 | return 0 |
| 3012 | else: |
| 3013 | return a[n+m]*prod(k for k in range(n+1,n+m+1)) |
| 3014 | |
| 3015 | class Integral(FormalPowerSeries): |
| 3016 | def __init__(self,a,c=0): |
| 3017 | """ |
| 3018 | Description and tests at FormalPowerSeries.integral |
| 3019 | sage: None # indirect doctest |
| 3020 | """ |
| 3021 | si = FormalPowerSeries.__init__ |
| 3022 | if c == 0: |
| 3023 | si(self,a.parent(),min_index=a.min_index+1) |
| 3024 | else: |
| 3025 | si(self,a.parent()) |
| 3026 | |
| 3027 | self.a = a |
| 3028 | self.c = c |
| 3029 | |
| 3030 | def coeffs(self,n): |
| 3031 | """ sage: None # indirect doctest """ |
| 3032 | a = self.a |
| 3033 | c = self.c |
| 3034 | |
| 3035 | if n == 0: |
| 3036 | return c |
| 3037 | if n < a.min_index+1: |
| 3038 | return 0 |
| 3039 | |
| 3040 | #This test can not be moved to the beginning because |
| 3041 | #it would break lazyness |
| 3042 | if a.min_index < 0: |
| 3043 | assert a[-1] == 0, "Coefficient at -1 must be 0, but it is: " + repr(a[-1]) |
| 3044 | return a[n-1]/Integer(n) |
| 3045 | |
| 3046 | class N(FormalPowerSeries): |
| 3047 | def __init__(self,a,*args,**kwargs): |
| 3048 | """ |
| 3049 | Description and tests at FormalPowerSeries.n |
| 3050 | sage: None # indirect doctest |
| 3051 | """ |
| 3052 | si = FormalPowerSeries.__init__ |
| 3053 | si(self,FormalPowerSeriesRing(n(0,*args,**kwargs).parent())) |
| 3054 | |
| 3055 | self.a = a |
| 3056 | self.args = args |
| 3057 | self.kwargs = kwargs |
| 3058 | |
| 3059 | def coeffs(self,k): |
| 3060 | """ sage: None # indirect doctest """ |
| 3061 | return n(self.a[k],*self.args,**self.kwargs) |
| 3062 | |
| 3063 | class Regit(FormalPowerSeries0): |
| 3064 | def __init__(self,a,t): |
| 3065 | """ |
| 3066 | Description and tests at FormalPowerSeries.regit |
| 3067 | sage: None # indirect doctest |
| 3068 | """ |
| 3069 | si = FormalPowerSeries.__init__ |
| 3070 | si(self,a.parent(),min_index=1) |
| 3071 | self.a = a |
| 3072 | self.t = t |
| 3073 | |
| 3074 | def coeffs(b,n): |
| 3075 | """ sage: None # indirect doctest """ |
| 3076 | a = b.a |
| 3077 | t = b.t |
| 3078 | if decidable0(a.K): |
| 3079 | assert a[1]!=1, a[1] |
| 3080 | assert a[1]!=0, a[1] |
| 3081 | |
| 3082 | #print "(" + repr(n) |
| 3083 | if n == 0: |
| 3084 | #print ")" |
| 3085 | return 0 |
| 3086 | if n == 1: |
| 3087 | #print ")" |
| 3088 | return a[1]**t |
| 3089 | res = a[n]*(b[1]**n)-b[1]*a[n] |
| 3090 | |
| 3091 | for m in range(2,n): |
| 3092 | res += a[m]*b.npow(m)[n] - b[m]*a.npow(m)[n] |
| 3093 | |
| 3094 | res /= a[1]**n - a[1] |
| 3095 | #print ")" |
| 3096 | return res |
| 3097 | |
| 3098 | class Inv(FormalPowerSeries0): |
| 3099 | def __init__(self,a): |
| 3100 | """ |
| 3101 | Description and tests at FormalPowerSeries.inv |
| 3102 | sage: None # indirect doctest |
| 3103 | """ |
| 3104 | si = FormalPowerSeries.__init__ |
| 3105 | si(self,a.parent(),min_index=1) |
| 3106 | self.a=a |
| 3107 | |
| 3108 | def coeffs(b,n): |
| 3109 | """ sage: None # indirect doctest """ |
| 3110 | a = b.a |
| 3111 | if n==0: |
| 3112 | return 0 |
| 3113 | if n==1: |
| 3114 | return 1/a[1] |
| 3115 | if n>1: |
| 3116 | return - sum([ b[k]*a.npow(k)[n] for k in range(1,n)])/a[1]**n |
| 3117 | |
| 3118 | class Julia(FormalPowerSeries01): |
| 3119 | def __init__(self,a): |
| 3120 | """ |
| 3121 | Description and tests at FormalPowerSeries.julia |
| 3122 | sage: None # indirect doctest |
| 3123 | """ |
| 3124 | si = FormalPowerSeries.__init__ |
| 3125 | si(self,a.parent(),min_index=1) |
| 3126 | self.a=a |
| 3127 | |
| 3128 | def coeffs(j,n): |
| 3129 | """ sage: None #indirect doctest """ |
| 3130 | |
| 3131 | a = j.a |
| 3132 | if n < 0: |
| 3133 | return 0 |
| 3134 | |
| 3135 | if n == 1: |
| 3136 | return 1 |
| 3137 | |
| 3138 | ap = a.diff() |
| 3139 | |
| 3140 | r = a.K(0) |
| 3141 | |
| 3142 | for k in range(1,n): |
| 3143 | r+=ap[n-k]*j[k] |
| 3144 | |
| 3145 | for k in range(1,n): |
| 3146 | r-=j[k]*a.npow(k)[n] |
| 3147 | |
| 3148 | return r/(a[1]**n-a[1]) |
| 3149 | |
| 3150 | class Schroeder(FormalPowerSeries01): |
| 3151 | def __init__(self,a): |
| 3152 | """ |
| 3153 | Description and tests at FormalPowerSeries.schroeder |
| 3154 | sage: None # indirect doctest |
| 3155 | """ |
| 3156 | si = FormalPowerSeries.__init__ |
| 3157 | si(self,a.parent(),min_index=1) |
| 3158 | self.a=a |
| 3159 | |
| 3160 | def coeffs(s,n): |
| 3161 | """ sage: None # indirect doctest """ |
| 3162 | a = s.a |
| 3163 | if decidable0(a.K): |
| 3164 | assert a[1] != 0, a[1] |
| 3165 | assert a[1] != 1, a[1] |
| 3166 | |
| 3167 | q = a[1] |
| 3168 | |
| 3169 | if n <= 0: |
| 3170 | return 0 |
| 3171 | if n == 1: |
| 3172 | return 1 |
| 3173 | return sum([s[m]*a.npow(m)[n] for m in range(1,n)])/(q - q**n) |
| 3174 | |
| 3175 | class InvSchroeder(FormalPowerSeries01): |
| 3176 | def __init__(self,a): |
| 3177 | """ |
| 3178 | Description and tests at FormalPowerSeries.inv_schroeder |
| 3179 | sage: None # indirect doctest |
| 3180 | """ |
| 3181 | si = FormalPowerSeries.__init__ |
| 3182 | si(self,a.parent(),min_index=1) |
| 3183 | self.a=a |
| 3184 | |
| 3185 | def coeffs(s,n): |
| 3186 | """sage: None #indirect doctest""" |
| 3187 | a = s.a |
| 3188 | q = a[1] |
| 3189 | |
| 3190 | if n <= 0: |
| 3191 | return 0 |
| 3192 | if n == 1: |
| 3193 | return 1 |
| 3194 | return sum([a[m]*s.npow(m)[n] for m in range(2,n+1)])/(q**n-q) |
| 3195 | |
| 3196 | class Regit01(FormalPowerSeries01): |
| 3197 | def __init__(self,a,t): |
| 3198 | """ |
| 3199 | Description and tests at FormalPowerSeries01.regit |
| 3200 | sage: None # indirect doctest |
| 3201 | """ |
| 3202 | si = FormalPowerSeries.__init__ |
| 3203 | si(self,a.parent(),min_index=1) |
| 3204 | self.a = a |
| 3205 | self.t = t |
| 3206 | |
| 3207 | def coeffs(self,n): |
| 3208 | """ sage: None # indirect doctest """ |
| 3209 | a = self.a |
| 3210 | t = self.t |
| 3211 | if decidable0(a.K): |
| 3212 | assert a[0] == 0, "The index of the lowest non-zero coefficient must be 1, but is " + repr(a.min_index) |
| 3213 | assert a[1] == 1, "The first coefficient must be 1, but is " + repr(a[1]) |
| 3214 | |
| 3215 | if n == 1: return 1 |
| 3216 | #def c(m): |
| 3217 | # return (-1)**(n-1-m)*binomial(t,m)*binomial(t-1-m,n-1-m) |
| 3218 | #res = sum([c(m)*a.nit(m)[n] for m in range(n)],a.K(0)) |
| 3219 | #return res |
| 3220 | |
| 3221 | r = a.K(0) |
| 3222 | for m in range(n): |
| 3223 | s = a.K(0) |
| 3224 | for k in range(m+1): |
| 3225 | s += binomial(m,k)*(-1)**(m-k)*a.nit(k)[n] |
| 3226 | s *= binomial(t,m) |
| 3227 | r += s |
| 3228 | return r |
| 3229 | |
| 3230 | class Julia01(FormalPowerSeries01): |
| 3231 | def __init__(self,a): |
| 3232 | """ |
| 3233 | Description and tests at FormalPowerSeries01.julia |
| 3234 | sage: None # indirect doctest |
| 3235 | """ |
| 3236 | P = a.parent() |
| 3237 | si = FormalPowerSeries.__init__ |
| 3238 | si(self,P,min_index=1) |
| 3239 | self.a=a |
| 3240 | |
| 3241 | def coeffs(self,n): |
| 3242 | """ sage: None # indirect doctest """ |
| 3243 | a = self.a |
| 3244 | r = a.K(0) |
| 3245 | P = a.parent() |
| 3246 | |
| 3247 | for m in range(n): |
| 3248 | s = a.K(0) |
| 3249 | for k in range(m+1): |
| 3250 | s += binomial(m,k)*(-1)**(m-k)*a.nit(k)[n] |
| 3251 | s *= P.Stirling1[m][1]/factorial(m) |
| 3252 | r += s |
| 3253 | |
| 3254 | return r |
| 3255 | |