[Openmcl-devel] Exposing math functions

[Steven Nunez]

Greetings all,
I've got need of a few functions that are in the standard C library, but not in common lisp. Since the functions are there, but not available, I set about wrapping them. Looking at the code for atanh, I was surprised that CCL's wrapping is so complicated. Now I get that the code, pasted below, has to cater for complex numbers, multiple architectures and the like, so I can likely simplify it for my use case in one architecture and the real domain.
So, my question is: is this best practice for getting good performance for a wrapper? In Allegro, where I don't have access to the code, things are much simpler, using their FFI, e.g.
  (ff:def-foreign-call (%erf "erf") ((x :double))
    :returning :double)
  (ff:def-foreign-call (%erff "erff") ((x :float))
    :returning :float)
  (defun erf (x)
    "Returns the error function value for x"
    (cond ((typep x 'double-float) (%erf x))
      ((typep x 'single-float) (%erff x))))

Whereas in the CCL code there's %setf-double-float, and target::with-stack-double-floats, all of which seem undocumented, but look like they're doing some magic to keep things fast and/or efficient. These functions are already in the database, so it's tempting to use the #_ macro to access them, but then I wonder why the CCL guys didn't do this.

Anyone got any ideas here? I guess I'm wondering as well whether or not I'll get any performance benefit from following the CCL pattern over a simpler approach.


#+windows-target
(progn
(defun %double-float-atanh! (n result)
  (declare (double-float n result))
  (with-stack-double-floats ((temp))
    (%setf-double-float TEMP (external-call "atanh" :double n :double))
    (%df-check-exception-1 'atanh n (%ffi-exception-status))
    (%setf-double-float result TEMP)))

#+32-bit-target
(defun %single-float-atanh! (n result)
  (declare (single-float n result)) 
  (target::with-stack-short-floats ((temp))
    #+arm-target (%set-fpscr-status 0)
    (%setf-short-float TEMP (external-call "atanhf" :float n :float))
    (%sf-check-exception-1 'atanh n (%ffi-exception-status))
    (%setf-short-float result TEMP)))

#+64-bit-target
(defun %single-float-atanh (n)
  (declare (single-float n)) 
  (let* ((result (external-call "atanhf" :float n :float)))
    (%sf-check-exception-1 'atanh n (%ffi-exception-status))
    result)))
(defun atan (y &optional (x nil x-p))
  "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
  (cond (x-p
         (cond ((or (typep x 'double-float)
                    (typep y 'double-float))
                (with-stack-double-floats ((dy y)
                                           (dx x))
                  (%df-atan2 dy dx)))
               (t
                (when (and (rationalp x) (rationalp y))
                  ;; rescale arguments so that the maximum absolute value is 1
                  (let ((x1 (abs x)) (y1 (abs y)))
                    (cond ((> y1 x1)
                           (setf x (/ x y1))
                           (setf y (signum y)))
                          ((not (zerop x))
                           (setf y (/ y x1))
                           (setf x (signum x))))))
                #+32-bit-target
                (target::with-stack-short-floats ((sy y)
                                                  (sx x))
                  (%sf-atan2! sy sx))
                #+64-bit-target
                (%sf-atan2 (%short-float y) (%short-float x)))))
        ((typep y 'double-float)
         (%double-float-atan! y (%make-dfloat)))
        ((typep y 'single-float)
         #+32-bit-target
         (%single-float-atan! y (%make-sfloat))
         #+64-bit-target
         (%single-float-atan y))
        ((typep y 'rational)
         (cond ((<= (abs y) most-positive-short-float)
                #+32-bit-target
                (target::with-stack-short-floats ((sy y))
                  (%single-float-atan! sy (%make-sfloat)))
                #+64-bit-target
                (%single-float-atan (%short-float y)))
               ((minusp y)
                #.(- single-float-half-pi))
               (t
                single-float-half-pi)))
        (t
         (let ((r (realpart y))
               (i (imagpart y)))
           (if (zerop i)
             (complex (atan r) i)
             (i* (%complex-atanh (complex (- i) r)) -1))))))

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.clozure.com/pipermail/openmcl-devel/attachments/20200618/eec8d2dd/attachment.htm>