tests/float.impure.lisp

   1 ;;;; This file is for floating-point-related tests which have side
   2 ;;;; effects (e.g. executing DEFUN).
   3
   4 ;;;; This software is part of the SBCL system. See the README file for
   5 ;;;; more information.
   6 ;;;;
   7 ;;;; While most of SBCL is derived from the CMU CL system, the test
   8 ;;;; files (like this one) were written from scratch after the fork
   9 ;;;; from CMU CL.
  10 ;;;;
  11 ;;;; This software is in the public domain and is provided with
  12 ;;;; absolutely no warranty. See the COPYING and CREDITS files for
  13 ;;;; more information.
  14
  15 (cl:in-package :cl-user)
  16
  17 ;;; Hannu Rummukainen reported a CMU CL bug on cmucl-imp@cons.org 26
  18 ;;; Jun 2000. This is the test case for it.
  19 ;;;
  20 ;;; The bug was listed as "39: .. Probably the same bug exists in
  21 ;;; SBCL" for a while until Martin Atzmueller showed that it's not
  22 ;;; present after all, presumably because the bug was introduced into
  23 ;;; CMU CL after the fork. But we'll test for it anyway, in case
  24 ;;; e.g. someone inadvertently ports the bad code.
  25 (defun point39 (x y)
  26   (make-array 2
  27               :element-type 'double-float
  28               :initial-contents (list x y)))
  29
  30 (declaim (inline point39-x point39-y))
  31 (defun point39-x (p)
  32   (declare (type (simple-array double-float (2)) p))
  33   (aref p 0))
  34 (defun point39-y (p)
  35   (declare (type (simple-array double-float (2)) p))
  36   (aref p 1))
  37 (defun order39 (points)
  38   (sort points  (lambda (p1 p2)
  39                   (let* ((y1 (point39-y p1))
  40                          (y2 (point39-y p2)))
  41                     (if (= y1 y2)
  42                         (< (point39-x p1)
  43                            (point39-x p2))
  44                         (< y1 y2))))))
  45 (defun test39 ()
  46   (order39 (make-array 4
  47                        :initial-contents (list (point39 0.0d0 0.0d0)
  48                                                (point39 1.0d0 1.0d0)
  49                                                (point39 2.0d0 2.0d0)
  50                                                (point39 3.0d0 3.0d0)))))
  51 (assert (equalp (test39)
  52                 #(#(0.0d0 0.0d0)
  53                   #(1.0d0 1.0d0)
  54                   #(2.0d0 2.0d0)
  55                   #(3.0d0 3.0d0))))
  56
  57 (defun complex-double-float-ppc (x y)
  58   (declare (type (complex double-float) x y))
  59   (declare (optimize speed))
  60   (+ x y))
  61 (compile 'complex-double-float-ppc)
  62 (assert (= (complex-double-float-ppc #c(0.0d0 1.0d0) #c(2.0d0 3.0d0))
  63            #c(2.0d0 4.0d0)))
  64
  65 (defun single-float-ppc (x)
  66   (declare (type (signed-byte 32) x) (optimize speed))
  67   (float x 1f0))
  68 (compile 'single-float-ppc)
  69 (assert (= (single-float-ppc -30) -30f0))
  70
  71 ;;; constant-folding irrational functions
  72 (declaim (inline df))
  73 (defun df (x)
  74   ;; do not remove the ECASE here: the bug this checks for indeed
  75   ;; depended on this configuration
  76   (ecase x (1 least-positive-double-float)))
  77 (macrolet ((test (fun)
  78              (let ((name (intern (format nil "TEST-CONSTANT-~A" fun))))
  79                `(progn
  80                   (defun ,name () (,fun (df 1)))
  81                   (,name)))))
  82   (test sqrt)
  83   (test log)
  84   (test sin)
  85   (test cos)
  86   (test tan)
  87   (test asin)
  88   (test acos)
  89   (test atan)
  90   (test sinh)
  91   (test cosh)
  92   (test tanh)
  93   (test asinh)
  94   (test acosh)
  95   (test atanh)
  96   (test exp))
  97
  98 ;;; Broken move-arg-double-float for non-rsp frame pointers on x86-64
  99 (defun test (y)
 100   (declare (optimize speed))
 101   (multiple-value-bind (x)
 102       (labels ((aux (x)
 103                  (declare (double-float x))
 104                  (etypecase y
 105                    (double-float
 106                     nil)
 107                    (fixnum
 108                     (aux x))
 109                    (complex
 110                     (format t "y=~s~%" y)))
 111                  (values x)))
 112         (aux 2.0d0))
 113     x))
 114
 115 (assert (= (test 1.0d0) 2.0d0))
 116
 117 (deftype myarraytype (&optional (length '*))
 118   `(simple-array double-float (,length)))
 119 (defun new-pu-label-from-pu-labels (array)
 120   (setf (aref (the myarraytype array) 0)
 121         sb-ext:double-float-positive-infinity))
 122
 123 ;;; bug 407
 124 ;;;
 125 ;;; FIXME: it may be that TYPE-ERROR is wrong, and we should
 126 ;;; instead signal an overflow or coerce into an infinity.
 127 (defun bug-407a ()
 128   (loop for n from (expt 2 1024) upto (+ 10 (expt 2 1024))
 129         do (handler-case
 130                (coerce n 'single-float)
 131              (simple-type-error ()
 132                (return-from bug-407a :type-error)))))
 133 (assert (eq :type-error (bug-407a)))
 134 (defun bug-407b ()
 135   (loop for n from (expt 2 1024) upto (+ 10 (expt 2 1024))
 136         do (handler-case
 137                (format t "~E~%" (coerce n 'single-float))
 138              (simple-type-error ()
 139                (return-from bug-407b :type-error)))))
 140 (assert (eq :type-error (bug-407b)))
 141
 142 ;; 1.0.29.44 introduces a ton of changes for complex floats
 143 ;; on x86-64. Huge test of doom to help catch weird corner
 144 ;; cases.
 145 ;; Abuse the framework to also test some float arithmetic
 146 ;; changes wrt constant arguments in 1.0.29.54.
 147 (defmacro def-compute (name real-type
 148                        &optional (complex-type `(complex ,real-type)))
 149   `(defun ,name (x y r)
 150      (declare (type ,complex-type x y)
 151               (type ,real-type r))
 152      (flet ((reflections (x)
 153               (values x
 154                       (conjugate x)
 155                       (complex (- (realpart x)) (imagpart x))
 156                       (- x)))
 157             (compute (x y r)
 158               (declare (type ,complex-type x y)
 159                        (type ,real-type r))
 160               (list (1+ x) (* 2 x) (/ x 2) (= 1 x)
 161                     (+ x y) (+ r x) (+ x r)
 162                     (- x y) (- r x) (- x r)
 163                     (* x y) (* x r) (* r x)
 164                     (unless (zerop y)
 165                       (/ x y))
 166                     (unless (zerop r)
 167                       (/ x r))
 168                     (unless (zerop x)
 169                       (/ r x))
 170                     (conjugate x) (conjugate r)
 171                     (abs r) (- r) (= 1 r)
 172                     (- x) (1+ r) (* 2 r) (/ r 2)
 173                     (complex r) (complex r r) (complex 0 r)
 174                     (= x y) (= r x) (= y r) (= x (complex 0 r))
 175                     (= r (realpart x)) (= (realpart x) r)
 176                     (> r (realpart x)) (< r (realpart x))
 177                     (> (realpart x) r) (< (realpart x) r)
 178                     (eql x y) (eql x (complex r)) (eql y (complex r))
 179                     (eql x (complex r r)) (eql y (complex 0 r))
 180                     (eql r (realpart x)) (eql (realpart x) r))))
 181        (declare (inline reflections))
 182        (multiple-value-bind (x1 x2 x3 x4) (reflections x)
 183          (multiple-value-bind (y1 y2 y3 y4) (reflections y)
 184            #.(let ((form '(list)))
 185                (dolist (x '(x1 x2 x3 x4) (reverse form))
 186                  (dolist (y '(y1 y2 y3 y4))
 187                    (push `(list ,x ,y r
 188                                 (append (compute ,x ,y r)
 189                                         (compute ,x ,y (- r))))
 190                          form)))))))))
 191
 192 (def-compute compute-number real number)
 193 (def-compute compute-single single-float)
 194 (def-compute compute-double double-float)
 195
 196 (labels ((equal-enough (x y)
 197            (cond ((eql x y))
 198                  ((or (complexp x)
 199                       (complexp y))
 200                   (or (eql (coerce x '(complex double-float))
 201                            (coerce y '(complex double-float)))
 202                       (and (equal-enough (realpart x) (realpart y))
 203                            (equal-enough (imagpart x) (imagpart y)))))
 204                  ((numberp x)
 205                   (or (eql (coerce x 'double-float) (coerce y 'double-float))
 206                       (< (abs (- x y))  1d-5))))))
 207   (let* ((reals     '(0 1 2))
 208          (complexes '#.(let ((reals '(0 1 2))
 209                              (cpx   '()))
 210                          (dolist (x reals (nreverse cpx))
 211                            (dolist (y reals)
 212                              (push (complex x y) cpx))))))
 213     (declare (notinline every))
 214     (dolist (r reals)
 215       (dolist (x complexes)
 216         (dolist (y complexes)
 217           (let ((value  (compute-number x y r))
 218                 (single (compute-single (coerce x '(complex single-float))
 219                                         (coerce y '(complex single-float))
 220                                         (coerce r 'single-float)))
 221                 (double (compute-double (coerce x '(complex double-float))
 222                                         (coerce y '(complex double-float))
 223                                         (coerce r 'double-float))))
 224             (assert (every (lambda (pos ref single double)
 225                              (declare (ignorable pos))
 226                              (every (lambda (ref single double)
 227                                       (or (and (equal-enough ref single)
 228                                                (equal-enough ref double))
 229                                           (and (not (numberp single)) ;; -ve 0s
 230                                                (equal-enough single double))))
 231                                     (fourth ref) (fourth single) (fourth double)))
 232                            '((0 0) (0 1) (0 2) (0 3)
 233                              (1 0) (1 1) (1 2) (1 3)
 234                              (2 0) (2 1) (2 2) (2 3)
 235                              (3 0) (3 1) (3 2) (3 3))
 236                            value single double))))))))
 237
 238 ;; The x86 port used not to reduce the arguments of transcendentals
 239 ;; correctly.
 240 ;; This test is valid only for x86: The x86 port uses the builtin x87
 241 ;; FPU instructions to implement the trigonometric functions; other
 242 ;; ports rely on the system's math library. These two differ in the
 243 ;; precision of pi used for the range reduction and so yield results
 244 ;; that can differ by arbitrarily large amounts for large inputs.
 245 ;; The test expects the x87 results.
 246 (with-test (:name (:range-reduction :x87)
 247             :skipped-on '(not :x86))
 248   (flet ((almost= (x y)
 249            (< (abs (- x y)) 1d-5)))
 250     (macrolet ((foo (op value)
 251                  `(let ((actual (,op ,value))
 252                         (expected (,op (mod ,value (* 2 pi)))))
 253                     (unless (almost= actual expected)
 254                       (error "Inaccurate result for ~a: expected ~a, got ~a"
 255                              (list ',op ,value) expected actual)))))
 256       (let ((big (* pi (expt 2d0 70)))
 257             (mid (coerce most-positive-fixnum 'double-float))
 258             (odd (* pi most-positive-fixnum)))
 259         (foo sin big)
 260         (foo sin mid)
 261         (foo sin odd)
 262         (foo sin (/ odd 2d0))
 263
 264         (foo cos big)
 265         (foo cos mid)
 266         (foo cos odd)
 267         (foo cos (/ odd 2d0))
 268
 269         (foo tan big)
 270         (foo tan mid)
 271         (foo tan odd)))))
 272
 273 ;; To test the range reduction of trigonometric functions we need a much
 274 ;; more accurate approximation of pi than CL:PI is. Calculating this is
 275 ;; more fun than copy-pasting a constant and Gauss-Legendre converges
 276 ;; extremely fast.
 277 (defun pi-gauss-legendre (n-bits)
 278   "Return a rational approximation to pi using the Gauss-Legendre
 279 algorithm. The calculations are done with integers, representing
 280 multiples of (expt 2 (- N-BITS)), and the result is an integral multiple
 281 of this number. The result is accurate to a few less than N-BITS many
 282 fractional bits."
 283   (let ((a (ash 1 n-bits))                     ; scaled 1
 284         (b (isqrt (expt 2 (1- (* n-bits 2))))) ; scaled (sqrt 1/2)
 285         (c (ash 1 (- n-bits 2)))               ; scaled 1/4
 286         (d 0))
 287     (loop
 288       (when (<= (- a b) 1)
 289         (return))
 290       (let ((a1 (ash (+ a b) -1)))
 291         (psetf a a1
 292                b (isqrt (* a b))
 293                c (- c (ash (expt (- a a1) 2) (- d n-bits)))
 294                d (1+ d))))
 295     (/ (round (expt (+ a b) 2) (* 4 c))
 296        (ash 1 n-bits))))
 297
 298 ;; Test that the range reduction of trigonometric functions is done
 299 ;; with a sufficiently accurate value of pi that the reduced argument
 300 ;; is correct to nearly double-float precision even for arguments of
 301 ;; very large absolute value.
 302 ;; This test is skipped on x86; as to why see the comment at the test
 303 ;; (:range-reduction :x87) above.
 304 (with-test (:name (:range-reduction :precise-pi)
 305             :skipped-on :x86
 306             :fails-on '(and :openbsd :x86-64))
 307   (let ((rational-pi-half (/ (pi-gauss-legendre 2200) 2)))
 308     (labels ((round-pi-half (x)
 309                "Return two values as if (ROUND X (/ PI 2)) was called
 310                 but where PI is precise enough that for all possible
 311                 double-float arguments the quotient is exact and the
 312                 remainder is exact to double-float precision."
 313                (declare (type double-float x))
 314                (multiple-value-bind (q r)
 315                    (round (rational x) rational-pi-half)
 316                  (values q (coerce r 'double-float))))
 317              (expected-val (op x)
 318                "Calculate (OP X) precisely by shifting the argument by
 319                 an integral multiple of (/ PI 2) into the range from
 320                 (- (/ PI 4)) to (/ PI 4) and applying the phase-shift
 321                 formulas for the trigonometric functions. PI here is
 322                 precise enough that the result is exact to double-float
 323                 precision."
 324                (labels ((precise-val (op q r)
 325                           (ecase op
 326                             (sin (let ((x (if (zerop (mod q 2))
 327                                               (sin r)
 328                                               (cos r))))
 329                                    (if (<= (mod q 4) 1)
 330                                        x
 331                                        (- x))))
 332                             (cos (precise-val 'sin (1+ q) r))
 333                             (tan (if (zerop (mod q 2))
 334                                      (tan r)
 335                                      (/ (- (tan r))))))))
 336                  (multiple-value-bind (q r)
 337                      (round-pi-half x)
 338                    (precise-val op q r))))
 339              (test (op x)
 340                (let ((actual (funcall op x))
 341                      (expected (expected-val op x)))
 342                  ;; Some of the test values are chosen to lie very near
 343                  ;; to an integral multiple of pi/2 (within a distance of
 344                  ;; between 1d-11 and 1d-8), making the absolute value of
 345                  ;; their sine or cosine this small, too. The absolute
 346                  ;; value of the tangent is then either similarly small or
 347                  ;; as large as the reciprocal of this value. Therefore we
 348                  ;; measure relative instead of absolute error.
 349                  (unless (or (= actual expected 0)
 350                              (and (= (signum actual) (signum expected))
 351                                   (< (abs (/ (- actual expected)
 352                                              (+ actual expected)))
 353                                      (* 8 double-float-epsilon))))
 354                    (error "Inaccurate result for ~a: expected ~a, got ~a"
 355                           (list op x) expected actual)))))
 356       (dolist (op '(sin cos tan))
 357         (dolist (val `(,(coerce most-positive-fixnum 'double-float)
 358                        ,@(loop for v = most-positive-double-float
 359                                then (expt v 4/5)
 360                                while (> v (expt 2 50))
 361                                collect v)
 362                        ;; The following values cover all eight combinations
 363                        ;; of values slightly below or above integral
 364                        ;; multiples of pi/2 with the integral factor
 365                        ;; congruent to 0, 1, 2 or 3 modulo 4.
 366                        5.526916451564098d71
 367                        4.913896894631919d229
 368                        7.60175752894437d69
 369                        3.8335637324151093d42
 370                        1.8178427396473695d155
 371                        9.41634760758887d89
 372                        4.2766818550391727d188
 373                        1.635888515419299d28))
 374           (test op val))))))