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 (defun halfangle (f a)
  57   (and (mtimesp a)
  58        (ratnump (cadr a))
  59        (equal (caddr (cadr a)) 2)
  60        (halfangleaux f (mul 2 a))))
  61
  62 (defun halfangleaux (f a) ;; f=function; a=twice argument
  63   (let ((sw (member f '(%cos %cot %coth %cosh) :test #'eq)))
  64     (cond ((member f '(%sin %cos) :test #'eq)
  65            (mul (halfangleaux-factor f a)
  66                 (power (div (add 1 (porm sw (take '(%cos) a))) 2) 1//2)))
  67           ((member f '(%tan %cot) :test #'eq)
  68            (div (add 1 (porm sw (take '(%cos) a))) (take '(%sin) a)))
  69           ((member f '(%sinh %cosh) :test #'eq)
  70            (mul (halfangleaux-factor f a)
  71                 (power (div (add (take '(%cosh) a) (porm sw 1)) 2) 1//2)))
  72           ((member f '(%tanh %coth) :test #'eq)
  73            (div (add (take '(%cosh) a) (porm sw 1)) (take '(%sinh) a)))
  74           ((member f '(%sec %csc %sech %csch) :test #'eq)
  75            (inv (halfangleaux (get f 'recip) a))))))
  76
  77 (defun halfangleaux-factor (f a)
  78   (cond
  79     ((member f '(%sin %cos))
  80      (let ((arg (div (if (eq f '%sin)
  81                          ($realpart a)
  82                          (add ($realpart a) '$%pi))
  83                      (mul 2 '$%pi))))
  84        (mul
  85          (power -1 (simplify (list '($floor) arg)))
  86          (sub 1
  87            (mul
  88              (add 1
  89                (power -1 (add (simplify (list '($floor) arg))
  90                               (simplify (list '($floor) (mul -1 arg))))))
  91                (simplify (list '($unit_step) (mul -1 ($imagpart a)))))))))
  92     ((member f '(%sinh %cosh))
  93      (let ((arg (div (add ($imagpart a) '$%pi) (mul 2 '$%pi)))
  94            (fac (if (eq f '%sinh)
  95                     (div (power (power a 2) (div 1 2)) a)
  96                     1)))
  97        (mul fac
  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) ($realpart a))))))))
 105     (t 1)))