; 


; 

; Die entscheidente Erkenntnis ist, dass sich
; Gleitkommazahlen als Paare ganzer Zahlen darstellen lassen. Wir
; wickeln um dieses Paar ganzer Zahlen rudimentäre Konstruktoren und
; Selektoren:
; 


(define make-fp
  (lambda (m e)
    (cons m e)))

(define fp-mantissa
  (lambda (f)
    (car f)))

(define fp-exponent
  (lambda (f)
    (cdr f)))

; Die volle Länge der fp-Mantisse würde sinnvoll ein Bit länger
; definiert als die Länge der Fractionfield des IEEE-Formats.
; 


(define fp-mantissa-len 53) ; mit Vorkommastelle

; Ein Paar ganze Zahlen $(m,e)$ stellt die reelle Zahl 
; $r=m2^e$ dar. Dabei ist sind $m$ und $e$ vorzeichenbehaftete
; ganze Zahlen.
; 


(define fp->number
  (lambda (f)
    (exact->inexact 
     (* (fp-mantissa f)
        (expt 2.0 (fp-exponent f))))))

; Sie $x$ nun die Zahl, die in ein Gleitkommazahl umgewandelt werden
; soll.  Der in der Vorlesung besprochene Algorithmus zur Umwandlung von
; {\em poitiven} Schemezahlen in Gleitkommazahlen besteht aus folgenden
; Schritten:
; \begin{enumerate}
; \item Setze die gebrochene Mantisse $r=x$. Setze den Exponenten $e=0$.
;       Die ganzzahligen Mantisse $m=0$, die wir erst aufbauen wollen,
;       soll ebenfalls $m=0$ gesetzt werden, auch wenn wir sie erst sehr viel
;       später brauchen.
; \item Normalisiere nun $r$, indem solange $r$ mit zwei multipliziert oder 
;       durch zwei geteilt wird,
;       bis $1\leq{}r<2$ und der Exponent entsprechend jeweils um 1 verringert
;       oder vergrössert wird. Also entweder
;       \[
;         (m,e) \longrightarrow (m/2,e+1)
;       \]
;       oder
;       \[
;         (m,e) \longrightarrow (2m,e-1).
;       \]
; 

(define number->fp-1 
  (lambda (s r m e)      ; r<>0
    (cond ((< r 1.0)
           (number->fp-1 s (* r 2.0) m (- e 1)))
          ((>= r 2.0)
           (number->fp-1 s (* r 0.5) m (+ e 1)))
          (else
           (number->fp-2 s r m e)))))
; 
; \item Nun wird die ganzzahlige Mantisse aufgebaut. Dies geschieht, indem
;   die einzelnen Bits von $r$ nacheinander über den Dezimalpunkt geschoben
;   werden, und jeweils in die ganzzahlige Mantisse eingebaut werden.
;   Es werden also die folgenden Schritte durchgeführt, solange, bis
;   $r=0$ ist.
;   \begin{description}
;   \item[Wenn $r\geq1$:] $(r,m,e)\longrightarrow(2(r-1),2m+1,e-1)$
;   \item[Wenn $r<1$:]    $(r,m,e)\longrightarrow(2r    ,2m  ,e-1)$
;   \end{description}
; 

(define number->fp-2
  (lambda (s r m e)      ; 1<r<=2
    (cond ((= 0.0 r)
           (make-fp (* s m) (+ e 1)))
          ((< r 1)
           (number->fp-2 s (* 2.0 r) (* m 2) (- e 1)))
          ((>= r 1)
           (number->fp-2 s (* 2.0 (- r 1.0)) 
                         (+ 1 (* m 2)) (- e 1))))))
; Mit Hilfe eine Wrapperfunktion lässt nun sich der allgemeine Fall auf den
; der positiven Schemezahlen zurückführen:
; 

(define number->fp
  (lambda (r)
    (if (= r 0.0)
        (make-fp 0 0)
        (let ((s (if (< r 0) -1 1 )))
          (number->fp-1 s (* r s) 0 0)))))
; Für das Produkt zweier durch Mantisse und Exponent dargestellter Zahlen gilt,
;   dass 
;   \begin{equation}
;       {r_1}{r_2}={m_1}{m_2}2^{e_1+e_2}.
;   \end{equation}
;   Wir erhalten also eine Darstellung des Produktes, wenn wir
;   $m={m_1}{m_2}$ und $e={e_1}+{e_2}$ setzen. Die Funktion
;   ; fp-multiply;  implementiert dieses Verfahren.
; 


(define fp-multiply
  (lambda (f1 f2)
    (make-fp
     (* (fp-mantissa f1) (fp-mantissa f2))
     (+ (fp-exponent f1) (fp-exponent f2)))))

; Analog wird bei der Addition die Addition auf den reellen Zahlen
;   auf Rechenoperationen mit der Darstellung zurückgeführt. Wir nehmen
;   dabei o.~B.~d.~A. an, dass $e_1>=e_2$, und erhalten:
;   \[
;        {r_1}+{r_2}={m_1}2^{e_1}+{m_2}2^{e_2}=
;                    ({m_1}2^{e_1-e_2}+{m_2})2^{e_2}  
;   \]
; 
;    Der Wert des Ausdrucks in der Klammer ist dabei eine ganze Zahl
;    (und genau deswegen wurde diese Klammerung so gewaehlt),
;    so dass wir eine Darstellung des Ergebnisses erhalten, wenn wir
;    $m={m_1}2^{e_1-e_2}+{m_2}$ und $e=e_2$ setzen. ; fp-add; 
;    implementiert dieses Verfahren.
; 
; 

(define fp-add
 (lambda (f1 f2)
   (let* ((result-mantissa
           (lambda (m1 m2 diff-e)
             (+ m2
                (* m1
                   (expt 2 diff-e)))))
          (m1 (fp-mantissa f1))
          (m2 (fp-mantissa f2))
          (e1 (fp-exponent f1))
          (e2 (fp-exponent f2)))
     (if (> e1 e2)
         (make-fp 
          (result-mantissa m1 m2 (- e1 e2))
          e2)
         (make-fp
          (result-mantissa m2 m1 (- e2 e1))
          e1)))))


; Das Fragment ; types.scm;  enthält einige Routinen zur Anfertigung
; getypter Werte.
; 


(define make-type
  (let ((type-id 0))
    (lambda (name)
      (set! type-id (+ 1 type-id))
      (cons type-id name))))

(define make-typed-object
  (lambda (type value)
    (cons type value)))

(define typed-object-value
  (lambda (type object)
    (if (eq? type (car object))
        (cdr object)
        (error "type mismatch"))))

(define typed-object-predicate
  (lambda (type)
    (lambda (value)
      (and (pair? value)
           (eq? type (car value))))))


; Der Betrag der Mantisse, Exponent und Vorzeichen werden in 
; IEEE-Kodierung in 3 separaten Bitfeldern verschiedener Länge
; gespeichert. Die Länge des Vorzeichenfeldes ist selbstverständlich 1.
; Für das hier relevante {\em Double\/} Format beträgt die Länge des
; Feldes für die Mantisse 52 Bit, für den Exponenten stehen 11 Bit zur
; Verfügung. Das Mantissenfeld wird hier {\em Nachkommafeld} nennen 
; (im Programm ; fractionfield; ), da
; in diesem Feld nicht die vollständige Mantisse abgespeichert wird,
; sondern lediglich der Nachkommateil einer auf das Intervall 
; $]1.0,2.0]$ normierten Mantisse.
; 

(define ieee754-exponent-field-bits
  11)
(define ieee754-fraction-field-bits
  52)

; Zwei Werte des Exponentenfeldes, der maximale (nur Einsen)
;   und der minimale (nur Nullen), sind im IEEE-Standard reserviert,
;   um besondere Zahlen zu signalisieren. Der maximale Werte 
;   zeigt NaN's und Undendlichkeiten an (die wir hier aussen vorlassen),
;   der minimale Wert zeigt {\em denormalisierte} Zahlen an.
; 

(define ieee754-nan-exponent-field
  (- (expt 2 ieee754-exponent-field-bits) 1))
(define ieee754-exponent-field-denormal 0)

; Mit Ausnahme der Spezialfälle erhält man den Exponenten einer
;   Darstellung $(r,e)$, welche die Zahl ${r}2^e$ kodiert aus dem Exponentenfeld,
;   indem vom Exponentenfeld einen {\em Bias} $B$ abzieht. Man spricht 
;   von einer {\em biased\/} oder {\em Exzessdarstellung\/} des Exponenten.
;   Der Bias ist bei IEEE i.~d.~R.\ die Hälfte des Wertebereichs des
;   Exponentenfeldes.
; 

(define ieee754-exponent-bias 
  (- (expt 2 (- ieee754-exponent-field-bits 1)) 1))

; Der Bias bezieht sich auf eine Darstellung, in der die Mantisse $r$
;   in der Regel zwischen 0 und 1 normiert ist. Unser Gleitkommatyp ; fp; 
;   hat jedoch als Mantisse nur ganze Zahlen. Um diese zu erhalten, 
;   multiplizieren wir später beim dekodieren der IEEE-Darstellung mit
;   $2^p$, dabei fest im Hinterkopf behaltend, dass
;   die mit den oben angeführten Beziehungen berechnete
;   Mantisse aus dem  Intervall $[0.0,2.0]$ (das schliesst jetzt den 
;   denormalisierten Fall mit ein) maximal $p$ Bits nach dem Komma hat.
;   Wir werden also dadurch eine ganze Zahl in $[0,2^p]$ erhalten.
;   Diese Renormierung der Mantisse muss durch eine Vergr\"osserung
;   des Exponenten um $p$ ausgeglichen werden.
;   Wir können also gleich mit einem entsprechend
;   grösseren Bias rechnen:
; 

(define ieee754-real-exponent-bias 
  (+ ieee754-exponent-bias ieee754-fraction-field-bits))        

; Nun ergeben sich die minimalen und die maximalen Exponenten
;   welche wir bei nicht denormalen Zahl erhalten können.
; 


(define ieee754-minimal-exponent
  (- 1 ieee754-real-exponent-bias))

(define ieee754-maximal-exponent
  (- (- ieee754-nan-exponent-field 1) 
     ieee754-real-exponent-bias))


; % Die von uns verwendete Mantisse wurde noch mit $2^p$ multipliziert, um eine
;   % ganze Zahl zu erhalten (für den Typ ; fp;  ist also die Mantisse keine
;   % Zahl in $[1,2[$ sondern eine Zahl im Intervall $[2^p,2^{p+1}[$.
;   Dementsprechend ist in unserer Darstellung der implizierte Teil der
;   Mantisse ist nicht 1, sondern $2^p$.
; 

(define ieee754-implied-mantissa-part
  (expt 2 ieee754-fraction-field-bits))

(define ieee754-maximal-mantissa
  (- (* 2 ieee754-implied-mantissa-part) 1))

; Die vorliegende Implementation interpretiert die Bitfelder 
;   der IEEE-Kodierung als positive ganze Zahlen und kapselt sie
;   in einem getypten Schemeobjekt.
; 


(define ieee754-type
  (make-type 'ieee754))

(define make-ieee754
  (lambda (S F E)
    (make-typed-object 
     ieee754-type
     (list S F E))))

(define ieee754?
  (typed-object-predicate ieee754-type))


(define ieee754-sign-field
  (lambda (i3e)
    (car (typed-object-value ieee754-type i3e))))

(define ieee754-fraction-field
  (lambda (i3e)
    (car (cdr (typed-object-value ieee754-type i3e)))))

(define ieee754-exponent-field
  (lambda (i3e)
    (car (cddr (typed-object-value ieee754-type i3e)))))


; Die folgende Routine implementiert schliesslich das Verfahren zur 
;   Dekodierung der  IEEE-Darstellung.
; 


(define ieee754->fp ; decode ieee754 representation
  (lambda (i3e)
    (let* ((sign 
            (if  
             (= 1 (ieee754-sign-field i3e))
             -1
             1))
           (fraction-field
            (ieee754-fraction-field i3e))
           (exponent-field
            (ieee754-exponent-field i3e)))
      (cond ((= exponent-field ieee754-nan-exponent-field)
             (error "ieee754->fp: not a number."))
            ((= exponent-field ieee754-exponent-field-denormal)
             (make-fp
              (* sign (fraction-field))
              (ieee754-exponent-denormal)))
            (else
              (make-fp
               (* sign (+ fraction-field ieee754-implied-mantissa-part))
               (- exponent-field ieee754-real-exponent-bias)))))))

; Die Kodierung von Mantissen-Exponenten-Paaren, wie wir sie bereits
; hatten, in IEEE-Format gestaltet sich nicht ganz so trivial. Dabei
; muss Rücksicht genommen werden auf die beschränkte Länge der Felder
; des IEEE-Formats. Wir müssen deshalb Mantissen-Exponenten-Paare erst
; 'normalisieren', d.~h.\ wir müssen unter Ausnutzung der Tatsache, dass
; $(m,e)$ und $(m*2,e-1)$ dieselbe relle Zahl darstellen, die Mantisse
; kurz genug machen, dass sie in das Fractionfeld der IEEE-Darstellung
; passt, bzw. auch vorne (auf dem $p+1$\/-ten Bit) eine 1 trägt, die wir
; zur Umwandlung in implizite Darstellung abschneiden können.  
; 
; Beim Versuch die Mantisse zu kürzen, könnten wir auf den Umstand stossen,
; dass die Mantisse, die wir verkleinern wollen, eine ungerade Zahl ist,
; also nicht durch zwei geteilt werden kann, so dass wir wieder eine
; ganze Zahl erhalten. Hier müssen wir einen Verlust an Genauigkeit in
; Kauf nehmen, indem wir trotzdem teilen, dabei jedoch das letzte Bit
; (also den Nachkommateil des Ergebnisses unter den Tisch fallen lassen.
; Dieses Verfahren wird von der folgenden Routine implementiert, die
; Tripel aus Vorzeichen, Mantisse und Exponent normalisiert und in
; IEEE-Format kodiert.
; 
; Die vorliegende Funktion ruft sich solange selbst rekursiv auf, bis 
; Mantisse und Exponent normalisiert sind.
; 

(define ieee754-sign-bit
  (lambda (sign)
    (if (> 0 sign) 1 0)))

; 

(define s/m/e->ieee754
  (lambda (sign mantissa exponent)
    (cond ((= 0 mantissa)
           (begin
             (make-ieee754 
              0 0 ieee754-exponent-field-denormal)))
          ((< exponent ieee754-minimal-exponent)
           (begin
             (s/m/e->ieee754
              sign
              (quotient mantissa 2) ; potential loss of precision
              (+ exponent 1))))
          ((> exponent ieee754-maximal-exponent)
           (begin
             (s/m/e->ieee754
              sign
              (* mantissa 2)
              (- exponent 1))))
          ((> mantissa ieee754-maximal-mantissa)
           (begin
             (if (< exponent ieee754-maximal-exponent)
                 (begin
                   (s/m/e->ieee754
                    sign
                    (quotient mantissa 2)
                    (+ exponent 1)))
                 (error "s/m/e->ieee754: overflow"))))
          (else         ; e in range 
           (begin
             (if (< mantissa ieee754-implied-mantissa-part) 
                 (if (> exponent ieee754-minimal-exponent)
                     (s/m/e->ieee754
                      sign
                      (* mantissa 2)
                      (- exponent 1))
                     (make-ieee754
                      (ieee754-sign-bit sign)
                      mantissa
                      ieee754-exponent-field-denormal))
                 (make-ieee754
                  (ieee754-sign-bit sign)
                  (- mantissa ieee754-implied-mantissa-part)
                  (+ exponent ieee754-real-exponent-bias))))))))

; Zum einen wollen wir unseren einfachen Gleitkommatyp ; fp; 
; in ; ieee754;  umwandeln. Dazu benötigen wir nun nur noch eine kleine
; Wrapperroutine.
; 

(define fp->ieee754
  (lambda (r)
    (let* ((signed-mantissa (fp-mantissa r))
           (sign 
            (if (< signed-mantissa 0)
                -1
                1))
           (mantissa (* signed-mantissa sign))
           (exponent (fp-exponent r)))
      (s/m/e->ieee754 sign mantissa exponent))))


; Zum anderen wollen wir Schemezahlen in den neuen ; ieee754;  Typ umwandeln.
; Wir erledigen das, indem wir über unseren vereinfachten Gleitkommatyp ; fp;  
; gehen.
; 

(define number->ieee754
  (lambda (x)
    (fp->ieee754 (number->fp (exact->inexact x)))))

(define ieee754->number
  (lambda (i3e)
    (fp->number (ieee754->fp i3e))))


; Die Rechenoperationen lassen sich ebenso als einfache Wrapper um die
;   Rechenoperationen für den Typ ; fp;  implementieren.
; 


(define ieee754+
  (lambda (x y)
    (fp->ieee754 (fp-add (ieee754->fp x)(ieee754->fp y)))))

(define ieee754*
  (lambda (x y)
    (fp->ieee754 (fp-multiply (ieee754->fp x)(ieee754->fp y)))))