[Openmcl-devel] Help on using FPC-compiler

[Diedrich Wolter]

Hi all,

to speed up my floating point calculation I'm using Randall Beer's 
excellent FPC compiler
(see http://vorlon.cwru.edu/~beer/). My code works fine in MCL, but the 
results are
screwed up using OpenMCL (0.14-031220). An example is attached.

Has anybody had a similar experience? Any pointers appreciated!

Best wishes,
   Diedrich


;; Standard Lisp to verify the results

(defun lin-solve (m11 m12 m21 m22 lx ly)
   "Solves the equation M * x = l, M being a 2*2 matrix. The solution x1 
and x2 is
    returned as multiple values."
   (declare (type double-float m11 m12 m21 m22 lx ly)
            (optimize (speed 3) (safety 0)))
   (let ((det (- (* m11 m22) (* m12 m21))))
     (declare (type double-float det)
              (dynamic-extent det))
     (if (<  (abs det) +eps+)
       nil
       (let ((det1 (/ 1.0d0 det)))
         (declare (type double-float det1)
                  (dynamic-extent det))
         (values (+ (* m22 det1 lx) (* (- m12) det1 ly))
                 (+ (* (- m21) det1 lx) (* m11 det1 ly)))))))


;; FPC version which *should* give the same results, just way faster

(defun lin-solve-fpc (m11 m12 m21 m22 lx ly)
   "Solves the equation M * x = l, M being a 2*2 matrix. The solution x 
and y is
    returned as multiple values."
   (declare (type double-float m11 m12 m21 m22 lx ly +eps+)
            (optimize (speed 3) (safety 0) (debug 0)))
   (ccl:with-temp-double-floats (det)
     (ccl::%set-double! det (- (* m11 m22) (* m12 m21)))
     (if (< (abs det) (the double-float +eps+))
       nil
       (ccl:with-temp-double-floats (x y)
         (ccl::%set-double! det (/ 1.0d0 det))
         (ccl:%set-double! x (+ (* m22 det lx) (* (- m12) det ly)))
         (ccl:%set-double! y (+ (* (- m21) det lx) (* m11 det ly)))
         (values (ccl:%copy-float x)
                 (ccl:%copy-float y))))))

;; Test function: should infinitely loop (and does that using MCL)
;; In OpenMCL differing results are printed out (see below)

(defun fpctest ()
   (declare (optimize (speed 0) (safety 3) (debug 3)))
   ;; (rnum) generates a random double-float
   (flet ((rnum ()
            (coerce (* 1e-4 (random 100000)) 'double-float)))
     (let ((i 0))
       (loop
         (let ((m11 (rnum))
               (m12 (rnum))
               (m21 (rnum))
               (m22 (rnum))
               (lx  (rnum))
               (ly  (rnum)))
           (multiple-value-bind (fpc-1 fpc-2) (lin-solve-fpc m11 m12 m21 
m22 lx ly)
             (multiple-value-bind (std-1 std-2) (lin-solve m11 m12 m21 
m22 lx ly)
               ;; Check, if results differ
               (when (or (and fpc-2 (not std-2))
                         (and std-2 (not fpc-2))
                         (and fpc-1 (not std-1))
                         (and std-1 (not fpc-1))
                         (not (zero? (- fpc-1 std-1)))
                         (not (zero? (- fpc-2 std-2))))
                 (format t
                         "~%~%m11 = ~a, m12 = ~a, m21 = ~a, m22 = ~a, lx 
= ~a, ly = ~a~%-> fpc-1 = ~a, fpc-2 = ~a,~%-> std-1 = ~a, std-2 = ~a"
                         m11 m12 m21 m22 lx ly fpc-1 fpc-2 std-1 
std-2)))))
         (incf i)
         (when (= 0 (mod i 100000))
           (print i))))))

;; Exemplary output:
;; m11 = 4.119199752807617D0, m12 = 1.4930999279022217D0, m21 = 
2.8643999099731445D0,
;; m22 = 9.78849983215332D0, lx = 6.079599857330322D0, ly = 
3.3651998043060303D0
;; -> fpc-1 = -5.370073764483212D0, fpc-2 = -5.370073764483212D0,
;; -> std-1 = 1.511642908277411D0, std-2 = -0.09855953645790666D0
;;
;; Observation: The FPC version always returns two identical values