;;; float-ieee.el --- encode and decode IEEE floats  -*- lexical-binding: t; -*-

;; Copyright (C) 2015  Leo Liu

;; Author: Leo Liu <sdl.web@gmail.com>
;; Keywords: extensions, data
;; Version: 0.1.0
;; Created: 2015-02-06

;; This program is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.

;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.

;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.

;;; Commentary:

;; Ref: http://en.wikipedia.org/wiki/IEEE_floating_point

;;; Code:

(require 'cl-lib)

;;;###autoload
(defun float-ieee-encode (type num)
  "Encode NUM as IEEE 754 binary or decimal floating point."
  (cl-ecase type
    ((binary16 :binary16 half :half)
     (float-ieee-encode-binary 16 (float num)))
    ((binary32 :binary32 single :single)
     (float-ieee-encode-binary 32 (float num)))
    ((binary64 :binary64 double :double)
     (float-ieee-encode-binary 64 (float num)))))

;;;###autoload
(defun float-ieee-decode (type bytes)
  "Decode BYTES encoded in TYPE as a floating point."
  (cl-ecase type
    ((binary :binary)
     (float-ieee-decode-binary (* (length bytes) 8) bytes))
    ((decimal :decimal) (error "Not yet implemented"))))

(eval-when-compile
  (if (ignore-errors (reverse []))
      (defalias 'float-ieee-reverse 'reverse) ;25.1
    (defsubst float-ieee-reverse (seq)
      (cl-etypecase seq
        (list (reverse seq))
        (array (let ((len (length seq))
                     (seq (copy-sequence seq)))
                 (dotimes (i (floor len 2))
                   (cl-rotatef (aref seq i) (aref seq (- len i 1))))
                 seq))))))

(defsubst float-ieee-minusp (n)
  "Return non-nil if N (including -0.0) is minus."
  ;; Use `copysign' to correctly handle -0.0.
  (< (copysign 1.0 n) 0))

(defsubst float-ieee-exponent-bits (precision)
  "Get the length of the exponent bits for PRECISION."
  (cl-ecase precision
    (16   5)
    (32   8)
    (64   11)
    (128  15)))

(defsubst float-ieee-exponent-bias (precision)
  (1- (expt 2 (1- (float-ieee-exponent-bits precision)))))

(defun float-ieee--bytes->bits (bytes)
  "Convert BYTES to a vector of 0's and 1's."
  (let (bits)
    (mapc (lambda (byte)
            (cl-check-type byte (integer 0 255))
            (cl-loop repeat 8 for x = byte then (lsh x -1) do
                     (push (logand x 1) bits)))
          (float-ieee-reverse bytes))
    (cl-coerce bits 'vector)))

(defun float-ieee--bits->bytes (bits)
  (cl-check-type bits vector)
  (cl-assert (zerop (mod (length bits) 8)) t)
  (let* ((len           (/ (length bits) 8))
         (result        (make-string len 0)))
    (cl-loop for i from 0 to (1- len) do
             (setf (aref result i)
                   (cl-loop repeat 8 for j from (* i 8)
                            sum (* (aref bits j)
                                   (expt 2 (- 7 (mod j 8)))))))
    result))

(cl-defun float-ieee--fill-exponent (bits E &key (start 0) (end (length bits)))
  (cl-loop for x = E then (lsh x -1)
           for i from (1- end) downto start
           do (setf (aref bits i) (logand x 1))))

(cl-defun float-ieee--fill-mantissa (bits M &key (start 0) (end (length bits)))
  (cl-loop for i from start to (1- end)
           for tmp = (- M (expt 0.5 (- i (1- start))))
           when (>= tmp 0) do (setf (aref bits i) 1 M tmp)))

(cl-defun float-ieee--sum-exponent (bits &key (start 0) (end (length bits)))
  (cl-loop for i from start to (1- end)
           for w = (expt 2 (- end start 1)) then (lsh w -1)
           sum (* w (aref bits i))))

(cl-defun float-ieee--sum-mantissa (bits &key (start 0) (end (length bits)))
  (cl-loop for i from start to (1- end)
           sum (* (aref bits i)
                  (expt 0.5 (- i (1- start))))))

(defun float-ieee-encode-binary (precision n)
  (let* ((EL    (float-ieee-exponent-bits precision))
         (bias  (float-ieee-exponent-bias precision))
         (E     (frexp n))
         (S     (pop E))
         (bits  (make-vector precision 0)))
    (when (float-ieee-minusp n)         ;sign
      (setf (aref bits 0) 1 S (abs S)))
    (cond ((zerop n) )
          ((isnan n) (cl-fill bits 1 :start 1))
          ((eql S 1.0e+INF) (cl-fill bits 1 :start 1 :end (1+ EL)))
          ((<= E (- 1 bias))            ;subnormals
           (float-ieee--fill-mantissa bits (* (expt 2 (+ E (1- bias))) S)
                                      :start (1+ EL)))
          (t (float-ieee--fill-exponent bits (+ (1- E) bias)
                                        :start 1 :end (1+ EL))
             (float-ieee--fill-mantissa bits (1- (* 2 S)) :start (1+ EL))))
    (float-ieee--bits->bytes bits)))

(defun float-ieee-decode-binary (precision bytes)
  (cl-assert (= (length bytes) (/ precision 8)) t)
  (let* ((bits  (float-ieee--bytes->bits bytes))
         (sign  (expt -1 (aref bits 0)))
         (EL    (float-ieee-exponent-bits precision))
         (E     (float-ieee--sum-exponent bits :start 1 :end (1+ EL)))
         (M     (float-ieee--sum-mantissa bits :start (1+ EL))))
    (cond ((= E (1- (expt 2 EL)))
           (if (zerop M) (* sign 1.0e+INF) 0.0e+NaN))
          ((zerop E)
           (* sign M (expt 2.0 (- 1 (float-ieee-exponent-bias precision)))))
          (t (* sign (1+ M)
                (expt 2 (- E (float-ieee-exponent-bias precision))))))))

(provide 'float-ieee)
;;; float-ieee.el ends here