Add some tests for dgemm.
[maxima.git] / tests / rtest16.mac
blob0ed1e990752b233336a59f6815b8785836619876
1 (kill (all), values);
2 [];
4 /* make sure things work after a reset().  ID: 1986726 and ID: 2787047
5  *
6  * display2d is a resettable option variable. We save the value of display2d
7  * and restore it after the reset. This allows to run the testsuite in both
8  * display modes.
9  */
10 (save:display2d, done);
11 done$
12 (reset(),0);
14 (display2d:save, done);
15 done$
17 /* From A. Reiner.  Fixed in buildq.lisp rev. 1.4.
18    Exposed a bug caused by the dynamical scope of VARLIST. */
19 buildq([foo:sin(baz+bar)],1);
22 buildq([foo:[sin(baz+bar),sin(baz-bar)]],+splice(foo));
23 sin(baz+bar)+sin(baz-bar)$
25 /* Verify that extended functionality of rhs/lhs works as advertised */
27 (kill (x, y, aa, bb, cc), infix("@@"), 0);
30 map (lhs, [aa < bb, aa <= bb, aa = bb, aa # bb, equal (aa, bb), notequal (aa, bb), aa >= bb, aa > bb]);
31 [aa, aa, aa, aa, aa, aa, aa, aa];
33 map (rhs, [aa < bb, aa <= bb, aa = bb, aa # bb, equal (aa, bb), notequal (aa, bb), aa >= bb, aa > bb]);
34 [bb, bb, bb, bb, bb, bb, bb, bb];
36 map (lhs, [foo(x) := 2*x, bar(y) ::= 3*y, '(aa : bb), '(aa :: bb), ?marrow(aa, bb)]);
37 ['(foo(x)), '(bar(y)), aa, aa, aa];
39 map (rhs, [foo(x) := 2*x, bar(y) ::= 3*y, '(aa : bb), '(aa :: bb), ?marrow(aa, bb)]);
40 [2*x, 3*y, bb, bb, bb];
42 [lhs (aa @@ bb), lhs (aa @@ bb @@ cc), rhs (aa @@ bb), rhs (aa @@ bb @@ cc)];
43 [aa, aa @@ bb, bb, cc];
45 map (lhs, [aa + bb, aa - bb, aa * bb, aa / bb, sin(aa), log(aa)]);
46 [aa + bb, aa - bb, aa * bb, aa / bb, sin(aa), log(aa)];
48 map (rhs, [aa + bb, aa - bb, aa * bb, aa / bb, sin(aa), log(aa)]);
49 [0, 0, 0, 0, 0, 0];
51 kill ("@@");
52 done;
54 /* Verify that grind treats nouns correctly. string calls MSIZE in src/grind.lisp.
55  */
57 kill (all);
58 done;
60 string ('(integrate(f(x), x) + integrate(g(x), x, minf, inf) + diff(u, x) + sum(h(x), x, 1, n)));
61 "sum(h(x),x,1,n)+integrate(g(x),x,minf,inf)+integrate(f(x),x)+diff(u,x)";
63 string ('integrate(f(x), x) + 'integrate(g(x), x, minf, inf) + 'diff(u, x) + 'sum(h(x), x, 1, n));
64 "'sum(h(x),x,1,n)+'integrate(g(x),x,minf,inf)+'integrate(f(x),x)+'diff(u,x,1)";
66 /* GREAT puts nounified atoms before others, it appears ... */
67 string (%a%a + %b%b + nounify(%c%c) + nounify(%d%d) + %e%e + %f%f);
68 "%d%d+%c%c+%f%f+%e%e+%b%b+%a%a";
70 string (sin(x) * cos(x) + tan(x));
71 "tan(x)+cos(x)*sin(x)";
73 string ('foo(x, y, z) / bar(a, b, c) + 'baz(%pi - 'quux(%e ^ mumble(%i))));
74 "'foo(x,y,z)/bar(a,b,c)+'baz(%pi-'quux(%e^mumble(%i)))";
76 /* It's conceivable that someday nounified arithmetic operators would be treated differently by grind.
77  * If/when that happens, revise this example accordingly.
78  */
79 string ('"+"(a, b, '"."(c, d), '"^"(e, f)));
80 "'?mplus(a,b,'?mnctimes(c,d),'?mexpt(e,f))";
82 /* Bug 626721 */
83 logarc(atan2(y,x));
84 -%i*log((%i*y+x)/sqrt(x^2+y^2));
85 rectform(ev(%,x=-1,y=1));
86 3*%pi/4; 
89  * Bug [ 1661490 ] An integral gives a wrong result.
90  */
91 (assume(a>0, b>0, sqrt(sqrt(b^2+a^2)-a)*(sqrt(b^2+a^2)+a)^(3/2)-b^2>0),0);
93 radcan(integrate(exp(-(a+%i*b)*x^2),x,minf,inf)/(sqrt(%pi)/sqrt(a+%i*b)));
97  * [ 1663704 ] integrate(sin(r*x)^7/x^4,x,0,inf) -> r^3*false
98  *
99  * Should return the integral instead of producing false.
100  */
101 integrate(sin(a*x)^7/x^4,x,0,inf);
102 'integrate(sin(a*x)^7/x^4,x,0,inf);
104 /* we have assumed a>0 */
105 integrate(%e^(-a*r)*sin(k*r),r,0,inf);
106 k/(k^2+a^2);
109  * Bug [ 1854888 ] hgfred([5],[5], 1) doesn't simplify
110  */
111 hgfred([5],[5],1);
115  * Bug [ 1858964 ] hgfred([7],[-1], x) --/--> error
116  */
117 hgfred([7],[-1],x);
118 und;
121  * Bug [ 1858939 ] hgfred([-1],[-2],x) --> error
122  */
123 hgfred([-1],[-2],x);
124 /* Because of revision 1.110 of hyp.lisp gen_laguerre simplifies
125    -gen_laguerre(1,-3,x)/2; */
126 1+x/2;
129  * Tests for the :: operator
130  */
131 a:b;
134 a::3;
140 p:concat('p,1);
143 p::5;
149 kill(all);
150 done$
152 /* Bug [ 1860250 ] erf(-inf) --> -erf(inf) */
154 erf(-inf);
157 erf(inf);
160 erf(-x) + erf(x);
163 erf(a-b) + erf(b-a);
166 /* Bug [ 1950653 ] bessel_j not simplified 
167  * A few additional related tests added too.
168  */
170 bessel_j(1/2,%pi),besselexpand:true;
172 bessel_y(1/2,%pi/2),besselexpand:true;
175 /* Bug [ 2149714 ] fpprintprec does not work correctly
178 fpprec:16;
181 block([fpprintprec:5], string(1.23b0));
182 "1.23b0";
184 block([fpprintprec:5], string(1.2345b0));
185 "1.2345b0";
187 block([fpprintprec:5], string(1.23456789b0));
188 "1.2346b0";
190 block([fpprintprec:25], string(1.2345678901234567890123456789b0));
191 "1.234567890123457b0";
193 /* verify that fpprintprec behavior matches its description */
195 block ([L1 : [["1.2E-10","1.2E-9","1.2E-8","1.2E-7","1.2E-6","1.2E-5","1.2E-4","0.0012","0.012","0.12","1.2","1.2E+1","1.2E+2",
196   "1.2E+3","1.2E+4","1.2E+5","1.2E+6","1.2E+7","1.2E+8","1.2E+9","1.2E+10"],
197  ["1.23E-10","1.23E-9","1.23E-8","1.23E-7","1.23E-6","1.23E-5","1.23E-4","0.00123","0.0123","0.123","1.23","12.3",
198   "1.23E+2","1.23E+3","1.23E+4","1.23E+5","1.23E+6","1.23E+7","1.23E+8","1.23E+9","1.23E+10"],
199  ["1.234E-10","1.234E-9","1.234E-8","1.234E-7","1.234E-6","1.234E-5","1.234E-4","0.001234","0.01234","0.1234","1.234",
200   "12.34","123.4","1.234E+3","1.234E+4","1.234E+5","1.234E+6","1.234E+7","1.234E+8","1.234E+9","1.234E+10"],
201  ["1.2344E-10","1.2344E-9","1.2344E-8","1.2344E-7","1.2344E-6","1.2344E-5","1.2344E-4","0.0012344","0.012344",
202   "0.12344","1.2344","12.344","123.44","1234.4","1.2344E+4","1.2344E+5","1.2344E+6","1.2344E+7","1.2344E+8",
203   "1.2344E+9","1.2344E+10"],
204  ["1.23443E-10","1.23443E-9","1.23443E-8","1.23443E-7","1.23443E-6","1.23443E-5","1.23443E-4","0.00123443","0.0123443",
205   "0.123443","1.23443","12.3443","123.443","1234.43","12344.3","1.23443E+5","1.23443E+6","1.23443E+7","1.23443E+8",
206   "1.23443E+9","1.23443E+10"],
207  ["1.234432E-10","1.234432E-9","1.234432E-8","1.234432E-7","1.234432E-6","1.234432E-5","1.234432E-4","0.001234432",
208   "0.01234432","0.1234432","1.234432","12.34432","123.4432","1234.432","12344.32","123443.2","1.234432E+6",
209   "1.234432E+7","1.234432E+8","1.234432E+9","1.234432E+10"],
210  ["1.2344321E-10","1.2344321E-9","1.2344321E-8","1.2344321E-7","1.2344321E-6","1.2344321E-5","1.2344321E-4",
211   "0.0012344321","0.012344321","0.12344321","1.2344321","12.344321","123.44321","1234.4321","12344.321","123443.21",
212   "1234432.1","1.2344321E+7","1.2344321E+8","1.2344321E+9","1.2344321E+10"],
213  ["1.23443211E-10","1.23443211E-9","1.23443211E-8","1.23443211E-7","1.23443211E-6","1.23443211E-5","1.23443211E-4",
214   "0.00123443211","0.0123443211","0.123443211","1.23443211","12.3443211","123.443211","1234.43211","12344.3211",
215   "123443.211","1234432.11","1.23443211E+7","1.23443211E+8","1.23443211E+9","1.23443211E+10"],
216  ["1.234432112E-10","1.234432112E-9","1.234432112E-8","1.234432112E-7","1.234432112E-6","1.234432112E-5",
217   "1.234432112E-4","0.001234432112","0.01234432112","0.1234432112","1.234432112","12.34432112","123.4432112",
218   "1234.432112","12344.32112","123443.2112","1234432.112","1.234432112E+7","1.234432112E+8","1.234432112E+9",
219   "1.234432112E+10"],
220  ["1.2344321123E-10","1.2344321123E-9","1.2344321123E-8","1.2344321123E-7","1.2344321123E-6","1.2344321123E-5",
221   "1.2344321123E-4","0.0012344321123","0.012344321123","0.12344321123","1.2344321123","12.344321123","123.44321123",
222   "1234.4321123","12344.321123","123443.21123","1234432.1123","1.2344321123E+7","1.2344321123E+8","1.2344321123E+9",
223   "1.2344321123E+10"],
224  ["1.23443211234E-10","1.23443211234E-9","1.23443211234E-8","1.23443211234E-7","1.23443211234E-6","1.23443211234E-5",
225   "1.23443211234E-4","0.00123443211234","0.0123443211234","0.123443211234","1.23443211234","12.3443211234",
226   "123.443211234","1234.43211234","12344.3211234","123443.211234","1234432.11234","1.23443211234E+7",
227   "1.23443211234E+8","1.23443211234E+9","1.23443211234E+10"],
228  ["1.234432112344E-10","1.234432112344E-9","1.234432112344E-8","1.234432112344E-7","1.234432112344E-6",
229   "1.234432112344E-5","1.234432112344E-4","0.001234432112344","0.01234432112344","0.1234432112344","1.234432112344",
230   "12.34432112344","123.4432112344","1234.432112344","12344.32112344","123443.2112344","1234432.112344",
231   "1.234432112344E+7","1.234432112344E+8","1.234432112344E+9","1.234432112344E+10"],
232  ["1.2344321123443E-10","1.2344321123443E-9","1.2344321123443E-8","1.2344321123443E-7","1.2344321123443E-6",
233   "1.2344321123443E-5","1.2344321123443E-4","0.0012344321123443","0.012344321123443","0.12344321123443",
234   "1.2344321123443","12.344321123443","123.44321123443","1234.4321123443","12344.321123443","123443.21123443",
235   "1234432.1123443","1.2344321123443E+7","1.2344321123443E+8","1.2344321123443E+9","1.2344321123443E+10"],
236  ["1.23443211234432E-10","1.23443211234432E-9","1.23443211234432E-8","1.23443211234432E-7","1.23443211234432E-6",
237   "1.23443211234432E-5","1.23443211234432E-4","0.00123443211234432","0.0123443211234432","0.123443211234432",
238   "1.23443211234432","12.3443211234432","123.443211234432","1234.43211234432","12344.3211234432","123443.211234432",
239   "1234432.11234432","1.23443211234432E+7","1.23443211234432E+8","1.23443211234432E+9","1.23443211234432E+10"],
240  ["1.234432112344321E-10","1.234432112344321E-9","1.234432112344321E-8","1.234432112344321E-7","1.234432112344321E-6",
241   "1.234432112344321E-5","1.234432112344321E-4","0.001234432112344321","0.01234432112344321","0.1234432112344321",
242   "1.234432112344321","12.34432112344321","123.4432112344321","1234.432112344321","12344.32112344321",
243   "123443.2112344321","1234432.112344321","1.234432112344321E+7","1.234432112344321E+8","1.234432112344321E+9",
244   "1.234432112344321E+10"]],
245  L2 : block ([foo : 1.2344321123443211234],
246              makelist (block ([fpprintprec : m], makelist ([m, string (foo*10^n)], n, -10, 10)), m, 2, 16))],
247  map (lambda ([s1, s2], if sequalignore (s1, second (s2)) then true else ['fpprintprec = first (s2), 'expected = s1, 'actual = second (s2)]), apply (append, L1), apply (append, L2)),
248  delete (true, %%));
251 /* verify that fpprintprec = 0 or  fpprintprec > 16 prints readably */
253 block ([foo: float (%pi), L],
254        L: makelist (block ([fpprintprec: m], makelist ([m, foo*10^n, string (foo*10^n)], n, -10, 10)), m, [0, 17, 18, 19, 20]),
255        map (lambda ([L1], if parse_string (L1[3]) = L1[2] then true else ['fpprintprec = first (L1), 'expected = L1[2], 'actual = L1[3]]), apply (append, L)),
256        delete (true, %%));
259 /* SF bug #3213: "fpprintprec do not round bfloat correctly." */
261 (reset (fpprec, bftrunc),
262  fpprintprec:4,
263  string (1.23456789b0));
264 "1.235b0";
266 /* default fpprintprec => print all digits */
268 (reset (fpprintprec), 0);
271 string (bfloat (1000000000000*(1 + 8/9)));
272 "1.888888888888889b12";
274 string (bfloat ((1 + 8/9)/1000000000000)), fpprec=50;
275 "1.8888888888888888888888888888888888888888888888889b-12";
277 string (bfloat (1000000000000*(1 + 8/9))), fpprec=10;
278 "1.888888889b12";
280 string (bfloat ((1 + 8/9)/1000000000000)), fpprec=2;
281 "1.9b-12";
283 /* expect to see bigfloat ending in 5 to round to even */
285 map (string, map (bfloat, 1 + [1, 2, 3, 4, 5, 6, 7]/8)), fpprintprec=3;
286 ["1.12b0", "1.25b0", "1.38b0", "1.5b0", "1.62b0", "1.75b0", "1.88b0"];
288 map (string, map (bfloat, 100 + makelist (i/64, i, 1, 63))), fpprintprec=8;
289 ["1.0001562b2", "1.0003125b2", "1.0004688b2", "1.000625b2", "1.0007812b2", "1.0009375b2", "1.0010938b2", "1.00125b2",
290  "1.0014062b2", "1.0015625b2", "1.0017188b2", "1.001875b2", "1.0020312b2", "1.0021875b2", "1.0023438b2", "1.0025b2",
291  "1.0026562b2", "1.0028125b2", "1.0029688b2", "1.003125b2", "1.0032812b2", "1.0034375b2", "1.0035938b2", "1.00375b2",
292  "1.0039062b2", "1.0040625b2", "1.0042188b2", "1.004375b2", "1.0045312b2", "1.0046875b2", "1.0048438b2", "1.005b2",
293  "1.0051562b2", "1.0053125b2", "1.0054688b2", "1.005625b2", "1.0057812b2", "1.0059375b2", "1.0060938b2", "1.00625b2",
294  "1.0064062b2", "1.0065625b2", "1.0067188b2", "1.006875b2", "1.0070312b2", "1.0071875b2", "1.0073438b2", "1.0075b2",
295  "1.0076562b2", "1.0078125b2", "1.0079688b2", "1.008125b2", "1.0082812b2", "1.0084375b2", "1.0085938b2", "1.00875b2",
296  "1.0089062b2", "1.0090625b2", "1.0092188b2", "1.009375b2", "1.0095312b2", "1.0096875b2", "1.0098438b2"];
298 /* bftrunc=false => don't strip trailing 0's */
300 map (string, map (bfloat, 1 + [1, 2, 3, 4, 5, 6, 7]/8)), fpprintprec=3, bftrunc=false;
301 ["1.12b0", "1.25b0", "1.38b0", "1.50b0", "1.62b0", "1.75b0", "1.88b0"];
303 map (string, map (bfloat, 100 + makelist (i/64, i, 1, 63))), fpprintprec=8, bftrunc=false;
304 ["1.0001562b2", "1.0003125b2", "1.0004688b2", "1.0006250b2", "1.0007812b2", "1.0009375b2", "1.0010938b2", "1.0012500b2",
305  "1.0014062b2", "1.0015625b2", "1.0017188b2", "1.0018750b2", "1.0020312b2", "1.0021875b2", "1.0023438b2", "1.0025000b2",
306  "1.0026562b2", "1.0028125b2", "1.0029688b2", "1.0031250b2", "1.0032812b2", "1.0034375b2", "1.0035938b2", "1.0037500b2",
307  "1.0039062b2", "1.0040625b2", "1.0042188b2", "1.0043750b2", "1.0045312b2", "1.0046875b2", "1.0048438b2", "1.0050000b2",
308  "1.0051562b2", "1.0053125b2", "1.0054688b2", "1.0056250b2", "1.0057812b2", "1.0059375b2", "1.0060938b2", "1.0062500b2",
309  "1.0064062b2", "1.0065625b2", "1.0067188b2", "1.0068750b2", "1.0070312b2", "1.0071875b2", "1.0073438b2", "1.0075000b2",
310  "1.0076562b2", "1.0078125b2", "1.0079688b2", "1.0081250b2", "1.0082812b2", "1.0084375b2", "1.0085938b2", "1.0087500b2",
311  "1.0089062b2", "1.0090625b2", "1.0092188b2", "1.0093750b2", "1.0095312b2", "1.0096875b2", "1.0098438b2"];
313 /* fpprintprec <= fpprec */
315 every (lambda ([s], s="5.6b2"), map (string, makelist (bfloat (5000/9), fpprec, fpprintprec, 40))), fpprintprec=2;
316 true;
318 every (lambda ([s], s="5.56b-4"), map (string, makelist (bfloat (5/9000), fpprec, fpprintprec, 40))), fpprintprec=3;
319 true;
321 every (lambda ([s], s="5.556b2"), map (string, makelist (bfloat (5000/9), fpprec, fpprintprec, 40))), fpprintprec=4;
322 true;
324 every (lambda ([s], s="5.5556b-4"), map (string, makelist (bfloat (5/9000), fpprec, fpprintprec, 40))), fpprintprec=5;
325 true;
327 every (lambda ([s], s="5.55556b2"), map (string, makelist (bfloat (5000/9), fpprec, fpprintprec, 40))), fpprintprec=6;
328 true;
330 every (lambda ([s], s="5.555556b-4"), map (string, makelist (bfloat (5/9000), fpprec, fpprintprec, 40))), fpprintprec=7;
331 true;
333 every (lambda ([s], s="5.5555556b2"), map (string, makelist (bfloat (5000/9), fpprec, fpprintprec, 40))), fpprintprec=8;
334 true;
336 /* fpprintprec >= fpprec */
338 /* bfloat generally produces a number which is not exactly the same as its argument.
339  * When bfloat returns a bigger number, it should round up when formatted,
340  * and when it's smaller, it should round down.
341  * When bfloat returns an exact result, in this case it should round up (to 6).
342  * Consider this when figuring out the correct output for these examples.
343  */
345 (fpprec : 2, 
346  mybf : bfloat(5/9000),
347  is (rationalize (mybf) >= (5*sum(10^i, i, 0, fpprec)/10^(fpprec + 1))/1000));
348 true;
350 every (lambda ([s], s="5.6b-4"), map (string, makelist (mybf, fpprintprec, fpprec, 40)));
351 true;
353 (fpprec : 3,
354  mybf : bfloat(5000/9),
355  is (rationalize (mybf) >= (5*sum(10^i, i, 0, fpprec)/10^(fpprec + 1))*1000));
356 true;
358 every (lambda ([s], s="5.56b2"), map (string, makelist (mybf, fpprintprec, fpprec, 40)));
359 true;
361 (fpprec : 4,
362  mybf : bfloat(5/9000),
363  is (rationalize (mybf) >= (5*sum(10^i, i, 0, fpprec)/10^(fpprec + 1))/1000));
364 true;
366 every (lambda ([s], s="5.556b-4"), map (string, makelist (mybf, fpprintprec, fpprec, 40)));
367 true;
369 (fpprec : 5,
370  mybf : bfloat(5000/9),
371  is (rationalize (mybf) >= (5*sum(10^i, i, 0, fpprec)/10^(fpprec + 1))*1000));
372 false;
374 every (lambda ([s], s="5.5555b2"), map (string, makelist (mybf, fpprintprec, fpprec, 40)));
375 true;
377 (fpprec : 10,
378  mybf : bfloat(5/9000),
379  is (rationalize (mybf) >= (5*sum(10^i, i, 0, fpprec)/10^(fpprec + 1))/1000));
380 true;
382 every (lambda ([s], s="5.555555556b-4"), map (string, makelist (mybf, fpprintprec, fpprec, 40)));
383 true;
385 (fpprec : 20,
386  mybf : bfloat(5000/9),
387  is (rationalize (mybf) >= (5*sum(10^i, i, 0, fpprec)/10^(fpprec + 1))*1000));
388 true;
390 every (lambda ([s], s="5.5555555555555555556b2"), map (string, makelist (mybf, fpprintprec, fpprec, 40)));
391 true;
393 (fpprec : 40,
394  mybf : bfloat(5/9000),
395  is (rationalize (mybf) >= (5*sum(10^i, i, 0, fpprec)/10^(fpprec + 1))/1000));
396 false;
398 every (lambda ([s], s="5.555555555555555555555555555555555555555b-4"), map (string, makelist (mybf, fpprintprec, fpprec, 40)));
399 true;
401 (reset (fpprec), 0);
405  * Bug 2142758: integrate(sqrt(2-2*x^2)*(sqrt(2)*x^2+sqrt(2))/(4-4*x^2),x,0,1)
406  */
407 integrate(sqrt(2-2*x^2)*(sqrt(2)*x^2+sqrt(2))/(4-4*x^2),x,0,1);
408 3*%pi/8;
411 integrate(sqrt(1-x^2)*(x^2+1)/(2-2*x^2),x,0,1);
412 3*%pi/8;
414 integrate(sqrt(1-x^2)*(x^2+1)/(1-x^2),x,0,1);
415 3*%pi/4;
418  * Bug [ 2208303 ] Problem with jacobi_dn and elliptic_kc
419  */
420 jacobi_dn(elliptic_kc(m)*t,m);
421 jacobi_dn(elliptic_kc(m)*t,m);
424  * Bug [ 2180110 ] GCL do not signal an overflow converting bigfloat to float
425  */
426 errcatch(float(2b400));
429 errcatch(float(bfloat(2^1024)));
433  * Bug [ 2055235 ] Plot leaves range with jacobi functions
435  * Actually jacobi_cn(100, .7) is computed inaccurately.  Just check that abs(jacobi_cn(100,.7)) < 1
436  */
438 is(abs(jacobi_cn(100.0, 0.7)) < 1);
439 true$
442  * Bug [ 1658067 ] jacobi_sn(elliptic_kc(1-m)*%i/2,m) isn't simplified
444  * This test (and Maxima) used to be wrong.  This is related to the
445  * jacobi_sc(elliptic_kc(m)/2,m) test below.
446  */
447 jacobi_sn(elliptic_kc(1-m)*%i/2,m);
448 %i/m^(1/4)$
450 jacobi_sc(elliptic_kc(m)+u,m);
451 -jacobi_cs(u,m)/sqrt(1-m)$
454  * Maxima used to get this wrong by returning the reciprocal instead
455  */
456 jacobi_sc(elliptic_kc(m)/2,m);
457 1/(1-m)^(1/4);
460  * Bug [ 2505945 ] - hgfred([2,-1/2],[3],-x^2);
462  * Shouldn't signal from diff about non-variable second arg.
464  * The expected value here is computed from
465  * factor(ratsimp(subst([z=-x^2],hgfred([2,-1/2],[3],z))))
466  */
467 factor(ratsimp(hgfred([2,-1/2],[3],-x^2)));
468 4*(3*x^4*sqrt(x^2+1)+x^2*sqrt(x^2+1)-2*sqrt(x^2+1)+2)/(15*x^4)$
471  * Bug 2534420: asinh(%i*2b0) causes error
472  */
473 is(abs(asinh(%i*2b0)-expand(bfloat(asinh(%i*2)))) < 3b-17);
474 true;
477  * Bug 2543079: bfloat(gamma(3/4)/gamma(1/4)) is wrong.
478  */
479 bfloat(gamma(3/4)/gamma(1/4));
480 3.379891200336424b-1;
483  * Bug 2582034 - hgfred([a/2,-a/2],[1/2],z) causes error
484  */
485 (assume(zn<0), done);
486 done;
488 hgfred([a/2,-a/2],[1/2],zn);
489 ((%i*sqrt(zn)+sqrt(1-zn))^a+(sqrt(1-zn)-%i*sqrt(zn))^a)/2$
492  * Bug 2618401 - bfloat produces incorrect answer
493  */
494 is(abs(bfloat((sqrt(2)+2)*%pi^(3/2)/(8*gamma(3/4)^2))-float((sqrt(2)+2)*%pi^(3/2)/(8*gamma(3/4)^2))) < 1d-15);
495 true;
497 /* (-1.0b0)^(1/3) vs (-1.0d0)^(1/3) - ID: 619927 */
498 (-1b0)^(1/3);
499 -1.0b0;
501 (-1.0)^(1/3);
502 -1.0;
504 domain:complex;
505 complex;
507 (-1b0)^(1/3);
508 1.0b0*(-1)^(1/3);
510 (-1.0)^(1/3);
511 1.0*(-1)^(1/3);
513 domain:real;
514 real;
516 (-1b0)^(1/3),domain:complex,m1pbranch:true;
517 1.0b0*(sqrt(3)*%i/2+1/2);
519 (-1.0)^(1/3),domain:complex,m1pbranch:true;
520  1.0*(sqrt(3)*%i/2+1/2);
524  * Bug [ 2688847 ] float of rats rounds incorrectly
525  */
526 float((2^60-1)/2^60)-1;
527 0.0;
528 float((2^1000-1)/2^1000)-1;
529 0.0;
532  * Bug [ 2687962 ] hgfred([-3/2,1],[-1/2],-t) division by zero
534  * Solution from functions.wolfram.com
535  */
536 ratsimp(hgfred([-3/2,1],[-1/2], t));
537 1+3*t-3*t^(3/2)*atanh(sqrt(t));
540  * Bug 2793827: internal error in integrate
541  */
542 (assume(n>0),declare(n,integer),0);
545 integrate((g32475^n*(g32475*n-n-1)/(g32475-1)^2+1/(g32475-1)^2)/(1-g32475)
546  -(g32475^(2*n+1)*(g32475*n-n-1)/(g32475-1)^2+g32475^(n+1)/(g32475-1)^2)
547   /(1-g32475),g32475,0,1);
548 (n^2+2*n+1)/2-1/2$
550 kill(all);
551 done;
554  * Bug 609464 : 1+%e,numer and %e^%e,numer
556  * The simplifier has been extended to handle %e like other constants.
557  * In addition functions with arguments which involve %e simplify
558  * accordingly.
559  */
561 %e,numer;
562 2.7182818284590451;
564 %e+1,numer;
565 3.7182818284590451;
567 %e^%e,numer;
568 15.154262241479262;
570 %e^x,numer;
571 %e^x;
573 sin(%e),numer;
574 0.41078129050290885;
576 sin(%e+1),numer;
577 -0.54525155669233449;
579 /* Do not simplify, when %e is the base of an expression and %enumer FALSE*/
581 sin(%e^(2*x+1)),numer;
582 sin(%e^(2*x+1));
584 sin(%e^(%e^(2*x+1))),numer;
585 sin(%e^(%e^(2*x+1)));
587 /* Additionally simplifications when %enumer TRUE */
589 %enumer:true;
590 true;
592 sin(%e^x),numer;
593 sin(2.7182818284590451^x);
595 sin(%e^(%e^(2*x+1))),numer;
596 sin(2.7182818284590451^(2.7182818284590451^(2*x+1)));
598 %enumer:false;
599 false;
602  * Bug ID: 2797885 - "problem with integration"
604  * integrate(exp(%i*x)*sin(x),x) generates a Lisp error.
606  * This is a special case for the integrand: exp(a*x)*sin(b*x),
607  * with a^2+b^2 equal to zero.
608  */
610 /* This is the general case for an integral with exp and sin or cos */
611 integrate(exp(a*x)*sin(b*x),x);
612 %e^(a*x)*(a*sin(b*x)-b*cos(b*x))/(b^2+a^2);
614 integrate(exp(a*x)*cos(b*x),x);
615 %e^(a*x)*(b*sin(b*x)+a*cos(b*x))/(b^2+a^2);
617 /* Now the special case with a=%i and b=1 */
618 expand(integrate(exp(%i*x)*sin(x),x));
619 %i*x/2-%e^(2*%i*x)/4;
621 expand(integrate(exp(x)*sin(%i*x),x));
622 %i*%e^(2*x)/4-%i*x/2;
624 expand(integrate(exp(%i*x)*cos(x),x));
625 x/2-%i*%e^(2*%i*x)/4;
627 expand(integrate(exp(x)*cos(%i*x),x));
628 %e^(2*x)/4+x/2;
630 /* Bug ID: 932076 -  ode2( 'diff(y,x)=%i*y+sin(x), y, x) => div by 0
632  * This bug is related to the Bug ID: 2797885 - "problem with integration"
633  */
635 ode2('diff(y,x)-%i*y-sin(x),y,x);
636 y = (%c-%i*(x-%i*%e^-(2*%i*x)/2)/2)*%e^(%i*x);
639  * Bug ID: 826623 "simplifer returns %i*%i"
640  * 
641  * Some examples to show simplification of expressions of the form
642  * (a*b*...)^q*(a*b*...)^r, where q+r=1
643  */
645 sqrt(-%i)*sqrt(-%i)*%i;
648 sqrt(a*b)*sqrt(a*b)*a*b;
649 a^2*b^2;
651 (a*b*c)^(3/4)*(a*b*c)^(1/4)*c;
652 a*b*c^2;
655  * Bug ID: 2792493 "hgfred([1],[-5.2],x);"
656  */
657 hgfred([1],[-5.2],x);
658 %f[1,1]([-6.2],[-5.2],-x)*%e^x$
660 /* BUG ID: 721575 2/sqrt(2) doesn\'t simplify */
661 2/sqrt(2);
662 sqrt(2);
664 (1/2)*sqrt(2);
665 1/sqrt(2);
667 sqrt(2)*(1/2);
668 1/sqrt(2);
670 /* BUG ID 2029041 a*sqrt(2)/2 unsimplified */
672 a*sqrt(2)/2;
673 a/sqrt(2);
675 /* BUG ID 1923119 1/sqrt(8)-sqrt(8)/8 */
677 1/sqrt(8)-sqrt(8)/8;
680 /* BUG ID 1927178 integrate(sin(t),t,%pi/4,3*%pi/4) */
682 integrate(sin(t),t,%pi/4,3*%pi/4);
683 sqrt(2);
685 /* BUG ID: 1480562 2*a*2^k isn't simplified to a*2^(k+1) */
687 2*a*2^k;
688 a*2^(k+1);
690 a*2^k*2;
691 a*2^(k+1);
693 /* Some examples to show simplification of expressions
694  * with floating point and bigfloat numbers after improvement
695  * of plusin
696  */
698 (4.0*x-4.0*x);
699 0.0;
700 (4.0*x-3.0*x);
701 1.0*x;
702 (4.0*x-3.0*x)/2;
703 0.5*x;
705 (4.0b0*x-4.0*x);
706 0.0b0;
707 (4.0b0*x-3.0*x);
708 1.0b0*x;
709 (4.0b0*x-3.0*x)/2;
710 0.5b0*x;
712 /* BUG ID: 1996354 unsimplifed result from expand */
714 expand((%e^(-2*sqrt(2))*(%e^(2*sqrt(2))+2*%e^sqrt(2)+1)^2)/16
715        +(%e^(-2*sqrt(2))*(%e^(2*sqrt(2))-2*%e^sqrt(2)+1)^2)/16
716        -(%e^(-2*sqrt(2))*(%e^(2*sqrt(2))-1)^2)/8);
719 /* BUG ID: 631216 - "horner([...],x)/FIX" 
720    horner now maps over lists, matrices and equations.
721  */
723 horner(x^2+x=a*x^2+b*x);
724 x*(x+1) = x*(a*x+b);
725 horner([x^2+x,x^3+x,x^4+x]);
726 [x*(x+1),x*(x^2+1),x*(x^3+1)];
728 /* BUG ID: 2699862 "derivative of polylogarithm"
729  * The noun form is not put on the property list, but NIL. The routine 
730  * sdiffgrad generates a noun form, when the derivative is not known.
731  */
733 diff(li[n](x),n);
734 'diff(li[n](x),n);
736 diff(li[n*x](x),x);
737 'diff(li[n*x](x),x);
739 diff(li[n](x),x,1,n,1);
740 'diff(li[n-1](x),n)/x;
742 /* Not reported as a bug, but the same problem for the function psi */
744 diff(psi[n](x),n);
745 'diff(psi[n](x),n);
747 diff(psi[n*x](x),x);
748 'diff(psi[n*x](x),x);
750 diff(psi[n](x),x,1,n,1);
751 'diff(psi[n+1](x),n);
753 /* BUG ID: 2824909 " exp(%i*%pi/4) not simplified" 
754  * Check the simplification of exp(%i*%pi/4) and exp(-%i*pi/4)
755  */
757 exp(%i*%pi/4);
758 1/sqrt(2)+%i/sqrt(2);
759 exp(-%i*%pi/4);
760 1/sqrt(2)-%i/sqrt(2);
763  * Bug ID: 2831259 - bfloat() underflow bug
764  */
765 fpprec:500;
766 500;
767 float(0b0);
768 0.0;
771  * BUG ID: 2835098 - SIGN-PREP strangeness
772  */
774 block ([?limitp : true], sign (foo (x)));
775 pnz;
777 integrate(sqrt(2*m*(E[n]-U(x))),x,-x[0],x[0])=(n-1/2)*%pi*hbar;
778 sqrt(2)*'integrate(sqrt(m*(E[n]-U(x))),x,-x[0],x[0]) = %pi*hbar*(n-1/2);
780 integrate(f(x),x,x[0],x[1]);
781 'integrate(f(x),x,x[0],x[1]);
784  * BUG ID: 2840566 - defint fails to determine if one of its limit is real
785  */
787 (assume(b>0,c>0),done);
788 done;
790 integrate(x,x,0,sqrt(b^2+(b-c)^2));
791 (c^2-2*b*c+2*b^2)/2;
794  * BUG ID: 2842060 - unsimplified result from integrate
795  */
797 /* The result for a general symbol x */
798 integrate(1/x/sqrt(x^2-1),x);
799 -asin(1/abs(x));
801 (assume(x>0), done);
802 done;
804 /* abs(x) simplifies to x for x>0 */
805 integrate(1/x/sqrt(x^2-1),x);
806 -asin(1/x);
808 (forget(x>0), done);
809 done;
811 /* 
812  * Bug ID: 2820202 - rootscontract(%i/2) 
813  */
814 rootscontract(%i/2);
815 %i/2;
817 /* Bug ID: 2872738 - sign(-(1/n)*(-1)^n)
818  * We got the error because of the simplification
819  *  (-1)^n*(-1) -> (-1)^(n+1) and not -(-1)^n
820  * The other case 
821     (-1)*(-1)^n simplifies already to -(-1)^n
822  * Adding tests for both cases.
823  */
824 kill(all);
825 done;
826 sign(-(1/n)*(-1)^n);
828 (-1)*(-1)^n;
829 -(-1)^n;
830 (-1)^n*(-1);
831 -(-1)^n;
833 /* Bug ID: 2835634 - logcontract broken
834  * Bug ID: 1467368 - logcontract returns unsimplified expr
835  */
836 logcontract(log(x)-log(2));
837 log(x/2);
838 /* Check that we do not break the following again */
839 logcontract(log(%e*k)-log(%e^-1*k));
841 log(%e^2),logexpand:false;
844 /* Bug ID: 2880923 - realpart --> floating-point-overflow
845  */
846 sign(exp(2009));
847 pos;
848 realpart(sqrt(4*%e^2009-3)-1);
849 sqrt(4*%e^2009-3)-1;
850 sqrt(4*exp(2009));
851 2*%e^(2009/2);
853 /* Bug ID: 640332 - Need to specdisrep more systematically
854    Add the examples of the bug report.
855  */
856 ratdisrep(diff(rat(x),rat(x)));
858 diff(x,rat(x));
860 outofpois(diff(intopois(sin(x)),x));
861 cos(x);
862 taylor(intopois(sin(x)),x,0,3);
863 x-x^3/6;
864 ratsimp(intopois(sin(x)));
865 sin(x);
867 /* Bug ID: 627759 - Ratdisrep of aggregates 
868  */
869 ratdisrep(rat(x=y));
870 x = y;
871 ratdisrep(rat([x=a,y=b]));
872 [x = a,y = b];
873 ratdisrep(rat(matrix([a,b],[c,d])));
874 matrix([a,b],[c,d]);
876 /* Bug ID: 711885 - Rootscontract with imaginaries fails
877  */
878 (oldvalue:radexpand, radexpand:false, done);
879 done;
881 rootscontract(((sqrt(3)*%i+1)^(3/2)-4*%i)/sqrt(sqrt(3)*%i+1));
882 ((sqrt(-3)+1)^(3/2)-4*%i)/sqrt(sqrt(-3)+1);
884 /* It is a problem of the simplifier. Show that it works */
885 sqrt(1/(1+sqrt(-3)));
886 1/sqrt(sqrt(-3)+1);
888 (radexpand:oldvalue, done);
889 done;
891 /* BUG ID: 767556 - Distributing operations over =
892  * The operators "." and "^^" distribute over equations.
893  */
895 x . (a=b);
896 x . a = x . b;
898 (a=b)^^x;
899 a^^x = b^^x;
901 /* A more complicated example */
902 x . ((2*a+b . c) = x . (y + z))^^w;
903 x . (b . c+2*a)^^w = x . (x . (z+y))^^w;
905 /* Bug ID: 2914176 - Conversion of rational to bfloat is inaccurate
907  * The difference should be 1/262144, but we don't check for that.
908  */
909 (oldfpprec:fpprec, fpprec:5, done);
910 done;
911 is(bfloat((2^20+1)/(2^20-1)) - 1b0 > 0);
912 true;
914 /* Related to the fix for 2914176.  Didn't handle the ratio 0/1 */
915 is(equal(0b0, 0));
916 true;
918 (fpprec:oldfpprec, done);
919 done;
921 /* Bug ID:2933882 - Power function: 0^a not fully implemented
922  * Show some simplifications of 0^a
923  */
925 assume(a>0);
926 [a>0];
927 0^a;
929 errcatch(0^-a);
931 0^(a+%i);
933 0^(1/2+%i);
935 errcatch(0^(-1/2+%));
937 errcatch(0^%i);
938 []; 
940 forget(a>0);
941 [a>0];
943 /* Bug ID: 2938078 - Crash on attached input
944  */
946 declare(n,integer, j,noninteger);
947 done;
948 assume(equal(x,n), equal(y,j), equal(z,i));
949 [equal(x, n), equal(y, j), equal(z,i)];
951 featurep(x,integer);
952 true;
953 featurep(x,noninteger);
954 false;
955 featurep(y,integer);
956 false;
957 featurep(y,noninteger);
958 true;
960 diff(z+1,z);
963 remove(n,integer, j,noninteger);
964 done;
965 forget(equal(x,n), equal(y,j), equal(z,i));
966 [equal(x, n), equal(y, j), equal(z, i)];
968 /* Bug ID:  2948800 - integrate((1-cos(2*x)^2)^2/x^4,x,0,inf) wrong
969  */
971 integrate((1-cos(2*x)^2)^2/x^4,x,0,inf);
972 8*%pi/3;
974 assume(a>0);
975 [a>0];
977 /* The more general type with an argument a*x and a positive */
978 integrate((1-cos(a*x)^2)^2/x^4,x,0,inf);
979 %pi*a^3/3;
981 forget(a>0);
982 [a>0];
984 /* Bug ID: 777564 - subtraction "-"(a,b) should work/FIX */
986 "-"();
989 "-"(a);
992 "-"(2*a);
993 -2*a;
995 "-"(a+b);
996 -b-a;
998 "-"(a+b+c);
999 -c-b-a;
1001 "-"(100,20,10);
1004 map("-",[a,x,100],[b,y,20]);
1005 [a-b,x-y,80];
1007 map("-",[a,x,100],[b,y,20],[c,z,10]);
1008 [-c-b+a,-z-y+x,70];
1010 /* Bug ID: 910270 - 1/+3*x parses as 1/(+3*x)
1011  * Show that the "+" operator can be used as a prefix operator too.
1012  */
1013 1/+3*x;
1014 x*1/3;
1015 1/+x/3;
1016 1/(3*x);
1017 a^+b*c;
1018 c*a^b;
1020 /* Bug ID: 2961822 - sinh(0.0b0) causes Maxima to abort
1021  */
1022 sinh(0.0b0);
1023 0.0b0;
1025 /* Bug ID: 1219846 - properties of translated functions
1026  * The property noun is already present
1027  */
1028 kill(f);
1029 done;
1030 f(x):=x;
1031 f(x):=x;
1032 properties(f);
1033 [function,noun];
1034 translate(f);
1035 [f];
1036 properties(f);
1037 [transfun,function,noun];
1038 kill(f);
1039 done;
1041 /* Bug ID: 2968344 - gamma_incomplete(1.0, 4.368265444147715e+19) fails
1042  */
1043 gamma_incomplete(1.0, 4.368265444147715e+19);
1044 0.0;
1046 /* Bug ID: 643254 - orderlessp([rat(x)], [rat(x)])
1047  */
1048 orderlessp([rat(x)],[rat(x)]);
1049 false;
1051 /* Bug ID: 781657 - binomial(x,x) => 1, but binomial(-1,-1) => 0
1052  * binomial(x,x) simplifies to 1 only if the sign of x is known not to be 
1053  * negative.
1054  */
1055 is(equal(binomial(x,x),1));
1056 'unknown;
1057 is(equal(binomial(x^2,x^2),1));
1058 true;
1060 /* Bug ID:856209 - inconsistency between facts() and facts(v)
1061  * Show that facts(expr) now works more general.
1062  */
1063 assume(z+a>0,b>z);
1064 [a+z>0,b>z];
1065 facts(a);
1066 [a+z>0];
1067 facts(b);
1068 [b>z];
1069 facts(z);
1070 [a+z>0,b>z];
1071 facts(a+z);
1072 [a+z>0];
1073 forget(z+a>0,b>z);
1074 [a+z>0,b>z];
1076 /* Bug ID: 840848 - trigreduce doesn't enter unknown functions
1077  */
1078 trigexpand(f(sin(2*x)));
1079 f(2*cos(x)*sin(x));
1080 trigreduce(%);
1081 f(sin(2*x));
1083 /* Bug ID: 2954472 - rectform with large floats gives bad answer
1085 is(abs(rectform(1e160/(1e160+%i))-1) < 1e-160);
1086 true;
1088 is(abs(rectform(1e160/(1e160+3/2*%i))-1) < 1.5e-160);
1089 true;
1091 /* Bug ID: 2953369 - Definite Integration of 1/(a-b*cos(x)) wrong
1093  * For simplicity we test the equivalent integrate(1/(1-r*cos(x)),x,0,%pi).
1094  */
1095 /* These assumes are to answer the questions integrate (from routine unitcir) will ask */
1096 (assume(r>0,r<1,abs(sqrt(1-r^2)-1)/r-1 < 0, sqrt(1-r^2)-r+1>0), 0);
1098 integrate(1/(1-r*cos(x)),x,0,%pi);
1099 %pi/sqrt(1-r^2);
1101 /* Bug ID: 2907727 - Incorrect Integral with option integrate_use_rootsof
1102  * :true
1103  */
1104 %rnum:0;
1106 integrate((d*x^2+2*c*x+3*b)/(g*r*x^3+d*x^2+c*x+b), x), integrate_use_rootsof:true;
1107 lsum((%r1^2*d+2*%r1*c+3*b)*log(x-%r1)/(3*%r1^2*g*r+2*%r1*d+c),%r1,
1108              rootsof(g*r*%r1^3+d*%r1^2+c*%r1+b,%r1));
1110 /* Bug ID: 2880797 - bad answer in integrate(sqrt(sin(t)^2+cos(t)^2),t,0,2*%pi)
1112  */
1113 integrate(sqrt(sin(t)^2+cos(t)^2),t,0,2*%pi);
1114 2*%pi;
1116 /* Bug ID: 2980551 - Inconsistent simplification of exp(x*%i*%pi)
1118  * These examples show consistent simplification for x an expression which
1119  * can contain float or bigfloat values
1120  */
1122 exp(2*%i*%pi);
1125 exp((2+x)*%i*%pi);
1126 exp(x*%i*%pi);
1128 exp(2*%i*%pi+x*%i*%pi);
1129 exp(x*%i*%pi);
1131 log(exp((2+x)^2*%i*%pi));
1132 (2+x)^2*%i*%pi;
1134 exp(2.0*%i*%pi);
1135 1.0;
1137 exp((2.0+x)*%i*%pi);
1138 exp(1.0*x*%i*%pi);
1140 exp(2.0*%i*%pi+x*%i*%pi);
1141 exp(1.0*x*%i*%pi);
1143 log(exp((2.0+x)^2*%i*%pi));
1144 (2.0+x)^2*%i*%pi;
1146 exp(2.0b0*%i*%pi);
1147 1.0b0;
1149 exp((2.0b0+x)*%i*%pi);
1150 exp(1.0b0*x*%i*%pi);
1152 exp(2.0b0*%i*%pi+x*%i*%pi);
1153 exp(1.0b0*x*%i*%pi);
1155 log(exp((2.0b0+x)^2*%i*%pi));
1156 (2.0b0+x)^2*%i*%pi;
1158 exp(3/2*%pi*%i);
1159 -%i;
1160 exp(1.5*%pi*%i);
1161 -1.0*%i;
1162 exp(1.5b0*%pi*%i);
1163 -1.0b0*%i;
1165 exp((3/2+x)*%pi*%i);
1166 -%i*exp(%i*%pi*x);
1167 exp((1.5+x)*%pi*%i);
1168 -%i*exp(1.0*%i*%pi*x);
1169 exp((1.5b0+x)*%pi*%i);
1170 -%i*exp(1.0b0*%i*%pi*x);
1172 /* Bug ID: 2781127 - bfpsi0 of complex
1174  * (The result was not in rectangular form but it should be.)
1175  */
1176 bfpsi0(4.5 + %i,15);
1177 2.43845477527606b-1*%i + 1.41875534014717b0;
1179 /* Bug ID: 2988190 - atan2(1b20,-1b0); badly wrong
1180  * It's really a bug in atan for x < -1, so test both.
1181  */
1182 (fpprec:16, atan(-1b20));
1183 -1.570796326794897b0;
1184 atan2(1b20,-1b0);
1185 1.570796326794897b0;
1188  * Bug ID: 2991924 - Incorrect integration of rational functions
1189  */
1190 integrate(1/(x^4-2),x,0,1) - integrate(1/(x^2-sqrt(2))/(x^2+sqrt(2)),x,0,1);
1193 integrate(1/(x^6-4),x,0,1) - integrate(1/(x^3-2)/(x^3+2),x,0,1);
1196 /* BUG ID:  2113751 - Incomprehensible behavior of coeff()
1197  */
1198 coeff(2*%e^x, x, 0);
1199 2*%e^x;
1201 /* For numerical tests */
1202 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
1203 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
1205 /* Bug ID: 2997276 - zeta(3),numer; gives Lisp error 
1206  * 
1207  * Also add a test for complex rational argument, which wasn't handled
1208  * correctly either.
1210  * Some Lisp implementations fail these tests because things like
1211  * (cl:expt 2d0 3) only gives single-float accuracy (but with
1212  * double-float precision).
1213  */
1214 closeto(zeta(3)-1.202056903159594,1e-15), numer:true;
1215 true;
1216 closeto(zeta(3+%i)-(1.10721440843141 - .1482908671781754*%i), 1e-15);
1217 true;
1219  * Reported on mailing list 2011-05-22 by Thomas Dean:
1221  * plot2d(abs(zeta(1/2+x*%i)),[x,0,36]) causes a Lisp error with
1222  * clisp.
1224  */
1225 closeto(abs(zeta(1/2+.5*%i)) - 1.06534921249378, 1e-14);
1226 true;
1228 /* Bug ID: 2997401 - float(log(200!)) produces an error
1230  */
1231 closeto(float(log(200!))-863.2319871924054, 1e-15);
1232 true;
1234 closeto(float(log((1+200!)/7))-861.2860770433501, 1e-15);
1235 true;
1237 /* Additional tests */
1238 closeto(float(log(-1))-float(%pi)*%i, 1e-15);
1239 true;
1241 closeto(float(log((1+200!)/(-7))) - (3.141592653589793*%i + 861.2860770433501), 1e-14);
1242 true;
1244 closeto(float(log((1+200!)+(1+199!)*%i))- (.004999958333958322*%i + 863.2319996922491), 1e-15);
1245 true;
1247 closeto(float(log((1+200!)/7+(1+199!)/11*%i)) - (.003181807444342708*%i + 869.9736929490153), 1e-15);
1248 true;
1250 /* Bug ID: 2306402 - scalarp bug
1251  * Bug ID: 1985748 - array and scalar declarations yield inconsistent results
1252  * Examples from the bug report to show consistent behavior of scalarp
1253  */
1255 declare(x,scalar);
1256 done;
1257 scalarp(foo(x));
1258 true;
1259 scalarp(foo(1));
1260 true;
1261 scalarp(foo(x,1));
1262 true;
1264 scalarp(x);
1265 true;
1266 scalarp(x[1]);
1267 true;
1268 array(x,5);
1270 scalarp(x);
1271 true;
1272 scalarp(x[1]);
1273 true;
1274 nonscalarp(x);
1275 false;
1277 kill(x);
1278 done;
1280 /* Bug ID:1723548 - gradef for variables: not used in diff
1281  * Show that the total differential of f works in expressions too.
1282  */
1283 depends(f,[x,y]);
1284 [f(x,y)];
1285 diff(f);
1286 'diff(f,y,1)*del(y)+'diff(f,x,1)*del(x)$
1287 diff(3*f);
1288 3*'diff(f,y,1)*del(y)+3*'diff(f,x,1)*del(x)$
1289 diff(a*f);
1290 a*'diff(f,y,1)*del(y)+a*'diff(f,x,1)*del(x)+f*del(a)$
1291 remove(f,dependency);
1292 done;
1294 /* Bug ID: 1089719 addrow creates strange matrix
1295  */
1296 m:matrix([0,0]);
1297 matrix([0,0]);
1298 m:addrow(m,m);
1299 matrix([0,0],[0,0])$
1300 m[1,1]:11;
1303 matrix([11,0],[0,0])$
1304 kill(m);
1305 done;
1307 /* Bug ID: 1663385 - declare multiplicative - wrong simplification
1308  */
1309 declare(f,additive,f,multiplicative);
1310 done;
1311 f(a*b+c);
1312 f(a)*f(b)+f(c);
1313 kill(f);
1314 done;
1316 /* Bug ID: 816808 - subst(in)part of rat -- internal errs
1317  */
1318 substpart(x,2/3,2);
1319 2/x;
1320 substinpart(4,2/3,2);
1321 1/2;
1323 /* Bug ID: 1117533 - letsimp complains about assignment to %pi
1324  */
1325 matchdeclare(a,true);
1326 done;
1327 (let(%pi*a,foo(a)),done);
1328 done;
1329 letsimp(%pi*x);
1330 foo(x);
1331 remlet(%pi*a);
1332 done;
1334 /* Bug ID: 2805600 depends() partially prevents diff() to work
1335  */
1336 depends(t,x);
1337 [t(x)];
1338 diff(f(t),z);
1340 remove(t,dependency);
1341 done;
1343 /* Bug ID: 1184718 - AT needs soime basic simplifications
1344  */
1345 'at(2,x=0);
1348 /* Bug ID: 2998227 - spurious at(0,A=0)
1349  */
1350 taylor(integrate(gamma(x+1),x,0,A),A,0,3),nouns;
1351 A-%gamma*A^2/2+(6*%gamma^2+%pi^2)*A^3/36;
1353 /* Bug ID: 3010829 - numerical evaluation of elliptic_ec fails for argument > 1
1354  */
1355 closeto(elliptic_ec(2.0)-(.5990701173677959*%i+0.599070117367796), 1.5e-15);
1356 true;
1358 /* Bug ID: 1929287 - 0.0 + [0] ---> [0]
1359  */
1360 0.0+[0];
1361 [0.0];
1362 0.0b0+[0];
1363 [0.0b0];
1364 0.0+matrix([0,1/2,1,x]);
1365 matrix([0.0,0.5,1.0,x]);
1366 0.0b0+matrix([0,1/2,1,x]);
1367 matrix([0.0b0,5.0b-1,1.0b0,x]);
1369 /* Bug ID: 2996106 - at(diff(f(x,y),x,1,y,1),[x=a,y=b]) is wrong
1370  */
1371 at(diff(f(x,y),x,1,y,1),[x=a,y=b]);
1372 'at('diff(f(x,y),x,1,y,1),[x = a,y = b]);
1374 /* Bug report ID: 2556133 - "at" should do parallel substitutions
1375  */
1376 errcatch(at(atan2(y^2+1,x),[y=%i,x=0]));
1378 errcatch(at(atan2(y^2+1,x),[x=0,y=%i]));
1381 /* Bug report ID: 2014941 - compositions of 'at'
1382  */
1383 at(at(diff(f(x),x),[x=b]),[b=y]);
1384 'at('diff(f(x),x,1),[x = y]);
1386 at(diff(f(x,y),x,1,y,1),[x=a,y=b]) - at(diff(f(x,y),x,1,y,1),[y=b,x=a]);
1389 /* Bug report ID: 1677217 - composistions of 'at'
1390  */
1391 depends(y,[x,z]);
1392 [y(x,z)];
1393 at(at(diff(y,x),x=a),z=b);
1394 'at('at('diff(y,x,1),x = a),z = b);
1395 remove(y,dependency);
1396 done;
1398 /* Bug report ID: 3023978 - integrate(x^x+x,x) is wrong
1399  */
1400 integrate(x^x+x,x);
1401 'integrate(exp(x*log(x)),x)+x^2/2;
1403 /* Bug report ID: 2465066 - unsimplified result from integrate
1404  */
1405 matchdeclare(x, symbolp);
1406 done;
1407 (tellsimpafter('integrate(f(x),x), g(x)),done);
1408 done;
1409 integrate(5*f(x) + 7,x);
1410 5*g(x)+7*x;
1411 kill(rules);
1412 done;
1414 /* Bug report ID: 2789110 - solve, tan and atan depend on order of variables
1415  */
1416 solve(tan(x - atan(a/b)) = 0, x);
1417 [x = atan(a/b)];
1418 solve(tan(x - atan(b/a)) = 0, x);
1419 [x = atan(b/a)];
1421 /* Bug report ID: 1961494 - translated functions & values list
1422  */
1423 (kill(all), f():= x:2, translate(f));
1424 [f];
1425 f();
1427 values;
1428 [x];
1429 kill(f,x);
1430 done;
1431 /* The value of x has been removed. */
1435 /* Bug report ID: 3025038 - gruntz needs logexpand:true
1436  */
1437 gruntz( (x + 2^x) / 3^x, x, inf),logexpand:false;
1440 /* Bug report ID: 2977217 - maxima can not integrate x*exp(-1/2*(x-m)^2)
1441  */
1442 integrate(x*exp(-1/2*(x-m)^2),x);
1443 %i*(sqrt(2)*%i*gamma_incomplete(1,(m-x)^2/2)*(m-x)^2/(x-m)^2
1444    -%i*gamma_incomplete(1/2,(m-x)^2/2)*m*(m-x)/abs(x-m))/sqrt(2);
1446 /* Bug report ID: 2996542 - log(x) integration is incorrect
1447  */
1448 assume(a>0);
1449 [a>0];
1450 integrate(log(x),x,0,a);
1451 a*log(a)-a;
1452 forget(a>0);
1453 [a>0];
1455 /* Bug report ID: 3062883 - diff does not recognize indirect dependencies 
1456  *                          in expressions
1457  */
1459 depends([a,b],x,x,t);
1460 [a(x),b(x),x(t)];
1461 diff(-a,t);
1462 -'diff(a,x,1)*'diff(x,t,1);
1463 diff(a*b,t);
1464 a*'diff(b,x,1)*'diff(x,t,1)+'diff(a,x,1)*b*'diff(x,t,1);
1465 remove([a,b,x],dependency);
1466 done;
1468 /* Bug report ID: 3080397 - laplace(unit_step(-t),t,s) generates an error.
1469  */
1470 laplace(unit_step(-t),t,s);
1473 /* Bug report ID: 3081820 - lbfgs causes error
1475  * Still generates an error, but a different error that maxima
1476  * signals.
1477  */
1478 block([V:0.75, a:24, b:68, e],
1479   C(r) := 2*%pi*b*r^2 + 4*a*%pi*r + 2*b*V/r + a*V/(%pi*r^2),
1480   load(lbfgs),
1481   /* This should signal an error that we catch */
1482   e : errcatch(lbfgs(C(r), [r], [1], 1e-4, [1,0])),
1483   [e, error]);
1484 [[], ["Evaluation of gradient at ~M failed.  Bad initial point?~%", [0.0]]];
1486 /* Bug report ID: 875089 - defint(f(x)=g(x),x,0,1) -> false = false
1488  * We distribute defint more early in the code of bags to get a correct result.
1489  */
1490  defint(f(x)=g(x),x,0,1);
1491  'integrate(f(x),x,0,1)='integrate(g(x),x,0,1);
1493 /* Bug report ID: 2796194 - error doing a Fourier transform */
1494 (assume(equal(x,0)),done);
1495 done;
1497 errcatch(integrate(%pi*exp(-2*%pi*t)*exp(2*%pi*x*t*%i),t,minf,inf));
1500 error;
1501 ["defint: integral is divergent."];
1503 (forget(equal(x,0)),done);
1504 done;
1506 /* Bug reported on the mailing list
1507  * <http://www.math.utexas.edu/pipermail/maxima/2010/023024.html>
1508  * integrate(cos(2*x)*cos(x),x) is wrong.
1510  * Add a few more test that are similar to test the part of
1511  * monstertrig that deals with trig1(m*x)*trig2(n*x) where trig1 and
1512  * trig2 are either sin or cos.
1513  */
1514 integrate(cos(2*x)*cos(x),x);
1515 sin(3*x)/6+sin(x)/2;
1517 integrate(sin(2*x)*sin(x),x);
1518 sin(x)/2-sin(3*x)/6;
1520 integrate(cos(2*x)*sin(x),x);
1521 cos(x)/2-cos(3*x)/6;
1523 integrate(sin(x)*cos(2*x),x);
1524 cos(x)/2-cos(3*x)/6;
1526 /* Bug ID: 3111568 - subsequent calls to gradef hide variable dependencies
1527  */
1528 gradef(f,x,g);
1530 gradef(f,y,h);
1532 dependencies;
1533 [f(y,x)];
1534 kill(f);
1535 done;
1537 /* Bug ID: 3118770 - %edispflag:true causes a bug
1538  */
1539 %edispflag:true;
1540 true;
1541 integrate(x/(%e)^(2*x), x, 0, 1);
1542 1/4-3/(4*%e^2);
1543 reset(%edispflag);
1544 [%edispflag];
1546 /* Bug ID: 3067098 - The command timer for a Lisp function
1547  * Check that the function trisplit does not go away, when we collect
1548  * timing statistics for this function and call later kill(all).
1549  */
1550 timer(?trisplit);
1551 [?trisplit];
1552 kill(all);
1553 done;
1554 rectform(1+%i);
1555 1+%i;
1557 /* Bug ID: 3133916 - scanmap(minfactorial,a!) infinite loop
1558  */
1559 scanmap(minfactorial, a!);
1562 /* Bug ID: 3131324 - simplification of sqrt
1563  */
1564 sqrt(x^3)/sqrt(x^3);
1567 /* Bug ID: 1285104 - trigsimp and trigreduce & square roots
1568  */
1570 trigreduce(sqrt(r^2*sin(x)^2+r^2*cos(x)^2));
1571 abs(r);
1572 trigreduce(sqrt(r^2*sin(x)^2+r^2*cos(x)^2)),radexpand:all;
1574 radexpand:false;
1575 false;
1576 trigreduce(sqrt(r^2*sin(x)^2+r^2*cos(x)^2));
1577 sqrt(r^2);
1578 reset(radexpand);
1579 [radexpand];
1581 /* Bug ID: 917283 - Comment syntax confused
1582  * Show that nested comments work as expected.
1583  */
1584 a/*/**/*/+b;
1585 a+b;
1586 a/*/**/*/+/*/**/*/b;
1587 a+b;
1589 /* Bug ID: 3138054 -  bfloat problem / FIX -
1590  */
1591 exp(gamma(1/3)),bfloat;
1592 1.45696199392313b1;
1594 /* Bug ID: 3288989 - Lisp functions and linear display
1595  * Show that we do not get a Lisp error.
1596  */
1597 grind(?cdr([a,b,c]));
1598 done;
1600 /* Bug ID: 3291590 - Problems with fast arrays
1601  */
1602 (a:make_array(hashed), done);
1603 done;
1604 a[100]:100;
1605 100;
1606 a[x]:sin(x);
1607 sin(x);
1608 a[x*y]:x^2+y;
1609 y+x^2;
1611 /* Cutting out these two examples.
1612  * The ordering of the lists is different depending on the underlying Lisp.
1614  * arrayinfo(a);
1615  * [hash_table,1,100,x,x*y];
1616  * listarray(a);
1617  * [100,sin(x),y+x^2];
1618  */
1620 (f:make_array(functional, 'factorial, hashed), done);
1621 done;
1622 f[10];
1623 3628800;
1625 (kill(f), 0);
1628 (a: make_array(fixnum, 2, 2), done);
1629 done;
1630 listarray(a);
1631 [0, 0, 0, 0];
1633 use_fast_arrays:true;
1634 true;
1636 (array(a, any, 2, 2), done);
1637 done;
1638 arrayinfo(a);
1639 [declared, 2, [2, 2]];
1640 (array(a, fixnum, 2, 2), done);
1641 done;
1642 arrayinfo(a);
1643 [declared, 2, [2, 2]];
1644 (array(a, flonum, 2, 2), done);
1645 done;
1646 arrayinfo(a);
1647 [declared, 2, [2, 2]];
1648 (array(a, hashed), done);
1649 done;
1650 arrayinfo(a);
1651 [hash_table, 1];
1653 reset(use_fast_arrays);
1654 [use_fast_arrays];
1655 kill(a);
1656 done;
1658 /* Bug ID: 3247367 - expand returns unsimplified
1659  */
1660 sqrt(2)+sqrt(2);
1661 2^(3/2);
1662 sqrt(2)+sqrt(2)+sqrt(2);
1663 3*sqrt(2);
1664 sqrt(2)+sqrt(2)+sqrt(2)+sqrt(2);
1665 2^(5/2);
1666 sqrt(2)+sqrt(2)+sqrt(2)+sqrt(2)+sqrt(2);
1667 5*sqrt(2);
1669 2*sqrt(2)+3*sqrt(2);
1670 5*sqrt(2);
1671 3*sqrt(2)+2*sqrt(2);
1672 5*sqrt(2);
1673 3*sqrt(2)+2*sqrt(2)+sqrt(2);
1674 3*2^(3/2);
1675 sqrt(2)+3*sqrt(2)+2*sqrt(2);
1676 3*2^(3/2);
1678 sqrt(1/2)+sqrt(1/2);
1679 sqrt(2);
1680 sqrt(1/2)+sqrt(1/2)+sqrt(1/2);
1681 3/sqrt(2);
1682 sqrt(1/2)+sqrt(1/2)+sqrt(1/2)+sqrt(1/2);
1683 2^(3/2);
1684 sqrt(1/2)+sqrt(1/2)+sqrt(1/2)+sqrt(1/2)+sqrt(1/2);
1685 5/sqrt(2);
1687 2*sqrt(1/2)+3*sqrt(1/2);
1688 5/sqrt(2);
1689 3*sqrt(1/2)+2*sqrt(1/2);
1690 5/sqrt(2);
1691 3*sqrt(1/2)+2*sqrt(1/2)+sqrt(1/2);
1692 3*sqrt(2);
1693 sqrt(1/2)+3*sqrt(1/2)+2*sqrt(1/2);
1694 3*sqrt(2);
1696 2^(1/5)+2^(1/5);
1697 2^(6/5);
1698 2^(1/5)+2^(1/5)+2^(1/5);
1699 3*2^(1/5);
1700 2^(1/5)+2^(1/5)+2^(1/5)+2^(1/5);
1701 2^(11/5);
1702 2^(1/5)+2^(1/5)+2^(1/5)+2^(1/5)+2^(1/5);
1703 5*2^(1/5);
1705 (1/2)^(1/5)+(1/2)^(1/5);
1706 2^(4/5);
1707 (1/2)^(1/5)+(1/2)^(1/5)+(1/2)^(1/5);
1708 3/2^(1/5);
1709 (1/2)^(1/5)+(1/2)^(1/5)+(1/2)^(1/5)+(1/2)^(1/5);
1710 2^(9/5);
1711 (1/2)^(1/5)+(1/2)^(1/5)+(1/2)^(1/5)+(1/2)^(1/5)+(1/2)^(1/5);
1712 5/2^(1/5);
1714 2*(1/2)^(1/5)+3*(1/2)^(1/5);
1715 5/2^(1/5);
1716 3*(1/2)^(1/5)+2*(1/2)^(1/5);
1717 5/2^(1/5);
1719 2^sin(x)+2^sin(x);
1720 2^(sin(x)+1);
1721 2^sin(x)+2^sin(x)+2^sin(x);
1722 3*2^sin(x);
1723 2^sin(x)+2^sin(x)+2^sin(x)+2^sin(x);
1724 2^(sin(x)+2);
1725 2^sin(x)+2^sin(x)+2^sin(x)+2^sin(x)+2^sin(x);
1726 5*2^sin(x);
1728 (1/2)^sin(x)+(1/2)^sin(x);
1729 2^(1-sin(x));
1730 (1/2)^sin(x)+(1/2)^sin(x)+(1/2)^sin(x);
1731 3/2^sin(x);
1732 (1/2)^sin(x)+(1/2)^sin(x)+(1/2)^sin(x)+(1/2)^sin(x);
1733 2^(2-sin(x));
1734 (1/2)^sin(x)+(1/2)^sin(x)+(1/2)^sin(x)+(1/2)^sin(x)+(1/2)^sin(x);
1735 5/2^sin(x);
1737 (1-sqrt(5))^3-4*(1-sqrt(5))^2+8, expand;
1739 1/sqrt(2)+1/sqrt(2)+1/sqrt(2);
1740 3/sqrt(2);
1741 2^(9/5)+2^(4/5);
1742 3*2^(4/5);
1743 3*sqrt(2)+2*sqrt(2);
1744 5*sqrt(2);
1745 2*sqrt(2)+3*sqrt(2);
1746 5*sqrt(2);
1747 (1-sqrt(5))^3, expand;
1748 16-8*sqrt(5);
1749 p : z^3-2^(3/2)*%i*z^2-4*z^2+2^(5/2)*%i*z+2*z;
1750 z^3-2^(3/2)*%i*z^2-4*z^2+2^(5/2)*%i*z+2*z;
1751 divide(p, (z-2-sqrt(2)*%i),z);
1752 [z^2+(-sqrt(2)*%i-2)*z,0];
1753 2^a + 3*2^(a+1);
1754 7*2^a;
1755 2^(3/5)+2^(-2/5);
1756 3/2^(2/5);
1757 2^(3/5+x)+2^(-2/5+x);
1758 3*2^(x-2/5);
1760 /* -----------------------------------------------------------------------------
1761  * Bug ID: 1439566 - zerobern & bernpoly
1762  * Show that the option variable zerobern does not change the results.
1763  * -------------------------------------------------------------------------- */
1764 zerobern:false;
1765 false;
1766 bernpoly(x,3);
1767 x^3-3*x^2/2+x/2$
1768 bernpoly(x,5);
1769 x^5-5*x^4/2+5*x^3/3-x/6$
1770 reset(zerobern);
1771 [zerobern];
1773 /* Show that bern no longer fails with zerobern:false.
1775  * The compared values are from A&S p810, Table 23.2.
1776  */
1777 bern(28);
1778 -23749461029/870$
1779 bern(40);
1780 -261082718496449122051/13530$
1781 euler(16);
1782 19391512145$
1784 zerobern:false;
1785 false$
1786 bern(17);
1787 -7709321041217/510$
1788 bern(24);
1789 596451111593912163277961/282$
1790 euler(10);
1791 370371188237525$
1793 reset(zerobern);
1794 [zerobern]$
1796 /* -----------------------------------------------------------------------------
1797  * Bug ID: 2905929 - gcdex
1798  * -------------------------------------------------------------------------- */
1799 q0[2] : 6;
1801 ratdisrep(gcdex(x-7, x-8));
1802 [1,-1,1]$
1804 is(equal(gcdex(z^2-1, 0, z), [1,0,z^2-1]));
1805 true;
1807 is(equal(gcdex(0, z^2-1, z), [0, 1, z^2-1]));
1808 true;
1810 /* Examples from the Maxima Manual */
1811 is(equal(gcdex(x^2+1, x^3+4), [-(x^2+4*x-1)/17,(x+4)/17,1]));
1812 true$
1813 is(equal(gcdex(x*(y+1), y^2-1, x), [0,1/(y^2-1),1]));
1814 true$
1816 kill(q0);
1817 done;
1819 /* -----------------------------------------------------------------------------
1820  * Bug ID: 3389830 - Error break in rtest15 with linear display
1822  * Show that we do not get an error with grind for a prefix and a postfix
1823  * Operator when displaying the definition of a function or a macro.
1824  * This is a test of the function msz-mdef in grind.lisp.
1825  * -------------------------------------------------------------------------- */
1826 (postfix("f"), prefix("g"), done);
1827 done$
1828 grind("f"(x):=x);
1829 done$
1830 grind("f"(x)::=x);
1831 done$
1832 grind("g"(x):=x);
1833 done$
1834 grind("g"(x)::=x);
1835 done$
1837 kill ("f", "g");
1838 done;
1840 /* -----------------------------------------------------------------------------
1841  * Bug ID: 3396631 - equal terms produce different results
1842  * Correcting a bug in plusin revision 23.08.2011
1843  * -------------------------------------------------------------------------- */
1844 5*sqrt(5)+2*sqrt(3)+6*sqrt(5);
1845 11*sqrt(5)+2*sqrt(3)$
1846 5*sqrt(5)+4*sqrt(3)+6*sqrt(5);
1847 11*sqrt(5)+4*sqrt(3)$
1848 5*sqrt(5)+3*sqrt(5)+5*sqrt(3)+3*sqrt(3)+2*sqrt(75)+2*sqrt(45);
1849 14*sqrt(5)+2*3^(5/2)$
1850 5*sqrt(5)+5*sqrt(3)+3*sqrt(3)+2*sqrt(75)+2*sqrt(45)+3*sqrt(5);
1851 14*sqrt(5)+2*3^(5/2)$
1852 5*sqrt(5)+5*sqrt(3)+2*sqrt(75)+2*sqrt(45)+3*sqrt(5)+3*sqrt(3);
1853 14*sqrt(5)+2*3^(5/2)$
1855 /* -----------------------------------------------------------------------------
1856  * Bug ID 3437268: expand doesn't fully expand
1857  * Check the expected results after revision 14.11.2011 of simp.lisp.
1858  * -------------------------------------------------------------------------- */
1859 -3*a/sqrt(2)+sqrt(2)*a+a/sqrt(2);
1861 expand(-(atan((sqrt(2)*sin(8))/cos(8))+3*%pi)/sqrt(2)
1862         +atan((sqrt(2)*sin(8))/cos(8))/sqrt(2)
1863         +sqrt(2)*%pi
1864         +%pi/sqrt(2));
1868  * SF Bug: elliptic_e error observed - ID: 3438911
1870  * Test the bug by using the identity:
1871  * elliptic_e(%pi,m) = 2*elliptic_ec(m)
1873  * There's also a bug in elliptic_f.  We test this using the identity:
1874  * elliptic_f(%pi,m) = 2*elliptic_kc(m)
1875  */
1877 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
1878 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
1880 closeto(elliptic_e(float(%pi), .1) - 2*elliptic_ec(.1), 1e-15);
1881 true;
1882 closeto(elliptic_e(float(%pi), .9) - 2*elliptic_ec(.9), 1e-15);
1883 true;
1884 closeto(elliptic_e(bfloat(%pi), .1b0) - 2*elliptic_ec(.1b0), 1e-16);
1885 true;
1886 closeto(elliptic_e(bfloat(%pi), .9b0) - 2*elliptic_ec(.9b0), 1e-16);
1887 true;
1889 closeto(elliptic_f(float(%pi), .1) - 2*elliptic_kc(.1), 1e-15);
1890 true;
1891 closeto(elliptic_f(float(%pi), .9) - 2*elliptic_kc(.9), 1e-15);
1892 true;
1893 closeto(elliptic_f(bfloat(%pi), .1b0) - 2*elliptic_kc(.1b0), 1e-16);
1894 true;
1895 closeto(elliptic_f(bfloat(%pi), .9b0) - 2*elliptic_kc(.9b0), 1e-16);
1896 true;
1899  * Bug 3526111 - float erf (%i) not working
1900  */
1902 closeto(float(erf(%i)) - 1.650425758797543*%i, 1e-15);
1903 true;
1906  * Bug 3529992: Shi (sinh integral) wrong branch, integrate inconsistent
1907  */
1908 closeto(float(expintegral_shi(1/2) - 0.50699674981966719583), 3e-16);
1909 true;
1911 /* integrate changes k[0] --> k(0) - ID: 3530767 */
1912 integrate(x * (x^2 + k[0])/(1 + x^2),x);
1913 ((k[0]-1)*log(x^2+1))/2+x^2/2$
1915 /* polarform error on simple case - ID: 3517034 */
1916 polarform((a+1)/2 - a/2 - 1/2);
1919 polarform((a+1)/2 - a/2 - 0.5);
1920 0.0$
1922 polarform((a+1)/2 - a/2 - 0.5b0);
1923 0.0b0$
1926 /* #2531 Integration with inf */
1927 errcatch(integrate((1+1/x)^(1/2),x,1,inf));
1930 error;
1931 ["defint: integral is divergent."];
1936  * ID: 3440046: elliptic_f(0.5,1) signals error
1938  * Add a few more tests for invalid values.
1939  */
1940 closeto(elliptic_f(0.5,1)-elliptic_f(1/2,1), 1e-15);
1941 true;
1943 closeto(elliptic_f(0.5b0,1) - bfloat(elliptic_f(1/2,1)), 1b-16);
1944 true;
1946 errcatch(elliptic_f(2.0,1));
1949 errcatch(elliptic_f(2b0,1));
1953  * Bug 3428734:  integrate(bessel_y(1,z),z)  with ?z : 1
1954  */
1956 (?z:1, 0);
1959 integrate(bessel_y(1,z),z);
1960 -bessel_y(0,z);
1962 integrate(bessel_j(1,z),z);
1963 -bessel_j(0,z);
1965 integrate(bessel_k(1,z),z);
1966 -bessel_k(0,z);
1968 integrate(bessel_i(1,z),z);
1969 bessel_i(0,z);
1972  * Bug 3381301: log(-1.0b0) has small realpart
1973  */
1974 realpart(log(-1b0));
1978  * Bug 3559064: elliptic_f(2,1) is wrong.
1979  */
1980 errcatch(elliptic_f(2,1));
1984  * Bug 2528: A variable should be real if it is both real and complex
1985  */
1986 (declare(foo, real), declare(foo, complex), 0);
1989 [realpart(foo), imagpart(foo)];
1990 [foo, 0]$
1992 (kill(foo), 0);
1995 map(lambda([x], featurep(x, 'irrational)),[42,%pi,%phi,%e]);
1996 [false, true, true, true]$
1998 /* #2501 %pi/8 is definitely not an integer */
1999 integrate(log(cot(x)-1),x,0,%pi/4);
2000 (%i*li[2]((%i+1)/2)-%i*li[2](-((%i-1)/2)))/2 -(%i*(2*li[2](%i+1)-2*li[2](1-%i))+%pi*log(2))/4$
2002 integrate(log(cos(x)),x,0,%pi/2);
2003 %i*%pi^2/24-(6*%pi*log(4)+%i*%pi^2)/24$
2005 map(lambda([x], featurep(x, noninteger)),[sqrt(5),%pi,%pi/3, %pi/8,log(42),99/2013]);
2006 [true,true,true,true,true,true]$
2008 map(lambda([x], featurep(x, noninteger)),[0,1,2013]);
2009 [false,false,false]$
2011 /* #2583 sign error for integrate(x^(8*%i-1),x) */
2012 block([domain : 'real], integrate(x^(8*%i-1),x));
2013 -%i*x^(8*%i)/8$
2015 /* #2602: some-bfloatp and some-floatp recursed wrongly on rat expressions */
2016 ?some\-bfloatp(rat(1/2));
2017 false$
2019 /* #2594: Error in trigreduce for complicated expressions */
2020 subst(0, x, trigreduce(product(cos(k*x), k, 1, 8)));
2023 /* SF bug #2818: Problem with trigreduce */
2024 trigreduce(sin(1/8*%pi)*sin(3/8*%pi)*sin(5/8*%pi)*sin(7/8*%pi));
2025 1/8;
2027 /* #2591: Risch gives lisp error */
2028 (risch(asinh((z^2-1)/z)/z,z), 0);
2031 /* #2682: Function zeta fails numerically for large numbers */
2032 closeto(zeta(40.0) - 1, 1e-12);
2033 true$
2035 /* #2675 (1/3): Integration yields noun form with subscripted variable */
2036 integrate(exp(-(1+%i)*x[1]),x[1],0,inf);
2037 1/2-%i/2$
2039 /* #2688: %e^^A returns element by element exponent */
2040 is(%e^^matrix([1,2],[3,4]) = %e^matrix([1,2],[3,4]));
2041 false$
2043 /* #2676: Integral incorrect when variable is subscripted */
2044 integrate(x[1]*exp(x[1]), x[1]);
2045 exp(x[1])*(x[1]-1)$
2047 /* #2726: Integrate produces wrong answer for Gaussian Moments */
2048 block([tmp, m2726],
2049    assume(m2726+1 > 0),
2050    tmp : integrate(exp(-x^2/2)/sqrt(2*%pi) * x^m2726, x, -1/4, 1/4),
2051    tmp : sign(subst(m2726 = 4, tmp)),
2052    forget(m2726+1 > 0),
2053    tmp);
2054 pos$
2056 /* # 2697: Inconsistent handling of Greek symbols */
2057 integrate(y(%theta)=sin(%theta),%theta,%theta[0], %theta[1]);
2058 'integrate(y(%theta),%theta,%theta[0],%theta[1]) = cos(%theta[0])-cos(%theta[1])$
2060 integrate(y(t)=sin(t),t,t[0], t[1]);
2061 'integrate(y(t),t,t[0],t[1]) = cos(t[0])-cos(t[1])$
2063 integrate(y(tau)=sin(tau),tau,tau[0], tau[1]);
2064 'integrate(y(tau),tau,tau[0],tau[1]) = cos(tau[0])-cos(tau[1])$
2066 integrate ([foo(x), bar(x)], x, x[1], x[2]);
2067 ['integrate (foo(x), x, x[1], x[2]), 'integrate (bar(x), x, x[1], x[2])];
2069 /* # 2738: Integrate encountered a Lisp error: The value 2 is not of type LIST. */
2071 (kill (x, y, I, J),
2072  x(t):=2*cos(t),
2073  y(t):=2*sin(t),
2074  I : (x(t)+y(t)^2)*sqrt(diff(x(t),t)^2+diff(y(t),t)^2),
2075  J : integrate (I, t));
2076 4*(t-tan(t)/(tan(t)^2+1))+4*sin(t)$
2078 trigsimp (diff (J, t) - I);
2080 /* bug #2980: Infinite recursion with (e: log(e), rectform(e)) */
2081 block ([e: log(e)], rectform(e));
2082 log(abs(e)) + %i*atan2(0, e)$
2084 /* bug #2159: integration_with_logabs ("integrate(tan(x),x);" etc. do not take "logabs" flag into account) */
2085 integrate([tan(x),csc(x),sec(x),cot(x),tanh(x),coth(x),csch(x)],x);
2086 [log(sec(x)),-log(csc(x)+cot(x)),log(sec(x)+tan(x)),log(sin(x)),log(cosh(x)),log(sinh(x)),log(tanh(x/2))]$
2087 integrate([tan(x),csc(x),sec(x),cot(x),tanh(x),coth(x),csch(x)],x),logabs;
2088 [log(abs(sec(x))),-log(abs(csc(x)+cot(x))),log(abs(sec(x)+tan(x))),log(abs(sin(x))),log(cosh(x)),log(abs(sinh(x))),log(abs(tanh(x/2)))]$
2090 /* Bug #3075: #3075 answer "3*false" from "integrate(3*asinh(x),x,-inf,inf)" */
2091 /* We can't check for the correct answer zero (yet), because Maxima can't solve integrate(asinh(x),x,-inf,inf). */
2092 is(integrate(3*asinh(x),x,-inf,inf)#3*false);
2093 true$
2096  * Bug #3056: exp(1b19) signals error that 1b19 doesn't have enough
2097  * precision to compute its integer part.  Add test for this and also
2098  * the original test from commit 576c7508.)
2099  */
2101 /* Don't really care what the answer is as long as we don't signal an error */
2102 is(errcatch(exp(1b19)) # []);
2103 true;
2105 /* Test from the commit 576c7508 */
2106 ceiling((207300647060*%e^-563501581931)/(403978495031*%e^-1098127402131));
2107 ceiling((207300647060*%e^534625820200)/403978495031);
2109 /* Bug #3098: numerical evaluation of li[3] */
2110 li[3](0.0);
2111 0.0;
2113 closeto(li[3](1.0) - zeta(3.0), 1e-15);
2114 true;
2116 closeto(li[3](0.5) - li[3](1/2), 1e-15);
2117 true;
2119 closeto(li[3](-2.0) - float(subst(z=-2.0, li[3](1/z) - log(-z)/6*(log(-z)^2+%pi^2))), 1e-15);
2120 true;
2122 closeto(abs(li[3](2.0) - expand(float(subst(z=2.0, li[3](1/z) - log(-z)/6*(log(-z)^2+%pi^2))))), 1e-15);
2123 true;
2125 li[2](-1);
2126 -%pi^2/12;
2128 li[2](1);
2129 %pi^2/6;
2131 li[2](1/2);
2132 %pi^2/12 - log(2)^2/2;
2134 closeto(li[2](0.5) - (%pi^2/12 - log(2)^2/2), 1.1103e-16), numer;
2135 true;
2137 closeto(li[2](2.0) - (%pi^2/4 - %i*%pi*log(2)), 2.513e-15), numer;
2138 true;
2140 closeto(li[2]((1-sqrt(5))/2) - (log((sqrt(5)-1)/2)^2/2-%pi^2/15), 1.1103e-16), numer;
2141 true;
2143 /* Catalan's constant: 0.915965594... */
2144 closeto(li[2](1.0*%i) - (-%pi^2/48 + %i*0.915965594177219015054603514932384110774149374281672134266), 1.3878e-15), numer;
2145 true;
2147 closeto(li[2](1.0-%i) - (%pi^2/16-%i*0.915965594177219015054603514932384110774149374281672134266 - %pi*%i*log(2)/4), 4.5e-16), numer;
2148 true;
2150 /* Make sure li[3](1/7),numer returns a float and not a bfloat */
2151 ?floatp(li[3](float(1/7)));
2152 true;
2154 ?floatp(realpart(li[2](1.0+1.0*%i)));
2155 true;
2157 closeto(li[3](0.5) - float((7*zeta(3))/8+log(2)^3/6-(%pi^2*log(2))/12), 1.1103e-16);
2158 true;
2160 closeto(float(li[3](exp(%pi*%i/3))) - float(zeta(3)/3 + %i*5*%pi^3/162), 5.6611e-16);
2161 true;
2164  * li[3] should not return a small imaginary part for this value.
2165  * (There are others, but the solution fixes those as well.
2166  */
2167 imagpart(li[3](-.862));
2171  * Likewise li[s](-1.5) was returning small imaginary parts.  In fact,
2172  * this was true for -2 < z < -1 because the log series was used to
2173  * compute the value.  Check that we return exactly zero now for
2174  * select values of s
2175  */
2177 makelist(imagpart(li[s](-1.5)), s, 4, 10);
2178 [0, 0, 0, 0, 0, 0, 0];  
2180 /* Bug 3582: numerical evaluation of li[2] */
2181 /* this gave an infinite loop */
2182 closeto(li[2]( 0.5 + %i* 0.6) - ( 0.7535498331267791 * %i + 0.4081870726952078 ), 1e-15);
2183 true;
2185 /* Bug 3112: zeta(n) for negative even n is inaccurate */
2186 zeta(-4.0);
2187 0.0;
2189 zeta(-6b0);
2190 0b0;
2192 /* Bug 3105: li[s](1.0) doesn't simplify */
2193 closeto(li[4](1.0)-li[4](1), 1.11022e-16);
2194 true;
2196 closeto(li[4](-1.0)-li[4](-1), 1.11022e-16);
2197 true;
2199 closeto(li[4](1b0)-li[4](1), 10^-fpprec);
2200 true;
2202 closeto(li[4](-1b0)-li[4](-1), 10^-fpprec);
2203 true;
2207  * More numeric test for polylog function li[s](z).
2209  * The reference values were obtained from
2210  * http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=PolyLog
2211  */
2212 fpprec:24;
2215 closeto(li[4](0.25) - 0.25411619074634353405967371315352,
2216   5.5512e-17);
2217 true;
2219 closeto(li[4](0.25b0) - 0.25411619074634353405967371315352b0,
2220   2.06796b-25);
2221 true;
2223 closeto(li[4](0.75) - 0.79222102797282777952948578955736,
2224   2.2205e-16);
2225 true;
2227 closeto(li[4](0.75b0) - 0.79222102797282777952948578955736b0,
2228   6.20386b-25);
2229 true;
2231 closeto(li[4](1.5) - (1.7347570807760620737768805117515 - 0.0349027048283367002627421237287*%i),
2232   2.2205e-16);
2233 true;
2235 closeto(li[4](1.5b0) - (1.7347570807760620737768805117515b0 - 0.0349027048283367002627421237287b0*%i),
2236   4.14398b-25);
2237 true;
2239 closeto(li[4](10.0) - (9.6140263862742968515251940747859 - 6.3921313179656069159740055708257*%i),
2240   5.618e-15);
2241 true;
2243 closeto(li[4](10.0b0) - (9.6140263862742968515251940747859b0 - 6.3921313179656069159740055708257b0*%i),
2244   1.9364b-23);
2245 true;
2247 closeto(li[4](-5.0) - (-4.1064679790949702621073505378164),
2248   8.8818e-16);
2249 true;
2251 closeto(li[4](-5b0) - (-4.1064679790949702621073505378164b0),
2252   4.9631b-24);
2253 true;
2255 closeto(li[4](2.0*%i) - (-0.2113943747614829764326347997923 + 1.9289340331586646502356787803360*%i),
2256   7.22e-16);
2257 true;
2259 closeto(li[4](2b0*%i) - (-0.2113943747614829764326347997923b0 + 1.9289340331586646502356787803360b0*%i),
2260   1.9788b-24);
2261 true;
2263 closeto(li[5](0.25) - 0.25202158817857420100669519623555,
2264   5.55112e-17);
2265 true;
2267 closeto(li[5](0.25b0) - 0.25202158817857420100669519623555b0,
2268   5.16988b-25);
2269 true;
2271 closeto(li[5](0.75) - 0.76973541059975738097269173152535,
2272   1.111e-16);
2273 true;
2275 closeto(li[5](0.75b0) - 0.76973541059975738097269173152535b0,
2276   2.06796b-25);
2277 true;
2279 closeto(li[5](1.5) - (1.5961739456813534102689069143338 - 0.0035379572466222227823786042761*%i),
2280   4.33681e-19);
2281 true;
2283 closeto(li[5](1.5b0) - (1.5961739456813534102689069143338b0 - 0.0035379572466222227823786042761b0*%i),
2284   4.13594b-25);
2285 true;
2287 closeto(li[5](10.0) - (11.2390407376112991620107110964539 - 3.6796065713019972004384472107791*%i),
2288   9.93014e-15);
2289 true;
2291 closeto(li[5](10.0b0) - (11.2390407376112991620107110964539b0 - 3.6796065713019972004384472107791b0*%i),
2292   2.18382b-23);
2293 true;
2295 closeto(li[5](-5.0) - (-4.4800824065112010228046981292686),
2296   1e-16);
2297 true;
2299 closeto(li[5](-5b0) - (-4.4800824065112010228046981292686b0),
2300   6.61745b-24);
2301 true;
2303 closeto(li[5](2.0*%i) - (-0.1139660114783041974283299769126 + 1.9734617121122917650272273558071*%i),
2304   6.98e-16);
2305 true;
2307 closeto(li[5](2b0*%i) - (-0.1139660114783041974283299769126b0 + 1.9734617121122917650272273558071b0*%i),
2308   6.42086b-25);
2309 true;
2311 closeto(li[20](5.0) - (5.0000238783147176199891129814336 - 2.182054400942151422e-13*%i),
2312   2.41797e-15);
2313 true;
2315 closeto(li[20](5b0) - (5.0000238783147176199891129814336b0 - 2.182054400942151422b-13*%i),
2316   6.76366b-24);
2317 true;
2320  * Bug 3966: li[s](1) = zeta(s)
2321  */
2322 makelist(li[s+1/2](1), s, 1, 5);
2323 [zeta(3/2), zeta(5/2), zeta(7/2), zeta(9/2), zeta(11/2)];
2325 li[7/3](1);
2326 zeta(7/3);
2328 li[-5/2](1);
2329 li[-5/2](1);
2331 li[-5/2](1.0);
2332 li[-5/2](1.0);
2334 /* Bug #3194: No simplification of "tan(x+n*%pi)" and "cot(x+n*%pi)" with "n" being a declared integer */
2335 /* Also test that the "modulus 1/2" check in "%piargs-tan/cot" works correctly. */
2336 declare(i, integer, ei, even, oi, odd);
2337 done$
2338 [tan(x+i*%pi),cot(x+i*%pi)];
2339 [tan(x),cot(x)]$
2340 [tan(x+ei*%pi),cot(x+ei*%pi)];
2341 [tan(x),cot(x)]$
2342 [tan(x+oi*%pi),cot(x+oi*%pi)];
2343 [tan(x),cot(x)]$
2344 [tan(x-3/2*%pi),tan(x-1/2*%pi),tan(x+1/2*%pi),tan(x+3/2*%pi)];
2345 [-cot(x),-cot(x),-cot(x),-cot(x)]$
2346 [cot(x-3/2*%pi),cot(x-1/2*%pi),cot(x+1/2*%pi),cot(x+3/2*%pi)];
2347 [-tan(x),-tan(x),-tan(x),-tan(x)]$
2348 kill(i, ei, oi);
2349 done$
2351 /* Bug #3148: sign can't figure out sign(a - b) but it knows sign(b - a) where a and b are exponentials */
2352 assume(t>0);
2353 [t>0]$
2354 sign(2^(500000*t)-2^(500007*t));
2355 neg$
2356 sign(2^(500007*t)-2^(500000*t));
2357 pos$
2358 forget(t>0);
2359 [t>0]$
2361 /* Bug #3246: Integrating u'(x) * f(u(x) + c) fails for f any inverse trigonometric/hyperbolic function */
2362 /* Also make sure that all the antiderivatives of (inverse) trigonometric/hyperbolic functions are correct. */
2363 integrate(map(lambda([f], diff(u(x), x) * f(u(x) + c)), [asin,acos,atan,acsc,asec,acot,asinh,acosh,atanh,acsch,asech,acoth]), x);
2364 [(u(x)+c)*asin(u(x)+c)+sqrt(1-(u(x)+c)^2),(u(x)+c)*acos(u(x)+c)-sqrt(1-(u(x)+c)^2),(u(x)+c)*atan(u(x)+c)-log((u(x)+c)^2+1)/2,log(sqrt(1-1/(u(x)+c)^2)+1)/2-log(1-sqrt(1-1/(u(x)+c)^2))/2+(u(x)+c)*acsc(u(x)+c),-log(sqrt(1-1/(u(x)+c)^2)+1)/2+log(1-sqrt(1-1/(u(x)+c)^2))/2+(u(x)+c)*asec(u(x)+c),log((u(x)+c)^2+1)/2+(u(x)+c)*acot(u(x)+c),(u(x)+c)*asinh(u(x)+c)-sqrt((u(x)+c)^2+1),(u(x)+c)*acosh(u(x)+c)-sqrt((u(x)+c)^2-1),log(1-(u(x)+c)^2)/2+(u(x)+c)*atanh(u(x)+c),log(sqrt(1/(u(x)+c)^2+1)+1)/2-log(sqrt(1/(u(x)+c)^2+1)-1)/2+(u(x)+c)*acsch(u(x)+c),(u(x)+c)*asech(u(x)+c)-atan(sqrt(1/(u(x)+c)^2-1)),log(1-(u(x)+c)^2)/2+(u(x)+c)*acoth(u(x)+c)]$
2365 radcan(exponentialize(map(lambda([f], f(x) - diff(integrate(f(x), x), x)), [sin,cos,tan,csc,sec,cot,sinh,cosh,tanh,csch,sech,coth,asin,acos,atan,acsc,asec,acot,asinh,acosh,atanh,acsch,asech,acoth]))), domain:complex;
2366 [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]$
2368 /* Bug #3258: integrate's answer contains diff operator with negative order */
2369 freeof(diff, integrate(%e^(c*x-r^2/(4*x))/x^(3/2),x));
2370 true$
2372 /* Bug #3220: create_list doesn't bind variables properly */
2373 errcatch (create_list (%i, %i, [1, 2, 3]));
2376 create_list (bfloat (1 - 10^-50) - 1, fpprec, [16, 100]);
2377 [0.0b0, -1.0b-50]$
2379 create_list (bfloat (3.14), fpprec, 1, 3);
2380 [3.0b0, 3.1b0, 3.14b0]$
2382 /* Bug #3356: sign(nz * nz) = nz */
2383 assume(apos > 0, aneg < 0, apz >= 0, anz <= 0, notequal(apn, 0), bpos > 0, bneg < 0, bpz >= 0, bnz <= 0, notequal(bpn, 0));
2384 [apos > 0, aneg < 0, apz >= 0, anz <= 0, notequal(apn, 0), bpos > 0, bneg < 0, bpz >= 0, bnz <= 0, notequal(bpn, 0)]$
2385 map('sign,
2386         [
2387                 apos * bpos, apos * bneg, apos * bpz, apos * bnz, apos * bpn, apos * bpnz,
2388                 aneg * bpos, aneg * bneg, aneg * bpz, aneg * bnz, aneg * bpn, aneg * bpnz,
2389                 apz * bpos, apz * bneg, apz * bpz, apz * bnz, apz * bpn, apz * bpnz,
2390                 anz * bpos, anz * bneg, anz * bpz, anz * bnz, anz * bpn, anz * bpnz,
2391                 apn * bpos, apn * bneg, apn * bpz, apn * bnz, apn * bpn, apn * bpnz,
2392                 apnz * bpos, apnz * bneg, apnz * bpz, apnz * bnz, apnz * bpn, apnz * bpnz
2393         ]
2396         'pos, 'neg, 'pz, 'nz, 'pn, 'pnz,
2397         'neg, 'pos, 'nz, 'pz, 'pn, 'pnz,
2398         'pz, 'nz, 'pz, 'nz, 'pnz, 'pnz,
2399         'nz, 'pz, 'nz, 'pz, 'pnz, 'pnz,
2400         'pn, 'pn, 'pnz, 'pnz, 'pn, 'pnz,
2401         'pnz, 'pnz, 'pnz, 'pnz, 'pnz, 'pnz
2403 forget(apos > 0, aneg < 0, apz >= 0, anz <= 0, notequal(apn, 0), bpos > 0, bneg < 0, bpz >= 0, bnz <= 0, notequal(bpn, 0));
2404 [apos > 0, aneg < 0, apz >= 0, anz <= 0, notequal(apn, 0), bpos > 0, bneg < 0, bpz >= 0, bnz <= 0, notequal(bpn, 0)]$
2406 /* test cases mentioned in bug report #3356 */
2408 (kill (a, b), assume(a<=0, b<=0), sign(a*b));
2411 (kill (q, r), assume(q<0, r<0), sign(q*r));
2412 pos;
2414 sign(a*r);
2417 sign(a*-2);
2420 sign(a*q*r);
2423 sign(a*b*q);
2426 (forget (a <= 0, b <= 0, q < 0, r < 0), 0);
2429 /* Bug #3403: Function/lambda parameters declared constant cause error */
2430 declare(c, constant);
2431 done$
2432 emptyp(errcatch(f(c) := c));
2433 false$
2434 emptyp(errcatch(lambda([c], c)));
2435 false$
2436 emptyp(errcatch(f(%e) := %e));
2437 true$
2438 emptyp(errcatch(lambda([%e], %e)));
2439 true$
2440 kill(c, f);
2441 done$
2443 /*        Bug #3009: factoring exponentials causes MAKE-ARRAY error */
2444 /* a.k.a. Bug #3146: max() runs out of memory when comparing exponential functions */
2445 sign(((-1)-%e^-(4.075321792706671E-4*d_1))*%e^-(4.075321792706671E-4*d_2)+%e^-(8.150643585413342E-4*d_2)-%e^-(4.075321792706671E-4*d_1)+%e^-(8.150643585413342E-4*d_1)+1);
2446 pnz$
2447 factor(a^1000000+a+1);
2448 a^1000000+a+1$
2449 max(2.0325-6.9825*%e^-(492380.0*t),2.103-7.053*%e^-(406810.0*t));
2450 'max(2.0325-6.9825*%e^-(492380.0*t),2.103-7.053*%e^-(406810.0*t))$
2451 sign(exp(-492380*t)-2793);
2452 pnz$
2454 /* Bug #3147: sign of expressions with exponents crashes */
2455 assume(t>0);
2456 [t>0]$
2457 sign(2^(500005*t)-2^(500001*t));
2458 pos$
2459 kill(t);
2460 done$
2462 /* Test factor_max_degree. */
2463 block([factor_max_degree : 0], [factor(x^2+2*x+1), factor(x^2-1), factor(x^2+3*x+2), factor(x^2000+x^1999)]);
2464 [(x+1)^2, (x-1)*(x+1), (x+1)*(x+2), x^1999*(x+1)]$
2465 block([factor_max_degree : 1, factor_max_degree_print_warning : false], [factor(x^2+2*x+1), factor(x^2-1), factor(x^2+3*x+2), factor(x^2000+x^1999)]);
2466 [(x+1)^2, x^2-1, x^2+3*x+2, x^1999*(x+1)]$
2467 block([factor_max_degree : 2], [factor(x^2+2*x+1), factor(x^2-1), factor(x^2+3*x+2), factor(x^2000+x^1999)]);
2468 [(x+1)^2, (x-1)*(x+1), (x+1)*(x+2), x^1999*(x+1)]$
2470 /*        Bug #2928: Loops forever: ev(ratexpand(1/(x^(2/3)+1)), algebraic:true); */
2471 /* a.k.a. Bug #2994: Infinite loop: ratsimp(x/(sqrt(x^(2/3)+1)),x),algebraic; */
2472 /* a.k.a. Bug #3419: Endless loop in BPROG (rat(1/(x^(2/3)+1)), algebraic) */
2474 /* examples present in previous version of rtest16 */
2476 string (ev (rat (1/(x^(1/3)+1)), algebraic));
2477 "((x^(1/3))^2-x^(1/3)+1)/(x+1)";
2479 string (ev (rat (1/(x^(2/3)+1)), algebraic));
2480 "(x^(1/3)*x-(x^(1/3))^2+1)/(x^2+1)";
2482 string (ev (rat(1/((x-1)^(2/3)+1)), algebraic));
2483 "((x-1)^(1/3)*x-((x-1)^(1/3))^2-(x-1)^(1/3)+1)/(x^2-2*x+2)";
2485 string (ev (rat(x/(sqrt(x^(2/3)+1))), algebraic));
2486 "(sqrt(x^(2/3)+1)*x^(1/3)*x^2+(-(sqrt(x^(2/3)+1)*(x^(1/3))^2)+sqrt(x^(2/3)+1))*x)/(x^2+1)";
2488 /* additional examples courtesy of M. Talon */
2490 string (ev (rat(1/(x^(2/5)+1)), algebraic));
2491 "(((x^(1/5))^3-x^(1/5))*x+(x^(1/5))^4-(x^(1/5))^2+1)/(x^2+1)";
2493 /* Maxima cannot yet handle this next one (incomplete factorization, gets one factor but misses the other)
2495 rat(1/(x^(3/5)+x^(1/3))),algebraic;
2496  */
2498 string (ev (rat(1/(x^(4/5)+1)), algebraic));
2499 "(x^(1/5)*x^3-(x^(1/5))^2*x^2+(x^(1/5))^3*x-(x^(1/5))^4+1)/(x^4+1)";
2501 /* Bug #3431: error system variable holds unsimplified list, causing errors to be repeated when trying to access it */
2502 errcatch(0^0);
2504 errcatch(print(error), true);
2505 [true]$
2507 /* SF bug #3458: "does not interpret addcol as it should but instead interprets same as addrow" */
2509 addcol (matrix (), [11, 22, 33]);
2510 matrix ([11], [22], [33]);
2512 /* additional tests for addcol and addrow */
2514 addcol (matrix (), [11, 22, 33], [111, 222, 333]);
2515 matrix ([11, 111], [22, 222], [33, 333]);
2517 addcol (matrix (), [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2518 matrix ([11, 111, 1111], [22, 222, 2222], [33, 333, 3333]);
2520 addcol (matrix ([1, 2, 3]), [11, 22, 33]);
2521 matrix ([1, 2, 3, 11, 22, 33]);
2523 addcol (matrix ([1, 2, 3]), [11, 22, 33], [111, 222, 333]);
2524 matrix ([1, 2, 3, 11, 22, 33, 111, 222, 333]);
2526 addcol (matrix ([1, 2, 3]), [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2527 matrix ([1, 2, 3, 11, 22, 33, 111, 222, 333, 1111, 2222, 3333]);
2529 addcol (matrix ([1], [2], [3]), [11, 22, 33]);
2530 matrix ([1, 11], [2, 22], [3, 33]);
2532 addcol (matrix ([1], [2], [3]), [11, 22, 33], [111, 222, 333]);
2533 matrix ([1, 11, 111], [2, 22, 222], [3, 33, 333]);
2535 addcol (matrix ([1], [2], [3]), [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2536 matrix ([1, 11, 111, 1111], [2, 22, 222, 2222], [3, 33, 333, 3333]);
2538 addrow (matrix (), [11, 22, 33]);
2539 matrix ([11, 22, 33]);
2541 addrow (matrix (), [11, 22, 33], [111, 222, 333]);
2542 matrix ([11, 22, 33], [111, 222, 333]);
2544 addrow (matrix (), [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2545 matrix ([11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2547 addrow (matrix ([1, 2, 3]), [11, 22, 33]);
2548 matrix ([1, 2, 3], [11, 22, 33]);
2550 addrow (matrix ([1, 2, 3]), [11, 22, 33], [111, 222, 333]);
2551 matrix ([1, 2, 3], [11, 22, 33], [111, 222, 333]);
2553 addrow (matrix ([1, 2, 3]), [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2554 matrix ([1, 2, 3], [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2557  * Seems like these should succeed, but they fail
2558  * with error about "incompatible structure", although similar calls succeed.
2559  * Just omit these tests unless addrow/addcol are ever revised for greater consistency.
2561 addrow (matrix ([1], [2], [3]), [11, 22, 33]);
2562 matrix ([1], [2], [3], [11], [22], [33]);
2564 addrow (matrix ([1], [2], [3]), [11, 22, 33], [111, 222, 333]);
2565 matrix ([1], [2], [3], [11], [22], [33], [111], [222], [333]);
2567 addrow (matrix ([1], [2], [3]), [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2568 matrix ([1], [2], [3], [11], [22], [33], [111], [222], [333], [1111], [2222], [3333]);
2570  */
2572 addcol (matrix (), matrix ([11], [22], [33]));
2573 matrix ([11], [22], [33]);
2575 addcol (matrix (), matrix ([11, 22, 33]));
2576 matrix ([11, 22, 33]);
2578 addcol (matrix (), matrix ([11], [22], [33]), matrix ([111], [222], [333]));
2579 matrix ([11, 111], [22, 222], [33, 333]);
2581 addcol (matrix (), matrix ([11, 22, 33]), matrix ([111, 222, 333]));
2582 matrix ([11, 22, 33, 111, 222, 333]);
2584 addcol (matrix (), matrix ([11], [22], [33]), matrix ([111], [222], [333]), matrix ([1111], [2222], [3333]));
2585 matrix ([11, 111, 1111], [22, 222, 2222], [33, 333, 3333]);
2587 addcol (matrix (), matrix ([11, 22, 33]), matrix ([111, 222, 333]), matrix ([1111, 2222, 3333]));
2588 matrix ([11, 22, 33, 111, 222, 333, 1111, 2222, 3333]);
2590 addcol (matrix ([1, 2, 3]), matrix ([11, 22, 33]));
2591 matrix ([1, 2, 3, 11, 22, 33]);
2593 addcol (matrix ([1, 2, 3]), matrix ([11, 22, 33]), matrix ([111, 222, 333]));
2594 matrix ([1, 2, 3, 11, 22, 33, 111, 222, 333]);
2596 addcol (matrix ([1, 2, 3]), matrix ([11, 22, 33]), matrix ([111, 222, 333]), matrix ([1111, 2222, 3333]));
2597 matrix ([1, 2, 3, 11, 22, 33, 111, 222, 333, 1111, 2222, 3333]);
2599 addrow (matrix (), matrix ([11], [22], [33]));
2600 matrix ([11], [22], [33]);
2602 addrow (matrix (), matrix ([11, 22, 33]));
2603 matrix ([11, 22, 33]);
2605 addrow (matrix (), matrix ([11], [22], [33]), matrix ([111], [222], [333]));
2606 matrix ([11], [22], [33], [111], [222], [333]);
2608 addrow (matrix (), matrix ([11, 22, 33]), matrix ([111, 222, 333]));
2609 matrix ([11, 22, 33], [111, 222, 333]);
2611 addrow (matrix (), matrix ([11], [22], [33]), matrix ([111], [222], [333]), matrix ([1111], [2222], [3333]));
2612 matrix ([11], [22], [33], [111], [222], [333], [1111], [2222], [3333]);
2614 addrow (matrix (), matrix ([11, 22, 33]), matrix ([111, 222, 333]), matrix ([1111, 2222, 3333]));
2615 matrix ([11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2617 addrow (matrix ([1, 2, 3]), [11, 22, 33]);
2618 matrix ([1, 2, 3], [11, 22, 33]);
2620 addrow (matrix ([1, 2, 3]), [11, 22, 33], [111, 222, 333]);
2621 matrix ([1, 2, 3], [11, 22, 33], [111, 222, 333]);
2623 addrow (matrix ([1, 2, 3]), [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2624 matrix ([1, 2, 3], [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]);
2626 /* addrow/addcol that should fail */
2628 errcatch (addrow (matrix ([1], [2], [3]), [11, 22, 33]));
2632  * Seems like these should fail, but they succeed.
2633  * Just omit these tests unless addrow/addcol are ever revised for greater consistency.
2635 errcatch (addrow (matrix ([1], [2], [3]), matrix ([11, 22, 33])));
2638 errcatch (addrow (matrix ([1, 2, 3]), matrix ([1], [2], [3])));
2641 errcatch (addcol (matrix ([1, 2, 3]), matrix ([1], [2], [3])));
2644 errcatch (addcol (matrix ([1], [2], [3]), matrix ([11, 22, 33])));
2647  */
2649 /* Bug #3532: Integrator doesn't use a new variable internally, causing facts on the original variable to be used for the substitution variable */
2650 is(integrate(cos(x) * abs(cos(x)), x, 0, %pi) = %pi / 2);
2651 false$
2653 /* Bug found on mailing list thread "Wrong integral, but antiderivative and limits are correct" */
2654 integrate(x^8/(5+x^9)^2, x, 0, inf);
2655 1/45$
2657 /* Test case for bug mentioned at https://sourceforge.net/p/maxima/mailman/message/8488105/ */
2658 is(abs(integrate(1/(x^4+1), x, 0, 1/2) - 0.49396) < 0.00001);
2659 true$
2661 /* Another bug found on mailing list thread "Wrong integral, but antiderivative and limits are correct" */
2662 freeof(?xz, residue(plog(x)/(x-ratsimp(rectform(%e^(%i*%pi/9)))), x, %e^(%i*%pi/9)));
2663 true$
2665 /* Bug #3562: "integrate(1/(1+tan(x)), x, 0, %pi/2) gives complex result, should be %pi/4" */
2666 integrate(1/(1+tan(x)), x, 0, %pi/2);
2667 %pi/4$
2669 /* Also from bug #3562 */
2670 limit(log((x^2+1)/x^2)-2*log((x-1)/x), x, 0, minus);
2673 /* #3787 fixnump checks in simpexpt */
2674 (declare(n,integer),0);
2677 (-1)^((2^62-1)*n);
2678 (-1)^n$
2680 (-1)^((2^62+1)*n);
2681 (-1)^n$
2683 (-1)^((2^62+2)*n);
2686 (remove(n,integer),0);
2690  * Check sech/csch don't overflow for large numbers.  Use relative
2691  *  error to test for accuracy.
2692  */
2694 closetorel(actual, expected, tol) := block([numer:true, relerr: abs(actual-expected)/expected], if (relerr < tol) then true else relerr);
2695 closetorel(actual, expected, tol) := block([numer:true, relerr: abs(actual-expected)/expected], if (relerr < tol) then true else relerr);
2697 closeto(sech(715.0), 6.033b-311);
2698 true;
2699 closeto(csch(715.0), 6.033b-311);
2700 true;
2704  * Check acsch doesn't overflow for small numbers.  But some lisps
2705  * like clisp don't support denormals, so we have to be careful.
2707  */
2709 if (1e-309 = 0.0) then
2710     closetorel(acsch(1e-300), acsch(bfloat(1e-300)), 1.1375b-16)
2711   else
2712     closetorel(acsch(1e-309), acsch(bfloat(1b-309)), 7.1559b-17);
2713 true;
2715 /* SF bug #3825: "apply('forget, facts()) gives Lisp error" */
2717 (kill (all), declare (F, increasing));
2718 done;
2720 facts ();
2721 [kind (F, increasing)];
2723 featurep (F, increasing);
2724 true;
2726 apply ('forget, facts ());
2727 [kind (F, increasing)];
2729 facts ();
2732 featurep (F, increasing);
2733 false;
2735 declare ([F, G], rational, [H, I], irrational);
2736 done;
2738 facts ();
2739 [kind (F, rational), kind (G, rational), kind (H, irrational), kind (I, irrational)];
2741 [featurep (F, rational),
2742  featurep (G, rational),
2743  featurep (H, irrational),
2744  featurep (I, irrational)];
2745 [true, true, true, true];
2747 forget (kind (G, rational), kind (I, irrational));
2748 [kind (G, rational), kind (I, irrational)];
2750 facts ();
2751 [kind (F, rational), kind (H, irrational)];
2753 forget (kind (F, rational), kind (H, irrational));
2754 [kind (F, rational), kind (H, irrational)];
2756 facts ();
2759 [featurep (F, rational),
2760  featurep (G, rational),
2761  featurep (H, irrational),
2762  featurep (I, irrational)];
2763 [false, false, false, false];
2765 /* Bug #3934: expand(1/(1+%i)^4) => (-4)^(-1) (unsimplified) */
2767 expand((1+%i)^-2);
2768 -%i/2;
2770 expand(1/(1+%i)^4);
2771 -1/4;
2773 expand(4/(%i+1)^4)+1;
2776 expand(1/(1+%i)^8-1/16);
2779 expand((sqrt(-1/2)+sqrt(1/2))^-8);
2782 /* Bug #3956: expand(1/((sqrt(2)-1)*(sqrt(2)+1))) => 1/1 (unsimplified) */
2784 expand(1/((sqrt(2)-1)*(sqrt(2)+1)));
2787 expand(1/(-(sqrt(3)-1)*(sqrt(3)+1)));
2788 -1/2;
2790 expand(1/(-(sqrt(n+1)-1)*(sqrt(n+1)+1)));
2791 -1/n;
2793 /* Bugs #4045 and #4062:
2794    - Different results for integration in Maxima 5.45.1 and 5.46.0
2795    - limit(li[3](x), x, inf) gives li[3](inf) */
2797 integrate((log(x)^2)/(1+x), x, 0, 1);
2798 3*zeta(3)/2;
2800 limit(li[2](x), x, -inf);
2801 minf;
2803 limit(li[2](x), x, inf);
2804 infinity;
2806 limit(li[3](x), x, -inf);
2807 minf;
2809 limit(li[3](x), x, inf);
2810 infinity;
2812 limit(li[n](x), x, -inf);
2813 'limit(li[n](x), x, -inf);
2815 limit(li[2](sin(x)), x, inf);
2816 'limit(li[2](sin(x)), x, inf);
2818 /* Bug #4048: "An incorrect limit" */
2820 makelist(limit(diff(log(tan(%pi/2*tanh(x))), x), x, inf), simpfun, ['identity, lambda([e], trigreduce(trigsimp(trigreduce(e)))), 'exponentialize, 'trigreduce, 'trigrat]);
2821 [2, 2, 2, 2, 2];
2823 /* Bug 4064: "Simple limit triggers Lisp error '1 is not of type LIST'" */
2825 [limit(tan(%pi*tanh(x)), x, inf), limit(tan(%pi*tanh(x)), x, minf)];
2826 [0, 0];
2828 limit(tan(%pi/2*tanh(x)+1), x, inf);
2829 -cot(1);
2831 [limit(cot(%pi*tanh(x)), x, inf), limit(cot(%pi*tanh(x)), x, minf)];
2832 [minf, inf];
2834 [limit(tan(%pi/2*tanh(x)), x, inf), limit(tan(%pi/2*tanh(x)), x, minf)];
2835 [inf, minf];
2837 [limit(cot(%pi/2*tanh(x)), x, inf), limit(cot(%pi/2*tanh(x)), x, minf)];
2838 [0, 0];
2840 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
2841 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
2843 closeto(float(li[2](%i))-(0.915965594177219*%i-0.20561675835602827), 1e-15);
2844 true;
2846 closeto(bfloat(li[2](%i))-(9.15965594177219015054603b-1*%i - 2.05616758356028304559052b-1), 10*10^(-fpprec));
2847 true;
2849 /* Bug #4317: "Maxima doesn't know that signum(real) can only be -1, 0 or +1" */
2851 is(signum(x)<=1);
2852 true;
2854 is(signum(x)>=-1);
2855 true;
2857 is(signum(x)<2);
2858 true;
2860 is(signum(x)>-2);
2861 true;
2863 is(equal(signum(x),1));
2864 unknown;
2866 is(equal(signum(x),-1));
2867 unknown;
2869 is(equal(signum(x),0));
2870 unknown;
2872 is(equal(signum(x),1/2));
2873 false;
2875 is(equal(signum(x),-1/2));
2876 false;
2878 declare(x,complex);
2879 done;
2881 is(signum(x)<=1);
2882 unknown;
2884 is(signum(x)>=-1);
2885 unknown;
2887 kill(x);
2888 done;
2890 /* Bug #4326: "sin(asec(x)), cos(acsc(x)) and many more are incorrectly simplified" */
2893         err : 0.0,
2894         for outer in [sin, cos, tan, csc, sec, cot] do
2895         (
2896                 for inner in [asin, acos, atan, acsc, asec, acot] do
2897                 (
2898                         a : outer(inner(x)),
2899                         b : block([exponentialize : true, logarc : true, triginverses : false], outer(inner(x))),
2900                         for x in [-0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] do
2901                         (
2902                                 err : err + abs(at(a, 'x = x) - at(b, 'x = x))
2903                         )
2904                 )
2905         ),
2906         is(err < 1e-10)
2908 true$
2909 kill(err);
2910 done$