tests/rtest_powerseries.mac

   1 /*
   2   test_powerseries.mac
   3
   4   Tests for the powerseries() function
   5 */
   6
   7 kill(all);
   8 done$
   9
  10 /*
  11   Basic identities
  12 */
  13 powerseries(a, x, 0);
  14 a$
  15
  16 powerseries(x, x, 0);
  17 x$
  18
  19 /* Check that we get an error if the variable is a non-symbol atom */
  20 errcatch(powerseries(x, 0, 0));
  21 []$
  22
  23 errcatch(powerseries(x, "foo", 0));
  24 []$
  25
  26 /* Check we distribute over bags */
  27 niceindices(powerseries([sin(x), cos(x)], x, 0));
  28 ''(niceindices([powerseries(sin(x), x, 0), powerseries(cos(x), x, 0)]))$
  29
  30 /* Check we distribute over sum() and product() */
  31 niceindices(powerseries(sum(f(x,i), i, 0, inf), x, 0));
  32 ''(niceindices(sum(powerseries(f(x,i), x, 0), i, 0, inf)));
  33
  34 niceindices(powerseries(product(f(x,i), i, 0, inf), x, 0));
  35 ''(niceindices(product(powerseries(f(x,i), x, 0), i, 0, inf)));
  36
  37
  38 /* Check that we can expand with respect to a non-atomic expression
  39    (or at least check that we do something vaguely sensible) */
  40 niceindices(powerseries(sin(x^2), x^2, 0) - powerseries(sin(x^2), x, 0));
  41 0$
  42
  43 /* Check that sums get expanded sanely (in SP1) */
  44 niceindices(powerseries(a + sin(x), x, 0) - powerseries(sin(x), x, 0));
  45 a$
  46
  47 /* An element of the transformation dictionary for special trig sums in SP1 */
  48 niceindices(powerseries (1+tan(x)^2, x, 0) - powerseries(sec(x)^2, x, 0));
  49 0$
  50
  51 /* Check that constant factors get pulled out (in SP1TIMES) */
  52 niceindices(powerseries(a*sin(x), x, 0) - a*powerseries(sin(x), x, 0));
  53 0$
  54
  55 /* ... even if they aren't atoms. */
  56 niceindices(powerseries((a+b)*sin(x), x, 0) - (a+b)*powerseries(sin(x), x, 0));
  57 0$
  58
  59 /* "Atomic" factors also get pulled out */
  60 niceindices(powerseries(x*sin(x), x, 0) - x*powerseries(sin(x), x, 0));
  61 0$
  62
  63 /* Check that we evaluate powerseries at inf and minf sensibly */
  64 powerseries(1/x, x, inf);
  65 1/x$
  66
  67 powerseries(1/x, x, minf);
  68 1/x$
  69
  70 /*
  71   Exponentiation of a sum is dealt with in SRBINEXPAND, which has to
  72   cope with various possibilities depending on how much we know about
  73   which terms are nonzero.
  74 */
  75 powerseries(x^y^z, x, 0);
  76 x^y^z$
  77
  78 (declare([int, posint], integer),
  79  assume(posint > 0),
  80  assume(pos > 0),
  81  assume(notequal(nz, 0)),
  82  assume(equal(zero, 0)), 0);
  83 0$
  84
  85 powerseries((x^y+1)^zero, x, 0);
  86 1$
  87
  88 powerseries((zero*x + 2)^y, x, 0);
  89 2^y$
  90
  91 powerseries((x + zero)^y, x, 0);
  92 x^y$
  93
  94 niceindices(powerseries((nz + 1)^n, nz, 0));
  95 'sum(nz^i/beta(i+1, n-i+1), i, 0, inf)/(n+1)$
  96
  97 niceindices(powerseries((nz + 1)^posint, nz, 0));
  98 'sum(nz^i/beta(i+1, posint-i+1), i, 0, posint)/(posint+1)$
  99
 100 niceindices(powerseries((nz + 1)^k, nz, 0));
 101 'sum(nz^i/beta(i+1, k-i+1), i, 0, inf)/(k+1)$
 102
 103 niceindices(powerseries((x + 1)^k, x, 0));
 104 1 + 'sum(x^i/beta(i+1, k-i+1), i, 1, inf)/(k+1)$
 105
 106 niceindices(powerseries((x + u)^posint, x, 0));
 107 ('sum(x^i*u^(posint-i)/beta(i+1, posint-i+1),
 108       i, 0, posint-1) / (posint+1)) + x^posint$
 109
 110 powerseries((x + y)^k, x, 0);
 111 'powerseries((x + y)^k, x, 0)$
 112
 113 /* Some more white-box tests, based on the structure of SRATEXPND */
 114
 115 /* If the numerator is monomial, factor it out and recurse */
 116 niceindices(powerseries(a*x^u/(1+x), x, 0) - a*x^u*powerseries(1/(1+x), x, 0));
 117 0$
 118
 119 /* If the denominator is free of x and the numerator is a polynomial,
 120    return the quotient */
 121 powerseries((x^2 + x + 1)/(a + b + c), x, 0);
 122 (x^2 + x + 1)/(a + b + c)$
 123
 124 /*
 125   If the denominator is free of y, but the numerator is sufficiently
 126   gnarly, we give up. (Obviously, it would be nice to do better here!)
 127 */
 128 powerseries((x+x^(1+y)+1)/(2+y), x, 0);
 129 powerseries((x+x^(1+y)+1)/(2+y), x, 0)$
 130
 131 /*
 132   Another gnarly example we can't deal with. We can't use partial
 133   fraction expansion because we don't know that k < 2 (and, even if we
 134   did, I don't know how we'd parametrise the result on k).
 135 */
 136 powerseries((x+1)^k/((x+2)*(x+3)), x, 0);
 137 powerseries((x+1)^k/((x+2)*(x+3)), x, 0)$
 138
 139 /* The only way of getting a monomial denominator in sratexpnd is if
 140    the numerator is a polynomial. If the denominator is a sum, we return the
 141    sum of quotients: */
 142 powerseries((x^2 + x + 1)/(a*x^2), x, 0);
 143 1/a + 1/(a*x) + 1/(a*x^2)$
 144
 145 /* If the denominator isn't a sum, we just return the quotient of the
 146    top and bottom. */
 147 powerseries((1+x)^n/(a*x^2), x, 0);
 148 (1+x)^n/(a*x^2)$
 149
 150 /* (I haven't yet managed to find an example that triggers the
 151     SRATSUBST clause and doesn't yield a noun form) */
 152
 153 /* Partial fraction expansion */
 154 niceindices(powerseries(1/((1-2*x)*(1-3*x)), x, 0));
 155 'sum((3^(i+1) - 2^(i+1))*x^i, i, 0, inf)$
 156
 157 /*
 158   Make sure we give up (rather than possibly recursing infinitely) in
 159   the general case at the bottom of SRATEXPND.
 160 */
 161 powerseries((x+2)*(x+3)^k/(x*(x+4)), x, 0);
 162 powerseries((x+2)*(x+3)^k/(x*(x+4)), x, 0)$
 163
 164 /*
 165   A complicated example of a partial fraction expansion. The expansion
 166   into partial fractions is 1 + 2/(x+3) - 6/(x+4), and the terms then
 167   expand into two infinite series and one constant term. This test
 168   makes sure that we coalesce the two series properly (in PSP2FOLDSUM).
 169 */
 170 factor(niceindices(powerseries((x+1)*(x+2)/((x+3)*(x+4)), x, 0)));
 171 (sum(3^(-i-1)*(-1)^i*(4^(i+1)-3^(i+2))*x^i/4^i,i,0,inf)+2)/2$
 172
 173 /*
 174   Make sure we successfully remove zero roots from the denominator
 175   before trying to expand partial fractions.
 176 */
 177 niceindices(powerseries((x+2)*(x+3)/(x^2*(x+1)^k), x, 0)*x^2 -
 178             powerseries((x+2)*(x+3)/(x+1)^k, x, 0));
 179 0$
 180
 181 /* Fixed in series.lisp rev 1.5, 04 Feb 2004.
 182    First case returned second result */
 183 niceindices(powerseries(1/sqrt(1 + x), x, 0));
 184 1 + 2*'sum(x^i/beta(1/2-i,i+1),i,1,inf);
 185 niceindices (powerseries(sqrt(1+x),x,0));
 186 1 + 2*('sum(x^i/beta(3/2-i,i+1),i,1,inf))/3;
 187
 188 /* sf bug 1730044 - powerseries(1+x^n,x,0) wrong; see also #2764 power series of 1 + x^n */
 189 powerseries(1+x^n, x, 0);
 190 1+x^n$
 191
 192 (gensumnum : 0, declare(n, integer), powerseries(1/(1+x^n), x, 0));
 193 'sum((-1)^i1*x^(i1*n),i1,0,inf)$
 194
 195 /* #2756 powerseries of constant plus power of linear  */
 196 (gensumnum : 0, powerseries(1 + (1-x)^(1/3),x,0));
 197 (3*sum(((-1)^i1*x^i1)/beta(4/3-i1,i1+1),i1,1,inf))/4+2$
 198
 199 /* Examples from SF bug report [ 1722156 ] powerseries((x+1)/(1-2*x^2), x, 0);
 200  * tnx Dan Gildea
 201  */
 202
 203 gensumnum : 0;
 204 0;
 205
 206 /* two simple roots */
 207 powerseries((1)/((1-2*x)*(1-3*x)), x, 0);
 208 'sum((3^(i1+1)+(-2*2^i1))*x^i1,i1,0,inf);
 209
 210 /* fibonacci */
 211 powerseries((1)/(1-x-x^2), x, 0);
 212 -'sum((-2*(sqrt(5)-1)^(-i2-1)*2^i2/sqrt(5)
 213  -2*(sqrt(5)+1)^(-i2-1)*(-2)^i2/sqrt(5))
 214  *x^i2,i2,0,inf);
 215
 216 /* 1 1 2 2 4 4 8 8 */
 217 factor(powerseries((1+x)/(1-2*x^2), x, 0));
 218 /*
 219 'sum( (-1*(1/-(1/sqrt(2)))*(1/sqrt(2)-1)*(1/-(4/sqrt(2)))*(1/-(1/sqrt(2)))^i3 +
 220       -1*(1/(1/sqrt(2)))*(-1/sqrt(2)-1)*(1/(4/sqrt(2)))*(1/(1/sqrt(2)))^i3 ) * x^i3,i3,0,inf);
 221 */
 222
 223 'sum(2^(i3/2-3/2)*((sqrt(2)-1)*(-1)^i3+sqrt(2)+1)*x^i3,i3,0,inf);
 224
 225 /* multiple root */
 226 powerseries((1+x)/(1-x)^2, x, 0);
 227 'sum(2*(i4+1)*x^i4-x^i4,i4,0,inf);
 228
 229 /* numerator higher order poly than denom */
 230 powerseries((1+x^3)/(1-x-x^2), x, 0);
 231 -2*x
 232  *'sum((-2*(sqrt(5)-1)^(-i5-1)*2^i5/sqrt(5)
 233    -2*(sqrt(5)+1)^(-i5-1)*(-2)^i5/sqrt(5))
 234    *x^i5,i5,0,inf)
 235   -x+1;
 236
 237 /* zero root in denom */
 238 powerseries((1)/((1-2*x)*(x)), x, 0);
 239 ('sum(2^i6*x^i6,i6,0,inf))/x;
 240
 241 /* one simple and one repeated root in denom */
 242 powerseries((1+x+x^2)/((1-2*x)*(1+x)^2), x, 0);
 243 'sum(7*2^i7*x^i7/9+(i7+1)*(-1)^i7*x^i7/3-(-1)^i7*x^i7/9,i7,0,inf);
 244
 245 /* gcd of exps is two */
 246 powerseries((1-x^2)/(1-4*x^2+x^4), x, 0);
 247 'sum((-1*(1/(2-sqrt(3)))*(sqrt(3)-1)*(1/(2*(2-sqrt(3))-4))*(1/(2-sqrt(3)))^i8 +
 248       -1*(1/(sqrt(3)+2))*(-sqrt(3)-1)*(1/(2*(sqrt(3)+2)-4))*(1/(sqrt(3)+2))^i8 )*x^(2*i8),
 249       i8, 0, inf);
 250
 251 /* #2750 powerseries(x^x,x,0) gives Lisp error */
 252 powerseries(x^x, x, 0);
 253 'powerseries(x^x, x, 0)$
 254
 255 sumcontract(intosum(powerseries(1+ (1-x)^(a),x,0) - powerseries((1-x)^(a),x,0)));
 256 1$
 257
 258 /*
 259   #2775: Don't expand a powerseries with log(ab)=log(a)+log(b) if a
 260          is negative
 261
 262   (and also #2767)
 263 */
 264 (gensumnum : 0, powerseries (log(2-x), x, 0));
 265 log(2) - sum (2^(-i2-1)*x^(i2+1)/(i2+1), i2, 0, inf)$
 266
 267 /*
 268   #2760: powerseries at minf
 269
 270   For an analytic function like this, the power series at minf should
 271   match the power series at inf.
 272 */
 273 block([inf_exp, minf_exp],
 274   gensumnum : 0,
 275   inf_exp: powerseries (1/(1+x), x, inf),
 276   gensumnum : 0,
 277   minf_exp: powerseries (1/(1+x), x, inf),
 278   inf_exp - minf_exp);
 279 0$
 280
 281 /* #2765: powerseries with non-integer order */
 282 powerseries (diff (f(x), x, biggles), x, 0);
 283 'powerseries (diff (f(x), x, biggles), x, 0)$
 284
 285 /* #2755: powerseries of natural exponential */
 286 niceindices (powerseries (%e^x, x, 1));
 287 sum (1/i!, i, 0, inf) * sum((x-1)^i/i!, i, 0, inf)$
 288
 289 /*
 290   #2760
 291
 292   Make sure that we give up rather than trying to substitute an
 293   arbitrary expression for an integration / differentiation variable
 294   after expansion.
 295 */
 296 powerseries(f(x), x, inf);
 297 'powerseries(f(x), x, inf)$
 298
 299 /*
 300   The test above gives up when it sees either an integral or a
 301   derivative wrt the expansion variable as it's trying to substitute
 302   answers back in. Make sure that we still allow either with respect
 303   to some other variable. (Expand somewhere other than zero, otherwise
 304   we skip the substitution step)
 305 */
 306 factor(powerseries(integrate(f(x), x)*y, y, 1));
 307 'integrate(f(x), x)*y$
 308
 309 niceindices(powerseries(log(sin(x)/x),x,0));
 310 /*('sum((-1)^i1*2^(2*i1)*bern(2*i1)*x^(2*i1)/(i1*(2*i1)!),i1,1,inf))/2$*/
 311 'sum((-1)^i*2^(2*i-1)*bern(2*i)*x^(2*i)/(i*(2*i)!),i,1,inf);
 312
 313 /*
 314   Check that we don't substitute blindly into the bound variables in
 315   at() expressions.
 316
 317   The powerseries expands the integral of f(x) by computing an
 318   antiderivative of the power series of f(x). Since we don't know
 319   anything about f, this contains terms like "at(diff(f(x), x, n),
 320   x=0)". If we then substituted the endpoints blindly, we'd get an
 321   error from substituting in a number as a differentiation variable.
 322 */
 323 op(niceindices(powerseries(integrate(f(x), x, 0, y), y, 0)));
 324 ''(nounify('sum))$
 325
 326 /*
 327   Check we get a noun form (rather than something exploding) when
 328   expanding a known power series and then trying to substitute in some
 329   arbitrary function. At the moment, we only support arguments that
 330   are monomials (see SP2SUB in series.lisp).
 331 */
 332 powerseries(sin(f(x)), x, 0);
 333 powerseries(sin(f(x)), x, 0)$
 334
 335 /* Basic support for expanding log(1+x) */
 336 factor(sumcontract(intosum(powerseries(log(x+1), x, 0) -
 337                            sum((-1)^(i+1)*x^i/i, i, 1, inf))));
 338 0$
 339
 340 /* Check that we can expand log(1+a*x^k) for various a and k. */
 341 niceindices(powerseries(log(1+a*x), x, 0));
 342 -sum((-1)^i * a^i * x^i / i, i, 1, inf)$
 343
 344 niceindices(powerseries(log(1-a*x), x, 0));
 345 -sum(a^i * x^i / i, i, 1, inf)$
 346
 347 /*
 348   A particular example of a log expansion that didn't work when you
 349   only allowed contents that SMONO accepted
 350 */
 351 niceindices(powerseries(log(1-z*exp(-4)), z, 0));
 352 -sum(exp(-4*i)*z^i/i, i, 1, inf)$
 353
 354 (kill(all), 0);
 355 0$