Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagesmc
Path: blob/master/src/sage/calculus/wester.py
8815 views
1
r"""
2
Further examples from Wester's paper
3
4
These are all the problems at
5
http://yacas.sourceforge.net/essaysmanual.html
6
7
They come from the 1994 paper "Review of CAS mathematical
8
capabilities", by Michael Wester, who put forward 123 problems that
9
a reasonable computer algebra system should be able to solve and
10
tested the then current versions of various commercial CAS on this
11
list. Sage can do most of the problems natively now, i.e., with no
12
explicit calls to Maxima or other systems.
13
14
::
15
16
sage: # (YES) factorial of 50, and factor it
17
sage: factorial(50)
18
30414093201713378043612608166064768844377641568960512000000000000
19
sage: factor(factorial(50))
20
2^47 * 3^22 * 5^12 * 7^8 * 11^4 * 13^3 * 17^2 * 19^2 * 23^2 * 29 * 31 * 37 * 41 * 43 * 47
21
22
::
23
24
sage: # (YES) 1/2+...+1/10 = 4861/2520
25
sage: sum(1/n for n in range(2,10+1)) == 4861/2520
26
True
27
28
::
29
30
sage: # (YES) Evaluate e^(Pi*Sqrt(163)) to 50 decimal digits
31
sage: a = e^(pi*sqrt(163)); a
32
e^(sqrt(163)*pi)
33
sage: print RealField(150)(a)
34
2.6253741264076874399999999999925007259719820e17
35
36
::
37
38
sage: # (YES) Evaluate the Bessel function J[2] numerically at z=1+I.
39
sage: bessel_J(2, 1+I).n()
40
0.0415798869439621 + 0.247397641513306*I
41
42
::
43
44
sage: # (YES) Obtain period of decimal fraction 1/7=0.(142857).
45
sage: a = 1/7
46
sage: print a
47
1/7
48
sage: print a.period()
49
6
50
51
::
52
53
sage: # (YES) Continued fraction of 3.1415926535
54
sage: a = 3.1415926535
55
sage: continued_fraction(a)
56
[3, 7, 15, 1, 292, 1, 1, 6, 2, 13, 4]
57
58
::
59
60
sage: # (YES) Sqrt(2*Sqrt(3)+4)=1+Sqrt(3).
61
sage: # The Maxima backend equality checker does this;
62
sage: # note the equality only holds for one choice of sign,
63
sage: # but Maxima always chooses the "positive" one
64
sage: a = sqrt(2*sqrt(3) + 4); b = 1 + sqrt(3)
65
sage: print float(a-b)
66
0.0
67
sage: print bool(a == b)
68
True
69
sage: # We can, of course, do this in a quadratic field
70
sage: k.<sqrt3> = QuadraticField(3)
71
sage: asqr = 2*sqrt3 + 4
72
sage: b = 1+sqrt3
73
sage: asqr == b^2
74
True
75
76
::
77
78
sage: # (NOT REALLY) Sqrt(14+3*Sqrt(3+2*Sqrt(5-12*Sqrt(3-2*Sqrt(2)))))=3+Sqrt(2).
79
sage: a = sqrt(14+3*sqrt(3+2*sqrt(5-12*sqrt(3-2*sqrt(2)))))
80
sage: b = 3+sqrt(2)
81
sage: a, b
82
(sqrt(3*sqrt(2*sqrt(-12*sqrt(-2*sqrt(2) + 3) + 5) + 3) + 14), sqrt(2) + 3)
83
sage: bool(a==b)
84
False
85
sage: abs(float(a-b)) < 1e-10
86
True
87
sage: # 2*Infinity-3=Infinity.
88
sage: 2*infinity-3 == infinity
89
True
90
91
::
92
93
sage: # (YES) Standard deviation of the sample (1, 2, 3, 4, 5).
94
sage: v = vector(RDF, 5, [1,2,3,4,5])
95
sage: v.standard_deviation()
96
1.58113883008
97
98
::
99
100
sage: # (NO) Hypothesis testing with t-distribution.
101
sage: # (NO) Hypothesis testing with chi^2 distribution
102
sage: # (But both are included in Scipy and R)
103
104
::
105
106
sage: # (YES) (x^2-4)/(x^2+4*x+4)=(x-2)/(x+2).
107
sage: R.<x> = QQ[]
108
sage: (x^2-4)/(x^2+4*x+4) == (x-2)/(x+2)
109
True
110
sage: restore('x')
111
112
::
113
114
sage: # (YES -- Maxima doesn't immediately consider them
115
sage: # equal, but simplification shows that they are)
116
sage: # (Exp(x)-1)/(Exp(x/2)+1)=Exp(x/2)-1.
117
sage: f = (exp(x)-1)/(exp(x/2)+1)
118
sage: g = exp(x/2)-1
119
sage: f
120
(e^x - 1)/(e^(1/2*x) + 1)
121
sage: g
122
e^(1/2*x) - 1
123
sage: f.simplify_radical()
124
e^(1/2*x) - 1
125
sage: g
126
e^(1/2*x) - 1
127
sage: f(x=10.0).n(53), g(x=10.0).n(53)
128
(147.413159102577, 147.413159102577)
129
sage: bool(f == g)
130
True
131
132
::
133
134
sage: # (YES) Expand (1+x)^20, take derivative and factorize.
135
sage: # first do it using algebraic polys
136
sage: R.<x> = QQ[]
137
sage: f = (1+x)^20; f
138
x^20 + 20*x^19 + 190*x^18 + 1140*x^17 + 4845*x^16 + 15504*x^15 + 38760*x^14 + 77520*x^13 + 125970*x^12 + 167960*x^11 + 184756*x^10 + 167960*x^9 + 125970*x^8 + 77520*x^7 + 38760*x^6 + 15504*x^5 + 4845*x^4 + 1140*x^3 + 190*x^2 + 20*x + 1
139
sage: deriv = f.derivative()
140
sage: deriv
141
20*x^19 + 380*x^18 + 3420*x^17 + 19380*x^16 + 77520*x^15 + 232560*x^14 + 542640*x^13 + 1007760*x^12 + 1511640*x^11 + 1847560*x^10 + 1847560*x^9 + 1511640*x^8 + 1007760*x^7 + 542640*x^6 + 232560*x^5 + 77520*x^4 + 19380*x^3 + 3420*x^2 + 380*x + 20
142
sage: deriv.factor()
143
(20) * (x + 1)^19
144
sage: restore('x')
145
sage: # next do it symbolically
146
sage: var('y')
147
y
148
sage: f = (1+y)^20; f
149
(y + 1)^20
150
sage: g = f.expand(); g
151
y^20 + 20*y^19 + 190*y^18 + 1140*y^17 + 4845*y^16 + 15504*y^15 + 38760*y^14 + 77520*y^13 + 125970*y^12 + 167960*y^11 + 184756*y^10 + 167960*y^9 + 125970*y^8 + 77520*y^7 + 38760*y^6 + 15504*y^5 + 4845*y^4 + 1140*y^3 + 190*y^2 + 20*y + 1
152
sage: deriv = g.derivative(); deriv
153
20*y^19 + 380*y^18 + 3420*y^17 + 19380*y^16 + 77520*y^15 + 232560*y^14 + 542640*y^13 + 1007760*y^12 + 1511640*y^11 + 1847560*y^10 + 1847560*y^9 + 1511640*y^8 + 1007760*y^7 + 542640*y^6 + 232560*y^5 + 77520*y^4 + 19380*y^3 + 3420*y^2 + 380*y + 20
154
sage: deriv.factor()
155
20*(y + 1)^19
156
157
::
158
159
sage: # (YES) Factorize x^100-1.
160
sage: factor(x^100-1)
161
(x^40 - x^30 + x^20 - x^10 + 1)*(x^20 + x^15 + x^10 + x^5 + 1)*(x^20 - x^15 + x^10 - x^5 + 1)*(x^8 - x^6 + x^4 - x^2 + 1)*(x^4 + x^3 + x^2 + x + 1)*(x^4 - x^3 + x^2 - x + 1)*(x^2 + 1)*(x + 1)*(x - 1)
162
sage: # Also, algebraically
163
sage: x = polygen(QQ)
164
sage: factor(x^100 - 1)
165
(x - 1) * (x + 1) * (x^2 + 1) * (x^4 - x^3 + x^2 - x + 1) * (x^4 + x^3 + x^2 + x + 1) * (x^8 - x^6 + x^4 - x^2 + 1) * (x^20 - x^15 + x^10 - x^5 + 1) * (x^20 + x^15 + x^10 + x^5 + 1) * (x^40 - x^30 + x^20 - x^10 + 1)
166
sage: restore('x')
167
168
::
169
170
sage: # (YES) Factorize x^4-3*x^2+1 in the field of rational numbers extended by roots of x^2-x-1.
171
sage: k.< a> = NumberField(x^2 - x -1)
172
sage: R.< y> = k[]
173
sage: f = y^4 - 3*y^2 + 1
174
sage: f
175
y^4 - 3*y^2 + 1
176
sage: factor(f)
177
(y - a) * (y - a + 1) * (y + a - 1) * (y + a)
178
179
::
180
181
sage: # (YES) Factorize x^4-3*x^2+1 mod 5.
182
sage: k.< x > = GF(5) [ ]
183
sage: f = x^4 - 3*x^2 + 1
184
sage: f.factor()
185
(x + 2)^2 * (x + 3)^2
186
sage: # Alternatively, from symbol x as follows:
187
sage: reset('x')
188
sage: f = x^4 - 3*x^2 + 1
189
sage: f.polynomial(GF(5)).factor()
190
(x + 2)^2 * (x + 3)^2
191
192
::
193
194
sage: # (YES) Partial fraction decomposition of (x^2+2*x+3)/(x^3+4*x^2+5*x+2)
195
sage: f = (x^2+2*x+3)/(x^3+4*x^2+5*x+2); f
196
(x^2 + 2*x + 3)/(x^3 + 4*x^2 + 5*x + 2)
197
sage: f.partial_fraction()
198
3/(x + 2) - 2/(x + 1) + 2/(x + 1)^2
199
200
::
201
202
sage: # (YES) Assuming x>=y, y>=z, z>=x, deduce x=z.
203
sage: forget()
204
sage: var('x,y,z')
205
(x, y, z)
206
sage: assume(x>=y, y>=z,z>=x)
207
sage: print bool(x==z)
208
True
209
210
::
211
212
sage: # (YES) Assuming x>y, y>0, deduce 2*x^2>2*y^2.
213
sage: forget()
214
sage: assume(x>y, y>0)
215
sage: print list(sorted(assumptions()))
216
[x > y, y > 0]
217
sage: print bool(2*x^2 > 2*y^2)
218
True
219
sage: forget()
220
sage: print assumptions()
221
[]
222
223
::
224
225
sage: # (NO) Solve the inequality Abs(x-1)>2.
226
sage: # Maxima doesn't solve inequalities
227
sage: # (but some Maxima packages do):
228
sage: eqn = abs(x-1) > 2
229
sage: print eqn
230
abs(x - 1) > 2
231
232
::
233
234
sage: # (NO) Solve the inequality (x-1)*...*(x-5)<0.
235
sage: eqn = prod(x-i for i in range(1,5 +1)) < 0
236
sage: # but don't know how to solve
237
sage: eqn
238
(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5) < 0
239
240
::
241
242
sage: # (YES) Cos(3*x)/Cos(x)=Cos(x)^2-3*Sin(x)^2 or similar equivalent combination.
243
sage: f = cos(3*x)/cos(x)
244
sage: g = cos(x)^2 - 3*sin(x)^2
245
sage: h = f-g
246
sage: print h.trig_simplify()
247
0
248
249
::
250
251
sage: # (YES) Cos(3*x)/Cos(x)=2*Cos(2*x)-1.
252
sage: f = cos(3*x)/cos(x)
253
sage: g = 2*cos(2*x) - 1
254
sage: h = f-g
255
sage: print h.trig_simplify()
256
0
257
258
::
259
260
sage: # (GOOD ENOUGH) Define rewrite rules to match Cos(3*x)/Cos(x)=Cos(x)^2-3*Sin(x)^2.
261
sage: # Sage has no notion of "rewrite rules", but
262
sage: # it can simplify both to the same thing.
263
sage: (cos(3*x)/cos(x)).simplify_full()
264
4*cos(x)^2 - 3
265
sage: (cos(x)^2-3*sin(x)^2).simplify_full()
266
4*cos(x)^2 - 3
267
268
::
269
270
sage: # (YES) Sqrt(997)-(997^3)^(1/6)=0
271
sage: a = sqrt(997) - (997^3)^(1/6)
272
sage: a.simplify()
273
0
274
sage: bool(a == 0)
275
True
276
277
::
278
279
sage: # (YES) Sqrt(99983)-99983^3^(1/6)=0
280
sage: a = sqrt(99983) - (99983^3)^(1/6)
281
sage: bool(a==0)
282
True
283
sage: float(a)
284
1.1368683772...e-13
285
sage: print 13*7691
286
99983
287
288
::
289
290
sage: # (YES) (2^(1/3) + 4^(1/3))^3 - 6*(2^(1/3) + 4^(1/3))-6 = 0
291
sage: a = (2^(1/3) + 4^(1/3))^3 - 6*(2^(1/3) + 4^(1/3)) - 6; a
292
(4^(1/3) + 2^(1/3))^3 - 6*4^(1/3) - 6*2^(1/3) - 6
293
sage: bool(a==0)
294
True
295
sage: abs(float(a)) < 1e-10
296
True
297
sage: ## or we can do it using number fields.
298
sage: reset('x')
299
sage: k.<b> = NumberField(x^3-2)
300
sage: a = (b + b^2)^3 - 6*(b + b^2) - 6
301
sage: print a
302
0
303
304
::
305
306
sage: # (NO, except numerically) Ln(Tan(x/2+Pi/4))-ArcSinh(Tan(x))=0
307
# Sage uses the Maxima convention when comparing symbolic expressions and
308
# returns True only when it can prove equality. Thus, in this case, we get
309
# False even though the equality holds.
310
sage: f = log(tan(x/2 + pi/4)) - arcsinh(tan(x))
311
sage: bool(f == 0)
312
False
313
sage: [abs(float(f(x=i/10))) < 1e-15 for i in range(1,5)]
314
[True, True, True, True]
315
sage: # Numerically, the expression Ln(Tan(x/2+Pi/4))-ArcSinh(Tan(x))=0 and its derivative at x=0 are zero.
316
sage: g = f.derivative()
317
sage: abs(float(f(x=0))) < 1e-10
318
True
319
sage: abs(float(g(x=0))) < 1e-10
320
True
321
sage: g
322
-sqrt(tan(x)^2 + 1) + 1/2*(tan(1/4*pi + 1/2*x)^2 + 1)/tan(1/4*pi + 1/2*x)
323
324
::
325
326
sage: # (NO) Ln((2*Sqrt(r) + 1)/Sqrt(4*r 4*Sqrt(r) 1))=0.
327
sage: var('r')
328
r
329
sage: f = log( (2*sqrt(r) + 1) / sqrt(4*r + 4*sqrt(r) + 1))
330
sage: f
331
log((2*sqrt(r) + 1)/sqrt(4*r + 4*sqrt(r) + 1))
332
sage: bool(f == 0)
333
False
334
sage: [abs(float(f(r=i))) < 1e-10 for i in [0.1,0.3,0.5]]
335
[True, True, True]
336
337
::
338
339
sage: # (NO)
340
sage: # (4*r+4*Sqrt(r)+1)^(Sqrt(r)/(2*Sqrt(r)+1))*(2*Sqrt(r)+1)^(2*Sqrt(r)+1)^(-1)-2*Sqrt(r)-1=0, assuming r>0.
341
sage: assume(r>0)
342
sage: f = (4*r+4*sqrt(r)+1)^(sqrt(r)/(2*sqrt(r)+1))*(2*sqrt(r)+1)^(2*sqrt(r)+1)^(-1)-2*sqrt(r)-1
343
sage: f
344
(4*r + 4*sqrt(r) + 1)^(sqrt(r)/(2*sqrt(r) + 1))*(2*sqrt(r) + 1)^(1/(2*sqrt(r) + 1)) - 2*sqrt(r) - 1
345
sage: bool(f == 0)
346
False
347
sage: [abs(float(f(r=i))) < 1e-10 for i in [0.1,0.3,0.5]]
348
[True, True, True]
349
350
::
351
352
sage: # (YES) Obtain real and imaginary parts of Ln(3+4*I).
353
sage: a = log(3+4*I); a
354
log(4*I + 3)
355
sage: a.real()
356
log(5)
357
sage: a.imag()
358
arctan(4/3)
359
360
::
361
362
sage: # (YES) Obtain real and imaginary parts of Tan(x+I*y)
363
sage: z = var('z')
364
sage: a = tan(z); a
365
tan(z)
366
sage: a.real()
367
tan(real_part(z))/(tan(imag_part(z))^2*tan(real_part(z))^2 + 1)
368
sage: a.imag()
369
tanh(imag_part(z))/(tan(imag_part(z))^2*tan(real_part(z))^2 + 1)
370
371
372
::
373
374
sage: # (YES) Simplify Ln(Exp(z)) to z for -Pi<Im(z)<=Pi.
375
sage: # Unfortunately (?), Maxima does this even without
376
sage: # any assumptions.
377
sage: # We *would* use assume(-pi < imag(z))
378
sage: # and assume(imag(z) <= pi)
379
sage: f = log(exp(z)); f
380
log(e^z)
381
sage: f.simplify()
382
z
383
sage: forget()
384
385
::
386
387
sage: # (YES) Assuming Re(x)>0, Re(y)>0, deduce x^(1/n)*y^(1/n)-(x*y)^(1/n)=0.
388
sage: # Maxima 5.26 has different behaviours depending on the current
389
sage: # domain.
390
sage: # To stick with the behaviour of previous versions, the domain is set
391
sage: # to 'real' in the following.
392
sage: # See Trac #10682 for further details.
393
sage: n = var('n')
394
sage: f = x^(1/n)*y^(1/n)-(x*y)^(1/n)
395
sage: assume(real(x) > 0, real(y) > 0)
396
sage: f.simplify()
397
x^(1/n)*y^(1/n) - (x*y)^(1/n)
398
sage: maxima = sage.calculus.calculus.maxima
399
sage: maxima.set('domain', 'real') # set domain to real
400
sage: f.simplify()
401
0
402
sage: maxima.set('domain', 'complex') # set domain back to its default value
403
sage: forget()
404
405
::
406
407
sage: # (YES) Transform equations, (x==2)/2+(1==1)=>x/2+1==2.
408
sage: eq1 = x == 2
409
sage: eq2 = SR(1) == SR(1)
410
sage: eq1/2 + eq2
411
1/2*x + 1 == 2
412
413
::
414
415
sage: # (SOMEWHAT) Solve Exp(x)=1 and get all solutions.
416
sage: # to_poly_solve in Maxima can do this.
417
sage: solve(exp(x) == 1, x)
418
[x == 0]
419
420
::
421
422
sage: # (SOMEWHAT) Solve Tan(x)=1 and get all solutions.
423
sage: # to_poly_solve in Maxima can do this.
424
sage: solve(tan(x) == 1, x)
425
[x == 1/4*pi]
426
427
::
428
429
sage: # (YES) Solve a degenerate 3x3 linear system.
430
sage: # x+y+z==6,2*x+y+2*z==10,x+3*y+z==10
431
sage: # First symbolically:
432
sage: solve([x+y+z==6, 2*x+y+2*z==10, x+3*y+z==10], x,y,z)
433
[[x == -r1 + 4, y == 2, z == r1]]
434
435
::
436
437
sage: # (YES) Invert a 2x2 symbolic matrix.
438
sage: # [[a,b],[1,a*b]]
439
sage: # Using multivariate poly ring -- much nicer
440
sage: R.<a,b> = QQ[]
441
sage: m = matrix(2,2,[a,b, 1, a*b])
442
sage: zz = m^(-1)
443
sage: print zz
444
[ a/(a^2 - 1) (-1)/(a^2 - 1)]
445
[(-1)/(a^2*b - b) a/(a^2*b - b)]
446
447
::
448
449
sage: # (YES) Compute and factor the determinant of the 4x4 Vandermonde matrix in a, b, c, d.
450
sage: var('a,b,c,d')
451
(a, b, c, d)
452
sage: m = matrix(SR, 4, 4, [[z^i for i in range(4)] for z in [a,b,c,d]])
453
sage: print m
454
[ 1 a a^2 a^3]
455
[ 1 b b^2 b^3]
456
[ 1 c c^2 c^3]
457
[ 1 d d^2 d^3]
458
sage: d = m.determinant()
459
sage: d.factor()
460
(a - b)*(a - c)*(a - d)*(b - c)*(b - d)*(c - d)
461
462
::
463
464
sage: # (YES) Compute and factor the determinant of the 4x4 Vandermonde matrix in a, b, c, d.
465
sage: # Do it instead in a multivariate ring
466
sage: R.<a,b,c,d> = QQ[]
467
sage: m = matrix(R, 4, 4, [[z^i for i in range(4)] for z in [a,b,c,d]])
468
sage: print m
469
[ 1 a a^2 a^3]
470
[ 1 b b^2 b^3]
471
[ 1 c c^2 c^3]
472
[ 1 d d^2 d^3]
473
sage: d = m.determinant()
474
sage: print d
475
a^3*b^2*c - a^2*b^3*c - a^3*b*c^2 + a*b^3*c^2 + a^2*b*c^3 - a*b^2*c^3 - a^3*b^2*d + a^2*b^3*d + a^3*c^2*d - b^3*c^2*d - a^2*c^3*d + b^2*c^3*d + a^3*b*d^2 - a*b^3*d^2 - a^3*c*d^2 + b^3*c*d^2 + a*c^3*d^2 - b*c^3*d^2 - a^2*b*d^3 + a*b^2*d^3 + a^2*c*d^3 - b^2*c*d^3 - a*c^2*d^3 + b*c^2*d^3
476
sage: print d.factor()
477
(-1) * (c - d) * (-b + c) * (b - d) * (-a + c) * (-a + b) * (a - d)
478
479
::
480
481
sage: # (YES) Find the eigenvalues of a 3x3 integer matrix.
482
sage: m = matrix(QQ, 3, [5,-3,-7, -2,1,2, 2,-3,-4])
483
sage: m.eigenspaces_left()
484
[
485
(3, Vector space of degree 3 and dimension 1 over Rational Field
486
User basis matrix:
487
[ 1 0 -1]),
488
(1, Vector space of degree 3 and dimension 1 over Rational Field
489
User basis matrix:
490
[ 1 1 -1]),
491
(-2, Vector space of degree 3 and dimension 1 over Rational Field
492
User basis matrix:
493
[0 1 1])
494
]
495
496
::
497
498
sage: # (YES) Verify some standard limits found by L'Hopital's rule:
499
sage: # Verify(Limit(x,Infinity) (1+1/x)^x, Exp(1));
500
sage: # Verify(Limit(x,0) (1-Cos(x))/x^2, 1/2);
501
sage: limit( (1+1/x)^x, x = oo)
502
e
503
sage: limit( (1-cos(x))/(x^2), x = 1/2)
504
-4*cos(1/2) + 4
505
506
::
507
508
sage: # (OK-ish) D(x)Abs(x)
509
sage: # Verify(D(x) Abs(x), Sign(x));
510
sage: diff(abs(x))
511
x/abs(x)
512
513
::
514
515
sage: # (YES) (Integrate(x)Abs(x))=Abs(x)*x/2
516
sage: integral(abs(x), x)
517
1/2*x*abs(x)
518
519
::
520
521
sage: # (YES) Compute derivative of Abs(x), piecewise defined.
522
sage: # Verify(D(x)if(x<0) (-x) else x,
523
sage: # Simplify(if(x<0) -1 else 1))
524
Piecewise defined function with 2 parts, [[(-10, 0), -1], [(0, 10), 1]]
525
sage: # (NOT really) Integrate Abs(x), piecewise defined.
526
sage: # Verify(Simplify(Integrate(x)
527
sage: # if(x<0) (-x) else x),
528
sage: # Simplify(if(x<0) (-x^2/2) else x^2/2));
529
sage: f = piecewise([ [[-10,0], -x], [[0,10], x]])
530
sage: f.integral(definite=True)
531
100
532
533
::
534
535
sage: # (YES) Taylor series of 1/Sqrt(1-v^2/c^2) at v=0.
536
sage: var('v,c')
537
(v, c)
538
sage: taylor(1/sqrt(1-v^2/c^2), v, 0, 7)
539
1/2*v^2/c^2 + 3/8*v^4/c^4 + 5/16*v^6/c^6 + 1
540
541
::
542
543
sage: # (OK-ish) (Taylor expansion of Sin(x))/(Taylor expansion of Cos(x)) = (Taylor expansion of Tan(x)).
544
sage: # TestYacas(Taylor(x,0,5)(Taylor(x,0,5)Sin(x))/
545
sage: # (Taylor(x,0,5)Cos(x)), Taylor(x,0,5)Tan(x));
546
sage: f = taylor(sin(x), x, 0, 8)
547
sage: g = taylor(cos(x), x, 0, 8)
548
sage: h = taylor(tan(x), x, 0, 8)
549
sage: f = f.power_series(QQ)
550
sage: g = g.power_series(QQ)
551
sage: h = h.power_series(QQ)
552
sage: f - g*h
553
O(x^8)
554
555
::
556
557
sage: # (YES) Taylor expansion of Ln(x)^a*Exp(-b*x) at x=1.
558
sage: a,b = var('a,b')
559
sage: taylor(log(x)^a*exp(-b*x), x, 1, 3)
560
-1/48*(a^3*(x - 1)^a + a^2*(6*b + 5)*(x - 1)^a + 8*b^3*(x - 1)^a + 2*(6*b^2 + 5*b + 3)*a*(x - 1)^a)*(x - 1)^3*e^(-b) + 1/24*(3*a^2*(x - 1)^a + a*(12*b + 5)*(x - 1)^a + 12*b^2*(x - 1)^a)*(x - 1)^2*e^(-b) - 1/2*(a*(x - 1)^a + 2*b*(x - 1)^a)*(x - 1)*e^(-b) + (x - 1)^a*e^(-b)
561
562
::
563
564
sage: # (YES) Taylor expansion of Ln(Sin(x)/x) at x=0.
565
sage: taylor(log(sin(x)/x), x, 0, 10)
566
-1/467775*x^10 - 1/37800*x^8 - 1/2835*x^6 - 1/180*x^4 - 1/6*x^2
567
568
::
569
570
sage: # (NO) Compute n-th term of the Taylor series of Ln(Sin(x)/x) at x=0.
571
sage: # need formal functions
572
573
::
574
575
sage: # (NO) Compute n-th term of the Taylor series of Exp(-x)*Sin(x) at x=0.
576
sage: # (Sort of, with some work)
577
sage: # Solve x=Sin(y)+Cos(y) for y as Taylor series in x at x=1.
578
sage: # TestYacas(InverseTaylor(y,0,4) Sin(y)+Cos(y),
579
sage: # (y-1)+(y-1)^2/2+2*(y-1)^3/3+(y-1)^4);
580
sage: # Note that InverseTaylor does not give the series in terms of x but in terms of y which is semantically
581
sage: # wrong. But other CAS do the same.
582
sage: f = sin(y) + cos(y)
583
sage: g = f.taylor(y, 0, 10)
584
sage: h = g.power_series(QQ)
585
sage: k = (h - 1).reversion()
586
sage: print k
587
y + 1/2*y^2 + 2/3*y^3 + y^4 + 17/10*y^5 + 37/12*y^6 + 41/7*y^7 + 23/2*y^8 + 1667/72*y^9 + 3803/80*y^10 + O(y^11)
588
589
::
590
591
sage: # (OK) Compute Legendre polynomials directly from Rodrigues's formula, P[n]=1/(2^n*n!) *(Deriv(x,n)(x^2-1)^n).
592
sage: # P(n,x) := Simplify( 1/(2*n)!! *
593
sage: # Deriv(x,n) (x^2-1)^n );
594
sage: # TestYacas(P(4,x), (35*x^4)/8+(-15*x^2)/4+3/8);
595
sage: P = lambda n, x: simplify(diff((x^2-1)^n,x,n) / (2^n * factorial(n)))
596
sage: P(4,x).expand()
597
35/8*x^4 - 15/4*x^2 + 3/8
598
599
::
600
601
sage: # (YES) Define the polynomial p=Sum(i,1,5,a[i]*x^i).
602
sage: # symbolically
603
sage: ps = sum(var('a%s'%i)*x^i for i in range(1,6)); ps
604
a5*x^5 + a4*x^4 + a3*x^3 + a2*x^2 + a1*x
605
sage: ps.parent()
606
Symbolic Ring
607
sage: # algebraically
608
sage: R = PolynomialRing(QQ,5,names='a')
609
sage: S.<x> = PolynomialRing(R)
610
sage: p = S(list(R.gens()))*x; p
611
a4*x^5 + a3*x^4 + a2*x^3 + a1*x^2 + a0*x
612
sage: p.parent()
613
Univariate Polynomial Ring in x over Multivariate Polynomial Ring in a0, a1, a2, a3, a4 over Rational Field
614
615
::
616
617
sage: # (YES) Convert the above to Horner's form.
618
sage: # Verify(Horner(p, x), ((((a[5]*x+a[4])*x
619
sage: # +a[3])*x+a[2])*x+a[1])*x);
620
sage: # We use the trick of evaluating the algebraic poly at a symbolic variable:
621
sage: restore('x')
622
sage: p(x)
623
((((a4*x + a3)*x + a2)*x + a1)*x + a0)*x
624
625
::
626
627
sage: # (NO) Convert the result of problem 127 to Fortran syntax.
628
sage: # CForm(Horner(p, x));
629
630
::
631
632
sage: # (YES) Verify that True And False=False.
633
sage: (True and False) == False
634
True
635
636
::
637
638
sage: # (YES) Prove x Or Not x.
639
sage: for x in [True, False]:
640
... print x or (not x)
641
True
642
True
643
644
::
645
646
sage: # (YES) Prove x Or y Or x And y=>x Or y.
647
sage: for x in [True, False]:
648
... for y in [True, False]:
649
... if x or y or x and y:
650
... if not (x or y):
651
... print "failed!"
652
"""
653
654