src/logarc.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12 (macsyma-module logarc)
  13
  14 ;;;  Logarc and Halfangles
  15
  16 (defmfun $logarc (exp)
  17   (cond ((atom exp) exp)
  18         ((arcp (caar exp)) (logarc (caar exp) ($logarc (cadr exp))))
  19         ((eq (caar exp) '%atan2)
  20          (logarc '%atan2 (list ($logarc (second exp)) ($logarc (third exp)))))
  21         (t (recur-apply #'$logarc exp))))
  22
  23 (defun logarc (f x)
  24   ;; Gives the logarithmic form of arc trig and hyperbolic functions
  25   (cond ((eq f '%acos)
  26          ;; -%i * log(x + %i*sqrt(1-x^2))
  27          (mul -1 '$%i (take '(%log) (add x (mul '$%i (root (sub 1 (power x 2)) 2))))))
  28         ((eq f '%asin)
  29          ;; -%i * log(sqrt(1-x^2)+%i*x)
  30          (mul -1 '$%i (take '(%log) (add (mul '$%i x) (root (sub 1 (power x 2)) 2)))))
  31         ((eq f '%atan)
  32          ;; (log(1 + %i*x) - log(1 - %i*x)) /(2 %i)
  33          (div (sub (take '(%log) (add 1 (mul '$%i x))) (take '(%log) (sub 1 (mul '$%i x))))
  34               (mul 2 '$%i)))
  35         ((eq f '%atan2)
  36          ;; atan2(y,x) = -%i*log((x + %i*y)/sqrt(x^2+y^2))
  37          (destructuring-bind (y x)
  38              x
  39            (mul -1 '$%i
  40                 (take '(%log) (div (add x (mul '$%i y))
  41                                    (root (add (mul x x) (mul y y)) 2))))))
  42         ((eq f '%asinh)
  43          ;; log(sqrt(x^2+1)+x)
  44          (take '(%log) (add x (root (add 1 (power x 2)) 2))))
  45         ((eq f '%acosh)
  46          ;; log(x+sqrt(x-1)*sqrt(x+1))
  47          (take '(%log) (add x (mul (root (add x -1) 2) (root (add x 1) 2)))))
  48         ((eq f '%atanh)
  49          ;; (log(x+1)-log(1-x))/2
  50          (div (sub (take '(%log) (add 1 x)) (take '(%log) (sub 1 x))) 2))
  51         ((member f '(%asec %acsc %acot %asech %acsch %acoth) :test #'eq)
  52          ;; asec(x) = acos(1/x), and etc.
  53          (logarc (zl-get (zl-get (get f '$inverse) 'recip) '$inverse) (inv x)))
  54         (t (merror "LOGARC: unrecognized argument: ~M" f))))
  55
  56 ;; Conditionally apply a logarc transformation to operators that either have the
  57 ;; arcp property or that are %atan2 expressions but are *not* members of the list `l`.
  58 ;; We could blend this functionality into $logarc, but I'm not sure there is much
  59 ;; demand for it.
  60 (defun partial-logarc (e l)
  61   (cond ((atom e) e)
  62         ((and (arcp (caar e)) (not (member (caar e) l)))
  63           (logarc (caar e) (partial-logarc (cadr e) l)))
  64         ((eq (caar e) '%atan2)
  65             (logarc '%atan2 (list (partial-logarc (second e) l)
  66                                 (partial-logarc (third e) l))))
  67         (t (recur-apply #'(lambda (q) (partial-logarc q l)) e))))
  68
  69 (defun halfangle (f a)
  70   (and (mtimesp a)
  71        (ratnump (cadr a))
  72        (equal (caddr (cadr a)) 2)
  73        (halfangleaux f (mul 2 a))))
  74
  75 (defun halfangleaux (f a) ;; f=function; a=twice argument
  76   (let ((sw (member f '(%cos %cot %coth %cosh) :test #'eq)))
  77     (cond ((member f '(%sin %cos) :test #'eq)
  78            (mul (halfangleaux-factor f a)
  79                 (power (div (add 1 (porm sw (take '(%cos) a))) 2) 1//2)))
  80           ((member f '(%tan %cot) :test #'eq)
  81            (div (add 1 (porm sw (take '(%cos) a))) (take '(%sin) a)))
  82           ((member f '(%sinh %cosh) :test #'eq)
  83            (mul (halfangleaux-factor f a)
  84                 (power (div (add (take '(%cosh) a) (porm sw 1)) 2) 1//2)))
  85           ((member f '(%tanh %coth) :test #'eq)
  86            (div (add (take '(%cosh) a) (porm sw 1)) (take '(%sinh) a)))
  87           ((member f '(%sec %csc %sech %csch) :test #'eq)
  88            (inv (halfangleaux (get f 'recip) a))))))
  89
  90 (defun halfangleaux-factor (f a)
  91   (cond
  92     ((member f '(%sin %cos))
  93      (let ((arg (div (if (eq f '%sin)
  94                          ($realpart a)
  95                          (add ($realpart a) '$%pi))
  96                      (mul 2 '$%pi))))
  97        (mul
  98          (power -1 (simplify (list '($floor) arg)))
  99          (sub 1
 100            (mul
 101              (add 1
 102                (power -1 (add (simplify (list '($floor) arg))
 103                               (simplify (list '($floor) (mul -1 arg))))))
 104                (simplify (list '($unit_step) (mul -1 ($imagpart a)))))))))
 105     ((member f '(%sinh %cosh))
 106      (let ((arg (div (add ($imagpart a) '$%pi) (mul 2 '$%pi)))
 107            (fac (if (eq f '%sinh)
 108                     (div (power (power a 2) (div 1 2)) a)
 109                     1)))
 110        (mul fac
 111          (power -1 (simplify (list '($floor) arg)))
 112          (sub 1
 113            (mul
 114              (add 1
 115                (power -1 (add (simplify (list '($floor) arg))
 116                               (simplify (list '($floor) (mul -1 arg))))))
 117                (simplify (list '($unit_step) ($realpart a))))))))
 118     (t 1)))