Merge branch 'master' of ssh://git.code.sf.net/p/maxima/code
[maxima.git] / tests / rtest_sqrt.mac
blob84411e60199406b272ae2716dab36bea1154a16f
1 /*******************************************************************************
2  *
3  *          Sqrt function
4  *
5  * The examples show what a Maxima user can expect from the sqrt function.
6  * Examples are taken from functions.wolfram.com.
7  *
8  ******************************************************************************/
10 /* ---------- Initialization ------------------------------------------------ */
12 kill(all);
13 done;
15 (reset(domain), reset(radexpand), done);
16 done;
18 /* ---------- Specific values ----------------------------------------------- */
20 sqrt([0, 0.0, 0.0b0]);
21 [0, 0.0, 0.0b0];
23 sqrt([-1, -1.0, -1.0b0]);
24 [%i, 1.0*%i, 1.0b0*%i];
26 /* Maxima does not simplify -%i -> -(-1)^(3/4) */
27 sqrt([%i, -%i]);
28 [(-1)^(1/4), sqrt(-%i)];
30 /* --- Polarform of specific values --- */
32 %emode:false;
33 false;
35 polarform(sqrt([-1, %i, -%i]));
36 [%e^(%i*%pi/2), %e^(%i*%pi/4), %e^(-%i*%pi/4)];
38 %emode:true;
39 true;
41 /* ---------- Values at infinities ------------------------------------------ */
43 /* Only limit can handle infinities. We check this. */
45 limit([sqrt(inf), sqrt(-inf), sqrt(minf), sqrt(-minf), sqrt(infinity)]);
46 [inf, infinity, infinity, inf, infinity];
48 limit(sqrt(x),x,inf);
49 inf;
50 limit(sqrt(x),x,-inf);
51 infinity;
52 limit(sqrt(x),x,minf);
53 infinity;
54 limit(sqrt(x),x,-minf);
55 inf;
57 /* ---------- Numerical evaluation ------------------------------------------ */
59 /* Check only for some values to be sure that it works */
61 (closeto(value,compare,tol):=
62   block(
63     [abse],
64     abse:abs(value-compare),if(abse<tol) then true else abse),
65     done);
66 done;
68 /* --- Float precision --- */
70 closeto(
71   sqrt(0.5),
72   0.707106781186547524400844362104849039284835937688474036588339868995366239,
73   1.0e-15);
74 true;
76 closeto(
77   sqrt(0.5+%i),
78   0.899453719973933636130613791812127819163772491812606935297702716148303720 +
79   0.555892970251421171992048047897569250691057814997037061577650731037508243*%i,
80   1.0e-15);
81 true;
83 closeto(
84   sqrt(0.5*%i),
85   0.500000000000000000000000000000000000000000000000000000000000000000000000 +
86   0.500000000000000000000000000000000000000000000000000000000000000000000000*%i,
87   1.0e-15);
88 true;
90 closeto(
91   sqrt(-0.5+%i),
92   0.555892970251421171992048047897569250691057814997037061577650731037508243 + 
93   0.899453719973933636130613791812127819163772491812606935297702716148303720*%i,
94   1.0e-15);
95 true;
97 closeto(
98   sqrt(-0.5),
99   0.707106781186547524400844362104849039284835937688474036588339868995366239*%i,
100   1.0e-15);
101 true;
103 closeto(
104   sqrt(-0.5-%i),
105   0.555892970251421171992048047897569250691057814997037061577650731037508243 - 
106   0.899453719973933636130613791812127819163772491812606935297702716148303720*%i,
107   1.0e-15);
108 true;
110 closeto(
111   sqrt(-0.5*%i),
112   0.500000000000000000000000000000000000000000000000000000000000000000000000 - 
113   0.500000000000000000000000000000000000000000000000000000000000000000000000*%i,
114   1.0e-15);
115 true;
117 closeto(
118   sqrt(0.5-%i),
119   0.899453719973933636130613791812127819163772491812606935297702716148303720 - 
120   0.555892970251421171992048047897569250691057814997037061577650731037508243*%i,
121   1.0e-15);
122 true;
124 /* --- Bigfloat precision with fpprec 72 --- */
126 (save_fpprec: fpprec, fpprec: 72);
129 closeto(
130   sqrt(0.5b0),
131   0.707106781186547524400844362104849039284835937688474036588339868995366239b0,
132   1.0b-72);
133 true;
135 closeto(
136   sqrt(0.5b0+%i),
137   0.899453719973933636130613791812127819163772491812606935297702716148303720b0 +
138   0.555892970251421171992048047897569250691057814997037061577650731037508243b0*%i,
139   1.0b-72);
140 true;
142 closeto(
143   sqrt(0.5b0*%i),
144   0.500000000000000000000000000000000000000000000000000000000000000000000000b0 +
145   0.500000000000000000000000000000000000000000000000000000000000000000000000b0*%i,
146   1.0b-72);
147 true;
149 closeto(
150   sqrt(-0.5b0+%i),
151   0.555892970251421171992048047897569250691057814997037061577650731037508243b0 + 
152   0.899453719973933636130613791812127819163772491812606935297702716148303720b0*%i,
153   1.0b-72);
154 true;
156 closeto(
157   sqrt(-0.5b0),
158   0.707106781186547524400844362104849039284835937688474036588339868995366239b0*%i,
159   1.0b-72);
160 true;
162 closeto(
163   sqrt(-0.5b0-%i),
164   0.555892970251421171992048047897569250691057814997037061577650731037508243b0 - 
165   0.899453719973933636130613791812127819163772491812606935297702716148303720b0*%i,
166   1.0b-72);
167 true;
169 closeto(
170   sqrt(-0.5b0*%i),
171   0.500000000000000000000000000000000000000000000000000000000000000000000000b0 - 
172   0.500000000000000000000000000000000000000000000000000000000000000000000000b0*%i,
173   1.0b-72);
174 true;
176 closeto(
177   sqrt(0.5b0-%i),
178   0.899453719973933636130613791812127819163772491812606935297702716148303720b0 - 
179   0.555892970251421171992048047897569250691057814997037061577650731037508243b0*%i,
180   1.0b-72);
181 true;
183 fpprec:save_fpprec;
184 ''save_fpprec;
186 /* ---------- Mirror symmetry ----------------------------------------------- */
188 /* We have no Mirror symmetry on the real negative axis. Maxima knows this. */
190 conjugate(sqrt(2+%i*y));
191 sqrt(2-%i*y);
193 conjugate(sqrt(-2+%i*y));
194 conjugate(sqrt(-2+%i*y));
196 conjugate(sqrt(x+%i*y));
197 conjugate(sqrt(x+%i*y));
199 assume(x>0);
200 [x>0];
202 conjugate(sqrt(x+%i*y));
203 sqrt(x-%i*y);
205 forget(x>0);
206 [x>0];
208 /* ---------- Series representations ---------------------------------------- */
210 /* Check taylor expansion */
212 taylor(sqrt(x), x, x0, 2);
213 sqrt(x0)+sqrt(x0)*(x-x0)/2/x0-sqrt(x0)*(x-x0)^2/8/x0^2;
215 /* Check taylor expansion as an argument to sqrt */
216 sqrt(taylor(x,x,x0,2));
217 sqrt(x0)+sqrt(x0)*(x-x0)/2/x0-sqrt(x0)*(x-x0)^2/8/x0^2;
219 taylor(sqrt(x), x, 1, 3);
220 1 + (x-1)/2 - (x-1)^2/8 + (x-1)^3/16;
222 taylor(sqrt(x+1), x, 0, 3);
223 1+x/2-x^2/8+x^3/16;
225 taylor(sqrt(x), x, inf, 2);
226 sqrt(x);
228 /* ---------- Differential equations ---------------------------------------- */
230 depends(y,x);
231 [y(x)];
232 %e_to_numlog:true; /* Simplify %e^log(x) -> x */
233 true;
235 ode2('diff(y,x)-y/(2*x), y, x);
236 y=%c*sqrt(x);
238 ode2(4*'diff(y,x,2)*x^2+y, y, x);
239 y=sqrt(x)*(%k1 + %k2*log(x)/2);
241 depends([g,h],x);
242 [g(x), h(x)];
244 ode2('diff(y,x)-'diff(g(x),x)/(2*g(x))*y, y,x);
245 y=%c*sqrt(g(x));
247 ode2('diff(y,x)-('diff(g(x),x)/(2*g(x))+'diff(h(x),x)/h(x))*y, y, x), expand;
248 y=%c*h(x)*sqrt(g(x));
250 remove([y,g,h], depends);
251 done;
252 %e_to_numlog:false;
253 false;
255 /* ---------- Transformations and argument simplifications ------------------ */
257 declare(z,complex);
258 done;
260 assume(xp>0, xn<0);
261 [xp>0, xn<0];
263 /* --- sqrt(-z) --- */
265 sqrt([-z,-x,-xp]);
266 [sqrt(-z), sqrt(-x), %i*sqrt(xp)];
268 sqrt([-1/2, -3/4, -1, -3/2, -2]);
269 [%i/sqrt(2), %i*sqrt(3)/2, %i, %i*sqrt(3)/sqrt(2), %i*sqrt(2)];
271 sqrt(-(xp+2));
272 sqrt(-(xp+2));
274 /* --- sqrt(%i*z) --- */
276 sqrt([%i*z,%i*x,%i*xp]);
277 [sqrt(%i*z), sqrt(%i*x), (-1)^(1/4)*sqrt(xp)];
279 sqrt([%i*1/2, %i*3/4, %i, %i*3/2, %i*2]);
280 [(-1)^(1/4)/sqrt(2), (-1)^(1/4)*sqrt(3)/2, (-1)^(1/4), 
281  (-1)^(1/4)*sqrt(3)/sqrt(2), (-1)^(1/4)*sqrt(2)];
283 sqrt([%i*(xp+2), %i*xp+%i*2]);
284 [(-1)^(1/4)*sqrt(xp+2),sqrt(%i*xp+%i*2)];
286 /* --- sqrt(-%i*z) --- */
288 sqrt([-%i*z,-%i*x,-%i*xp]);
289 [sqrt(-%i*z), sqrt(-%i*x), sqrt(-%i)*sqrt(xp)];
291 sqrt([-%i*1/2, -%i*3/4, -%i, -%i*3/2, -%i*2]);
292 [sqrt(-%i)/sqrt(2), sqrt(-%i)*sqrt(3)/2, sqrt(-%i), 
293  sqrt(-%i)*sqrt(3)/sqrt(2), sqrt(-%i)*sqrt(2)];
295 sqrt([-%i*(xp+2), -%i*xp-%i*2]);
296 [sqrt(-%i)*sqrt(xp+2),sqrt(-%i*xp-%i*2)];
298 /* --- sqrt(1/z) --- */
300 sqrt(1/xp) - 1/sqrt(xp);
302 sqrt(1/-xp) + 1/sqrt(-xp);
304 sqrt(1/-xn) - 1/sqrt(-xn);
307 /* --- sqrt(-1/z) --- */
309 sqrt([-1/x, -1/xp, -1/xn, -1/(xn-2), -1/(xp+2)]);
310 [%i/sqrt(x), %i/sqrt(xp), 1/sqrt(-xn), 1/sqrt(-xn+2), %i/sqrt(xp+2)];
312 /* --- sqrt(%i/z) --- */
314 sqrt([%i/x, %i/xp, %i/xn, %i/(xn-2), %i/(xp+2)]);
315 [sqrt(%i/x), (-1)^(1/4)/sqrt(xp), sqrt(-%i)/sqrt(-xn), sqrt(-%i)/sqrt(-xn+2), 
316  (-1)^(1/4)/sqrt(xp+2)];
318  /* --- sqrt(-%i/z) --- */
320 sqrt([-%i/x, -%i/xp, -%i/xn, -%i/(xn-2), -%i/(xp+2)]);
321 [sqrt(-%i/x), sqrt(-%i)/sqrt(xp), (-1)^(1/4)/sqrt(-xn), (-1)^(1/4)/sqrt(-xn+2), 
322  sqrt(-%i)/sqrt(xp+2)];
324 /* --- Half-angle formula --- */
326 sqrt(z/2);
327 sqrt(z)/sqrt(2);
329 /* --- Multiple arguments for products --- */
331 assume(a>0);
332 [a>0];
334 sqrt(a*z);
335 sqrt(a)*sqrt(z);
337 sqrt(b*z); /* Do not simplify if the sign is not known to be positive */
338 sqrt(b*z);
340 radcan(sqrt(z-z^2));
341 sqrt(z)*sqrt(1-z);
343 radcan(sqrt(-z^2-z));
344 sqrt(z)*sqrt(-z-1);
346 /* --- Multiple arguments for quotients --- */
348 radcan(sqrt(z)/sqrt(z+1));
349 sqrt(z)/sqrt(z+1);
351 radcan(sqrt(z/(z-1)));
352 sqrt(z)/sqrt(z-1);
354 /* radcan always simplifies sqrt(z1/z2) -> sqrt(z1)/sqrt(z2). 
355  * That is not correct.
356  */
358 declare([z1,z2],complex);
359 done;
361 radcan(sqrt(z1/z2));
362 sqrt(z1)/sqrt(z2);
364 /* More examples of this type */
366 radcan(sqrt(z/(1-z)));
367 sqrt(z)*sqrt(1/(1-z));
369 radcan(sqrt(-z/(1+z)));
370 sqrt(z)*sqrt(1/(-1+-z));
372 /* --- Power of arguments --- */
374 sqrt([z^2, x^2, xp^2, xn^2]);
375 [sqrt(z^2), abs(x), xp, -xn];
377 (domain:complex, done);
378 done;
380 res:sqrt(-z^2); /* Wrong: Maxima simplifies to %i*sqrt(z) */
381 sqrt(-z^2);
383 res,z=%i; /* We insert %i and get a wrong sign */
386 (reset(domain), done);
387 done;
389 /* --- Exponent of arguments --- */
391 /* Maxima simplifes immediately sqrt(exp(z)) -> exp(z/2) 
392  * That is wrong for imagpart(z) > %pi or < -%pi
393  * We give an example.
394  */
396 sqrt(exp(2*%pi*%i)); /* That is the expected result */
399 sqrt(exp(%i*y)); /* This expression simplifies wrongly */
400 sqrt(exp(%i*y));
402 %,y=2*%pi; /* We get a wrong sign */
405 /* ---------- Nested sqrt functions ----------------------------------------- */
407 /* The following 19 examples should simplify or evaluate to zero for all real
408  * or complex numbers. We check this. This way we can detect simplification
409  * errors for the sqrt function.
410  */
412 (domain:complex, done);
413 done;
415 /* --- Example 1 --- */
417 /* Wrong simplification because of sqrt(1/z) -> 1/sqrt(z) */
419 expr(z):=sqrt(-1-z)*sqrt(1/(1+z))*sqrt(1/z)*sqrt(-z)+1;
420 expr(z):=sqrt(-1-z)*sqrt(1/(1+z))*sqrt(1/z)*sqrt(-z)+1;
422 /* Check correct simplification of expression */
423 res:expr(z);
424 sqrt(-1-z)*sqrt(1/(1+z))*sqrt(1/z)*sqrt(-z)+1;
426 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
427 true;
428 closeto(float(rectform(res)),0.0,1e-15), z=-1/2; /* wrong -> 2 */
429 true;
430 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
431 true;
432 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
433 true;
434 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
435 true;
436 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
437 true;
438 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
439 true;
440 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
441 true;
443 /* --- Example 2 --- */
445 /* Wrong simplification because of sqrt(-1/z) -> 1/sqrt(-z) */
447 expr(z):=sqrt((z-1)/z)*sqrt(z/(1-z))-sqrt(-1/z)*sqrt(z);
448 expr(z):=sqrt((z-1)/z)*sqrt(z/(1-z))-sqrt(-1/z)*sqrt(z);
450 /* Check correct simplification of equation */
451 res:expr(z);
452 sqrt((z-1)/z)*sqrt(z/(1-z))-sqrt(-1/z)*sqrt(z);
454 closeto(float(rectform(res)),0.0,1e-15), z=1/2; /* Wrong -> 2 */
455 true;
456 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
457 true;
458 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
459 true;
460 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
461 true;
462 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
463 true;
464 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
465 true;
466 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
467 true;
468 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
469 true;
471 /* --- Example 3 --- */
473 /* Wrong simplification because of sqrt(z^2)/sqrt(-z^2) -> -%i */
475 expr(z):=sqrt(sqrt(z^2+1)-1)/sqrt(1-sqrt(z^2+1))-sqrt(z^2)/sqrt(-z^2);
476 expr(z):=sqrt(sqrt(z^2+1)-1)/sqrt(1-sqrt(z^2+1))-sqrt(z^2)/sqrt(-z^2);
478 /* Check correct simplification of equation */
479 res:expr(z);
480 sqrt(sqrt(z^2+1)-1)/sqrt(1-sqrt(z^2+1))-sqrt(z^2)/sqrt(-z^2);
482 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
483 true;
484 closeto(float(rectform(res)),0.0,1e-15), z=-1/2; 
485 true;
486 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;      /* Wrong -> 2 */
487 true;
488 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;     /* Wrong -> 2 */
489 true;
490 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;  /* Wrong -> 2 */
491 true;
492 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;  
493 true; 
494 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i; /* Wrong -> 2 */
495 true;
496 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
497 true;
499 /* --- Example 4 --- */
501 /* Wrong simplification sqrt(-z^2)/sqrt(z^2) -> %i */
503 expr(z):=sqrt(sqrt(1-z^2)-1)/sqrt(1-sqrt(1-z^2))-sqrt(-z^2)/sqrt(z^2);
504 expr(z):=sqrt(sqrt(1-z^2)-1)/sqrt(1-sqrt(1-z^2))-sqrt(-z^2)/sqrt(z^2);
506 /* Check correct simplification of equation */
507 res:expr(z);
508 sqrt(sqrt(1-z^2)-1)/sqrt(1-sqrt(1-z^2))-sqrt(-z^2)/sqrt(z^2);
510 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
511 true;
512 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
513 true;
514 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;     /* Wrong -> 2 */
515 true;
516 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;    /* Wrong -> 2 */
517 true;
518 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i; /* Wrong -> 2 */
519 true;
520 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
521 true;
522 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i; /* Wrong -> 2 */
523 true;
524 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
525 true;
527 /* --- Example 5 --- */
529 /* Wrong simplification sqrt(-1/z) -> 1/sqrt(-z) */
531 expr(z):=sqrt(sqrt(z^2+1)-z)/sqrt(z-sqrt(z^2+1))+sqrt(z)*sqrt(-1/z);
532 expr(z):=sqrt(sqrt(z^2+1)-z)/sqrt(z-sqrt(z^2+1))+sqrt(z)*sqrt(-1/z);
534 /* Check correct simplification of equation */
535 res:expr(z);
536 sqrt(sqrt(z^2+1)-z)/sqrt(z-sqrt(z^2+1))+sqrt(z)*sqrt(-1/z);
538 closeto(float(rectform(res)),0.0,1e-15), z=1/2; /* Wrong -> 2 */
539 true;
540 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
541 true;
542 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
543 true;
544 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
545 true;
546 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
547 true;
548 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
549 true;
550 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
551 true;
552 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
553 true;
555 /* --- Example 6 --- */
557 /* Wrong simplification sqrt(1/z) -> 1/sqrt(z) */
559 expr(z):=sqrt(z+sqrt(z^2+1))/sqrt(-z-sqrt(z^2+1))+sqrt(-z)*sqrt(1/z);
560 expr(z):=sqrt(z+sqrt(z^2+1))/sqrt(-z-sqrt(z^2+1))+sqrt(-z)*sqrt(1/z);
562 /* Check correct simplification of equation */
563 res:expr(z);
564 sqrt(z+sqrt(z^2+1))/sqrt(-z-sqrt(z^2+1))+sqrt(-z)*sqrt(1/z);
566 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
567 true;
568 closeto(float(rectform(res)),0.0,1e-15), z=-1/2; /* Wrong -> 2 */
569 true;
570 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
571 true;
572 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
573 true;
574 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
575 true;
576 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
577 true;
578 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
579 true;
580 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
581 true;
583 /* --- Example 7 --- */
585 /* This expression does not simplify at all. All numeric values are correct.
586  * rootscontract gives the result 2*%i. That is wrong. Because we have no
587  * problem with the simplification of the expression, the error is in 
588  * rootscontract for this example */
590 expr(z):=sqrt(sqrt(1-z^2)-%i*z)/sqrt(%i*z-sqrt(1-z^2))+sqrt(%i*z)*sqrt(%i/z);
591 expr(z):=sqrt(sqrt(1-z^2)-%i*z)/sqrt(%i*z-sqrt(1-z^2))+sqrt(%i*z)*sqrt(%i/z);
593 /* Check correct simplification of equation */
594 res:expr(z);
595 sqrt(sqrt(1-z^2)-%i*z)/sqrt(%i*z-sqrt(1-z^2))+sqrt(%i*z)*sqrt(%i/z);
597 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
598 true;
599 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
600 true;
601 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
602 true;
603 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
604 true;
605 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
606 true;
607 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
608 true;
609 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
610 true;
611 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
612 true;
614 /* --- Example 8 --- */
616 /* This expression does not simplify at all. All numeric values are correct.
617  * rootscontract gives the result 2*%i. That is wrong. Because we have no
618  * problem with the simplification of the expression, the error is in 
619  * rootscontract for this example. See also example 7.
620  */
622 expr(z):=sqrt(%i*z+sqrt(1-z^2))/sqrt(-%i*z-sqrt(1-z^2))+sqrt(-%i*z)*sqrt(-%i/z);
623 expr(z):=sqrt(%i*z+sqrt(1-z^2))/sqrt(-%i*z-sqrt(1-z^2))+sqrt(-%i*z)*sqrt(-%i/z);
625 /* Check correct simplification of equation */
626 res:expr(z);
627 sqrt(%i*z+sqrt(1-z^2))/sqrt(-%i*z-sqrt(1-z^2))+sqrt(-%i*z)*sqrt(-%i/z);
629 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
630 true;
631 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
632 true;
633 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
634 true;
635 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
636 true;
637 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
638 true;
639 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
640 true;
641 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
642 true;
643 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
644 true;
646 /* --- Example 9 --- */
648 /* This examples simplifies to zero immediately. This is by accident.
649  * We know that the applied simplifications are not correct.
650  */
652 expr(z):=sqrt(z^2/(1-sqrt(1-z^2)))*sqrt((1-sqrt(1-z^2))/z^2)-1;
653 expr(z):=sqrt(z^2/(1-sqrt(1-z^2)))*sqrt((1-sqrt(1-z^2))/z^2)-1;
655 /* Check correct simplification of equation */
656 res:expr(z);
657 sqrt(z^2/(1-sqrt(1-z^2)))*sqrt((1-sqrt(1-z^2))/z^2)-1;
659 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
660 true;
661 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
662 true;
663 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
664 true;
665 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
666 true;
667 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
668 true;
669 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
670 true;
671 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
672 true;
673 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
674 true;
676 /* --- Example 10 --- */
678 /* Wrong simplification sqrt(1/z) -> 1/sqrt(z) */
680 expr(z):=sqrt(z/(sqrt(z^2+1)-z))*sqrt((sqrt(z^2+1)-z)/z)-sqrt(1/(z^2+1))*sqrt(z^2+1)-sqrt(1/z)*sqrt(z)+1;
681 expr(z):=sqrt(z/(sqrt(z^2+1)-z))*sqrt((sqrt(z^2+1)-z)/z)-sqrt(1/(z^2+1))*sqrt(z^2+1)-sqrt(1/z)*sqrt(z)+1;
683 /* Check correct simplification of equation */
684 res:expr(z);
685 sqrt(z/(sqrt(z^2+1)-z))*sqrt((sqrt(z^2+1)-z)/z)-sqrt(1/(z^2+1))*sqrt(z^2+1)-sqrt(1/z)*sqrt(z)+1;
687 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
688 true;
689 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
690 true;
691 closeto(float(rectform(res)),0.0,1e-15), z=%i/2; /* Wrong -> 2 */
692 true;
693 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2; /* Wrong -> 2 */
694 true;
695 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
696 true;
697 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
698 true;
699 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
700 true;
701 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
702 true;
704 /* --- Example 11 --- */
706 /* This example is correct:
707  * No simplification of this expression. All numeric results are correct.
708  * radcan does not simplify too. That is not wrong.
709  * rootscontract simplifies correctly to zero.
710  */
712 expr(z):=sqrt(sqrt(z)-1)*sqrt(sqrt(z)+1)-sqrt(z-1);
713 expr(z):=sqrt(sqrt(z)-1)*sqrt(sqrt(z)+1)-sqrt(z-1);
715 /* Check correct simplification of equation */
716 res:expr(z);
717 sqrt(sqrt(z)-1)*sqrt(sqrt(z)+1)-sqrt(z-1);
719 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
720 true;
721 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
722 true;
723 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
724 true;
725 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
726 true;
727 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
728 true;
729 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
730 true;
731 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
732 true;
733 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
734 true;
736 /* --- Example 12 --- */
738 /* This example is correct:
739  * No simplification of this expression. All numeric results are correct.
740  * radcan does not simplify too. That is not wrong.
741  * rootscontract simplifies correctly to zero. See also example 11.
742  */
744 expr(z):=sqrt(sqrt(z+1)-1)*sqrt((sqrt(z+1)+1))-sqrt(z);
745 expr(z):=sqrt(sqrt(z+1)-1)*sqrt((sqrt(z+1)+1))-sqrt(z);
747 /* Check correct simplification of equation */
748 res:expr(z);
749 sqrt(sqrt(z+1)-1)*sqrt((sqrt(z+1)+1))-sqrt(z);
751 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
752 true;
753 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
754 true;
755 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
756 true;
757 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
758 true;
759 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
760 true;
761 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
762 true;
763 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
764 true;
765 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
766 true;
768 /* --- Example 13 --- */
770 /* This example does not simplify. All numeric results are correct.
771  * radcan gives a modified result, which is correct too.
772  * rootscontract gives a wrong expression.
773  */
775 expr(z):=sqrt(2)*sqrt(1-sqrt(1-z))*sqrt(sqrt(1-z)+%i*sqrt(z))/(%i*(1-sqrt(1-z))+sqrt(z))-1;
776 expr(z):=sqrt(2)*sqrt(1-sqrt(1-z))*sqrt(sqrt(1-z)+%i*sqrt(z))/(%i*(1-sqrt(1-z))+sqrt(z))-1;
778 /* Check correct simplification of equation */
779 res:expr(z);
780 sqrt(2)*sqrt(1-sqrt(1-z))*sqrt(sqrt(1-z)+%i*sqrt(z))/(%i*(1-sqrt(1-z))+sqrt(z))-1;
782 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
783 true;
784 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
785 true;
786 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
787 true;
788 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
789 true;
790 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
791 true;
792 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
793 true;
794 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
795 true;
796 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
797 true;
799 /* --- Example 14 --- */
801 /* This example does not simplify. All numeric results are correct.
802  * radcan gives the correct result zero.
803  * rootscontract gives a wrong expression.
804  */ 
806 expr(z):=sqrt(1-z)*(1+z-%i*sqrt(1-z^2))/(sqrt(z+1)*(1-z+%i*sqrt(1-z^2)))+%i;
807 expr(z):=sqrt(1-z)*(1+z-%i*sqrt(1-z^2))/(sqrt(z+1)*(1-z+%i*sqrt(1-z^2)))+%i;
809 /* Check correct simplification of equation */
810 res:expr(z);
811 sqrt(1-z)*(1+z-%i*sqrt(1-z^2))/(sqrt(z+1)*(1-z+%i*sqrt(1-z^2)))+%i;
813 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
814 true;
815 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
816 true;
817 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
818 true;
819 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
820 true;
821 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
822 true;
823 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
824 true;
825 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
826 true;
827 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
828 true;
830 /* --- Example 15 --- */
832 /* This example does not simplify. All numeric results are correct.
833  * radcan gives the correct result zero.
834  * rootscontract gives a wrong expression.
835  */ 
837 expr(z):=sqrt(1-z)*(1+z+%i*sqrt(1-z^2))/(sqrt(z+1)*(1-z-%i*sqrt(1-z^2)))-%i;
838 expr(z):=sqrt(1-z)*(1+z+%i*sqrt(1-z^2))/(sqrt(z+1)*(1-z-%i*sqrt(1-z^2)))-%i;
840 /* Check correct simplification of equation */
841 res:expr(z);
842 sqrt(1-z)*(1+z+%i*sqrt(1-z^2))/(sqrt(z+1)*(1-z-%i*sqrt(1-z^2)))-%i;
844 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
845 true;
846 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
847 true;
848 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
849 true;
850 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
851 true;
852 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
853 true;
854 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
855 true;
856 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
857 true;
858 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
859 true;
861 /* --- Example 16 --- */
863 /* This example does not simplify. All numeric results are correct.
864  * radcan gives a changed expression, which is correct too.
865  * rootscontract gives only a small change.
866  */ 
868 expr(z):=1/sqrt(1-sqrt(z))+1/sqrt(1+sqrt(z))-sqrt(2)*sqrt(sqrt(1-z)+1)/sqrt(1-z);
869 expr(z):=1/sqrt(1-sqrt(z))+1/sqrt(1+sqrt(z))-sqrt(2)*sqrt(sqrt(1-z)+1)/sqrt(1-z);
871 /* Check correct simplification of equation */
872 res:expr(z);
873 1/sqrt(1-sqrt(z))+1/sqrt(1+sqrt(z))-sqrt(2)*sqrt(sqrt(1-z)+1)/sqrt(1-z);
875 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
876 true;
877 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
878 true;
879 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
880 true;
881 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
882 true;
883 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
884 true;
885 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
886 true;
887 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
888 true;
889 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
890 true;
892 /* --- Example 17 --- */
894 /* This example does not simplify. All numeric results are correct.
895  * radcan gives a changed expression, which is correct too.
896  * rootscontract gives only a small change.
897  */ 
899 expr(z):=1/sqrt(1-sqrt(z))-1/sqrt(1+sqrt(z))-sqrt(2)*sqrt(1-sqrt(1-z))/sqrt(1-z);
900 expr(z):=1/sqrt(1-sqrt(z))-1/sqrt(1+sqrt(z))-sqrt(2)*sqrt(1-sqrt(1-z))/sqrt(1-z);
902 /* Check correct simplification of equation */
903 res:expr(z);
904 1/sqrt(1-sqrt(z))-1/sqrt(1+sqrt(z))-sqrt(2)*sqrt(1-sqrt(1-z))/sqrt(1-z);
906 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
907 true;
908 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
909 true;
910 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
911 true;
912 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
913 true;
914 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
915 true;
916 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
917 true;
918 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
919 true;
920 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
921 true;
923 /* --- Example 18 --- */
925 /* All numeric values are correct. radcan and rootscontract does not
926  * change the expression.
927  */
929 expr(z):=sqrt(1+sqrt(z))+sqrt(1-sqrt(z))-sqrt(2)*sqrt(1+sqrt(1-z));
930 expr(z):=sqrt(1+sqrt(z))+sqrt(1-sqrt(z))-sqrt(2)*sqrt(1+sqrt(1-z));
932 /* Check correct simplification of equation */
933 res:expr(z);
934 sqrt(1+sqrt(z))+sqrt(1-sqrt(z))-sqrt(2)*sqrt(1+sqrt(1-z));
936 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
937 true;
938 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
939 true;
940 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
941 true;
942 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
943 true;
944 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
945 true;
946 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
947 true;
948 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
949 true;
950 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
951 true;
953 /* --- Example 19 --- */
955 /* All numeric values are correct. radcan and rootscontract does not
956  * change the expression.
957  */
959 expr(z):=sqrt(1+sqrt(z))-sqrt(1-sqrt(z))-sqrt(2)*sqrt(1-sqrt(1-z));
960 expr(z):=sqrt(1+sqrt(z))-sqrt(1-sqrt(z))-sqrt(2)*sqrt(1-sqrt(1-z));
962 /* Check correct simplification of equation */
963 res:expr(z);
964 sqrt(1+sqrt(z))-sqrt(1-sqrt(z))-sqrt(2)*sqrt(1-sqrt(1-z));
966 closeto(float(rectform(res)),0.0,1e-15), z=1/2;
967 true;
968 closeto(float(rectform(res)),0.0,1e-15), z=-1/2;
969 true;
970 closeto(float(rectform(res)),0.0,1e-15), z=%i/2;
971 true;
972 closeto(float(rectform(res)),0.0,1e-15), z=-%i/2;
973 true;
974 closeto(float(rectform(res)),0.0,1e-15), z=1/2+%i;
975 true;
976 closeto(float(rectform(res)),0.0,1e-15), z=-1/2+%i;
977 true;
978 closeto(float(rectform(res)),0.0,1e-15), z=-1/2-%i;
979 true;
980 closeto(float(rectform(res)),0.0,1e-15), z=1/2-%i;
981 true;
983 (reset(domain), done);
984 done;
986 /* ---------- Complex characteristics --------------------------------------- */
988 /* --- Real part --- */
990 realpart(sqrt(x+%i*y));
991 (x^2+y^2)^(1/4)*cos(1/2*atan2(y,x));
993 /* --- Imaginary part --- */
995 imagpart(sqrt(x+%i*y));
996 (x^2+y^2)^(1/4)*sin(1/2*atan2(y,x));
998 /* --- Absolute value --- */
1000 abs(sqrt(x+%i*y));
1001 abs(sqrt(%i*y+x))$
1003 /* --- Argument --- */
1005 carg(sqrt(x+%i*y));
1006 1/2*atan2(y,x);
1008 /* --- Conjugate value --- */
1010 conjugate(sqrt(xp+%i*y));
1011 sqrt(xp-%i*y);
1013 rectform(conjugate(sqrt(xp+%i*y)));
1014 (xp^2+y^2)^(1/4)*cos(1/2*atan(y/xp))-%i*(xp^2+y^2)^(1/4)*sin(1/2*atan(y/xp));
1016 /* --- Signum value --- */
1018 /* Maxima has no signum function for complex expressions */
1020 signum(sqrt(z));
1021 'signum(sqrt(z));
1023 signum(sqrt(xp));
1026 signum(sqrt(xn));
1027 'signum(sqrt(xn));
1029 /* ---------- Differentiation ----------------------------------------------- */
1031 diff(sqrt(z),z);
1032 1/(2*sqrt(z));
1034 diff(sqrt(z),z,2);
1035 -1/(4*z^(3/2));
1037 /* ---------- Integration --------------------------------------------------- */
1039 /* --- Indefinite integration --- */
1041 integrate(sqrt(z),z);
1042 2/3*z^(3/2);
1044 /* --- Definite integration --- */
1046 integrate(sqrt(t),t,0,1);
1047 2/3;
1049 integrate(sqrt(t)*log(t),t,0,1);
1050 -4/9;
1052 /* Maxima does not know the general result: -> 2/3*2F1(3/2,-a;5/2;-1) */
1053 integrate(sqrt(t)*(t+1)^a,t,0,1);
1054 'integrate(sqrt(t)*(t+1)^a,t,0,1);
1056 /* Maxima can evaluate the integral for a an integer or a half integral value */
1057 integrate(sqrt(t)*(t+1)^2,t,0,1);
1058 184/105;
1060 integrate(sqrt(t)*(t+1)^-2,t,0,1);
1061 %pi/4-1/2$
1063 integrate(sqrt(t)*(t+1)^(1/2),t,0,1);
1064 (-log(sqrt(2)+1)+log(1-sqrt(2))+3*2^(3/2))/8-log(-1)/8$
1066 integrate(sqrt(t)*(t+1)^(-1/2),t,0,1);
1067 (-log(sqrt(2)+1)/2)+log(sqrt(2)-1)/2+sqrt(2);
1069 integrate(sqrt(t)*(t+1)^(3/2),t,0,1);
1070 (-3*log(sqrt(2)+1)+3*log(1-sqrt(2))+25*2^(3/2))/48-log(-1)/16$
1072 integrate(sqrt(t)*(t+1)^(-3/2),t,0,1);
1073 log(sqrt(2)+1)-log(sqrt(2)-1)-sqrt(2);
1075 /* ---------- Integral transforms ------------------------------------------- */
1077 /* --- Laplace transform --- */
1079 laplace(sqrt(t),t,s);
1080 sqrt(%pi)/(2*s^(3/2));
1082 /* Maxima does not know the inverse laplace transform 
1083  *    -> -1/(2*t^(3/2)*sqrt(%pi)
1084  */
1085 ilt(sqrt(s),s,t);
1086 'ilt(sqrt(s),s,t);
1088 /* ---------- Limits -------------------------------------------------------- */
1090 assume(b>1);
1091 [b>1];
1093 limit(sqrt(z)/b^z,z,inf);
1096 forget(b>1);
1097 [b>1];
1099 limit(sqrt(z)*log(z),z,0);
1102 limit(log(z)/sqrt(z),z,inf);
1105 /* ---------- Representation through hypergeometric functions --------------- */
1107 hgfred([-1/2],[],1-z); 
1108 sqrt(z);
1110 hgfred([],[], 1/2*log(z)), %e_to_numlog:true;
1111 sqrt(z);
1113 hgfred([-1/2,b],[b],1-z);
1114 sqrt(z);
1116 (reset(domain), done);
1117 done;
1119 /* bug reported to mailing list 2018-11-01: "radcan(sqrt(-1/(sqrt(-3)+1))*sqrt(2))" */
1121 radcan (sqrt(-1/(sqrt(-3)+1))*sqrt(2));
1122 sqrt(2)/sqrt(-sqrt(3)*%i - 1);
1124 radcan (sqrt(-1/(sqrt(-3)+1)));
1125 1/sqrt(-sqrt(3)*%i - 1);
1127 radcan (sqrt(1/(sqrt(-3)+1))*sqrt(2));
1128 sqrt(2)/sqrt((sqrt(3)*%i + 1));