src/sqrtdenest.lisp

   1 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   2 ;;;
   3 ;;; sqrtdenest - denest an expression containing square roots of square roots
   4 ;;;
   5 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   6 ;;; This library is free software; you can redistribute it and/or modify it
   7 ;;; under the terms of the GNU General Public License as published by the
   8 ;;; Free Software Foundation; either version 2 of the License, or (at
   9 ;;; your option) any later version.
  10 ;;;
  11 ;;; This library is distributed in the hope that it will be useful, but
  12 ;;; WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ;;; Library General Public License for more details.
  15 ;;;
  16 ;;; You should have received a copy of the GNU General Public License along
  17 ;;; with this library; if not, write to the Free Software
  18 ;;; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  19 ;;;
  20 ;;; Copyright (C) 2016 David Billinghurst
  21 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  22
  23 (in-package :maxima)
  24
  25 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  26 ;;;
  27 ;;; SQRTDENEST simplifies square roots of square roots.
  28 ;;;
  29 ;;; This implementation uses the same algorithm as the maxima language
  30 ;;; SQRTDENEST in share package sqdnst.mac.  It only handles simple cases.
  31 ;;;
  32 ;;; Reference:
  33 ;;;
  34 ;;;  Jeffrey, David J. and Rich, Albert D. (1999)
  35 ;;;  Simplifying Square Roots of Square Roots by Denesting,
  36 ;;;  in Computer Algebra Systems (Ed. M. J. Wester)
  37 ;;;  <http://www.cybertester.com/data/denest.pdf>
  38 ;;;
  39 ;;; Further reading:
  40 ;;;
  41 ;;; A. Borodin and R. Fagin and J. Hopcroft and M. Tompa. (1985)
  42 ;;; Decreasing the Nesting Depth of expressions Involving Square Roots.
  43 ;;; J. Symbolic Computation, 1:169-188
  44 ;;; <http://www.almaden.ibm.com/cs/people/fagin/symb85.pdf>
  45 ;;;
  46 ;;; Landau, Susan (1992) A Note on Zippel Denesting.
  47 ;;; Journal of Symbolic Computation 13:41-45
  48 ;;;
  49 ;;; Landau, Susan (1992) Simplification of Nested Radicals.
  50 ;;; SIAM Journal on Computing, 21: 85-110.
  51 ;;;
  52 ;;; Landau, Susan (1994). How to Tangle with a Nested Radical.
  53 ;;; Math. Intell. 16:49-55
  54 ;;;
  55 ;;; Zippel, Richard (1985) Simplification of Expressions involving Radicals.
  56 ;;; Journal of Symbolic Computation, 1:189-210.
  57 ;;;
  58 ;;; Python sympy sqrtdenest function
  59 ;;; http://docs.sympy.org/latest/_modules/sympy/simplify/sqrtdenest.html
  60 ;;; which references Borodin et al (1985), Jeffrey and Rich (2009)
  61 ;;;
  62 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  63 (defmfun $sqrtdenest (e)
  64   "Denest square roots in expression e"
  65   (sqrtdenest1 e))
  66
  67 (defun sqrtdenest1 (e)
  68   "Denest square roots in expression e"
  69   (cond
  70    ((mapatom e) e)
  71    ((eq (mop e) 'mexpt) (sqrtdenest2 e))
  72    (t `((,(mop e)) ,@(mapcar #'sqrtdenest1 (rest e))))))
  73
  74 ;;; Simple denesting of square roots.  Uses algorithm in maxima share package
  75 ;;; sqdnst.mac.  Should give identical results.
  76 ;;;
  77 ;;; Ref: Jeffrey and Rich (1999), Section 4.5.2
  78 ;;;
  79 ;;; Let X,Y real with X > Y > 0
  80 ;;;
  81 ;;; Consider
  82 ;;;     sqrt(X+Y) = sqrt(A) + sqrt(B)         (4.9)
  83 ;;;
  84 ;;; Squaring both sides
  85 ;;;     X + Y = A + B + 2 sqrt(AB)           (4.10)
  86 ;;;
  87 ;;; One way to satify this is to set X=A+B and Y^2=4AB
  88 ;;;
  89 ;;; Also have
  90 ;;;     sqrt(X-Y) = sqrt(A) - sqrt(B)          (4.11)
  91 ;;;
  92 ;;; Solve for A and B in terms of X and Y to derive:
  93 ;;;
  94 ;;; Theorem 4.12: Let X,Y real with X > Y > 0
  95 ;;; sqrt(X +/- Y) = sqrt(X/2+sqrt(X^2-Y^2)/2) +/- sqrt(X/2-sqrt(X^2-Y^2)/2)
  96 ;;;
  97 ;;; Apply this below by testing discriminant D = sqrt(1-(Y/X)^2)
  98 ;;;   If $numberp(D) then theorem 4.12 denests the square-root as
  99 ;;;   sqrt(X +/- Y) =  sqrt(X*(1+D)/2) +/- *sqrt(X*(1-D)/2)
 100 ;;; Note: X>0 and D<1 so both sqrt terms on RHS are real.
 101 (defun sqrtdenest2 (e)
 102   "Denest square roots in maxima expression e of form a^b"
 103   (let ((a (simplify (sqrtdenest1 (second e))))
 104         (b (simplify (sqrtdenest1 (third e))))
 105         x y D)
 106     (cond ((and ($evenp ($denom b))
 107                 (not (mapatom a))
 108                 (eq (mop a) 'mplus)
 109                 (progn
 110                   (setq x ; maximum of args(a)
 111                         (simplify `(($max) ,@(rest a))))
 112                   (setq y (sub a x)) ; a-x
 113                   (setq D ; sqrt(1-(y/x)^2)
 114                         (root (sub 1 (power (div y x) 2)) 2))
 115                   ($numberp D)))
 116            ;; (sqrt(x*(1+D)/2)+signum(y)*sqrt(x*(1-D)/2))^(2*b)
 117            (power (add (root (mul 1//2 x (add 1 D)) 2)
 118                        (mul (take '(%signum) y)
 119                             (root (mul 1//2 x (sub 1 D)) 2)))
 120                   (mul 2 b)))
 121           ;; Didn't denest at this level, but may have at a lower level.
 122           (t (power a b)))))