| blob | 89c7f452e3673e527a5e046d28648070579605f4 |
1 /* vim: set ft=maxima softtabstop=4 shiftwidth=4 expandtab: */
3 /*
4 Most of the code in this Maxima package is a rather direct port of Sympy's
5 denesting code:
6 http://docs.sympy.org/1.0/_modules/sympy/simplify/sqrtdenest.html
7 http://docs.sympy.org/1.0/_modules/sympy/simplify/radsimp.html
8 The code is licensed under the GPL-compatible "3-clause BSD license".
9 */
11 /**************************************************************************
12 Copyright (c) 2006-2017 SymPy Development Team (original code)
13 Copyright (c) 2017 Gilles Schintgen (Maxima port and extensions)
15 All rights reserved.
17 Redistribution and use in source and binary forms, with or without
18 modification, are permitted provided that the following conditions are met:
20 a. Redistributions of source code must retain the above copyright notice,
21 this list of conditions and the following disclaimer.
22 b. Redistributions in binary form must reproduce the above copyright
23 notice, this list of conditions and the following disclaimer in the
24 documentation and/or other materials provided with the distribution.
25 c. Neither the name of SymPy nor the names of its contributors
26 may be used to endorse or promote products derived from this software
27 without specific prior written permission.
30 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
31 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
32 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
33 ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
34 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
35 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
36 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
37 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
38 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
39 OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
40 DAMAGE.
41 **************************************************************************/
43 /*
44 RADDENEST
46 This package denests different classes of expressions with nested
47 roots. Most of the algorithms are specific to square roots, but some
48 handle higher roots.
50 The implemented algorithms are documented in the following articles:
52 [1] Allan Borodin, Ronald Fagin, John E. Hopcroft, and Martin Tompa:
53 "Decreasing the Nesting Depth of Expressions Involving Square Roots"
54 J. Symbolic Computation (1985) 1, 169-188
55 http://researcher.watson.ibm.com/researcher/files/us-fagin/symb85.pdf
57 [2] David J. Jeffrey and Albert D. Rich:
58 "Symplifying Square Roots of Square Roots by Denesting"
59 in "Computer Algebra Systems: A Practical Guide",
60 M.J. Wester, Ed., Wiley 1999
61 http://www.cybertester.com/data/denest.pdf
63 [3] Thomas J. Osler:
64 "Cardan Polynomials and the Reduction of Radicals"
65 Mathematics Magazine 74(1), Feb. 2001
66 http://www.rowan.edu/open/depts/math/osler/mathmag026-032.pdf
68 [4] Mascha Honsbeek:
69 "Radical Extensions and Galois Groups", Chapter 3
70 (PhD Thesis)
71 http://www.math.kun.nl/~bosma/students/honsbeek/M_Honsbeek_thesis.pdf
73 As such, the raddenest function should be able to denest some or all
74 of the denestable expressions of the following types:
76 * square roots of sums of (multiple) unnested square roots
77 (_sqrtdenest_rec)
79 * expressions of the form sqrt(a + b*sqrt(r)) that denest using square
80 roots can be denested numerically or symbolically
81 (_sqrt_symbolic_denest)
83 * expressions of the form sqrt(a + b*sqrt(r)) where a, b, and r are
84 linear combinations of square roots of positive rationals
85 (_sqrt_biquadratic_denest)
87 * expressions of the form sqrt(a + b*sqrt(r)) that can be denested
88 using fourth roots
90 * n-th roots of a+b*sqrt(r)
91 (_rad_denest_cardano)
93 * square roots of a sum of two cube roots
94 (_rad_denest_ramanujan)
97 INTERPRETATION OF RADICAL EXPRESSIONS
99 In general a^(1/n) has n different values.
101 However Maxima always simplifies such numeric expressions for positive
102 real 'a' and for negative real 'a' in case of a square root:
103 sqrt(-4)=2*%i, 8^(1/3)=2
104 This is independent of the global 'domain' variable.
106 With domain=complex Maxima no longer simplifies a^(1/n) for negative
107 real 'a'. It will only factor out (-1)^(1/n) which is then treated
108 according to the global variable m1pbranch. It also no longer simplifies
109 sqrt(x^2).
111 Raddenest should give results that are consistent with the above.
113 Example:
114 (2-sqrt(5))^(1/3) = 1/2-sqrt(5)/2
115 This denesting is only performed when domain=real and the cube root
116 is interpreted as a real-valued function over the reals.
119 IMPLEMENTATION NOTES
121 Functions whose names begin with an underscore are internal helper
122 functions. As such the only user-level function is the main command
123 "raddenest".
125 raddenest and _raddenest0 are relatively simple wrapper functions that
126 apply the various algorithm functions to the given expression and its
127 subexpressions.
129 POSSIBLE IMPROVEMENTS AND FURTHER READING
131 [5] Richard Zippel:
132 "Simplification of Expressions Involving Radicals"
133 J. Symbolic Computation (1985) 1, 189-210
134 http://www.sciencedirect.com/science/article/pii/S0747717185800146
136 [6] Susan Landau:
137 "Simplification of Nested Radicals"
138 1993
139 https://www.researchgate.net/publication/2629046_Simplification_of_Nested_Radicals
141 [7] Johannes Blömer:
142 "Simplifying Expressions Involving Radicals"
143 PhD Thesis, 1993
144 http://www.cs.nyu.edu/exact/pap/rootBounds/sumOfSqrts/bloemerThesis.pdf
145 */
147 /*
148 Helper functions
149 */
151 _sqrtp(expr):=block(
152 /* Returns true iff expr is a literal square root,
153 i.e. radicand^(1/2).
154 */
155 [inflag:true],
156 not atom(expr) and op(expr)="^" and args(expr)[2]=1/2
157 )$
158 _sqrtpowerp(expr):=block(
159 /* Returns true iff expr can be interpreted as a square root.
160 Maxima automatically simplifies sqrt(125) to 5^(3/2).
161 Negative exponents are rejected.
162 */
163 [inflag:true],
164 not atom(expr) and op(expr)="^" and denom(args(expr)[2])=2
165 and num(args(expr)[2])>0
166 )$
168 _algebraicp(p):=block(
169 /* Returns true if p is comprised only of rationals or
170 square(!) roots of rationals and algebraic operations. */
171 [inflag:true],
172 if ratnump(p) then return(true)
173 elseif atom(p) then return(false)
174 elseif op(p)="^" then
175 ((args(p)[2]=-1 or denom(args(p)[2])=2)
176 and _algebraicp(args(p)[1]))
177 else
178 (operatorp(p,["+","*"])
179 and apply("and",map('_algebraicp,args(p))))
180 )$
182 _sqrt_depth(p):=block(
183 /* Returns the maximum nesting depth of any square root argument of p,
184 where p is an algebraic expression as defined by _algebraicp.
186 Neither of these square roots contains any other square roots
187 so the depth is 1:
188 (%i1) _sqrt_depth(1 + sqrt(2)*(1 + sqrt(3)));
189 (%o1) 1
190 The sqrt(3) is contained within a square root so the depth is 2:
191 (%i2) _sqrt_depth(1 + sqrt(2)*sqrt(1 + sqrt(3)));
192 (%o2) 2
193 */
194 [inflag:true],
195 if atom(p) or ratnump(p) then return(0)
196 elseif (op(p)="+" or op(p)="*")
197 then apply('max,maplist('_sqrt_depth,p))
198 elseif _sqrtpowerp(p)
199 then _sqrt_depth(args(p)[1])+1
200 else 0
201 )$
203 _complexity(p):=block(
204 /* Returns a complexity measure of a given expression p.
205 It's similar to _sqrt_depth but increases with any nesting
206 of radicals (not only square roots) and penalize sums
207 of radicals by adding up individual complexities.
208 See [2] for more details.
209 It is used as a secondary recursion termination criterion
210 (after _sqrt_depth) in some algorithms.
211 */
212 [inflag:true,e],
213 if atom(p) or ratnump(p) then return(0)
214 elseif (op(p)="+")
215 then map('_complexity,p)
216 elseif (op(p)="*")
217 then apply('max,maplist('_complexity,p))
218 elseif op(p)="^"
219 and ratnump(e:args(p)[2])
220 and denom(e)#1
221 then _complexity(args(p)[1])+1
222 else 0
223 )$
225 _sqrtcontract0(prod):= block(
226 /* Contracts _all_ square roots in a given term/quotient.
228 Maxima's built-in commands sometimes have difficulties
229 handling quotients with square roots in the denominator.
230 Forcing them to contract allows for easier simplification.
232 (%i1) sqrt(2)/(7*sqrt(27));
233 (%o1) sqrt(2)/(7*3^(3/2))
234 (%i2) _sqrtcontract(%);
235 (%o2) sqrt(6)/63
236 */
237 [inflag:true,rootsconmode:false,r,e,rad:1],
238 if atom(prod) or op(prod)#"*" then return(prod),
239 for f in args(prod) do (
240 if not atom(f) and op(f)="^" then (
241 [r,e]: args(f),
242 if integerp(r) and ratnump(e) and e<0 and denom(e)=2
243 then rad: rad*r
244 )
245 ),
246 rootscontract(num(prod)*sqrt(rad))/rootscontract(denom(prod)*sqrt(rad))
247 )$
248 _sqrtcontract(expr):= scanmap('_sqrtcontract0,expr)$
251 /* _subsets(n)
252 Returns all possible subsets of the set (0, 1, ..., n-1) except the
253 empty set, listed in reversed lexicographical order according to binary
254 representation, so that the case of the fourth root is treated last.
255 Example:
256 (%i1) _subsets(3);
257 (%o2) [[1,0,0],[0,1,0],[1,1,0],[0,0,1],[1,0,1],[0,1,1],[1,1,1]]
258 */
259 _ss(n):= if n=1 then [[0],[1]] else
260 block([r:_ss(n-1)],
261 apply('append, map(lambda([z], [cons(0,z),cons(1,z)]),r)));
262 _subsets(n):=rest(_ss(n));
265 _splitcoef(p):= block(
266 /* Splits a given term p in two factors:
267 a rational coefficient and the remaining factors.
268 The return value is a list.
269 Note: 7^(3/2) is handled as 7*sqrt(7),
270 contrary to maxima's default.
271 Also: 15/sqrt(2) -> 15/2*sqrt(2)
272 */
273 [inflag:true,c,base,expt,r,spl],
274 /* single rational number */
275 if ratnump(p) then return([p,1]),
276 p: _sqrtcontract(p),
277 /* single power that can be read as a sqrt
278 e.g. 5^(13/2) -> [5^13, sqrt(5)]
279 2^(-3/2) -> [1/4, sqrt(2)] */
280 if not atom(p) and op(p)="^" and ratnump(base:args(p)[1])
281 and ratnump(expt:args(p)[2]) and denom(expt)=2
282 then (
283 if expt>0 then
284 return([base^quotient(num(expt),2), sqrt(base)])
285 else
286 return([base^quotient(num(expt)-2,2), sqrt(base)])
287 )
288 /* product with ratnum as first factor.
289 check if second factor is a power -> recurse.
290 this handles e.g. 3*5^(3/2) */
291 elseif not atom(p) and op(p)="*" and ratnump(c:first(args(p)))
292 then (
293 r: rest(p),
294 if not atom(r) and op(r)="^" then (
295 spl:_splitcoef(r),
296 return([c*spl[1], spl[2]])
297 )
298 else return([c,r])
299 )
300 else return([1,p])
301 )$
303 _sqrt_match(p):=block(
304 /* Returns [a, b, r] for _sqrt_match(a + b*sqrt(r)) where, in addition to
305 matching, sqrt(r) also has then maximal _sqrt_depth among addends of p.
306 If p can't be matched, the return value is false.
308 Example:
309 (%o1) _sqrt_match(2*sqrt(sqrt(5)+1)+sqrt(2)*sqrt(3)+sqrt(2)+1);
310 (%o2) [sqrt(2)*sqrt(3)+sqrt(2)+1,2,sqrt(5)+1]
311 */
312 [inflag: true, p:expand(p), res,pargs,v,nmax,depth,
313 ar,i,r,b,bv,rv,a1,b1,x2,x2args,a],
314 if numberp(p) then res:[p,0,0]
315 elseif not atom(p) and op(p)="+" then (
316 pargs:args(p),
317 /* handle "trivial" case without nested radicals: */
318 if apply("and",maplist(lambda([x],ratnump(x)and x>0),pargs^2))
319 then return(reverse(_split_surds(p))),
320 /* to make the process canonical, the argument is included in the tuple
321 so when the max is selected, it will be the largest arg having a
322 given depth */
323 v: makelist([_sqrt_depth(pargs[i]),pargs[i],i], i, 1, length(pargs)),
324 nmax:v[1],
325 for da in v do (
326 if da[1]>nmax[1]
327 or (da[1]=nmax[1] and ordergreatp(da[2],nmax[2]))
328 then nmax:da
329 ),
330 if nmax[1]=0 then res:false
331 else (
332 /* select r (term with max. depth)*/
333 [depth, ar, i]: nmax,
334 r: pargs[i],
335 pargs: delete(r,pargs),
336 v: delete(nmax,v),
337 b: 1,
338 /* if product, keep only factor with max. depth for r,
339 other factors will be part of coefficient b */
340 if not atom(r) and op(r)="*" then (
341 bv:[],
342 rv:[],
343 for x in args(r) do (
344 if _sqrt_depth(x)<depth then bv:endcons(x,bv)
345 else rv:endcons(x,rv)
346 ),
347 b: apply("*",bv),
348 r: apply("*",rv)
349 ),
350 /* collect terms containing r */
351 a1:[],
352 b1:[b],
353 for x in v do (
354 if x[1]<depth then a1:endcons(x[2],a1)
355 else (
356 x2: x[2],
357 if x2=r then b1:endcons(1,b1)
358 else (
359 if not atom(x2) and op(x2)="*" then (
360 x2args: args(x2),
361 if member(r,x2args) then (
362 x2args: delete(r,x2args,1),
363 x2args: apply("*",x2args),
364 push(x2args,b1)
365 )
366 else a1:endcons(x[2],a1)
367 )
368 else a1:endcons(x[2],a1)
369 )
370 )
371 ), /* for */
372 a: apply("+",a1),
373 b: apply("+",b1),
374 res: [a, b, r^2]
375 )
376 )
377 else (
378 [b, r]: _split_coef(p),
379 if _sqrtp(r) then res:[0,b,r^2]
380 else res:false
381 ),
382 return(res)
383 )$
385 _split_gcd(a):= block(
386 /* Splits the list of integers a into two lists of integers, a1 and a2,
387 such that a1 have some common divisor and the integers in a2
388 are not divisible by g=gcd(a1).
389 The result is returned as [g, a1, a2].
391 Example:
393 (%i1) _split_gcd(55, 35, 22, 14, 77, 10);
394 (%o2) [5, [55, 35, 10], [22, 14, 77]]
395 */
396 [g:first(a), g1, b1:[], b2:[]],
397 push(g,b1),
398 for x in rest(a) do (
399 g1: gcd(g,x),
400 if g1=1 then push(x,b2)
401 else (g:g1, push(x,b1))
402 ),
403 [g,reverse(b1),reverse(b2)]
404 )$
406 _split_surds(expr):= block(
407 /* Splits an expression with terms whose squares are rationals
408 into a sum of terms whose surds squared have gcd equal to g
409 and a sum of terms with surds squared prime with g.
411 Example:
413 (%i1) _split_surds(3*sqrt(3)+sqrt(5)/7+sqrt(6)+sqrt(10)+sqrt(15))
414 (%o1) [3, sqrt(2) + sqrt(5) + 3, sqrt(5)/7 + sqrt(10)]
415 */
416 [inflag:false,ar,coeff_muls,surds:[],g,b1,b2,g2,g1,b1n:[],c,s,s1,a1v:[],a2v:[],a,b],
417 ar: args(rootscontract(expr)),
418 coeff_muls: maplist('_splitcoef,ar),
419 for cs in coeff_muls do
420 if _sqrtp(cs[2]) then push(cs[2]^2,surds),
421 [g, b1, b2] : _split_gcd(sort(surds)),
422 g2: g,
423 if b2=[] and length(b1)>=2 then (
424 for x in b1 do if x#g then b1n:endcons(x/g, b1n),
425 /* only a common factor has been factored; split again */
426 [g1, b1n, b2] : _split_gcd(b1n),
427 g2: g*g1
428 ),
429 for cs in coeff_muls do (
430 [c, s]: cs,
431 if _sqrtp(s) then (
432 s1: first(args(s)),
433 if member(s1,b1)
434 then a1v:endcons(c*sqrt(s1/g2),a1v)
435 else a2v:endcons(c*s,a2v)
436 )
437 else a2v:endcons(c*s,a2v)
438 ),
439 a: apply("+",a1v),
440 b: apply("+",a2v),
441 [g2,a,b]
442 )$
444 /*
445 Denesting functions
446 */
448 raddenest(expr, [max_iter]):= block(
449 /* Denests sqrts in an expression that contain other square roots
450 if possible, otherwise returns the expr unchanged.
451 */
452 [algebraic:true, rootsconmode:false, z,res],
453 if max_iter=[] then max_iter: 3
454 else max_iter: max_iter[1],
455 res: for i:1 thru max_iter do (
456 z: _sqrtcontract(_raddenest0(expr)),
457 if expr=z then
458 return(expr),
459 expr: z
460 ),
461 if res#done then return(res)
462 else return(expr)
463 )$
465 _raddenest0(expr):= block(
466 /* Returns expr after (recursively) denesting its arguments. */
467 [inflag:true,algebraic:true,ar,ret,radicand,ex],
468 if _sqrtpowerp(expr) then (
469 /* expr could be of the form ( )^(n/2) */
470 radicand: expand(ratsimp(args(expr)[1])),
471 ex: 2*args(expr)[2],
472 if not atom(radicand) and op(radicand)="+" then (
473 ar: args(radicand),
474 if length(ar)>2 and apply("and",maplist('integerp,ar^2)) then (
475 ret: catch(_sqrtdenest_rec(sqrt(radicand))),
476 if ret#"raddenestStopIteration" then
477 return(ret^ex)
478 ),
479 if length(ar)=2 and ratnump(ar[1]^3) and ratnump(ar[2]^3) then (
480 ret: _rad_denest_ramanujan(sqrt(radicand)),
481 if ret#false then
482 return(ret^ex)
483 ),
484 radicand: maplist('_raddenest0,ar),
485 radicand: expand(apply("+",radicand))
486 ),
487 return(_sqrtdenest1(sqrt(radicand),true)^ex)
488 )
489 elseif not atom(expr) and op(expr)="^"
490 /* try Cardano polynomials for (a+b*sqrt(r))^(m/n) */
491 and ratnump(ex:args(expr)[2])
492 and not atom(radicand:args(expr)[1])
493 and op(radicand)="+"
494 and length(ar:args(radicand))=2
495 and apply("and",maplist(lambda([x],ratnump(x)and x>0),ar^2))
496 then return(_rad_denest_cardano(expr))
497 elseif not atom(expr) and op(expr)="^" and evenp(denom(ex:args(expr)[2]))
498 /* Some fourth roots denest in two steps
499 e.g. (19601-3465*2^(5/2))^(7/4)
500 -> (99-35*2^(3/2))^(7/2)
501 -> (5*sqrt(2)-7)^7
502 This is only useful if the global nesting depth
503 indeed decreases, but this is not checked by the
504 current code.
505 (Same behaviour as sqdnst.mac) */
506 then (
507 radicand: expand(ratsimp(args(expr)[1])),
508 ex: 2*args(expr)[2],
509 return(_raddenest0(sqrt(radicand))^ex)
510 )
511 elseif not atom(expr) then (
512 inflag:false,
513 if op(expr)="/"
514 then return(_raddenest0(args(expr)[1])/_raddenest0(args(expr)[2]))
515 else return(map('_raddenest0,expr))
516 ),
517 return(expr)
518 )$
520 _sqrtdenest_rec(expr):= block(
521 /* Helper that denests the square root of three or more surds.
523 It returns the denested expression; if it cannot be denested it
524 throws raddenestStopIteration
526 Algorithm: expr.base is in the extension Q_m = Q(sqrt(r_1),..,sqrt(r_k));
527 split expr.base = a + b*sqrt(r_k), where `a` and `b` are on
528 Q_(m-1) = Q(sqrt(r_1),..,sqrt(r_(k-1))); then a^2 - b^2*r_k is
529 on Q_(m-1); denest sqrt(a^2 - b^2*r_k) and so on.
531 See [1], section 6.
533 Examples
534 ========
536 (%i1) _sqrtdenest_rec(sqrt(-72*sqrt(2) + 158*sqrt(5) + 498));
537 (%o1) (-sqrt(10))+9*sqrt(5)+sqrt(2)+9
538 (%i2) w:-6*sqrt(55)-6*sqrt(35)-2*sqrt(22)-2*sqrt(14)
539 +2*sqrt(77)+6*sqrt(10)+65$
540 (%i3) _sqrtdenest_rec(sqrt(w));
541 (%o3) (-sqrt(11))-sqrt(7)+3*sqrt(5)+sqrt(2)
542 */
543 [algebraic:true,inflag:true, g,a,b,c2,a1,b1,c2_1,c_1,d_1,c,ac,d,r,f],
544 if atom(expr) or op(expr)#"^" then
545 return(raddenest(expr)),
546 if (r:args(expr)[1]) < 0 then
547 return(%i*_sqrtdenest_rec(sqrt(-r))),
548 [g,a,b]: _split_surds(r),
549 a: a*sqrt(g),
550 if a<b then [a,b]:[b,a],
551 c2: rootscontract(expand(a^2-b^2)),
552 if not atom(c2) and length(args(c2)) > 2 then (
553 [g,a1,b1]: _split_surds(c2),
554 a1: a1*sqrt(g),
555 if a1<b1 then [a1,b1]:[b1,a1],
556 c2_1: rootscontract(expand(a1^2-b1^2)),
557 c_1: _sqrtdenest_rec(sqrt(c2_1)),
558 d_1: _sqrtdenest_rec(sqrt(a1+c_1)),
559 f: _sqrtcontract(expand(ratsimp(b1/d_1))),
560 c: expand(d_1*sqrt(2)/2+f*sqrt(2)/2)
561 )
562 else c: _sqrtdenest1(sqrt(c2),true),
563 if _sqrt_depth(c)>1 then
564 throw("raddenestStopIteration"),
565 ac: rootscontract(a+c),
566 /* check for end of recursion */
567 if not atom(ac) and length(args(ac)) >= length(args(expr))
568 and _complexity(ac)>=_complexity(args(expr)[1])
569 then throw("raddenestStopIteration"),
570 d: raddenest(sqrt(ac)),
571 if _sqrt_depth(d) > 1
572 then throw("raddenestStopIteration"),
573 r: rootscontract(expand(ratsimp(b/d))),
574 r: rootscontract(ratsimp(expand(d*sqrt(2)+r*sqrt(2))/2)),
575 return(expand(r))
576 )$
578 _sqrtdenest1(expr, denester):= block(
579 /* Return denested expr after denesting with simpler
580 methods or, that failing, using the denester.
581 'denester' is a boolean indicating if _denester
582 should be called.
583 */
584 [inflag:true, algebraic: true, a,val,b,r,d2,z,dr2,dr,res,av0],
585 if not _sqrtp(expr) then return(expr),
586 a: args(expr)[1], /*radicand*/
587 if atom(a) then return(expr),
588 val: _sqrt_match(a),
589 if val=false then return(expr),
590 [a,b,r]: val,
591 /* try a quick numeric denesting */
592 d2: rootscontract(expand(a^2-b^2*r)),
593 if ratnump(d2) then(
594 if d2>0 then (
595 z: _sqrt_numeric_denest(a,b,r,d2),
596 if z#false then return(z)
597 )
598 else (
599 /* fourth root case, e.g.
600 sqrtdenest(sqrt(3+2*sqrt(3)))
601 -> sqrt(2)*3^(1/4)/2+sqrt(2)*3^(3/4)/2 */
602 dr2: expand(-d2*r),
603 dr: sqrt(dr2),
604 if ratnump(dr) then (
605 z: _sqrt_numeric_denest(expand(b*r),a,r,dr2),
606 if z#false then return(expand(ratsimp(z/r^(1/4))))
607 )
608 )
609 )
610 else (
611 z: _sqrt_symbolic_denest(a,b,r),
612 if z#false then return(z)
613 ),
614 if not denester or not _algebraicp(expr) then
615 return(expr),
616 res: _sqrt_biquadratic_denest(expr,a,b,r,d2),
617 if res#false then return(res),
618 /* now call the denester */
619 av0: [a,b,r,d2],
620 z: _denester([expr^2], av0, 0, _sqrt_depth(expr))[1],
621 /* denester has side-effect on av0! */
622 if av0[2]=false then return(expr),
623 if z#false then (
624 if _sqrt_depth(z)=_sqrt_depth(expr)
625 and _complexity(z)>=_complexity(expr)
626 then return(expr)
627 else return(z)
628 ),
629 return(expr)
630 )$
632 _sqrt_symbolic_denest(a,b,r):= block(
633 /* Given an expression, sqrt(a + b*sqrt(r)), return the denested
634 expression or false.
636 Algorithm 1:
637 Check if the standard denesting formula with d2=a^2-b^2*r
638 can the applied, i.e. check if d2=a^2-b^2*r is a square.
639 For a+b*sqrt(r) = (c+d*sqrt(r))^2, a=c^2+b^2*r and
640 a^2-b^2*r = d^4*r^2-2*c*d^2*r+c^4, quadratic in r.
641 -> Compute the discriminant.
642 a and r must be positive; hence denesting depends on the
643 facts database.
645 Algorithm 2 (Sympy):
646 If r = ra + rb*sqrt(rr), try replacing sqrt(rr) in 'a' with
647 (y^2 - ra)/rb, and if the result is a quadratic, ca*y^2 + cb*y + cc, and
648 (cb + b)^2 - 4*ca*cc is 0, then sqrt(a + b*sqrt(r)) can be rewritten as
649 sqrt(ca*(sqrt(r) + (cb + b)/(2*ca))^2).
651 Examples (alg. 1):
652 ==================
653 (%i2) raddenest(sqrt(9*y+6*x^2*sqrt(y)+x^4));
654 (%o2) sqrt(9*y + 6*x^2 sqrt(y) + x^4)
655 (%i3) assume(y>0)$
656 (%i4) raddenest(sqrt(9*y+6*x^2*sqrt(y)+x^4));
657 (%o4) 3* sqrt(y) + x^2
659 Examples (alg. 2):
660 ==================
662 (%i1) [a,b,r]: [16-2*sqrt(29), 2, -10*sqrt(29)+55]$
663 (%i2) _sqrt_symbolic_denest(a, b, r);
664 (%o2) sqrt(-2*sqrt(29) + 11) + sqrt(5)
666 If the expression is numeric, it will be simplified:
668 (%i3) w: sqrt(sqrt(sqrt(3) + 1) + 1) + 1 + sqrt(2)$
669 (%i4) [a,b,r]: _sqrt_match(expand(w^2))$
670 (%i5) _sqrt_symbolic_denest(a,b,r);
671 (%o5) sqrt(sqrt(sqrt(3)+1)+1)+sqrt(2)+1
673 Otherwise, it will be simplified depending on the value of
674 the global "domain" variable:
676 (%i6) assume(x >= 0)$
677 (%i7) w: sqrt(expand((sqrt(1+sqrt(x))-1)^2));
678 (%o6) sqrt(sqrt(x)-2*sqrt(1+sqrt(x))+2)
679 (%i7) [a,b,r]: _sqrt_match(w^2)$
680 (%i8) _sqrt_symbolic_denest(a,b,r), domain:'real;
681 (%o8) abs(sqrt(1+sqrt(x))-1)
682 (%i9) _sqrt_symbolic_denest(a,b,r), domain:'complex;
683 (%o9) sqrt((sqrt(1+sqrt(x))-1)^2)
684 */
685 [algebraic:true, rval,ra,rb,rr,y,newd2,newa,ca,cb,cc,z,cont,d2,d,s],
686 supcontext(cont),
687 y: gensym(),
688 assume(y>0),
689 /* make sure a and r are positive */
690 if not (is(r>=0)=true and is(ratsubst(y,r,a)>=0)=true)
691 then go(tryalg2),
692 d2: ratsubst(y,r,a^2-b^2*r),
693 if hipow(d2,y)=2 then (
694 [ca,cb,cc]: [ratcoef(d2,y,2),ratcoef(d2,y,1),ratcoef(d2,y,0)],
695 if is(equal(cb^2-4*ca*cc,0)) then (
696 /* Discriminant is 0, polynomial d2 can be factored
697 as a perfect square. Extract a sqrt. */
698 d: ratsimp(sqrt(ca)*(r+cb/(2*ca))),
699 /* If b is of constant sign (>=0 or <=0), a denesting
700 should not be rejected simply because b could be zero
701 for isolated values of the involved variables. */
702 assume(not equal(b,0)),
703 s: signum(b),
704 if not (integerp(s) and s#0 and ratnump(sqrt(ca)))
705 then go(tryalg2),
706 z: sqrt(a/2+d/2)+signum(b)*sqrt(a/2-d/2),
707 z: ratsimp(z),
708 killcontext(cont),
709 return(sqrt(z^2))
710 )
711 ),
712 tryalg2,
713 rval: _sqrt_match(r),
714 if rval=false then (
715 killcontext(cont),
716 return(false)
717 ),
718 [ra,rb,rr]: rval,
719 if rb#0 then (
720 newa: subst((y^2-ra)/rb,sqrt(rr),a),
721 if hipow(newa,y)=2 then (
722 [ca,cb,cc]: [ratcoef(newa,y,2),ratcoef(newa,y,1),ratcoef(newa,y,0)],
723 cb: cb+b,
724 if is(equal(cb^2-4*ca*cc,0))=true then (
725 z: sqrt(ca*(sqrt(r)+cb/(2*ca))^2),
726 if constantp(z) then z: rootscontract(expand(z)),
727 killcontext(cont),
728 return(z)
729 )
730 )
731 ),
732 killcontext(cont),
733 return(false)
734 )$
736 _sqrt_numeric_denest(a,b,r,d2):=block(
737 /* Helper that denest expr = a + b*sqrt(r),
738 with d2 = a^2 - b^2*r > 0
739 or returns false if not denested.
740 */
741 [algebraic:true, depthr,d,vad,vad1],
742 depthr: _sqrt_depth(r),
743 d: sqrt(d2),
744 vad: a+d,
745 /* _sqrt_depth(res) <= _sqrt_depth(vad) + 1
746 _sqrt_depth(expr) = depthr + 2
747 There is denesting if _sqrt_depth(vad)+1 < depthr + 2;
748 If vad^2 is a rational there is a fourth root */
749 if _sqrt_depth(vad) < depthr + 1 or ratnump(vad^2) then (
750 vad1: ratsimp(1/vad),
751 return(expand(sqrt(vad/2)+signum(b)*sqrt(b^2*r*vad1/2)))
752 ) else return(false)
753 )$
755 _sqrt_biquadratic_denest(expr,a,b,r,d2):= block(
756 /* denest expr = sqrt(a + b*sqrt(r))
757 where a, b, r are linear combinations of square roots of
758 positive rationals on the rationals (SQRR) and r > 0, b != 0,
759 d2 = a^2 - b^2*r > 0
761 If it cannot denest it returns false.
763 ALGORITHM
764 Search for a solution A of type SQRR of the biquadratic equation
765 4*A^4 - 4*a*A^2 + b^2*r = 0 (1)
766 sqd = sqrt(a^2 - b^2*r)
767 Choosing the sqrt to be positive, the possible solutions are
768 A = sqrt(a/2 +/- sqd/2)
769 Since a, b, r are SQRR, then a^2 - b^2*r is a SQRR,
770 so if sqd can be denested, it is done by
771 _sqrtdenest_rec, and the result is a SQRR.
772 Similarly for A.
774 Examples of solutions (in both cases a and sqd are positive):
776 Example of expr with solution sqrt(a/2 + sqd/2) but not
777 solution sqrt(a/2 - sqd/2):
778 expr: sqrt(-sqrt(15) - sqrt(2)*sqrt(-sqrt(5) + 5) - sqrt(3) + 8);
779 a: -sqrt(15) - sqrt(3) + 8;
780 sqd = -2*sqrt(5) - 2 + 4*sqrt(3)
782 Example of expr with solution sqrt(a/2 - sqd/2) but not
783 solution sqrt(a/2 + sqd/2):
784 w: 2 + sqrt(2) + sqrt(3) + (1 + sqrt(3))*sqrt(2 + sqrt(2) + 5*sqrt(3))
785 expr: sqrt(expand(w^2))
786 a: 4*sqrt(6) + 8*sqrt(2) + 47 + 28*sqrt(3)
787 sqd: 29 + 20*sqrt(3)
789 Define B = b/(2*A); eq.(1) implies a = A^2 + B^2*r; then
790 expr^2 = a + b*sqrt(r) = (A + B*sqrt(r))^2
792 Examples
793 ========
794 (%i1) z: sqrt((2*sqrt(2) + 4)*sqrt(2 + sqrt(2)) + 5*sqrt(2) + 8)$
795 [a,b,r] : _sqrt_match(z^2)$
796 d2: a^2 - b^2*r$
797 _sqrt_biquadratic_denest(z, a, b, r, d2);
798 (%o1) sqrt(2) + sqrt(sqrt(2) + 2) + 2
799 */
800 [inflag:true, algebraic:true, res,y2,sqd,x1,x2,A,B,z],
801 if r<=0 or d2<0 or b=0 or _sqrt_depth(args(expr)[1])<2 then
802 return(false),
803 /* check that all terms in the linear combinations a, b, r
804 have integer square */
805 res: for x in [a,b,r] do (
806 if atom(x) or op(x)#"+" then x: [x],
807 res: for y in args(x) do (
808 y2: y^2,
809 if not integerp(y2) or y2<=0 then
810 return(false)
811 ),
812 if res#done then return(false)
813 ),
814 if res#done then return(false),
815 /*sqd: expand(raddenest(sqrt(ratsimp(d2)))), */
816 sqd: raddenest(sqrt(d2)),
817 if _sqrt_depth(sqd)>1 then
818 return(false),
819 [x1,x2]: [a/2+sqd/2, a/2-sqd/2],
820 /* look for a solution A with depth 1 */
821 res: for x in [x1,x2] do (
822 A: raddenest(sqrt(x)),
823 if not _sqrt_depth(A)>1 then (
824 B: b/(2*A),
825 B: ratsimp(rootscontract(expand(ratsimp(B)))),
826 z: A+B*sqrt(r),
827 if z<0 then z: -z,
828 return(expand(z))
829 )
830 ),
831 if res#done then return(res)
832 else return(false)
833 )$
835 _denester(nested,av0,h,max_depth_level):= block(
836 /* Denests a list of expressions that contain nested square roots.
838 nested: list of radicands (no top level sqrt!)
839 av0: [a,b,r,d2] (a+b*sqrt(r) and d2=a^2-b^2*r)
840 h: recursion counter, starts at 0
841 max_depth_level: nesting depth
843 Algorithm based on [1].
845 It is assumed that all of the elements of 'nested' share the same
846 bottom-level radicand. (This is stated in the paper, on page 177, in
847 the paragraph immediately preceding the algorithm.)
849 When evaluating all of the arguments in parallel, the bottom-level
850 radicand only needs to be denested once. This means that calling
851 _denester with x arguments results in a recursive invocation with x+1
852 arguments; hence _denester has polynomial complexity.
854 However, if the arguments were evaluated separately, each call would
855 result in two recursive invocations, and the algorithm would have
856 exponential complexity.
858 This is discussed in the paper in the middle paragraph of page 179.
860 The return value is [denested, subset] where subset corresponds
861 to the subset of terms that led to the denesting.
862 Otherwise [false,false] is returned.
863 */
864 [algebraic:true,p,sqp,R,vals,nested2,d,f,df,v,vad,sqvad,sqvad1,s2,FR,s,res],
865 if h>max_depth_level
866 then return([false, false]),
867 if av0[2]=false
868 then return([false, false]),
869 if (av0[1]=false and apply("and", maplist('numberp,nested))) then (
870 /* no arguments are nested */
871 res: for f in _subsets(length(nested)) do (
872 /* test subset 'f' of nested */
873 p:[],
874 for i:1 thru length(nested) do (
875 if f[i]=1 then push(nested[i],p)
876 ),
877 p: expand(apply("*",p)),
878 if apply("+",f)>1 and last(f)=1 then
879 p: -p,
880 sqp: sqrt(p),
881 if ratnump(sqp) then (
882 /* Got a perfect square so return its sqrt */
883 return([sqp,f])
884 )
885 ),
886 if res#done then return(res)
887 else /* ... otherwise, return the radicand
888 from the previous invocation. */
889 return([sqrt(last(nested)), makelist(0,length(nested))])
890 )
891 else (
892 R: false,
893 if av0[1]#false then (
894 vals: [[av0[1],av0[2]]],
895 R: av0[3],
896 nested2: [av0[4],R],
897 av0[1]: false
898 )
899 else (
900 vals: maplist('_sqrt_match,nested),
901 vals: sublist(vals,lambda([x],x#false)),
902 res: for v in vals do (
903 /* if b=0, r is not defined */
904 if v[3]#0 then (
905 if R#false then (
906 if R#v[3] then (
907 av0[2]: false,
908 return([false,false])
909 )
910 )
911 else R: v[3]
912 )
913 ),
914 if res#done then return(res),
915 if R=false then
916 /* return the radicand of the previous invocation */
917 return([sqrt(last(nested)), makelist(0,length(nested))]),
918 nested2: makelist(expand(v[1]^2)-expand(R*v[2]^2),v,vals),
919 nested2: endcons(R,nested2)
920 ),
921 [d,f]: _denester(nested2,av0,h+1,max_depth_level),
922 if f=false then return([false,false]),
923 if not member(1,makelist(f[i],i,1,length(nested))) then (
924 v: last(vals),
925 return([sqrt(v[1]+expand(v[2]*d)), f])
926 )
927 else (
928 p: [],
929 for i:1 thru length(nested) do (
930 if f[i]=1 then push(nested[i],p)
931 ),
932 p: apply("*",p),
933 v: _sqrt_match(p),
934 if member(1,f)
935 and sublist_indices(f,lambda([x],x=1))[1] < length(nested)
936 and f[length(nested)]=1
937 then (
938 v[1]: -v[1],
939 v[2]: -v[2]
940 ),
941 if f[length(nested)+1]=0 then (
942 /* solution denests with square roots */
943 vad: expand(v[1]+d),
944 if vad<=0 then
945 /* return the radicand from the previous invocation */
946 return([sqrt(last(nested)), makelist(0,length(nested))]),
947 if not (_sqrt_depth(vad) <= _sqrt_depth(R)+1
948 or nump((vad)^2)) then (
949 av0[2]: false,
950 return([false,false])
951 ),
952 sqvad: _sqrtdenest1(sqrt(vad), false),
953 if not (_sqrt_depth(sqvad) <= _sqrt_depth(R)+1) then (
954 av0[2]: false,
955 return([false,false])
956 ),
957 sqvad1: ratsimp(1/sqvad),
958 res: _sqrtcontract(expand(sqvad*sqrt(2)/2 + v[2]*sqrt(R)*sqvad1*sqrt(2)/2)),
959 return([res,f])
960 )
961 else (
962 /* Solution requires a fourth root */
963 s2: expand(v[2]*R) + d,
964 if s2 <= 0 then
965 return([sqrt(last(nested)), makelist(0,length(nested))]),
966 FR: expand(R)^(1/4),
967 s: sqrt(s2),
968 return([expand(s/(sqrt(2)*FR)+v[1]*FR/(sqrt(2)*s)), f])
969 )
970 )
971 )
972 )$
974 _cardano_poly(n,c,x):= block(
975 /* "Cardano" polynomials as defined in [3] */
976 [poly,poly1,poly2],
977 poly2: rat('x), /* degree n-2 */
978 poly1: rat('x^2-2*c), /* degree n-1 */
979 if n=1 then return(poly2)
980 elseif n=2 then return(poly1),
981 for i:3 thru n do (
982 poly: 'x*poly1-c*poly2,
983 poly2: poly1,
984 poly1: poly
985 ),
986 poly
987 )$
989 _rad_denest_cardano(expr):= block(
990 /* Denests expr=(a+b*sqrt(r))^(m/n) using
991 Cardano polynomials. See [3] for details. */
992 [a,b,r,m,n,divs,d,c,p,s,em,res,sig,x,
993 ratprint:false,bftorat:false,fpprec:50,ratepsilon:1e-40],
994 [a,b,r]: _sqrt_match(args(expr)[1]),
995 m: num(args(expr)[2]),
996 n: denom(args(expr)[2]),
997 /* keep sign; reinterpret as a+sqrt(b) as in ref. article */
998 sig: signum(b),
999 b: b^2*r,
1000 divs: reverse(rest(listify(divisors(n)))),
1001 /* check each divisor of the radical's index */
1002 res:while divs#[] do (
1003 d: pop(divs),
1004 c:(a^2-b)^(1/d),
1005 if ratnump(c) then (
1006 /* expr can be denested if p has a rational root
1007 equal to s. In order to avoid excessive numeric
1008 evaluations of the polynomial a single rational
1009 candidate is checked. This candidate is a rational
1010 approximation of s. The given values for fpprec
1011 and ratepsilon are heuristic in nature. */
1012 p: _cardano_poly(d,c,x)-2*a,
1013 p: horner(p*denom(p)),
1014 s: rat(bfloat((a+sqrt(b))^(1/d)+(a-sqrt(b))^(1/d))),
1015 if subst(s,x,p)=0 then
1016 return(ratsimp((s+sig*sqrt(s^2-4*c))/2)^(d/n))
1017 )
1018 ),
1019 if res#done then return(expand(res^m))
1020 else return(expr)
1021 )$
1023 _rad_denest_ramanujan(expr):= block(
1024 /* Denest expr=sqrt(alpha^(1/3)+beta^(1/3)) using the results
1025 of [4], Chapter 3. If expr can't be denested, return false.
1026 alpha and beta must be rationals; parameter checking must
1027 be done by the calling function. */
1028 [inflag:true,alpha,beta,c,Fba,t,l:[],s,tt,res],
1029 [alpha,beta]: (args(args(expr)[1]))^3,
1030 /* If domain=complex, Maxima interprets a^(1/3) with a>0
1031 as the real branch of the cube root, but leaves the
1032 expression untouched if a<0. */
1033 if domain='complex and (alpha<0 or beta<0) then return(false),
1034 /* The perfect cube case can be trivially denested
1035 but this code segment should, in fact, not be reached
1036 since the simplifier already handles this case. */
1037 if ratnump(c:(beta/alpha)^(1/3)) then
1038 return(sqrt(alpha^(1/3)*sqrt(1+c))),
1039 Fba: t^4+4*t^3+8*beta/alpha*t-4*beta/alpha,
1040 Fba: multthru(denom(4*beta/alpha),Fba),
1041 /* apply the rational root theorem to the polynomial Fba */
1042 for n in divisors(coeff(Fba,t,0)) do
1043 for m in divisors(coeff(Fba,t,4)) do (
1044 push(n/m,l),
1045 push(-n/m,l)
1046 ),
1047 l: setify(l),
1048 s: for i in l do (
1049 if subst(i,t,Fba)=0 then return(i)
1050 ),
1051 if s=done then return(false),
1052 tt: beta-s^3*alpha,
1053 /* Sign: either alpha^(1/3)+beta^(1/3) is positive, then
1054 so is its sqrt, or alpha^(1/3)+beta^(1/3) is negative,
1055 then its sqrt is %i times a positive real.
1056 Both cases are treated correctly by the formula below. */
1057 res: sqrt(1/tt)*abs(-1/2*s^2*(alpha^2)^(1/3)
1058 +s*(alpha*beta)^(1/3)+(beta^2)^(1/3)),
1059 return(ratsimp(res))
1060 )$