home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / Misc / CLISP-1.LHA / CLISP960530-sr.lha / src / sfloat.d < prev    next >
Encoding:
Text File  |  1996-04-15  |  32.1 KB  |  772 lines
  1. # Grundfunktionen für Short-Floats
  2.  
  3. # Entpacken eines Short-Float:
  4. # SF_decode(obj, zero_statement, sign=,exp=,mant=);
  5. # zerlegt ein Short-Float obj.
  6. # Ist obj=0.0, wird zero_statement ausgeführt.
  7. # Sonst: signean sign = Vorzeichen (0 = +, -1 = -),
  8. #        sintWL exp = Exponent (vorzeichenbehaftet),
  9. #        uintL mant = Mantisse (>= 2^SF_mant_len, < 2^(SF_mant_len+1))
  10.   #if defined(GNU) && defined(MC680X0) && !defined(NO_ASM) && (SF_exp_shift==16) && (SF_exp_len==8)
  11.     #define SF_uexp(x)  \
  12.       ({var register uint32 __x = as_oint(x);            \
  13.         var register uint8 __uexp;                       \
  14.         __asm__("swap %0" : "=d" (__uexp) : "0" (__x) ); \
  15.         __uexp;                                          \
  16.        })
  17.   #else
  18.     #define SF_uexp(x)  ((as_oint(x) >> SF_exp_shift) & (bit(SF_exp_len)-1))
  19.   #endif
  20.   #define SF_decode(obj, zero_statement, sign_zuweisung,exp_zuweisung,mant_zuweisung)  \
  21.     { var reg1 object _x = (obj);                                         \
  22.       var reg2 uintBWL uexp = SF_uexp(_x);                                \
  23.       if (uexp==0)                                                        \
  24.         { zero_statement } # e=0 -> Zahl 0.0                              \
  25.         else                                                              \
  26.         { exp_zuweisung (sintWL)((uintWL)uexp - SF_exp_mid); # Exponent   \
  27.           unused (sign_zuweisung R_sign(_x));                # Vorzeichen \
  28.           mant_zuweisung (SF_mant_len == 16                  # Mantisse   \
  29.                           ? highlow32( bit(SF_mant_len-16), (as_oint(_x) >> SF_mant_shift) & (bit(SF_mant_len)-1) ) \
  30.                           : (bit(SF_mant_len) | ((uint32)(as_oint(_x) >> SF_mant_shift) & (bit(SF_mant_len)-1)) ) \
  31.                          );                                               \
  32.     }   }
  33.  
  34. # Einpacken eines Short-Float:
  35. # encode_SF(sign,exp,mant, ergebnis=);
  36. # liefert ein Short-Float.
  37. # > signean sign: Vorzeichen, 0 für +, -1 für negativ.
  38. # > sintWL exp: Exponent
  39. # > uintL mant: Mantisse, sollte >= 2^SF_mant_len und < 2^(SF_mant_len+1) sein.
  40. # < object ergebnis: ein Short-Float
  41. # Der Exponent wird auf Überlauf/Unterlauf getestet.
  42.   #define encode_SF(sign,exp,mant, erg_zuweisung)  \
  43.     { if ((exp) < (sintWL)(SF_exp_low-SF_exp_mid))  \
  44.         { if (underflow_allowed())                  \
  45.             { fehler_underflow(); }                 \
  46.             else                                    \
  47.             { erg_zuweisung SF_0; }                 \
  48.         }                                           \
  49.       else                                          \
  50.       if ((exp) > (sintWL)(SF_exp_high-SF_exp_mid)) \
  51.         { fehler_overflow(); }                      \
  52.       else                                          \
  53.       erg_zuweisung as_object                       \
  54.         (   ((oint)SF_type << oint_type_shift)      \
  55.           | ((soint)(sign) & wbit(vorz_bit_o))      \
  56.           | ((oint)((uint8)((exp)+SF_exp_mid) & (bit(SF_exp_len)-1)) << SF_exp_shift) \
  57.           | ((oint)((uint32)((mant) & (bit(SF_mant_len)-1))) << SF_mant_shift) \
  58.         );                                          \
  59.     }
  60.  
  61. # SF_zerop(x) stellt fest, ob ein Short-Float x = 0.0 ist.
  62.   #define SF_zerop(x)  (eq(x,SF_0))
  63.  
  64. # Liefert zu einem Short-Float x : (ftruncate x), ein SF.
  65. # SF_ftruncate_SF(x)
  66. # x wird zur 0 hin zur nächsten ganzen Zahl gerundet.
  67.   local object SF_ftruncate_SF (object x);
  68. # Methode:
  69. # x = 0.0 oder e<=0 -> Ergebnis 0.0
  70. # 1<=e<=16 -> letzte (17-e) Bits der Mantisse auf 0 setzen,
  71. #             Exponent und Vorzeichen beibehalten
  72. # e>=17 -> Ergebnis x
  73.   local object SF_ftruncate_SF(x)
  74.     var reg2 object x;
  75.     { var reg1 uintBWL uexp = SF_uexp(x); # e + SF_exp_mid
  76.       if (uexp <= SF_exp_mid) # 0.0 oder e<=0 ?
  77.         { return SF_0; }
  78.         else
  79.         { if (uexp > SF_exp_mid+SF_mant_len) # e > 16 ?
  80.             { return x; }
  81.             else
  82.             { return as_object
  83.                 ( as_oint(x) & # Bitmaske: Bits 16-e..0 gelöscht, alle anderen gesetzt
  84.                   ~(wbit(SF_mant_len+SF_mant_shift + 1+SF_exp_mid-uexp) - wbit(SF_mant_shift))
  85.                 );
  86.     }   }   }
  87.  
  88. # Liefert zu einem Short-Float x : (futruncate x), ein SF.
  89. # SF_futruncate_SF(x)
  90. # x wird von der 0 weg zur nächsten ganzen Zahl gerundet.
  91.   local object SF_futruncate_SF (object x);
  92. # Methode:
  93. # x = 0.0 -> Ergebnis 0.0
  94. # e<=0 -> Ergebnis 1.0 oder -1.0, je nach Vorzeichen von x.
  95. # 1<=e<=16 -> Greife die letzten (17-e) Bits von x heraus.
  96. #             Sind sie alle =0 -> Ergebnis x.
  97. #             Sonst setze sie alle und erhöhe dann die letzte Stelle um 1.
  98. #             Kein Überlauf der 16 Bit -> fertig.
  99. #             Sonst (Ergebnis eine Zweierpotenz): Mantisse := .1000...000,
  100. #               e:=e+1. (Test auf Überlauf wegen e<=17 überflüssig)
  101. # e>=17 -> Ergebnis x.
  102.   local object SF_futruncate_SF(x)
  103.     var reg2 object x;
  104.     { var reg1 uintBWL uexp = SF_uexp(x); # e + SF_exp_mid
  105.       if (uexp==0) # 0.0 ?
  106.         { return x; }
  107.       if (uexp <= SF_exp_mid) # e<=0 ?
  108.         { # Exponent auf 1, Mantisse auf .1000...000 setzen.
  109.           return as_object
  110.                  ( (as_oint(x) & ~(wbitm(oint_data_len+oint_data_shift)-wbit(oint_data_shift)))
  111.                   | ((oint)(SF_exp_mid+1) << SF_exp_shift)
  112.                   | ((oint)0 << SF_mant_shift)
  113.                  );
  114.         }
  115.         else
  116.         { if (uexp > SF_exp_mid+SF_mant_len) # e > 16 ?
  117.             { return x; }
  118.             else
  119.             { var reg1 oint mask = # Bitmaske: Bits 16-e..0 gesetzt, alle anderen gelöscht
  120.                 wbit(SF_mant_len+SF_mant_shift + 1+SF_exp_mid-uexp) - wbit(SF_mant_shift);
  121.               if ((as_oint(x) & mask)==0) # alle diese Bits =0 ?
  122.                 { return x; }
  123.               return as_object
  124.                      ((as_oint(x) | mask) # alle diese Bits setzen
  125.                       + wbit(SF_mant_shift) # letzte Stelle erhöhen, dabei evtl. Exponenten incrementieren
  126.                      );
  127.     }   }   }
  128.  
  129. # Liefert zu einem Short-Float x : (fround x), ein SF.
  130. # SF_fround_SF(x)
  131. # x wird zur nächsten ganzen Zahl gerundet.
  132.   local object SF_fround_SF (object x);
  133. # Methode:
  134. # x = 0.0 oder e<0 -> Ergebnis 0.0
  135. # 0<=e<=16 -> letzte (17-e) Bits der Mantisse wegrunden,
  136. #             Exponent und Vorzeichen beibehalten.
  137. # e>16 -> Ergebnis x
  138.   local object SF_fround_SF(x)
  139.     var reg2 object x;
  140.     { var reg1 uintBWL uexp = SF_uexp(x); # e + SF_exp_mid
  141.       if (uexp < SF_exp_mid) # x = 0.0 oder e<0 ?
  142.         { return SF_0; }
  143.         else
  144.         { if (uexp > SF_exp_mid+SF_mant_len) # e > 16 ?
  145.             { return x; }
  146.             else
  147.             if (uexp > SF_exp_mid+1) # e>1 ?
  148.               { var reg4 oint bitmask = # Bitmaske: Bit 16-e gesetzt, alle anderen gelöscht
  149.                   wbit(SF_mant_len+SF_mant_shift + SF_exp_mid-uexp);
  150.                 var reg3 oint mask = # Bitmaske: Bits 15-e..0 gesetzt, alle anderen gelöscht
  151.                   bitmask - wbit(SF_mant_shift);
  152.                 if ( ((as_oint(x) & bitmask) ==0) # Bit 16-e =0 -> abrunden
  153.                      || ( ((as_oint(x) & mask) ==0) # Bit 16-e =1 und Bits 15-e..0 >0 -> aufrunden
  154.                           # round-to-even, je nach Bit 17-e :
  155.                           && ((as_oint(x) & (bitmask<<1)) ==0)
  156.                    )    )
  157.                   # abrunden
  158.                   { mask |= bitmask; # Bitmaske: Bits 16-e..0 gesetzt, alle anderen gelöscht
  159.                     return as_object( as_oint(x) & ~mask );
  160.                   }
  161.                   else
  162.                   # aufrunden
  163.                   { return as_object
  164.                            ((as_oint(x) | mask) # alle diese Bits 15-e..0 setzen (Bit 16-e schon gesetzt)
  165.                             + wbit(SF_mant_shift) # letzte Stelle erhöhen, dabei evtl. Exponenten incrementieren
  166.                            );
  167.                   }
  168.               }
  169.             elif (uexp == SF_exp_mid+1) # e=1 ?
  170.               # Wie bei 1 < e <= 16, nur daß Bit 17-e stets gesetzt ist.
  171.               { if ((as_oint(x) & wbit(SF_mant_len+SF_mant_shift-1)) ==0) # Bit 16-e =0 -> abrunden
  172.                   # abrunden
  173.                   { return as_object( as_oint(x) & ~(wbit(SF_mant_len+SF_mant_shift)-wbit(SF_mant_shift)) ); }
  174.                   else
  175.                   # aufrunden
  176.                   { return as_object
  177.                            ((as_oint(x) | (wbit(SF_mant_len+SF_mant_shift)-wbit(SF_mant_shift))) # alle diese Bits 16-e..0 setzen
  178.                             + wbit(SF_mant_shift) # letzte Stelle erhöhen, dabei evtl. Exponenten incrementieren
  179.                            );
  180.                   }
  181.               }
  182.             else # e=0 ?
  183.               # Wie bei 1 < e <= 16, nur daß Bit 16-e stets gesetzt
  184.               # und Bit 17-e stets gelöscht ist.
  185.               { if ((as_oint(x) & (wbit(SF_mant_len+SF_mant_shift)-wbit(SF_mant_shift))) ==0)
  186.                   # abrunden von +-0.5 zu 0.0
  187.                   { return SF_0; }
  188.                   else
  189.                   # aufrunden
  190.                   { return as_object
  191.                            ((as_oint(x) | (wbit(SF_mant_len+SF_mant_shift)-wbit(SF_mant_shift))) # alle Bits 15-e..0 setzen
  192.                             + wbit(SF_mant_shift) # letzte Stelle erhöhen, dabei Exponenten incrementieren
  193.                            );
  194.               }   }
  195.     }   }
  196.  
  197. # Liefert zu einem Short-Float x : (- x), ein SF.
  198. # SF_minus_SF(x)
  199.   local object SF_minus_SF (object x);
  200. # Methode:
  201. # Falls x=0.0, fertig. Sonst Vorzeichenbit umdrehen.
  202.   local object SF_minus_SF(x)
  203.     var reg1 object x;
  204.     { return (eq(x,SF_0) ? x : as_object(as_oint(x) ^ wbit(vorz_bit_o)) ); }
  205.  
  206. # SF_SF_comp(x,y) vergleicht zwei Short-Floats x und y.
  207. # Ergebnis: 0 falls x=y, +1 falls x>y, -1 falls x<y.
  208.   local signean SF_SF_comp (object x, object y);
  209. # Methode:
  210. # x und y haben verschiedenes Vorzeichen ->
  211. #    x < 0 -> x < y
  212. #    x >= 0 -> x > y
  213. # x und y haben gleiches Vorzeichen ->
  214. #    x >=0 -> vergleiche x und y (die rechten 24 Bits)
  215. #    x <0 -> vergleiche y und x (die rechten 24 Bits)
  216.   local signean SF_SF_comp(x,y)
  217.     var reg1 object x;
  218.     var reg2 object y;
  219.     { if (!R_minusp(y))
  220.         # y>=0
  221.         { if (!R_minusp(x))
  222.             # y>=0, x>=0
  223.             { if (as_oint(x) < as_oint(y)) return signean_minus; # x<y
  224.               if (as_oint(x) > as_oint(y)) return signean_plus; # x>y
  225.               return signean_null;
  226.             }
  227.             else
  228.             # y>=0, x<0
  229.             { return signean_minus; } # x<y
  230.         }
  231.         else
  232.         { if (!R_minusp(x))
  233.             # y<0, x>=0
  234.             { return signean_plus; } # x>y
  235.             else
  236.             # y<0, x<0
  237.             { if (as_oint(x) > as_oint(y)) return signean_minus; # |x|>|y| -> x<y
  238.               if (as_oint(x) < as_oint(y)) return signean_plus; # |x|<|y| -> x>y
  239.               return signean_null;
  240.             }
  241.         }
  242.     }
  243.  
  244. # Liefert zu zwei Short-Float x und y : (+ x y), ein SF.
  245. # SF_SF_plus_SF(x,y)
  246.   local object SF_SF_plus_SF (object x, object y);
  247. # Methode (nach [Knuth, II, Seminumerical Algorithms, Abschnitt 4.2.1., S.200]):
  248. # x1=0.0 -> Ergebnis x2.
  249. # x2=0.0 -> Ergebnis x1.
  250. # Falls e1<e2, vertausche x1 und x2.
  251. # Also e1 >= e2.
  252. # Falls e1 - e2 >= 16 + 3, Ergebnis x1.
  253. # Schiebe beide Mantissen um 3 Bits nach links (Vorbereitung der Rundung:
  254. #   Bei e1-e2=0,1 ist keine Rundung nötig, bei e1-e2>1 ist der Exponent des
  255. #   Ergebnisses =e1-1, =e1 oder =e1+1. Brauche daher 1 Schutzbit und zwei
  256. #   Rundungsbits: 00 exakt, 01 1.Hälfte, 10 exakte Mitte, 11 2.Hälfte.)
  257. # Schiebe die Mantisse von x2 um e0-e1 Bits nach rechts. (Dabei die Rundung
  258. # ausführen: Bit 0 ist das logische Oder der Bits 0,-1,-2,...)
  259. # Falls x1,x2 selbes Vorzeichen haben: Addiere dieses zur Mantisse von x1.
  260. # Falls x1,x2 verschiedenes Vorzeichen haben: Subtrahiere dieses von der
  261. #   Mantisse von x1. <0 -> (Es war e1=e2) Vertausche die Vorzeichen, negiere.
  262. #                    =0 -> Ergebnis 0.0
  263. # Exponent ist e1.
  264. # Normalisiere, fertig.
  265.   local object SF_SF_plus_SF(x1,x2)
  266.     var reg7 object x1;
  267.     var reg8 object x2;
  268.     { # x1,x2 entpacken:
  269.       var reg9 signean sign1;
  270.       var reg5 sintWL exp1;
  271.       var reg1 uintL mant1;
  272.       var reg9 signean sign2;
  273.       var reg10 sintWL exp2;
  274.       var reg4 uintL mant2;
  275.       SF_decode(x1, { return x2; }, sign1=,exp1=,mant1=);
  276.       SF_decode(x2, { return x1; }, sign2=,exp2=,mant2=);
  277.       if (exp1 < exp2)
  278.         { swap(reg9 object,  x1   ,x2   );
  279.           swap(reg9 signean, sign1,sign2);
  280.           swap(reg9 sintWL,  exp1 ,exp2 );
  281.           swap(reg9 uintL,   mant1,mant2);
  282.         }
  283.       # Nun ist exp1>=exp2.
  284.      {var reg3 uintL expdiff = exp1 - exp2; # Exponentendifferenz
  285.       if (expdiff >= SF_mant_len+3) # >= 16+3 ?
  286.         { return x1; }
  287.       mant1 = mant1 << 3; mant2 = mant2 << 3;
  288.       # Nun 2^(SF_mant_len+3) <= mant1,mant2 < 2^(SF_mant_len+4).
  289.       {var reg2 uintL mant2_last = mant2 & (bit(expdiff)-1); # letzte expdiff Bits von mant2
  290.        mant2 = mant2 >> expdiff; if (!(mant2_last==0)) { mant2 |= bit(0); }
  291.       }
  292.       # mant2 = um expdiff Bits nach rechts geschobene und gerundete Mantisse
  293.       # von x2.
  294.       if (!same_sign_p(x1,x2))
  295.         # verschiedene Vorzeichen -> Mantissen subtrahieren
  296.         { if (mant1 > mant2) { mant1 = mant1 - mant2; goto norm_2; }
  297.           if (mant1 == mant2) # Ergebnis 0 ?
  298.             { return SF_0; }
  299.           # negatives Subtraktionsergebnis
  300.           mant1 = mant2 - mant1; sign1 = sign2; goto norm_2;
  301.         }
  302.         else
  303.         # gleiche Vorzeichen -> Mantissen addieren
  304.         { mant1 = mant1 + mant2; }
  305.       # mant1 = Ergebnis-Mantisse >0, sign1 = Ergebnis-Vorzeichen,
  306.       # exp1 = Ergebnis-Exponent.
  307.       # Außerdem: Bei expdiff=0,1 sind die zwei letzten Bits von mant1 Null,
  308.       # bei expdiff>=2 ist mant1 >= 2^(SF_mant_len+2).
  309.       # Stets ist mant1 < 2^(SF_mant_len+5). (Daher werden die 2 Rundungsbits
  310.       # nachher um höchstens eine Position nach links geschoben werden.)
  311.       # [Knuth, S.201, leicht modifiziert:
  312.       #   N1. m>=1 -> goto N4.
  313.       #   N2. [Hier m<1] m>=1/2 -> goto N5.
  314.       #       N3. m:=2*m, e:=e-1, goto N2.
  315.       #   N4. [Hier 1<=m<2] m:=m/2, e:=e+1.
  316.       #   N5. [Hier 1/2<=m<1] Runde m auf 17 Bits hinterm Komma.
  317.       #       Falls hierdurch m=1 geworden, setze m:=m/2, e:=e+1.
  318.       # ]
  319.       # Bei uns ist m=mant1/2^(SF_mant_len+4),
  320.       # ab Schritt N5 ist m=mant1/2^(SF_mant_len+1).
  321.       norm_1: # [Knuth, S.201, Schritt N1]
  322.       if (mant1 >= bit(SF_mant_len+4)) goto norm_4;
  323.       norm_2: # [Knuth, S.201, Schritt N2]
  324.               # Hier ist mant1 < 2^(SF_mant_len+4)
  325.       if (mant1 >= bit(SF_mant_len+3)) goto norm_5;
  326.       # [Knuth, S.201, Schritt N3]
  327.       mant1 = mant1 << 1; exp1 = exp1-1; # Mantisse links schieben
  328.       goto norm_2;
  329.       norm_4: # [Knuth, S.201, Schritt N4]
  330.               # Hier ist 2^(SF_mant_len+4) <= mant1 < 2^(SF_mant_len+5)
  331.       exp1 = exp1+1;
  332.       mant1 = (mant1>>1) | (mant1 & bit(0)); # Mantisse rechts schieben
  333.       norm_5: # [Knuth, S.201, Schritt N5]
  334.               # Hier ist 2^(SF_mant_len+3) <= mant1 < 2^(SF_mant_len+4)
  335.       # Auf SF_mant_len echte Mantissenbits runden, d.h. rechte 3 Bits
  336.       # wegrunden, und dabei mant1 um 3 Bits nach rechts schieben:
  337.       {var reg2 uintL rounding_bits = mant1 & (bit(3)-1);
  338.        mant1 = mant1 >> 3;
  339.        if ( (rounding_bits < bit(2)) # 000,001,010,011 werden abgerundet
  340.             || ( (rounding_bits == bit(2)) # 100 (genau halbzahlig)
  341.                  && ((mant1 & bit(0)) ==0) # -> round-to-even
  342.           )    )
  343.          # abrunden
  344.          {}
  345.          else
  346.          # aufrunden
  347.          { mant1 = mant1+1;
  348.            if (mant1 >= bit(SF_mant_len+1))
  349.              # Bei Überlauf während der Rundung nochmals rechts schieben
  350.              # (Runden ist hier überflüssig):
  351.              { mant1 = mant1>>1; exp1 = exp1+1; } # Mantisse rechts schieben
  352.          }
  353.       }# Runden fertig
  354.       encode_SF(sign1,exp1,mant1, return);
  355.     }}
  356.  
  357. # Liefert zu zwei Short-Float x und y : (- x y), ein SF.
  358. # SF_SF_minus_SF(x,y)
  359.   local object SF_SF_minus_SF (object x, object y);
  360. # Methode:
  361. # (- x1 x2) = (+ x1 (- x2))
  362.   local object SF_SF_minus_SF(x1,x2)
  363.     var reg2 object x1;
  364.     var reg1 object x2;
  365.     { if (eq(x2,SF_0))
  366.         { return x1; }
  367.         else
  368.         { return SF_SF_plus_SF(x1, as_object(as_oint(x2) ^ wbit(vorz_bit_o)) ); }
  369.     }
  370.  
  371. # Liefert zu zwei Short-Float x und y : (* x y), ein SF.
  372. # SF_SF_mal_SF(x,y)
  373.   local object SF_SF_mal_SF (object x, object y);
  374. # Methode:
  375. # Falls x1=0.0 oder x2=0.0 -> Ergebnis 0.0
  376. # Sonst: Ergebnis-Vorzeichen = VZ von x1 xor VZ von x2.
  377. #        Ergebnis-Exponent = Summe der Exponenten von x1 und x2.
  378. #        Ergebnis-Mantisse = Produkt der Mantissen von x1 und x2, gerundet:
  379. #          2^-17 * (2^16 + m1)  *  2^-17 * (2^16 + m2)
  380. #          = 2^-34 * (2^32 + 2^16*m1 + 2^16*m2 + m1*m2),
  381. #          die Klammer ist >=2^32, <=(2^17-1)^2<2^34 .
  382. #          Falls die Klammer >=2^33 ist, um 17 Bit nach rechts schieben und
  383. #            runden: Falls Bit 16 Null, abrunden; falls Bit 16 Eins und
  384. #            Bits 15..0 alle Null, round-to-even; sonst aufrunden.
  385. #          Falls die Klammer <2^33 ist, um 16 Bit nach rechts schieben und
  386. #            runden: Falls Bit 15 Null, abrunden; falls Bit 15 Eins und
  387. #            Bits 14..0 alle Null, round-to-even; sonst aufrunden. Nach
  388. #            Aufrunden: Falls =2^17, um 1 Bit nach rechts schieben. Sonst
  389. #            Exponenten um 1 erniedrigen.
  390.   local object SF_SF_mal_SF(x1,x2)
  391.     var reg7 object x1;
  392.     var reg8 object x2;
  393.     { # x1,x2 entpacken:
  394.       var reg6 signean sign1;
  395.       var reg3 sintWL exp1;
  396.       var reg4 uintL mant1;
  397.       var reg10 signean sign2;
  398.       var reg9 sintWL exp2;
  399.       var reg5 uintL mant2;
  400.       SF_decode(x1, { return x1; }, sign1=,exp1=,mant1=);
  401.       SF_decode(x2, { return x2; }, sign2=,exp2=,mant2=);
  402.       exp1 = exp1 + exp2; # Summe der Exponenten
  403.       sign1 = sign1 ^ sign2; # Ergebnis-Vorzeichen
  404.      {var reg1 uintL manthi;
  405.       var reg2 uintL mantlo;
  406.       # Mantissen mant1 und mant2 multiplizieren:
  407.       #if (SF_mant_len<16)
  408.       mantlo = mulu16(mant1,mant2);
  409.       manthi = mantlo >> SF_mant_len;
  410.       mantlo = mantlo & (bit(SF_mant_len)-1);
  411.       #elif (SF_mant_len==16)
  412.       manthi = mulu16(low16(mant1),low16(mant2));
  413.       mantlo = low16(manthi);
  414.       manthi = (uint32)(high16(manthi)) + (uint32)(low16(mant1)) + mant2;
  415.       #else # (SF_mant_len>16)
  416.       mulu24(mant1,mant2, manthi=,mantlo=);
  417.       manthi = (manthi << (32-SF_mant_len)) | (mantlo >> SF_mant_len);
  418.       mantlo = mantlo & (bit(SF_mant_len)-1);
  419.       #endif
  420.       # Nun ist 2^SF_mant_len * manthi + mantlo = mant1 * mant2.
  421.       if (manthi >= bit(SF_mant_len+1))
  422.         # mant1*mant2 >= 2^(2*SF_mant_len+1)
  423.         { if ( ((manthi & bit(0)) ==0) # Bit SF_mant_len =0 -> abrunden
  424.                || ( (mantlo ==0) # Bit SF_mant_len =1 und Bits SF_mant_len-1..0 >0 -> aufrunden
  425.                     # round-to-even, je nach Bit SF_mant_len+1 :
  426.                     && ((manthi & bit(1)) ==0)
  427.              )    )
  428.             # abrunden
  429.             { manthi = manthi >> 1; goto ab; }
  430.             else
  431.             # aufrunden
  432.             { manthi = manthi >> 1; goto auf; }
  433.         }
  434.         else
  435.         # mant1*mant2 < 2^(2*SF_mant_len+1)
  436.         { exp1 = exp1-1; # Exponenten decrementieren
  437.           if ( ((mantlo & bit(SF_mant_len-1)) ==0) # Bit SF_mant_len-1 =0 -> abrunden
  438.                || ( ((mantlo & (bit(SF_mant_len-1)-1)) ==0) # Bit SF_mant_len-1 =1 und Bits SF_mant_len-2..0 >0 -> aufrunden
  439.                     # round-to-even, je nach Bit SF_mant_len :
  440.                     && ((manthi & bit(0)) ==0)
  441.              )    )
  442.             # abrunden
  443.             goto ab;
  444.             else
  445.             # aufrunden
  446.             goto auf;
  447.         }
  448.       auf:
  449.       manthi = manthi+1;
  450.       # Hier ist 2^SF_mant_len <= manthi <= 2^(SF_mant_len+1)
  451.       if (manthi >= bit(SF_mant_len+1)) # rounding overflow?
  452.         { manthi = manthi>>1; exp1 = exp1+1; } # Shift nach rechts
  453.       ab:
  454.       # Runden fertig, 2^SF_mant_len <= manthi < 2^(SF_mant_len+1)
  455.       encode_SF(sign1,exp1,manthi, return);
  456.     }}
  457.  
  458. # Liefert zu zwei Short-Float x und y : (/ x y), ein SF.
  459. # SF_SF_durch_SF(x,y)
  460.   local object SF_SF_durch_SF (object x, object y);
  461. # Methode:
  462. # x2 = 0.0 -> Error
  463. # x1 = 0.0 -> Ergebnis 0.0
  464. # Sonst:
  465. # Ergebnis-Vorzeichen = xor der beiden Vorzeichen von x1 und x2
  466. # Ergebnis-Exponent = Differenz der beiden Exponenten von x1 und x2
  467. # Ergebnis-Mantisse = Mantisse mant1 / Mantisse mant2, gerundet.
  468. #   mant1/mant2 > 1/2, mant1/mant2 < 2;
  469. #   nach Rundung mant1/mant2 >=1/2, <=2*mant1<2.
  470. #   Bei mant1/mant2 >=1 brauche 16 Nachkommabits,
  471. #   bei mant1/mant2 <1 brauche 17 Nachkommabits.
  472. #   Fürs Runden: brauche ein Rundungsbit (Rest gibt an, ob exakt).
  473. #   Brauche daher insgesamt 18 Nachkommabits von mant1/mant2.
  474. #   Dividiere daher (als Unsigned Integers) 2^18*(2^17*mant1) durch (2^17*mant2).
  475. #   Falls der Quotient >=2^18 ist, runde die letzten zwei Bits weg und
  476. #     erhöhe den Exponenten um 1.
  477. #   Falls der Quotient <2^18 ist, runde das letzte Bit weg. Bei rounding
  478. #     overflow schiebe um ein weiteres Bit nach rechts, incr. Exponenten.
  479.   local object SF_SF_durch_SF(x1,x2)
  480.     var reg8 object x1;
  481.     var reg9 object x2;
  482.     { # x1,x2 entpacken:
  483.       var reg7 signean sign1;
  484.       var reg3 sintWL exp1;
  485.       var reg5 uintL mant1;
  486.       var reg10 signean sign2;
  487.       var reg10 sintWL exp2;
  488.       var reg6 uintL mant2;
  489.       SF_decode(x2, { divide_0(); }, sign2=,exp2=,mant2=);
  490.       SF_decode(x1, { return x1; }, sign1=,exp1=,mant1=);
  491.       exp1 = exp1 - exp2; # Differenz der Exponenten
  492.       sign1 = sign1 ^ sign2; # Ergebnis-Vorzeichen
  493.       # Dividiere 2^18*mant1 durch mant2 oder (äquivalent)
  494.       # 2^i*2^18*mant1 durch 2^i*mant2 für irgendein i mit 0 <= i <= 32-17 :
  495.      {var reg1 uintL mant;
  496.       var reg4 uintL rest;
  497.       # wähle i = 32-(SF_mant_len+1), also i+(SF_mant_len+2) = 33.
  498.       divu_6432_3232(mant1<<1,0, mant2<<(32-(SF_mant_len+1)), mant=,rest=);
  499.       if (mant >= bit(SF_mant_len+2))
  500.         # Quotient >=2^18 -> 2 Bits wegrunden
  501.         { var reg2 uintL rounding_bits = mant & (bit(2)-1);
  502.           exp1 += 1; # Exponenten incrementieren
  503.           mant = mant >> 2;
  504.           if ( (rounding_bits < bit(1)) # 00,01 werden abgerundet
  505.                || ( (rounding_bits == bit(1)) # 10
  506.                     && (rest == 0) # und genau halbzahlig
  507.                     && ((mant & bit(0)) ==0) # -> round-to-even
  508.              )    )
  509.             # abrunden
  510.             {}
  511.             else
  512.             # aufrunden
  513.             { mant += 1; }
  514.         }
  515.         else
  516.         # Quotient <2^18 -> 1 Bit wegrunden
  517.         { var reg2 uintL rounding_bit = mant & bit(0);
  518.           mant = mant >> 1;
  519.           if ( (rounding_bit == 0) # 0 wird abgerundet
  520.                || ( (rest == 0) # genau halbzahlig
  521.                     && ((mant & bit(0)) ==0) # -> round-to-even
  522.              )    )
  523.             # abrunden
  524.             {}
  525.             else
  526.             # aufrunden
  527.             { mant += 1;
  528.               if (mant >= bit(SF_mant_len+1)) # rounding overflow?
  529.                 { mant = mant>>1; exp1 = exp1+1; }
  530.         }   }
  531.       encode_SF(sign1,exp1,mant, return);
  532.     }}
  533.  
  534. # Liefert zu einem Short-Float x>=0 : (sqrt x), ein SF.
  535. # SF_sqrt_SF(x)
  536.   local object SF_sqrt_SF (object x);
  537. # Methode:
  538. # x = 0.0 -> Ergebnis 0.0
  539. # Ergebnis-Vorzeichen := positiv,
  540. # Ergebnis-Exponent := ceiling(e/2),
  541. # Ergebnis-Mantisse:
  542. #   Bilde aus [1,m15,...,m0,(19 Nullbits)] bei geradem e,
  543. #         aus [0,1,m15,...,m0,(18 Nullbits)] bei ungeradem e
  544. #   die Ganzzahl-Wurzel, eine 18-Bit-Zahl mit einer führenden 1.
  545. #   Runde das letzte Bit weg:
  546. #     Bit 0 = 0 -> abrunden,
  547. #     Bit 0 = 1 und Wurzel exakt -> round-to-even,
  548. #     Bit 0 = 1 und Rest >0 -> aufrunden.
  549. #   Dabei um ein Bit nach rechts schieben.
  550. #   Bei Aufrundung auf 2^17 (rounding overflow) Mantisse um 1 Bit nach rechts
  551. #     schieben und Exponent incrementieren.
  552.   local object SF_sqrt_SF(x)
  553.     var reg3 object x;
  554.     { # x entpacken:
  555.       var reg2 sintWL exp;
  556.       var reg1 uint32 mant;
  557.       SF_decode(x, { return x; },_EMA_,exp=,mant=);
  558.       # Um die 64-Bit-Ganzzahl-Wurzel ausnutzen zu können, fügen wir beim
  559.       # Radikanden 46 bzw. 47 statt 18 bzw. 19 Nullbits an.
  560.       if (exp & bit(0))
  561.         # e ungerade
  562.         { mant = mant << (31-(SF_mant_len+1)); exp = exp+1; }
  563.         else
  564.         # e gerade
  565.         { mant = mant << (32-(SF_mant_len+1)); }
  566.       exp = exp >> 1; # exp := exp/2
  567.      {var reg4 boolean exactp;
  568.       isqrt_64_32(mant,0, mant=,exactp=); # mant := isqrt(mant*2^32), eine 32-Bit-Zahl
  569.       # Die hinteren 31-SF_mant_len Bits wegrunden:
  570.       if ( ((mant & bit(30-SF_mant_len)) ==0) # Bit 14 =0 -> abrunden
  571.            || ( ((mant & (bit(30-SF_mant_len)-1)) ==0) # Bit 14 =1 und Bits 13..0 >0 -> aufrunden
  572.                 && exactp                   # Bit 14 =1 und Bits 13..0 =0, aber Rest -> aufrunden
  573.                 # round-to-even, je nach Bit 15 :
  574.                 && ((mant & bit(31-SF_mant_len)) ==0)
  575.          )    )
  576.         # abrunden
  577.         { mant = mant >> (31-SF_mant_len); }
  578.         else
  579.         # aufrunden
  580.         { mant = mant >> (31-SF_mant_len);
  581.           mant += 1;
  582.           if (mant >= bit(SF_mant_len+1)) # rounding overflow?
  583.             { mant = mant>>1; exp = exp+1; }
  584.         }
  585.       encode_SF(0,exp,mant, return);
  586.     }}
  587.  
  588. # SF_to_I(x) wandelt ein Short-Float x, das eine ganze Zahl darstellt,
  589. # in ein Integer um.
  590. # kann GC auslösen
  591.   local object SF_to_I (object x);
  592. # Methode:
  593. # Falls x=0.0, Ergebnis 0.
  594. # Sonst (ASH Vorzeichen*Mantisse (e-17)).
  595.   local object SF_to_I(x)
  596.     var reg4 object x;
  597.     { # x entpacken:
  598.       var reg3 signean sign;
  599.       var reg2 sintWL exp;
  600.       var reg1 uint32 mant;
  601.       SF_decode(x, { return Fixnum_0; }, sign=,exp=,mant=);
  602.       exp = exp-(SF_mant_len+1);
  603.       return I_I_ash_I( (sign==0
  604.                           ? posfixnum(mant) # mant als Fixnum >0
  605.                           : negfixnum(-(oint)mant) # -mant als Fixnum <0
  606.                         ),
  607.                         L_to_FN(exp)
  608.                       );
  609.     }
  610.  
  611. # I_to_SF(x) wandelt ein Integer x in ein Short-Float um und rundet dabei.
  612. # kann GC auslösen
  613.   local object I_to_SF (object x);
  614. # Methode:
  615. # x=0 -> Ergebnis 0.0
  616. # Merke Vorzeichen von x.
  617. # x:=(abs x)
  618. # Exponent:=(integer-length x)
  619. #   Greife die 18 höchstwertigen Bits heraus (angeführt von einer 1).
  620. #   Runde das letzte Bit weg:
  621. #     Bit 0 = 0 -> abrunden,
  622. #     Bit 0 = 1 und Rest =0 -> round-to-even,
  623. #     Bit 0 = 1 und Rest >0 -> aufrunden.
  624. #   Dabei um ein Bit nach rechts schieben.
  625. #   Bei Aufrundung auf 2^17 (rounding overflow) Mantisse um 1 Bit nach rechts
  626. #     schieben und Exponent incrementieren.
  627.   local object I_to_SF(x)
  628.     var reg7 object x;
  629.     { if (eq(x,Fixnum_0)) { return SF_0; }
  630.      {var reg8 signean sign = R_sign(x); # Vorzeichen
  631.       if (!(sign==0)) { x = I_minus_I(x); } # bei x<0: x := (- x)
  632.       {   var reg9 uintL exp = I_integer_length(x); # (integer-length x)
  633.           # NDS zu x>0 bilden:
  634.        {  var reg2 uintD* MSDptr;
  635.           var reg5 uintC len;
  636.           I_to_NDS_nocopy(x, MSDptr=,len=,_EMA_);
  637.           # MSDptr/len/LSDptr ist die NDS zu x, len>0.
  638.           # Führende Digits holen: Brauche SF_mant_len+1 Bits, dazu intDsize
  639.           # Bits (die NDS kann mit bis zu intDsize Nullbits anfangen).
  640.           # Dann werden diese Bits um (exp mod intDsize) nach rechts geschoben.
  641.         { var reg4 uintD msd = *MSDptr++; # erstes Digit
  642.           var reg1 uint32 msdd = 0; # weitere min(len-1,32/intDsize) Digits
  643.           #define NEXT_DIGIT(i)  \
  644.             { if (--len == 0) goto ok;                            \
  645.               msdd |= (uint32)(*MSDptr++) << (32-(i+1)*intDsize); \
  646.             }
  647.           DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
  648.           #undef NEXT_DIGIT
  649.           --len; ok:
  650.           # Die NDS besteht aus msd, msdd, und len weiteren Digits.
  651.           # Das höchste in 2^32*msd+msdd gesetzte Bit ist Bit Nummer
  652.           # 31 + (exp mod intDsize).
  653.          {var reg6 uintL shiftcount = exp % intDsize;
  654.           var reg3 uint32 mant = # führende 32 Bits
  655.             (shiftcount==0
  656.              ? msdd
  657.              : (((uint32)msd << (32-shiftcount)) | (msdd >> shiftcount))
  658.             );
  659.           # Das höchste in mant gesetzte Bit ist Bit Nummer 31.
  660.           if ( ((mant & bit(30-SF_mant_len)) ==0) # Bit 14 =0 -> abrunden
  661.                || ( ((mant & (bit(30-SF_mant_len)-1)) ==0) # Bit 14 =1 und Bits 13..0 =0
  662.                     && ((msdd & (bit(shiftcount)-1)) ==0) # und weitere Bits aus msdd =0
  663.                     && (!test_loop_up(MSDptr,len)) # und alle weiteren Digits =0
  664.                     # round-to-even, je nach Bit 15 :
  665.                     && ((mant & bit(31-SF_mant_len)) ==0)
  666.              )    )
  667.             # abrunden
  668.             { mant = mant >> (31-SF_mant_len); }
  669.             else
  670.             # aufrunden
  671.             { mant = mant >> (31-SF_mant_len);
  672.               mant += 1;
  673.               if (mant >= bit(SF_mant_len+1)) # rounding overflow?
  674.                 { mant = mant>>1; exp = exp+1; }
  675.             }
  676.           encode_SF(sign,(sintL)exp,mant, return);
  677.     }}}}}}
  678.  
  679. # RA_to_SF(x) wandelt eine rationale Zahl x in ein Short-Float um
  680. # und rundet dabei.
  681. # kann GC auslösen
  682.   local object RA_to_SF (object x);
  683. # Methode:
  684. # x ganz -> klar.
  685. # x = +/- a/b mit Integers a,b>0:
  686. #   Seien n,m so gewählt, daß
  687. #     2^(n-1) <= a < 2^n, 2^(m-1) <= b < 2^m.
  688. #   Dann ist 2^(n-m-1) < a/b < 2^(n-m+1).
  689. #   Berechne n=(integer-length a) und m=(integer-length b) und
  690. #   floor(2^(-n+m+18)*a/b) :
  691. #   Bei n-m>=18 dividiere a durch (ash b (n-m-18)),
  692. #   bei n-m<18 dividiere (ash a (-n+m+18)) durch b.
  693. #   Der erste Wert ist >=2^17, <2^19.
  694. #   Falls er >=2^18 ist, runde 2 Bits weg,
  695. #   falls er <2^18 ist, runde 1 Bit weg.
  696.   local object RA_to_SF(x)
  697.     var reg3 object x;
  698.     { if (RA_integerp(x)) { return I_to_SF(x); }
  699.       # x Ratio
  700.       pushSTACK(TheRatio(x)->rt_den); # b
  701.       x = TheRatio(x)->rt_num; # +/- a
  702.      {var reg7 signean sign = R_sign(x); # Vorzeichen
  703.       if (!(sign==0)) { x = I_minus_I(x); } # Betrag nehmen, liefert a
  704.       pushSTACK(x);
  705.       # Stackaufbau: b, a.
  706.       {var reg4 sintL lendiff = I_integer_length(x) # (integer-length a)
  707.                                 - I_integer_length(STACK_1); # (integer-length b)
  708.        if (lendiff > SF_exp_high-SF_exp_mid) # Exponent >= n-m > Obergrenze ?
  709.          { fehler_overflow(); } # -> Overflow
  710.        if (lendiff < SF_exp_low-SF_exp_mid-2) # Exponent <= n-m+2 < Untergrenze ?
  711.          { if (underflow_allowed())
  712.              { fehler_underflow(); } # -> Underflow
  713.              else
  714.              { skipSTACK(2); return SF_0; }
  715.          }
  716.        { var reg5 object zaehler;
  717.          var reg6 object nenner;
  718.          if (lendiff >= SF_mant_len+2)
  719.            # n-m-18>=0
  720.            { nenner = I_I_ash_I(STACK_1,fixnum((uint32)(lendiff - (SF_mant_len+2)))); # (ash b n-m-18)
  721.              zaehler = popSTACK(); # a
  722.              skipSTACK(1);
  723.            }
  724.            else
  725.            { zaehler = I_I_ash_I(popSTACK(),fixnum((uint32)((SF_mant_len+2) - lendiff))); # (ash a -n+m+18)
  726.              nenner = popSTACK(); # b
  727.            }
  728.          # Division zaehler/nenner durchführen:
  729.          I_I_divide_I_I(zaehler,nenner);
  730.          # Stackaufbau: q, r.
  731.          # 2^17 <= q < 2^19, also ist q Fixnum.
  732.         {var reg1 uint32 mant = posfixnum_to_L(STACK_1);
  733.          if (mant >= bit(SF_mant_len+2))
  734.            # 2^18 <= q < 2^19, schiebe um 2 Bits nach rechts
  735.            { var reg2 uintL rounding_bits = mant & (bit(2)-1);
  736.              lendiff = lendiff+1; # Exponent := n-m+1
  737.              mant = mant >> 2;
  738.              if ( (rounding_bits < bit(1)) # 00,01 werden abgerundet
  739.                   || ( (rounding_bits == bit(1)) # 10
  740.                        && (eq(STACK_0,Fixnum_0)) # und genau halbzahlig (r=0)
  741.                        && ((mant & bit(0)) ==0) # -> round-to-even
  742.                 )    )
  743.                # abrunden
  744.                goto ab;
  745.                else
  746.                # aufrunden
  747.                goto auf;
  748.            }
  749.            else
  750.            { var reg2 uintL rounding_bit = mant & bit(0);
  751.              mant = mant >> 1;
  752.              if ( (rounding_bit == 0) # 0 wird abgerundet
  753.                   || ( (eq(STACK_0,Fixnum_0)) # genau halbzahlig (r=0)
  754.                        && ((mant & bit(0)) ==0) # -> round-to-even
  755.                 )    )
  756.                # abrunden
  757.                goto ab;
  758.                else
  759.                # aufrunden
  760.                goto auf;
  761.            }
  762.          auf:
  763.          mant += 1;
  764.          if (mant >= bit(SF_mant_len+1)) # rounding overflow?
  765.            { mant = mant>>1; lendiff = lendiff+1; }
  766.          ab:
  767.          skipSTACK(2);
  768.          # Fertig.
  769.          encode_SF(sign,lendiff,mant, return);
  770.     }}}}}
  771.  
  772.