external/oct/branch-test.lisp

   1 ;;;; -*- Mode: lisp -*-
   2 ;;;;
   3 ;;;; Copyright (c) 2007 Raymond Toy
   4 ;;;;
   5 ;;;; Permission is hereby granted, free of charge, to any person
   6 ;;;; obtaining a copy of this software and associated documentation
   7 ;;;; files (the "Software"), to deal in the Software without
   8 ;;;; restriction, including without limitation the rights to use,
   9 ;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
  10 ;;;; copies of the Software, and to permit persons to whom the
  11 ;;;; Software is furnished to do so, subject to the following
  12 ;;;; conditions:
  13 ;;;;
  14 ;;;; The above copyright notice and this permission notice shall be
  15 ;;;; included in all copies or substantial portions of the Software.
  16 ;;;;
  17 ;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  18 ;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  19 ;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  20 ;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  21 ;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  22 ;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  23 ;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  24 ;;;; OTHER DEALINGS IN THE SOFTWARE.
  25
  26 ;;; Some simple tests to see that we're computing the branch cuts
  27 ;;; correctly.
  28 ;;;
  29 ;;; NOTE: the tests assume that the functions for double-float are
  30 ;;; computing the values correctly for the branch cuts.  We need to
  31 ;;; fix this.
  32
  33 (in-package #:qd)
  34
  35 (defun check-signs (fun arg real-sign imag-sign)
  36   (let* ((z (funcall fun arg))
  37          (x (realpart z))
  38          (y (imagpart z)))
  39     (unless (and (= (float-sign x) real-sign)
  40                  (= (float-sign y) imag-sign))
  41       (format t "Sign of result doesn't match expected signs~%~
  42                  ~& fun = ~A~
  43                  ~& arg = ~A~
  44                  ~& res = ~A~
  45                  ~& expected = ~A ~A~%"
  46               fun arg z real-sign imag-sign))))
  47
  48 (defun get-signs (z)
  49   (values (float-sign (realpart z))
  50           (float-sign (imagpart z))))
  51
  52 ;; asin branch cut is the real axis |x| > 1.  For x < -1, it is
  53 ;; continuous with quadrant II; for x > 1, continuous with quadrant
  54 ;; IV.
  55 (defun test-asin ()
  56   ;; Check x < -1
  57   (multiple-value-bind (tr ti)
  58       (get-signs (asin #c(-2d0 +1d-20)))
  59     (check-signs #'asin -2d0 tr ti)
  60     (check-signs #'asin -2w0 tr ti)
  61     (check-signs #'asin #q-2 tr ti)
  62     (check-signs #'asin #c(-2d0 0d0) tr ti)
  63     (check-signs #'asin #c(-2w0 0w0) tr ti)
  64     (check-signs #'asin #q(-2 0) tr ti)
  65     (check-signs #'asin #c(-2d0 -0d0) tr (- ti))
  66     (check-signs #'asin #c(-2w0 -0w0) tr (- ti))
  67     (check-signs #'asin #q(-2 #q-0q0) tr (- ti))
  68     )
  69
  70   ;; Check x > 1
  71   (multiple-value-bind (tr ti)
  72       (get-signs (asin #c(2d0 -1d-20)))
  73     (check-signs #'asin 2d0 tr ti)
  74     (check-signs #'asin 2w0 tr ti)
  75     (check-signs #'asin #q2 tr ti)
  76     (check-signs #'asin #c(2d0 -0d0) tr ti)
  77     (check-signs #'asin #c(2w0 -0w0) tr ti)
  78     (check-signs #'asin #q(2 #q-0q0) tr ti)))
  79
  80 ;; acos branch cut is the real axis, |x| > 1.  For x < -1, it is
  81 ;; continuous with quadrant II; for x > 1, quadrant IV.
  82 (defun test-acos ()
  83   ;; Check x < -1
  84   (multiple-value-bind (tr ti)
  85       (get-signs (acos #c(-2d0 +1d-20)))
  86     (check-signs #'acos -2d0 tr ti)
  87     (check-signs #'acos -2w0 tr ti)
  88     (check-signs #'acos #q-2 tr ti))
  89
  90   ;; Check x > 1
  91   (multiple-value-bind (tr ti)
  92       (get-signs (acos #c(2d0 -1d-20)))
  93     (check-signs #'acos 2d0 tr ti)
  94     (check-signs #'acos 2w0 tr ti)
  95     (check-signs #'acos #q2 tr ti)))
  96
  97
  98 ;; atan branch cut is the imaginary axis, |y| > 1.  For y < -1, it is
  99 ;; continuous with quadrant IV; for x > 1, quadrant II.
 100 (defun test-atan ()
 101   ;; Check y < -1
 102   (multiple-value-bind (tr ti)
 103       (get-signs (atan #c(1d-20 -2d0)))
 104     (check-signs #'atan #c(0d0 -2d0) tr ti)
 105     (check-signs #'atan #c(0w0 -2w0) tr ti)
 106     (check-signs #'atan #q(#q0 #q-2) tr ti))
 107
 108   ;; Check y > 1
 109   (multiple-value-bind (tr ti)
 110       (get-signs (atan #c(-1d-20 2d0)))
 111     (check-signs #'atan #c(-0d0 2d0) tr ti)
 112     (check-signs #'atan #c(-0w0 2w0) tr ti)
 113     (check-signs #'atan #q(#q-0 2) tr ti)))
 114
 115
 116 (defun test-atanh ()
 117   ;; Check x < -1
 118   (multiple-value-bind (tr ti)
 119       (get-signs (atanh #c(-2d0 -1d-20)))
 120     (check-signs #'atanh -2d0 tr ti)
 121     (check-signs #'atanh -2w0 tr ti)
 122     (check-signs #'atanh #q-2 tr ti))
 123
 124   ;; Check x > 1
 125   (multiple-value-bind (tr ti)
 126       (get-signs (atanh #c(2d0 1d-20)))
 127     (check-signs #'atanh 2d0 tr ti)
 128     (check-signs #'atanh 2w0 tr ti)
 129     (check-signs #'atanh #q2 tr ti)))
 130
 131
 132