home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / Misc / CLISP-1.LHA / CLISP960530-sr.lha / src / dfloat.d < prev    next >
Encoding:
Text File  |  1996-04-15  |  73.1 KB  |  1,682 lines

  1. # Grundfunktionen für Double-Floats
  2.  
  3. # Entpacken eines Double-Float:
  4. #ifdef intQsize
  5. # DF_decode(obj, zero_statement, sign=,exp=,mant=);
  6. # zerlegt ein Double-Float obj.
  7. # Ist obj=0.0, wird zero_statement ausgeführt.
  8. # Sonst: signean sign = Vorzeichen (0 = +, -1 = -),
  9. #        sintWL exp = Exponent (vorzeichenbehaftet),
  10. #        uintQ mant = Mantisse (>= 2^DF_mant_len, < 2^(DF_mant_len+1))
  11.   #define float_value_semhi  float_value
  12.   #define DF_uexp(x)  (((x) >> DF_mant_len) & (bit(DF_exp_len)-1))
  13.   #define DF_decode(obj, zero_statement, sign_zuweisung,exp_zuweisung,mant_zuweisung)  \
  14.     { var reg1 dfloat _x = TheDfloat(obj)->float_value;                    \
  15.       var reg2 uintWL uexp = DF_uexp(_x);                                  \
  16.       if (uexp==0)                                                         \
  17.         { zero_statement } # e=0 -> Zahl 0.0                               \
  18.         else                                                               \
  19.         { exp_zuweisung (sintWL)(uexp - DF_exp_mid); # Exponent            \
  20.           unused (sign_zuweisung ((sint64)_x >> 63));  # Vorzeichen        \
  21.           mant_zuweisung (bit(DF_mant_len) | (_x & (bit(DF_mant_len)-1))); \
  22.     }   }
  23. #else
  24. # DF_decode2(obj, zero_statement, sign=,exp=,manthi=,mantlo=);
  25. # zerlegt ein Double-Float obj.
  26. # Ist obj=0.0, wird zero_statement ausgeführt.
  27. # Sonst: signean sign = Vorzeichen (0 = +, -1 = -),
  28. #        sintWL exp = Exponent (vorzeichenbehaftet),
  29. #        uintL manthi,mantlo = Mantisse 2^32*manthi+mantlo
  30. #                              (>= 2^DF_mant_len, < 2^(DF_mant_len+1))
  31.   #define float_value_semhi  float_value.semhi
  32.   #define DF_uexp(semhi)  (((semhi) >> (DF_mant_len-32)) & (bit(DF_exp_len)-1))
  33.   #define DF_decode2(obj, zero_statement, sign_zuweisung,exp_zuweisung,manthi_zuweisung,mantlo_zuweisung)  \
  34.     { var reg1 uint32 semhi = TheDfloat(obj)->float_value.semhi;                \
  35.       var reg3 uint32 mlo = TheDfloat(obj)->float_value.mlo;                    \
  36.       var reg2 uintWL uexp = DF_uexp(semhi);                                    \
  37.       if (uexp==0)                                                              \
  38.         { zero_statement } # e=0 -> Zahl 0.0                                    \
  39.         else                                                                    \
  40.         { exp_zuweisung (sintWL)(uexp - DF_exp_mid);             # Exponent     \
  41.           unused (sign_zuweisung sign_of_sint32((sint32)(semhi))); # Vorzeichen \
  42.           manthi_zuweisung (bit(DF_mant_len-32) | (semhi & (bit(DF_mant_len-32)-1))); \
  43.           mantlo_zuweisung mlo;                                                 \
  44.     }   }
  45. #endif
  46.  
  47. # Einpacken eines Double-Float:
  48. #ifdef intQsize
  49. # encode_DF(sign,exp,mant, ergebnis=);
  50. # liefert ein Double-Float.
  51. # > signean sign: Vorzeichen, 0 für +, -1 für negativ.
  52. # > sintWL exp: Exponent
  53. # > uintQ mant: Mantisse, sollte >= 2^DF_mant_len und < 2^(DF_mant_len+1) sein.
  54. # < object ergebnis: ein Double-Float
  55. # Der Exponent wird auf Überlauf/Unterlauf getestet.
  56. # kann GC auslösen
  57.   #define encode_DF(sign,exp,mant, erg_zuweisung)  \
  58.     { if ((exp) < (sintWL)(DF_exp_low-DF_exp_mid))                  \
  59.         { if (underflow_allowed())                                  \
  60.             { fehler_underflow(); }                                 \
  61.             else                                                    \
  62.             { erg_zuweisung DF_0; }                                 \
  63.         }                                                           \
  64.       else                                                          \
  65.       if ((exp) > (sintWL)(DF_exp_high-DF_exp_mid))                 \
  66.         { fehler_overflow(); }                                      \
  67.       else                                                          \
  68.       erg_zuweisung allocate_dfloat                                 \
  69.         (  ((sint64)(sign) & bit(63))                  # Vorzeichen \
  70.          | ((uint64)((exp)+DF_exp_mid) << DF_mant_len) # Exponent   \
  71.          | ((uint64)(mant) & (bit(DF_mant_len)-1))     # Mantisse   \
  72.         );                                                          \
  73.     }
  74. #else
  75. # encode2_DF(sign,exp,manthi,mantlo, ergebnis=);
  76. # liefert ein Double-Float.
  77. # > signean sign: Vorzeichen, 0 für +, -1 für negativ.
  78. # > sintWL exp: Exponent
  79. # > uintL manthi,mantlo: Mantisse 2^32*manthi+mantlo,
  80. #                        sollte >= 2^DF_mant_len und < 2^(DF_mant_len+1) sein.
  81. # < object ergebnis: ein Double-Float
  82. # Der Exponent wird auf Überlauf/Unterlauf getestet.
  83. # kann GC auslösen
  84.   #define encode2_DF(sign,exp,manthi,mantlo, erg_zuweisung)  \
  85.     { if ((exp) < (sintWL)(DF_exp_low-DF_exp_mid))                       \
  86.         { if (underflow_allowed())                                       \
  87.             { fehler_underflow(); }                                      \
  88.             else                                                         \
  89.             { erg_zuweisung DF_0; }                                      \
  90.         }                                                                \
  91.       else                                                               \
  92.       if ((exp) > (sintWL)(DF_exp_high-DF_exp_mid))                      \
  93.         { fehler_overflow(); }                                           \
  94.       else                                                               \
  95.       erg_zuweisung allocate_dfloat                                      \
  96.         (  ((sint32)(sign) & bit(31))                       # Vorzeichen \
  97.          | ((uint32)((exp)+DF_exp_mid) << (DF_mant_len-32)) # Exponent   \
  98.          | ((uint32)(manthi) & (bit(DF_mant_len-32)-1))     # Mantisse   \
  99.          , mantlo                                                        \
  100.         );                                                               \
  101.     }
  102. #endif
  103.  
  104. #ifdef FAST_DOUBLE
  105. # Auspacken eines Double:
  106.   #define DF_to_double(obj)  (TheDfloat(obj)->representation.machine_double)
  107. # Überprüfen und Einpacken eines von den 'double'-Routinen gelieferten
  108. # IEEE-Floats.
  109. # Klassifikation:
  110. #   1 <= e <= 2046 : normalisierte Zahl
  111. #   e=0, m/=0: subnormale Zahl
  112. #   e=0, m=0: vorzeichenbehaftete 0.0
  113. #   e=2047, m=0: vorzeichenbehaftete Infinity
  114. #   e=2047, m/=0: NaN
  115. # Angabe der möglicherweise auftretenden Sonderfälle:
  116. #   maybe_overflow: Operation läuft über, liefert IEEE-Infinity
  117. #   maybe_subnormal: Ergebnis sehr klein, liefert IEEE-subnormale Zahl
  118. #   maybe_underflow: Ergebnis sehr klein und /=0, liefert IEEE-Null
  119. #   maybe_divide_0: Ergebnis unbestimmt, liefert IEEE-Infinity
  120. #   maybe_nan: Ergebnis unbestimmt, liefert IEEE-NaN
  121. #ifdef intQsize
  122.   #define double_to_DF(expr,ergebnis_zuweisung,maybe_overflow,maybe_subnormal,maybe_underflow,maybe_divide_0,maybe_nan)  \
  123.     { var dfloatjanus _erg; _erg.machine_double = (expr);                \
  124.       if ((_erg.explicit_ & ((uint64)bit(DF_exp_len+DF_mant_len)-bit(DF_mant_len))) == 0) # e=0 ? \
  125.         { if ((maybe_underflow                                           \
  126.                || (maybe_subnormal && !((_erg.explicit_ << 1) == 0))      \
  127.               )                                                          \
  128.               && underflow_allowed()                                     \
  129.              )                                                           \
  130.             { fehler_underflow(); } # subnormal oder noch kleiner -> Underflow \
  131.             else                                                         \
  132.             { ergebnis_zuweisung DF_0; } # +/- 0.0 -> 0.0                \
  133.         }                                                                \
  134.       elif ((maybe_overflow || maybe_divide_0)                           \
  135.             && (((~_erg.explicit_) & ((uint64)bit(DF_exp_len+DF_mant_len)-bit(DF_mant_len))) == 0) # e=2047 ? \
  136.            )                                                             \
  137.         { if (maybe_nan && !((_erg.explicit_<<(64-DF_mant_len)) == 0))    \
  138.             { divide_0(); } # NaN, also Singularität -> "Division durch 0" \
  139.           else # Infinity                                                \
  140.           if (!maybe_overflow || maybe_divide_0)                         \
  141.             { divide_0(); } # Infinity, Division durch 0                 \
  142.             else                                                         \
  143.             { fehler_overflow(); } # Infinity, Overflow                  \
  144.         }                                                                \
  145.       else                                                               \
  146.         { ergebnis_zuweisung allocate_dfloat(_erg.explicit_); }           \
  147.     }
  148. #else
  149.   #define double_to_DF(expr,ergebnis_zuweisung,maybe_overflow,maybe_subnormal,maybe_underflow,maybe_divide_0,maybe_nan)  \
  150.     { var dfloatjanus _erg; _erg.machine_double = (expr);                 \
  151.       if ((_erg.explicit_.semhi & ((uint32)bit(DF_exp_len+DF_mant_len-32)-bit(DF_mant_len-32))) == 0) # e=0 ? \
  152.         { if ((maybe_underflow                                            \
  153.                || (maybe_subnormal                                        \
  154.                    && !(((_erg.explicit_.semhi << 1) == 0) && (_erg.explicit_.mlo == 0)) \
  155.               )   )                                                       \
  156.               && underflow_allowed()                                      \
  157.              )                                                            \
  158.             { fehler_underflow(); } # subnormal oder noch kleiner -> Underflow \
  159.             else                                                          \
  160.             { ergebnis_zuweisung DF_0; } # +/- 0.0 -> 0.0                 \
  161.         }                                                                 \
  162.       elif ((maybe_overflow || maybe_divide_0)                            \
  163.             && (((~_erg.explicit_.semhi) & ((uint32)bit(DF_exp_len+DF_mant_len-32)-bit(DF_mant_len-32))) == 0) # e=2047 ? \
  164.            )                                                              \
  165.         { if (maybe_nan && !(((_erg.explicit_.semhi<<(64-DF_mant_len)) == 0) && (_erg.explicit_.mlo==0))) \
  166.             { divide_0(); } # NaN, also Singularität -> "Division durch 0" \
  167.           else # Infinity                                                 \
  168.           if (!maybe_overflow || maybe_divide_0)                          \
  169.             { divide_0(); } # Infinity, Division durch 0                  \
  170.             else                                                          \
  171.             { fehler_overflow(); } # Infinity, Overflow                   \
  172.         }                                                                 \
  173.       else                                                                \
  174.         { ergebnis_zuweisung allocate_dfloat(_erg.explicit_.semhi,_erg.explicit_.mlo); }  \
  175.     }
  176. #endif
  177. #endif
  178.  
  179. # DF_zerop(x) stellt fest, ob ein Double-Float x = 0.0 ist.
  180.   # define DF_zerop(x)  (DF_uexp(TheDfloat(x)->float_value_semhi) == 0)
  181.   #define DF_zerop(x)  (TheDfloat(x)->float_value_semhi == 0)
  182.  
  183. # Liefert zu einem Double-Float x : (ftruncate x), ein DF.
  184. # DF_ftruncate_DF(x)
  185. # x wird zur 0 hin zur nächsten ganzen Zahl gerundet.
  186. # kann GC auslösen
  187.   local object DF_ftruncate_DF (object x);
  188. # Methode:
  189. # x = 0.0 oder e<=0 -> Ergebnis 0.0
  190. # 1<=e<=52 -> letzte (53-e) Bits der Mantisse auf 0 setzen,
  191. #             Exponent und Vorzeichen beibehalten
  192. # e>=53 -> Ergebnis x
  193. #ifdef intQsize
  194.   local object DF_ftruncate_DF(x)
  195.     var reg3 object x;
  196.     { var reg2 dfloat x_ = TheDfloat(x)->float_value;
  197.       var reg1 uintWL uexp = DF_uexp(x_); # e + DF_exp_mid
  198.       if (uexp <= DF_exp_mid) # 0.0 oder e<=0 ?
  199.         { return DF_0; }
  200.         else
  201.         { if (uexp > DF_exp_mid+DF_mant_len) # e > 52 ?
  202.             { return x; }
  203.             else
  204.             # 1<=e<=52
  205.             { return allocate_dfloat
  206.                 ( x_ & # Bitmaske: Bits 52-e..0 gelöscht, alle anderen gesetzt
  207.                   ~(bit(DF_mant_len+1+DF_exp_mid-uexp)-1)
  208.                 );
  209.     }   }   }
  210. #else
  211.   local object DF_ftruncate_DF(x)
  212.     var reg3 object x;
  213.     { var reg2 uint32 semhi = TheDfloat(x)->float_value.semhi;
  214.       var reg4 uint32 mlo = TheDfloat(x)->float_value.mlo;
  215.       var reg1 uintWL uexp = DF_uexp(semhi); # e + DF_exp_mid
  216.       if (uexp <= DF_exp_mid) # 0.0 oder e<=0 ?
  217.         { return DF_0; }
  218.         else
  219.         { if (uexp > DF_exp_mid+DF_mant_len) # e > 52 ?
  220.             { return x; }
  221.             else
  222.             # 1<=e<=52
  223.             if (uexp > DF_exp_mid+DF_mant_len+1-32) # e > 21 ?
  224.               { return allocate_dfloat
  225.                   ( semhi,
  226.                     mlo & # Bitmaske: Bits 52-e..0 gelöscht, alle anderen gesetzt
  227.                     ~(bit(DF_mant_len+1+DF_exp_mid-uexp)-1)
  228.                   );
  229.               }
  230.               else
  231.               { return allocate_dfloat
  232.                   ( semhi & # Bitmaske: Bits 20-e..0 gelöscht, alle anderen gesetzt
  233.                     ~(bit(DF_mant_len+1+DF_exp_mid-32-uexp)-1),
  234.                     0
  235.                   );
  236.     }   }     }
  237. #endif
  238.  
  239. # Liefert zu einem Double-Float x : (futruncate x), ein DF.
  240. # DF_futruncate_DF(x)
  241. # x wird von der 0 weg zur nächsten ganzen Zahl gerundet.
  242. # kann GC auslösen
  243.   local object DF_futruncate_DF (object x);
  244. # Methode:
  245. # x = 0.0 -> Ergebnis 0.0
  246. # e<=0 -> Ergebnis 1.0 oder -1.0, je nach Vorzeichen von x.
  247. # 1<=e<=52 -> Greife die letzten (53-e) Bits von x heraus.
  248. #             Sind sie alle =0 -> Ergebnis x.
  249. #             Sonst setze sie alle und erhöhe dann die letzte Stelle um 1.
  250. #             Kein Überlauf der 52 Bit -> fertig.
  251. #             Sonst (Ergebnis eine Zweierpotenz): Mantisse := .1000...000,
  252. #               e:=e+1. (Test auf Überlauf wegen e<=53 überflüssig)
  253. # e>=53 -> Ergebnis x.
  254. #ifdef intQsize
  255.   local object DF_futruncate_DF(x)
  256.     var reg3 object x;
  257.     { var reg2 dfloat x_ = TheDfloat(x)->float_value;
  258.       var reg1 uintWL uexp = DF_uexp(x_); # e + DF_exp_mid
  259.       if (uexp==0) # 0.0 ?
  260.         { return x; }
  261.       if (uexp <= DF_exp_mid) # e<=0 ?
  262.         { # Exponent auf 1, Mantisse auf .1000...000 setzen.
  263.           return ((x_ & bit(63))==0 ? DF_1 : DF_minus1);
  264.         }
  265.         else
  266.         { if (uexp > DF_exp_mid+DF_mant_len) # e > 52 ?
  267.             { return x; }
  268.             else
  269.             { var reg1 uint64 mask = # Bitmaske: Bits 52-e..0 gesetzt, alle anderen gelöscht
  270.                 bit(DF_mant_len+1+DF_exp_mid-uexp)-1;
  271.               if ((x_ & mask)==0) # alle diese Bits =0 ?
  272.                 { return x; }
  273.               return allocate_dfloat
  274.                 ((x_ | mask) # alle diese Bits setzen
  275.                  + 1 # letzte Stelle erhöhen, dabei evtl. Exponenten incrementieren
  276.                 );
  277.     }   }   }
  278. #else
  279.   local object DF_futruncate_DF(x)
  280.     var reg3 object x;
  281.     { var reg2 uint32 semhi = TheDfloat(x)->float_value.semhi;
  282.       var reg4 uint32 mlo = TheDfloat(x)->float_value.mlo;
  283.       var reg1 uintWL uexp = DF_uexp(semhi); # e + DF_exp_mid
  284.       if (uexp==0) # 0.0 ?
  285.         { return x; }
  286.       if (uexp <= DF_exp_mid) # e<=0 ?
  287.         { # Exponent auf 1, Mantisse auf .1000...000 setzen.
  288.           return ((semhi & bit(31))==0 ? DF_1 : DF_minus1);
  289.         }
  290.         else
  291.         { if (uexp > DF_exp_mid+DF_mant_len) # e > 52 ?
  292.             { return x; }
  293.             else
  294.             if (uexp > DF_exp_mid+DF_mant_len+1-32) # e > 21 ?
  295.               { var reg1 uint32 mask = # Bitmaske: Bits 52-e..0 gesetzt, alle anderen gelöscht
  296.                   bit(DF_mant_len+1+DF_exp_mid-uexp)-1;
  297.                 if ((mlo & mask)==0) # alle diese Bits =0 ?
  298.                   { return x; }
  299.                 mlo = (mlo | mask) # alle diese Bits setzen
  300.                       + 1; # letzte Stelle erhöhen,
  301.                 if (mlo==0) { semhi += 1; } # dabei evtl. Exponenten incrementieren
  302.                 return allocate_dfloat(semhi,mlo);
  303.               }
  304.               else
  305.               { var reg1 uint32 mask = # Bitmaske: Bits 20-e..0 gesetzt, alle anderen gelöscht
  306.                   bit(DF_mant_len+1+DF_exp_mid-32-uexp)-1;
  307.                 if ((mlo==0) && ((semhi & mask)==0)) # alle diese Bits und mlo =0 ?
  308.                   { return x; }
  309.                 return allocate_dfloat
  310.                   ((semhi | mask) # alle diese Bits setzen
  311.                    + 1, # letzte Stelle erhöhen, dabei evtl. Exponenten incrementieren
  312.                    0
  313.                   );
  314.     }   }     }
  315. #endif
  316.  
  317. # Liefert zu einem Double-Float x : (fround x), ein DF.
  318. # DF_fround_DF(x)
  319. # x wird zur nächsten ganzen Zahl gerundet.
  320. # kann GC auslösen
  321.   local object DF_fround_DF (object x);
  322. # Methode:
  323. # x = 0.0 oder e<0 -> Ergebnis 0.0
  324. # 0<=e<=52 -> letzte (53-e) Bits der Mantisse wegrunden,
  325. #             Exponent und Vorzeichen beibehalten.
  326. # e>52 -> Ergebnis x
  327. #ifdef intQsize
  328.   local object DF_fround_DF(x)
  329.     var reg3 object x;
  330.     { var reg2 dfloat x_ = TheDfloat(x)->float_value;
  331.       var reg1 uintWL uexp = DF_uexp(x_); # e + DF_exp_mid
  332.       if (uexp < DF_exp_mid) # x = 0.0 oder e<0 ?
  333.         { return DF_0; }
  334.         else
  335.         { if (uexp > DF_exp_mid+DF_mant_len) # e > 52 ?
  336.             { return x; }
  337.             else
  338.             if (uexp > DF_exp_mid+1) # e>1 ?
  339.               { var reg4 uint64 bitmask = # Bitmaske: Bit 52-e gesetzt, alle anderen gelöscht
  340.                   bit(DF_mant_len+DF_exp_mid-uexp);
  341.                 var reg3 uint64 mask = # Bitmaske: Bits 51-e..0 gesetzt, alle anderen gelöscht
  342.                   bitmask-1;
  343.                 if ( ((x_ & bitmask) ==0) # Bit 52-e =0 -> abrunden
  344.                      || ( ((x_ & mask) ==0) # Bit 52-e =1 und Bits 51-e..0 >0 -> aufrunden
  345.                           # round-to-even, je nach Bit 53-e :
  346.                           && ((x_ & (bitmask<<1)) ==0)
  347.                    )    )
  348.                   # abrunden
  349.                   { mask |= bitmask; # Bitmaske: Bits 52-e..0 gesetzt, alle anderen gelöscht
  350.                     return allocate_dfloat( x_ & ~mask );
  351.                   }
  352.                   else
  353.                   # aufrunden
  354.                   { return allocate_dfloat
  355.                       ((x_ | mask) # alle diese Bits 51-e..0 setzen (Bit 52-e schon gesetzt)
  356.                        + 1 # letzte Stelle erhöhen, dabei evtl. Exponenten incrementieren
  357.                       );
  358.                   }
  359.               }
  360.             elif (uexp == DF_exp_mid+1) # e=1 ?
  361.               # Wie bei 1 < e <= 52, nur daß Bit 53-e stets gesetzt ist.
  362.               { if ((x_ & bit(DF_mant_len-1)) ==0) # Bit 52-e =0 -> abrunden
  363.                   # abrunden
  364.                   { return allocate_dfloat( x_ & ~(bit(DF_mant_len)-1) ); }
  365.                   else
  366.                   # aufrunden
  367.                   { return allocate_dfloat
  368.                       ((x_ | (bit(DF_mant_len)-1)) # alle diese Bits 52-e..0 setzen
  369.                        + 1 # letzte Stelle erhöhen, dabei evtl. Exponenten incrementieren
  370.                       );
  371.                   }
  372.               }
  373.             else # e=0 ?
  374.               # Wie bei 1 < e <= 52, nur daß Bit 52-e stets gesetzt
  375.               # und Bit 53-e stets gelöscht ist.
  376.               { if ((x_ & (bit(DF_mant_len)-1)) ==0)
  377.                   # abrunden von +-0.5 zu 0.0
  378.                   { return DF_0; }
  379.                   else
  380.                   # aufrunden
  381.                   { return allocate_dfloat
  382.                       ((x_ | (bit(DF_mant_len)-1)) # alle Bits 51-e..0 setzen
  383.                        + 1 # letzte Stelle erhöhen, dabei Exponenten incrementieren
  384.                       );
  385.               }   }
  386.     }   }
  387. #else
  388.   local object DF_fround_DF(x)
  389.     var reg3 object x;
  390.     { var reg2 uint32 semhi = TheDfloat(x)->float_value.semhi;
  391.       var reg4 uint32 mlo = TheDfloat(x)->float_value.mlo;
  392.       var reg1 uintWL uexp = DF_uexp(semhi); # e + DF_exp_mid
  393.       if (uexp < DF_exp_mid) # x = 0.0 oder e<0 ?
  394.         { return DF_0; }
  395.         else
  396.         { if (uexp > DF_exp_mid+DF_mant_len) # e > 52 ?
  397.             { return x; }
  398.             else
  399.             if (uexp > DF_exp_mid+1) # e>1 ?
  400.               { if (uexp > DF_exp_mid+DF_mant_len-32) # e > 20 ?
  401.                   { var reg4 uint32 bitmask = # Bitmaske: Bit 52-e gesetzt, alle anderen gelöscht
  402.                       bit(DF_mant_len+DF_exp_mid-uexp);
  403.                     var reg3 uint32 mask = # Bitmaske: Bits 51-e..0 gesetzt, alle anderen gelöscht
  404.                       bitmask-1;
  405.                     if ( ((mlo & bitmask) ==0) # Bit 52-e =0 -> abrunden
  406.                          || ( ((mlo & mask) ==0) # Bit 52-e =1 und Bits 51-e..0 >0 -> aufrunden
  407.                               # round-to-even, je nach Bit 53-e :
  408.                               && ( ((bitmask<<1) == 0) # e=21 ?
  409.                                     ? ((semhi & bit(0)) ==0)
  410.                                     : ((mlo & (bitmask<<1)) ==0)
  411.                        )    )    )
  412.                       # abrunden
  413.                       { mask |= bitmask; # Bitmaske: Bits 52-e..0 gesetzt, alle anderen gelöscht
  414.                         return allocate_dfloat(semhi, mlo & ~mask );
  415.                       }
  416.                       else
  417.                       # aufrunden
  418.                       { mlo = (mlo | mask) # alle diese Bits 51-e..0 setzen (Bit 52-e schon gesetzt)
  419.                               + 1; # letzte Stelle erhöhen,
  420.                         if (mlo==0) { semhi += 1; } # dabei evtl. Exponenten incrementieren
  421.                         return allocate_dfloat(semhi,mlo);
  422.                       }
  423.                   }
  424.                   else
  425.                   { var reg4 uint32 bitmask = # Bitmaske: Bit 20-e gesetzt, alle anderen gelöscht
  426.                       bit(DF_mant_len+DF_exp_mid-32-uexp);
  427.                     var reg3 uint32 mask = # Bitmaske: Bits 19-e..0 gesetzt, alle anderen gelöscht
  428.                       bitmask-1;
  429.                     if ( ((semhi & bitmask) ==0) # Bit 52-e =0 -> abrunden
  430.                          || ( (mlo==0) && ((semhi & mask) ==0) # Bit 52-e =1 und Bits 51-e..0 >0 -> aufrunden
  431.                               # round-to-even, je nach Bit 53-e :
  432.                               && ((semhi & (bitmask<<1)) ==0)
  433.                        )    )
  434.                       # abrunden
  435.                       { mask |= bitmask; # Bitmaske: Bits 20-e..0 gesetzt, alle anderen gelöscht
  436.                         return allocate_dfloat( semhi & ~mask, 0 );
  437.                       }
  438.                       else
  439.                       # aufrunden
  440.                       { return allocate_dfloat
  441.                           ((semhi | mask) # alle diese Bits 19-e..0 setzen (Bit 20-e schon gesetzt)
  442.                            + 1, # letzte Stelle erhöhen, dabei evtl. Exponenten incrementieren
  443.                            0
  444.                           );
  445.                       }
  446.                   }
  447.               }
  448.             elif (uexp == DF_exp_mid+1) # e=1 ?
  449.               # Wie bei 1 < e <= 20, nur daß Bit 53-e stets gesetzt ist.
  450.               { if ((semhi & bit(DF_mant_len-32-1)) ==0) # Bit 52-e =0 -> abrunden
  451.                   # abrunden
  452.                   { return allocate_dfloat( semhi & ~(bit(DF_mant_len-32)-1) , 0 ); }
  453.                   else
  454.                   # aufrunden
  455.                   { return allocate_dfloat
  456.                       ((semhi | (bit(DF_mant_len-32)-1)) # alle diese Bits 52-e..0 setzen
  457.                        + 1, # letzte Stelle erhöhen, dabei evtl. Exponenten incrementieren
  458.                        0
  459.                       );
  460.                   }
  461.               }
  462.             else # e=0 ?
  463.               # Wie bei 1 < e <= 20, nur daß Bit 52-e stets gesetzt
  464.               # und Bit 53-e stets gelöscht ist.
  465.               { if ((mlo==0) && ((semhi & (bit(DF_mant_len-32)-1)) ==0))
  466.                   # abrunden von +-0.5 zu 0.0
  467.                   { return DF_0; }
  468.                   else
  469.                   # aufrunden
  470.                   { return allocate_dfloat
  471.                       ((semhi | (bit(DF_mant_len-32)-1)) # alle Bits 51-e..0 setzen
  472.                        + 1, # letzte Stelle erhöhen, dabei Exponenten incrementieren
  473.                        0
  474.                       );
  475.               }   }
  476.     }   }
  477. #endif
  478.  
  479. # Liefert zu einem Double-Float x : (- x), ein DF.
  480. # DF_minus_DF(x)
  481. # kann GC auslösen
  482.   local object DF_minus_DF (object x);
  483. # Methode:
  484. # Falls x=0.0, fertig. Sonst Vorzeichenbit umdrehen.
  485. #ifdef intQsize
  486.   local object DF_minus_DF(x)
  487.     var reg2 object x;
  488.     { var reg1 dfloat x_ = TheDfloat(x)->float_value;
  489.       return (DF_uexp(x_) == 0
  490.               ? x
  491.               : allocate_dfloat( x_ ^ bit(63) )
  492.              );
  493.     }
  494. #else
  495.   local object DF_minus_DF(x)
  496.     var reg2 object x;
  497.     { var reg1 uint32 semhi = TheDfloat(x)->float_value.semhi;
  498.       var reg3 uint32 mlo = TheDfloat(x)->float_value.mlo;
  499.       return (DF_uexp(semhi) == 0
  500.               ? x
  501.               : allocate_dfloat( semhi ^ bit(31), mlo )
  502.              );
  503.     }
  504. #endif
  505.  
  506. # DF_DF_comp(x,y) vergleicht zwei Double-Floats x und y.
  507. # Ergebnis: 0 falls x=y, +1 falls x>y, -1 falls x<y.
  508.   local signean DF_DF_comp (object x, object y);
  509. # Methode:
  510. # x und y haben verschiedenes Vorzeichen ->
  511. #    x < 0 -> x < y
  512. #    x >= 0 -> x > y
  513. # x und y haben gleiches Vorzeichen ->
  514. #    x >=0 -> vergleiche x und y (die rechten 53 Bits)
  515. #    x <0 -> vergleiche y und x (die rechten 53 Bits)
  516. #ifdef intQsize
  517.   local signean DF_DF_comp(x,y)
  518.     var reg3 object x;
  519.     var reg4 object y;
  520.     { var reg1 dfloat x_ = TheDfloat(x)->float_value;
  521.       var reg2 dfloat y_ = TheDfloat(y)->float_value;
  522.       if ((sint64)y_ >= 0)
  523.         # y>=0
  524.         { if ((sint64)x_ >= 0)
  525.             # y>=0, x>=0
  526.             { if (x_ < y_) return signean_minus; # x<y
  527.               if (x_ > y_) return signean_plus; # x>y
  528.               return signean_null;
  529.             }
  530.             else
  531.             # y>=0, x<0
  532.             { return signean_minus; } # x<y
  533.         }
  534.         else
  535.         { if ((sint64)x_ >= 0)
  536.             # y<0, x>=0
  537.             { return signean_plus; } # x>y
  538.             else
  539.             # y<0, x<0
  540.             { if (x_ > y_) return signean_minus; # |x|>|y| -> x<y
  541.               if (x_ < y_) return signean_plus; # |x|<|y| -> x>y
  542.               return signean_null;
  543.             }
  544.         }
  545.     }
  546. #else
  547.   local signean DF_DF_comp(x,y)
  548.     var reg3 object x;
  549.     var reg4 object y;
  550.     { var reg1 uint32 x_semhi = TheDfloat(x)->float_value.semhi;
  551.       var reg2 uint32 y_semhi = TheDfloat(y)->float_value.semhi;
  552.       var reg5 uint32 x_mlo = TheDfloat(x)->float_value.mlo;
  553.       var reg6 uint32 y_mlo = TheDfloat(y)->float_value.mlo;
  554.       if ((sint32)y_semhi >= 0)
  555.         # y>=0
  556.         { if ((sint32)x_semhi >= 0)
  557.             # y>=0, x>=0
  558.             { if (x_semhi < y_semhi) return signean_minus; # x<y
  559.               if (x_semhi > y_semhi) return signean_plus; # x>y
  560.               if (x_mlo < y_mlo) return signean_minus; # x<y
  561.               if (x_mlo > y_mlo) return signean_plus; # x>y
  562.               return signean_null;
  563.             }
  564.             else
  565.             # y>=0, x<0
  566.             { return signean_minus; } # x<y
  567.         }
  568.         else
  569.         { if ((sint32)x_semhi >= 0)
  570.             # y<0, x>=0
  571.             { return signean_plus; } # x>y
  572.             else
  573.             # y<0, x<0
  574.             { if (x_semhi > y_semhi) return signean_minus; # |x|>|y| -> x<y
  575.               if (x_semhi < y_semhi) return signean_plus; # |x|<|y| -> x>y
  576.               if (x_mlo > y_mlo) return signean_minus; # |x|>|y| -> x<y
  577.               if (x_mlo < y_mlo) return signean_plus; # |x|<|y| -> x>y
  578.               return signean_null;
  579.             }
  580.         }
  581.     }
  582. #endif
  583.  
  584. # Liefert zu zwei Double-Float x und y : (+ x y), ein DF.
  585. # DF_DF_plus_DF(x,y)
  586. # kann GC auslösen
  587.   local object DF_DF_plus_DF (object x, object y);
  588. # Methode (nach [Knuth, II, Seminumerical Algorithms, Abschnitt 4.2.1., S.200]):
  589. # x1=0.0 -> Ergebnis x2.
  590. # x2=0.0 -> Ergebnis x1.
  591. # Falls e1<e2, vertausche x1 und x2.
  592. # Also e1 >= e2.
  593. # Falls e1 - e2 >= 52 + 3, Ergebnis x1.
  594. # Schiebe beide Mantissen um 3 Bits nach links (Vorbereitung der Rundung:
  595. #   Bei e1-e2=0,1 ist keine Rundung nötig, bei e1-e2>1 ist der Exponent des
  596. #   Ergebnisses =e1-1, =e1 oder =e1+1. Brauche daher 1 Schutzbit und zwei
  597. #   Rundungsbits: 00 exakt, 01 1.Hälfte, 10 exakte Mitte, 11 2.Hälfte.)
  598. # Schiebe die Mantisse von x2 um e0-e1 Bits nach rechts. (Dabei die Rundung
  599. # ausführen: Bit 0 ist das logische Oder der Bits 0,-1,-2,...)
  600. # Falls x1,x2 selbes Vorzeichen haben: Addiere dieses zur Mantisse von x1.
  601. # Falls x1,x2 verschiedenes Vorzeichen haben: Subtrahiere dieses von der
  602. #   Mantisse von x1. <0 -> (Es war e1=e2) Vertausche die Vorzeichen, negiere.
  603. #                    =0 -> Ergebnis 0.0
  604. # Exponent ist e1.
  605. # Normalisiere, fertig.
  606. #ifdef FAST_DOUBLE
  607.   local object DF_DF_plus_DF(x1,x2)
  608.     var reg1 object x1;
  609.     var reg2 object x2;
  610.     { double_to_DF(DF_to_double(x1) + DF_to_double(x2), return ,
  611.                    TRUE, TRUE, # Overflow und subnormale Zahl abfangen
  612.                    FALSE, # kein Underflow mit Ergebnis +/- 0.0 möglich
  613.                           # (nach Definition der subnormalen Zahlen)
  614.                    FALSE, FALSE # keine Singularität, kein NaN als Ergebnis möglich
  615.                   );
  616.     }
  617. #else
  618. #ifdef intQsize
  619.   local object DF_DF_plus_DF(x1,x2)
  620.     var reg7 object x1;
  621.     var reg8 object x2;
  622.     { # x1,x2 entpacken:
  623.       var reg9 signean sign1;
  624.       var reg5 sintWL exp1;
  625.       var reg1 uint64 mant1;
  626.       var reg9 signean sign2;
  627.       var reg10 sintWL exp2;
  628.       var reg4 uint64 mant2;
  629.       DF_decode(x1, { return x2; }, sign1=,exp1=,mant1=);
  630.       DF_decode(x2, { return x1; }, sign2=,exp2=,mant2=);
  631.       if (exp1 < exp2)
  632.         { swap(reg9 object,  x1   ,x2   );
  633.           swap(reg9 signean, sign1,sign2);
  634.           swap(reg9 sintWL,  exp1 ,exp2 );
  635.           swap(reg9 uint64,   mant1,mant2);
  636.         }
  637.       # Nun ist exp1>=exp2.
  638.      {var reg3 uintL expdiff = exp1 - exp2; # Exponentendifferenz
  639.       if (expdiff >= DF_mant_len+3) # >= 52+3 ?
  640.         { return x1; }
  641.       mant1 = mant1 << 3; mant2 = mant2 << 3;
  642.       # Nun 2^(DF_mant_len+3) <= mant1,mant2 < 2^(DF_mant_len+4).
  643.       {var reg2 uint64 mant2_last = mant2 & (bit(expdiff)-1); # letzte expdiff Bits von mant2
  644.        mant2 = mant2 >> expdiff; if (!(mant2_last==0)) { mant2 |= bit(0); }
  645.       }
  646.       # mant2 = um expdiff Bits nach rechts geschobene und gerundete Mantisse
  647.       # von x2.
  648.       if (!(sign1==sign2))
  649.         # verschiedene Vorzeichen -> Mantissen subtrahieren
  650.         { if (mant1 > mant2) { mant1 = mant1 - mant2; goto norm_2; }
  651.           if (mant1 == mant2) # Ergebnis 0 ?
  652.             { return DF_0; }
  653.           # negatives Subtraktionsergebnis
  654.           mant1 = mant2 - mant1; sign1 = sign2; goto norm_2;
  655.         }
  656.         else
  657.         # gleiche Vorzeichen -> Mantissen addieren
  658.         { mant1 = mant1 + mant2; }
  659.       # mant1 = Ergebnis-Mantisse >0, sign1 = Ergebnis-Vorzeichen,
  660.       # exp1 = Ergebnis-Exponent.
  661.       # Außerdem: Bei expdiff=0,1 sind die zwei letzten Bits von mant1 Null,
  662.       # bei expdiff>=2 ist mant1 >= 2^(DF_mant_len+2).
  663.       # Stets ist mant1 < 2^(DF_mant_len+5). (Daher werden die 2 Rundungsbits
  664.       # nachher um höchstens eine Position nach links geschoben werden.)
  665.       # [Knuth, S.201, leicht modifiziert:
  666.       #   N1. m>=1 -> goto N4.
  667.       #   N2. [Hier m<1] m>=1/2 -> goto N5.
  668.       #       N3. m:=2*m, e:=e-1, goto N2.
  669.       #   N4. [Hier 1<=m<2] m:=m/2, e:=e+1.
  670.       #   N5. [Hier 1/2<=m<1] Runde m auf 53 Bits hinterm Komma.
  671.       #       Falls hierdurch m=1 geworden, setze m:=m/2, e:=e+1.
  672.       # ]
  673.       # Bei uns ist m=mant1/2^(DF_mant_len+4),
  674.       # ab Schritt N5 ist m=mant1/2^(DF_mant_len+1).
  675.       norm_1: # [Knuth, S.201, Schritt N1]
  676.       if (mant1 >= bit(DF_mant_len+4)) goto norm_4;
  677.       norm_2: # [Knuth, S.201, Schritt N2]
  678.               # Hier ist mant1 < 2^(DF_mant_len+4)
  679.       if (mant1 >= bit(DF_mant_len+3)) goto norm_5;
  680.       # [Knuth, S.201, Schritt N3]
  681.       mant1 = mant1 << 1; exp1 = exp1-1; # Mantisse links schieben
  682.       goto norm_2;
  683.       norm_4: # [Knuth, S.201, Schritt N4]
  684.               # Hier ist 2^(DF_mant_len+4) <= mant1 < 2^(DF_mant_len+5)
  685.       exp1 = exp1+1;
  686.       mant1 = (mant1>>1) | (mant1 & bit(0)); # Mantisse rechts schieben
  687.       norm_5: # [Knuth, S.201, Schritt N5]
  688.               # Hier ist 2^(DF_mant_len+3) <= mant1 < 2^(DF_mant_len+4)
  689.       # Auf DF_mant_len echte Mantissenbits runden, d.h. rechte 3 Bits
  690.       # wegrunden, und dabei mant1 um 3 Bits nach rechts schieben:
  691.       {var reg2 uint64 rounding_bits = mant1 & (bit(3)-1);
  692.        mant1 = mant1 >> 3;
  693.        if ( (rounding_bits < bit(2)) # 000,001,010,011 werden abgerundet
  694.             || ( (rounding_bits == bit(2)) # 100 (genau halbzahlig)
  695.                  && ((mant1 & bit(0)) ==0) # -> round-to-even
  696.           )    )
  697.          # abrunden
  698.          {}
  699.          else
  700.          # aufrunden
  701.          { mant1 = mant1+1;
  702.            if (mant1 >= bit(DF_mant_len+1))
  703.              # Bei Überlauf während der Rundung nochmals rechts schieben
  704.              # (Runden ist hier überflüssig):
  705.              { mant1 = mant1>>1; exp1 = exp1+1; } # Mantisse rechts schieben
  706.          }
  707.       }# Runden fertig
  708.       encode_DF(sign1,exp1,mant1, return);
  709.     }}
  710. #else
  711.   local object DF_DF_plus_DF(x1,x2)
  712.     var reg7 object x1;
  713.     var reg8 object x2;
  714.     { # x1,x2 entpacken:
  715.       var reg9 signean sign1;
  716.       var reg5 sintWL exp1;
  717.       var reg1 uintL manthi1;
  718.       var reg1 uintL mantlo1;
  719.       var reg9 signean sign2;
  720.       var reg10 sintWL exp2;
  721.       var reg4 uintL manthi2;
  722.       var reg4 uintL mantlo2;
  723.       DF_decode2(x1, { return x2; }, sign1=,exp1=,manthi1=,mantlo1=);
  724.       DF_decode2(x2, { return x1; }, sign2=,exp2=,manthi2=,mantlo2=);
  725.       if (exp1 < exp2)
  726.         { swap(reg9 object,  x1   ,x2   );
  727.           swap(reg9 signean, sign1,sign2);
  728.           swap(reg9 sintWL,  exp1 ,exp2 );
  729.           swap(reg9 uintL,   manthi1,manthi2);
  730.           swap(reg9 uintL,   mantlo1,mantlo2);
  731.         }
  732.       # Nun ist exp1>=exp2.
  733.      {var reg3 uintL expdiff = exp1 - exp2; # Exponentendifferenz
  734.       if (expdiff >= DF_mant_len+3) # >= 52+3 ?
  735.         { return x1; }
  736.       manthi1 = (manthi1 << 3) | (mantlo1 >> (32-3)); mantlo1 = mantlo1 << 3;
  737.       manthi2 = (manthi2 << 3) | (mantlo2 >> (32-3)); mantlo2 = mantlo2 << 3;
  738.       # Nun 2^(DF_mant_len+3) <= mant1,mant2 < 2^(DF_mant_len+4).
  739.       if (expdiff<32)
  740.         {if (!(expdiff==0))
  741.            {var reg2 uintL mant2_last = mantlo2 & (bit(expdiff)-1); # letzte expdiff Bits von mant2
  742.             mantlo2 = (mantlo2 >> expdiff) | (manthi2 << (32-expdiff));
  743.             manthi2 = manthi2 >> expdiff;
  744.             if (!(mant2_last==0)) { mantlo2 |= bit(0); }
  745.         }  }
  746.         else
  747.         {var reg2 uintL mant2_last = (manthi2 & (bit(expdiff-32)-1)) | mantlo2; # letzte expdiff Bits von mant2
  748.          mantlo2 = manthi2 >> (expdiff-32); manthi2 = 0;
  749.          if (!(mant2_last==0)) { mantlo2 |= bit(0); }
  750.         }
  751.       # mant2 = um expdiff Bits nach rechts geschobene und gerundete Mantisse
  752.       # von x2.
  753.       if (!(sign1==sign2))
  754.         # verschiedene Vorzeichen -> Mantissen subtrahieren
  755.         { if (manthi1 > manthi2)
  756.             { manthi1 = manthi1 - manthi2;
  757.               if (mantlo1 < mantlo2) { manthi1 -= 1; }
  758.               mantlo1 = mantlo1 - mantlo2;
  759.               goto norm_2;
  760.             }
  761.           if (manthi1 == manthi2)
  762.             { if (mantlo1 > mantlo2)
  763.                 { manthi1 = 0; mantlo1 = mantlo1 - mantlo2; goto norm_2; }
  764.               if (mantlo1 == mantlo2) # Ergebnis 0 ?
  765.                 { return DF_0; }
  766.             }
  767.           # Hier ((manthi1 < manthi2) || ((manthi1 == manthi2) && (mantlo1 < mantlo2))).
  768.           # negatives Subtraktionsergebnis
  769.           manthi1 = manthi2 - manthi1;
  770.           if (mantlo2 < mantlo1) { manthi1 -= 1; }
  771.           mantlo1 = mantlo2 - mantlo1;
  772.           sign1 = sign2;
  773.           goto norm_2;
  774.         }
  775.         else
  776.         # gleiche Vorzeichen -> Mantissen addieren
  777.         { manthi1 = manthi1 + manthi2;
  778.           if ((mantlo1 = mantlo1 + mantlo2) < mantlo2) { manthi1 += 1; }
  779.         }
  780.       # mant1 = Ergebnis-Mantisse >0, sign1 = Ergebnis-Vorzeichen,
  781.       # exp1 = Ergebnis-Exponent.
  782.       # Außerdem: Bei expdiff=0,1 sind die zwei letzten Bits von mant1 Null,
  783.       # bei expdiff>=2 ist mant1 >= 2^(DF_mant_len+2).
  784.       # Stets ist mant1 < 2^(DF_mant_len+5). (Daher werden die 2 Rundungsbits
  785.       # nachher um höchstens eine Position nach links geschoben werden.)
  786.       # [Knuth, S.201, leicht modifiziert:
  787.       #   N1. m>=1 -> goto N4.
  788.       #   N2. [Hier m<1] m>=1/2 -> goto N5.
  789.       #       N3. m:=2*m, e:=e-1, goto N2.
  790.       #   N4. [Hier 1<=m<2] m:=m/2, e:=e+1.
  791.       #   N5. [Hier 1/2<=m<1] Runde m auf 53 Bits hinterm Komma.
  792.       #       Falls hierdurch m=1 geworden, setze m:=m/2, e:=e+1.
  793.       # ]
  794.       # Bei uns ist m=mant1/2^(DF_mant_len+4),
  795.       # ab Schritt N5 ist m=mant1/2^(DF_mant_len+1).
  796.       norm_1: # [Knuth, S.201, Schritt N1]
  797.       if (manthi1 >= bit(DF_mant_len-32+4)) goto norm_4;
  798.       norm_2: # [Knuth, S.201, Schritt N2]
  799.               # Hier ist mant1 < 2^(DF_mant_len+4)
  800.       if (manthi1 >= bit(DF_mant_len-32+3)) goto norm_5;
  801.       # [Knuth, S.201, Schritt N3]
  802.       manthi1 = (manthi1 << 1) | (mantlo1 >> 31); # Mantisse links schieben
  803.       mantlo1 = mantlo1 << 1;
  804.       exp1 = exp1-1;
  805.       goto norm_2;
  806.       norm_4: # [Knuth, S.201, Schritt N4]
  807.               # Hier ist 2^(DF_mant_len+4) <= mant1 < 2^(DF_mant_len+5)
  808.       exp1 = exp1+1;
  809.       mantlo1 = (mantlo1 >> 1) | (manthi1 << 31) | (mantlo1 & bit(0)); # Mantisse rechts schieben
  810.       manthi1 = (manthi1 >> 1);
  811.       norm_5: # [Knuth, S.201, Schritt N5]
  812.               # Hier ist 2^(DF_mant_len+3) <= mant1 < 2^(DF_mant_len+4)
  813.       # Auf DF_mant_len echte Mantissenbits runden, d.h. rechte 3 Bits
  814.       # wegrunden, und dabei mant1 um 3 Bits nach rechts schieben:
  815.       {var reg2 uintL rounding_bits = mantlo1 & (bit(3)-1);
  816.        mantlo1 = (mantlo1 >> 3) | (manthi1 << (32-3)); manthi1 = manthi1 >> 3;
  817.        if ( (rounding_bits < bit(2)) # 000,001,010,011 werden abgerundet
  818.             || ( (rounding_bits == bit(2)) # 100 (genau halbzahlig)
  819.                  && ((mantlo1 & bit(0)) ==0) # -> round-to-even
  820.           )    )
  821.          # abrunden
  822.          {}
  823.          else
  824.          # aufrunden
  825.          { mantlo1 = mantlo1+1;
  826.            if (mantlo1==0)
  827.              { manthi1 = manthi1+1;
  828.                if (manthi1 >= bit(DF_mant_len-32+1))
  829.                  # Bei Überlauf während der Rundung nochmals rechts schieben
  830.                  # (Runden ist hier überflüssig):
  831.                  { manthi1 = manthi1>>1; exp1 = exp1+1; } # Mantisse rechts schieben
  832.          }   }
  833.       }# Runden fertig
  834.       encode2_DF(sign1,exp1,manthi1,mantlo1, return);
  835.     }}
  836. #endif
  837. #endif
  838.  
  839. # Liefert zu zwei Double-Float x und y : (- x y), ein DF.
  840. # DF_DF_minus_DF(x,y)
  841. # kann GC auslösen
  842.   local object DF_DF_minus_DF (object x, object y);
  843. # Methode:
  844. # (- x1 x2) = (+ x1 (- x2))
  845. #ifdef FAST_DOUBLE
  846.   local object DF_DF_minus_DF(x1,x2)
  847.     var reg1 object x1;
  848.     var reg2 object x2;
  849.     { double_to_DF(DF_to_double(x1) - DF_to_double(x2), return ,
  850.                    TRUE, TRUE, # Overflow und subnormale Zahl abfangen
  851.                    FALSE, # kein Underflow mit Ergebnis +/- 0.0 möglich
  852.                           # (nach Definition der subnormalen Zahlen)
  853.                    FALSE, FALSE # keine Singularität, kein NaN als Ergebnis möglich
  854.                   );
  855.     }
  856. #else
  857. #ifdef intQsize
  858.   local object DF_DF_minus_DF(x1,x2)
  859.     var reg3 object x1;
  860.     var reg1 object x2;
  861.     { var reg2 dfloat x2_ = TheDfloat(x2)->float_value;
  862.       if (DF_uexp(x2_) == 0)
  863.         { return x1; }
  864.         else
  865.         { pushSTACK(x1);
  866.           x2 = allocate_dfloat(x2_ ^ bit(63));
  867.           return DF_DF_plus_DF(popSTACK(),x2);
  868.     }   }
  869. #else
  870.   local object DF_DF_minus_DF(x1,x2)
  871.     var reg3 object x1;
  872.     var reg1 object x2;
  873.     { var reg2 uint32 x2_semhi = TheDfloat(x2)->float_value.semhi;
  874.       var reg4 uint32 x2_mlo = TheDfloat(x2)->float_value.mlo;
  875.       if (DF_uexp(x2_semhi) == 0)
  876.         { return x1; }
  877.         else
  878.         { pushSTACK(x1);
  879.           x2 = allocate_dfloat(x2_semhi ^ bit(31), x2_mlo);
  880.           return DF_DF_plus_DF(popSTACK(),x2);
  881.     }   }
  882. #endif
  883. #endif
  884.  
  885. # Liefert zu zwei Double-Float x und y : (* x y), ein DF.
  886. # DF_DF_mal_DF(x,y)
  887. # kann GC auslösen
  888.   local object DF_DF_mal_DF (object x, object y);
  889. # Methode:
  890. # Falls x1=0.0 oder x2=0.0 -> Ergebnis 0.0
  891. # Sonst: Ergebnis-Vorzeichen = VZ von x1 xor VZ von x2.
  892. #        Ergebnis-Exponent = Summe der Exponenten von x1 und x2.
  893. #        Ergebnis-Mantisse = Produkt der Mantissen von x1 und x2, gerundet:
  894. #          2^-53 * mant1  *  2^-53 * mant2  =  2^-106 * (mant1*mant2),
  895. #          die Klammer ist >=2^104, <=(2^53-1)^2<2^106 .
  896. #          Falls die Klammer >=2^105 ist, um 53 Bit nach rechts schieben und
  897. #            runden: Falls Bit 52 Null, abrunden; falls Bit 52 Eins und
  898. #            Bits 51..0 alle Null, round-to-even; sonst aufrunden.
  899. #          Falls die Klammer <2^105 ist, um 52 Bit nach rechts schieben und
  900. #            runden: Falls Bit 51 Null, abrunden; falls Bit 51 Eins und
  901. #            Bits 50..0 alle Null, round-to-even; sonst aufrunden. Nach
  902. #            Aufrunden: Falls =2^53, um 1 Bit nach rechts schieben. Sonst
  903. #            Exponenten um 1 erniedrigen.
  904. #ifdef FAST_DOUBLE
  905.   local object DF_DF_mal_DF(x1,x2)
  906.     var reg1 object x1;
  907.     var reg2 object x2;
  908.     { double_to_DF(DF_to_double(x1) * DF_to_double(x2), return ,
  909.                    TRUE, TRUE, # Overflow und subnormale Zahl abfangen
  910.                    !(DF_zerop(x1) || DF_zerop(x2)), # ein Ergebnis +/- 0.0
  911.                                # ist genau dann in Wirklichkeit ein Underflow
  912.                    FALSE, FALSE # keine Singularität, kein NaN als Ergebnis möglich
  913.                   );
  914.     }
  915. #else
  916.   local object DF_DF_mal_DF(x1,x2)
  917.     var reg7 object x1;
  918.     var reg8 object x2;
  919.     { # x1,x2 entpacken:
  920.       var reg6 signean sign1;
  921.       var reg3 sintWL exp1;
  922.       var reg4 uintL manthi1;
  923.       var reg4 uintL mantlo1;
  924.       var reg10 signean sign2;
  925.       var reg9 sintWL exp2;
  926.       var reg5 uintL manthi2;
  927.       var reg5 uintL mantlo2;
  928.       #ifdef intQsize
  929.       { var reg1 uint64 mant1;
  930.         DF_decode(x1, { return x1; }, sign1=,exp1=,mant1=);
  931.         manthi1 = (uint32)(mant1>>32); mantlo1 = (uint32)mant1;
  932.       }
  933.       { var reg1 uint64 mant2;
  934.         DF_decode(x2, { return x2; }, sign2=,exp2=,mant2=);
  935.         manthi2 = (uint32)(mant2>>32); mantlo2 = (uint32)mant2;
  936.       }
  937.       #else
  938.       DF_decode2(x1, { return x1; }, sign1=,exp1=,manthi1=,mantlo1=);
  939.       DF_decode2(x2, { return x2; }, sign2=,exp2=,manthi2=,mantlo2=);
  940.       #endif
  941.       exp1 = exp1 + exp2; # Summe der Exponenten
  942.       sign1 = sign1 ^ sign2; # Ergebnis-Vorzeichen
  943.      {# Mantissen mant1 und mant2 multiplizieren (64x64-Bit-Multiplikation):
  944.       var uintD mant1 [64/intDsize];
  945.       var uintD mant2 [64/intDsize];
  946.       var uintD mant [128/intDsize];
  947.       #if (intDsize==32) || (intDsize==16) || (intDsize==8)
  948.       set_32_Dptr(mant1,manthi1); set_32_Dptr(&mant1[32/intDsize],mantlo1);
  949.       set_32_Dptr(mant2,manthi2); set_32_Dptr(&mant2[32/intDsize],mantlo2);
  950.       #else
  951.       {var reg1 uintD* ptr;
  952.        ptr = &mant1[64/intDsize];
  953.        doconsttimes(32/intDsize, { *--ptr = (uintD)mantlo1; mantlo1 = mantlo1>>intDsize; } );
  954.        doconsttimes(32/intDsize, { *--ptr = (uintD)manthi1; manthi1 = manthi1>>intDsize; } );
  955.       }
  956.       {var reg1 uintD* ptr;
  957.        ptr = &mant2[64/intDsize];
  958.        doconsttimes(32/intDsize, { *--ptr = (uintD)mantlo2; mantlo2 = mantlo2>>intDsize; } );
  959.        doconsttimes(32/intDsize, { *--ptr = (uintD)manthi2; manthi2 = manthi2>>intDsize; } );
  960.       }
  961.       #endif
  962.       mulu_2loop_down(&mant1[64/intDsize],64/intDsize,
  963.                       &mant2[64/intDsize],64/intDsize,
  964.                       &mant[128/intDsize]
  965.                      );
  966.       {
  967.         #ifdef intQsize
  968.         var reg1 uint64 manterg;
  969.         #else
  970.         var reg1 uintL manthi;
  971.         var reg2 uintL mantlo;
  972.         #endif
  973.         # Produkt mant = mant1 * mant2 ist >= 2^104, < 2^106. Bit 105 abtesten:
  974.         #define mant_bit(k)  (mant[128/intDsize - 1 - floor(k,intDsize)] & bit((k)%intDsize))
  975.         if (mant_bit(2*DF_mant_len+1))
  976.           # mant>=2^(2*DF_mant_len+1), um DF_mant_len+1 Bits nach rechts schieben:
  977.           { # Bits 105..53 holen:
  978.             #if defined(intQsize) # && (intDsize==32)
  979.               manterg = ((uint64)mant[0] << 43) | ((uint64)mant[1] << 11) | ((uint64)mant[2] >> 21); # Bits 116..53
  980.               #define mantrest() ((mant[2] & (bit(21)-1)) || mant[3])
  981.             #elif (intDsize==32)
  982.               manthi = ((uint32)mant[0] << 11) | ((uint32)mant[1] >> 21); # Bits 116..85
  983.               mantlo = ((uint32)mant[1] << 11) | ((uint32)mant[2] >> 21); # Bits 84..53
  984.               #define mantrest() ((mant[2] & (bit(21)-1)) || mant[3])
  985.             #elif (intDsize==16)
  986.               manthi = # ((uint32)mant[0] << 27) | ((uint32)mant[1] << 11) | ((uint32)mant[2] >> 5); # Bits 116..85
  987.                        (highlow32_at(&mant[0])<<11) | ((uint32)mant[2] >> 5); # Bits 116..85
  988.               mantlo = # ((uint32)mant[2] << 27) | ((uint32)mant[3] << 11) | ((uint32)mant[4] >> 5); # Bits 84..53
  989.                        (highlow32_at(&mant[2])<<11) | ((uint32)mant[4] >> 5); # Bits 84..53
  990.               #define mantrest() ((mant[4] & (bit(5)-1)) || mant[5] || mant[6] || mant[7])
  991.             #elif (intDsize==8)
  992.               manthi = ((uint32)mant[1] << 27) | ((uint32)mant[2] << 19) | ((uint32)mant[3] << 11) | ((uint32)mant[4] << 3) | ((uint32)mant[5] >> 5); # Bits 116..85
  993.               mantlo = ((uint32)mant[5] << 27) | ((uint32)mant[6] << 19) | ((uint32)mant[7] << 11) | ((uint32)mant[8] << 3) | ((uint32)mant[9] >> 5); # Bits 84..53
  994.               #define mantrest() ((mant[9] & (bit(5)-1)) || mant[10] || mant[11] || mant[12] || mant[13] || mant[14] || mant[15])
  995.             #endif
  996.             if ( (mant_bit(DF_mant_len) ==0) # Bit DF_mant_len =0 -> abrunden
  997.                  || ( !mantrest() # Bit DF_mant_len =1 und Bits DF_mant_len-1..0 >0 -> aufrunden
  998.                       # round-to-even, je nach Bit DF_mant_len+1 :
  999.                       && (mant_bit(DF_mant_len+1) ==0)
  1000.                )    )
  1001.               # abrunden
  1002.               goto ab;
  1003.               else
  1004.               # aufrunden
  1005.               goto auf;
  1006.             #undef mantrest
  1007.           }
  1008.           else
  1009.           # mant<2^(2*DF_mant_len+1), um DF_mant_len Bits nach rechts schieben:
  1010.           { exp1 = exp1-1; # Exponenten decrementieren
  1011.             # Bits 104..52 holen:
  1012.             #if defined(intQsize) # && (intDsize==32)
  1013.               manterg = ((uint64)mant[0] << 44) | ((uint64)mant[1] << 12) | ((uint64)mant[2] >> 20); # Bits 115..52
  1014.               #define mantrest() ((mant[2] & (bit(20)-1)) || mant[3])
  1015.             #elif (intDsize==32)
  1016.               manthi = ((uint32)mant[0] << 12) | ((uint32)mant[1] >> 20); # Bits 115..84
  1017.               mantlo = ((uint32)mant[1] << 12) | ((uint32)mant[2] >> 20); # Bits 83..52
  1018.               #define mantrest() ((mant[2] & (bit(20)-1)) || mant[3])
  1019.             #elif (intDsize==16)
  1020.               manthi = # ((uint32)mant[0] << 28) | ((uint32)mant[1] << 12) | ((uint32)mant[2] >> 4); # Bits 115..84
  1021.                        (highlow32_at(&mant[0])<<12) | ((uint32)mant[2] >> 4); # Bits 115..84
  1022.               mantlo = # ((uint32)mant[2] << 28) | ((uint32)mant[3] << 12) | ((uint32)mant[4] >> 4); # Bits 83..52
  1023.                        (highlow32_at(&mant[2])<<12) | ((uint32)mant[4] >> 4); # Bits 83..52
  1024.               #define mantrest() ((mant[4] & (bit(4)-1)) || mant[5] || mant[6] || mant[7])
  1025.             #elif (intDsize==8)
  1026.               manthi = ((uint32)mant[1] << 28) | ((uint32)mant[2] << 20) | ((uint32)mant[3] << 12) | ((uint32)mant[4] << 4) | ((uint32)mant[5] >> 4); # Bits 115..84
  1027.               mantlo = ((uint32)mant[5] << 28) | ((uint32)mant[6] << 20) | ((uint32)mant[7] << 12) | ((uint32)mant[8] << 4) | ((uint32)mant[9] >> 4); # Bits 83..52
  1028.               #define mantrest() ((mant[9] & (bit(4)-1)) || mant[10] || mant[11] || mant[12] || mant[13] || mant[14] || mant[15])
  1029.             #endif
  1030.             if ( (mant_bit(DF_mant_len-1) ==0) # Bit DF_mant_len-1 =0 -> abrunden
  1031.                  || ( !mantrest() # Bit DF_mant_len-1 =1 und Bits DF_mant_len-2..0 >0 -> aufrunden
  1032.                       # round-to-even, je nach Bit DF_mant_len :
  1033.                       && (mant_bit(DF_mant_len) ==0)
  1034.                )    )
  1035.               # abrunden
  1036.               goto ab;
  1037.               else
  1038.               # aufrunden
  1039.               goto auf;
  1040.             #undef mantrest
  1041.           }
  1042.         #undef mant_bit
  1043.         auf:
  1044.         #ifdef intQsize
  1045.         manterg = manterg+1;
  1046.         # Hier ist 2^DF_mant_len <= manterg <= 2^(DF_mant_len+1)
  1047.         if (manterg >= bit(DF_mant_len+1)) # rounding overflow?
  1048.           { manterg = manterg>>1; exp1 = exp1+1; } # Shift nach rechts
  1049.         #else
  1050.         mantlo = mantlo+1;
  1051.         if (mantlo==0)
  1052.           { manthi = manthi+1;
  1053.             # Hier ist 2^(DF_mant_len-32) <= manthi <= 2^(DF_mant_len-32+1)
  1054.             if (manthi >= bit(DF_mant_len-32+1)) # rounding overflow?
  1055.               { manthi = manthi>>1; exp1 = exp1+1; } # Shift nach rechts
  1056.           }
  1057.         #endif
  1058.         ab:
  1059.         # Runden fertig, 2^DF_mant_len <= manterg < 2^(DF_mant_len+1)
  1060.         #ifdef intQsize
  1061.         encode_DF(sign1,exp1,manterg, return);
  1062.         #else
  1063.         encode2_DF(sign1,exp1,manthi,mantlo, return);
  1064.         #endif
  1065.     }}}
  1066. #endif
  1067.  
  1068. # Liefert zu zwei Double-Float x und y : (/ x y), ein DF.
  1069. # DF_DF_durch_DF(x,y)
  1070. # kann GC auslösen
  1071.   local object DF_DF_durch_DF (object x, object y);
  1072. # Methode:
  1073. # x2 = 0.0 -> Error
  1074. # x1 = 0.0 -> Ergebnis 0.0
  1075. # Sonst:
  1076. # Ergebnis-Vorzeichen = xor der beiden Vorzeichen von x1 und x2
  1077. # Ergebnis-Exponent = Differenz der beiden Exponenten von x1 und x2
  1078. # Ergebnis-Mantisse = Mantisse mant1 / Mantisse mant2, gerundet.
  1079. #   mant1/mant2 > 1/2, mant1/mant2 < 2;
  1080. #   nach Rundung mant1/mant2 >=1/2, <=2*mant1<2.
  1081. #   Bei mant1/mant2 >=1 brauche 52 Nachkommabits,
  1082. #   bei mant1/mant2 <1 brauche 53 Nachkommabits.
  1083. #   Fürs Runden: brauche ein Rundungsbit (Rest gibt an, ob exakt).
  1084. #   Brauche daher insgesamt 54 Nachkommabits von mant1/mant2.
  1085. #   Dividiere daher (als Unsigned Integers) 2^54*(2^53*mant1) durch (2^53*mant2).
  1086. #   Falls der Quotient >=2^54 ist, runde die letzten zwei Bits weg und
  1087. #     erhöhe den Exponenten um 1.
  1088. #   Falls der Quotient <2^54 ist, runde das letzte Bit weg. Bei rounding
  1089. #     overflow schiebe um ein weiteres Bit nach rechts, incr. Exponenten.
  1090. #if defined(FAST_DOUBLE) && !defined(I80Z86)
  1091.   local object DF_DF_durch_DF(x1,x2)
  1092.     var reg1 object x1;
  1093.     var reg2 object x2;
  1094.     { double_to_DF(DF_to_double(x1) / DF_to_double(x2), return ,
  1095.                    TRUE, TRUE, # Overflow und subnormale Zahl abfangen
  1096.                    !DF_zerop(x1), # ein Ergebnis +/- 0.0
  1097.                                # ist genau dann in Wirklichkeit ein Underflow
  1098.                    DF_zerop(x2), # Division durch Null abfangen
  1099.                    FALSE # kein NaN als Ergebnis möglich
  1100.                   );
  1101.     }
  1102. #else
  1103.   local object DF_DF_durch_DF(x1,x2)
  1104.     var reg8 object x1;
  1105.     var reg9 object x2;
  1106.     { # x1,x2 entpacken:
  1107.       var reg7 signean sign1;
  1108.       var reg3 sintWL exp1;
  1109.       var reg5 uintL manthi1;
  1110.       var reg5 uintL mantlo1;
  1111.       var reg10 signean sign2;
  1112.       var reg10 sintWL exp2;
  1113.       var reg6 uintL manthi2;
  1114.       var reg6 uintL mantlo2;
  1115.       #ifdef intQsize
  1116.       var reg5 uint64 mant1;
  1117.       var reg6 uint64 mant2;
  1118.       DF_decode(x2, { divide_0(); }, sign2=,exp2=,mant2=);
  1119.       DF_decode(x1, { return x1; }, sign1=,exp1=,mant1=);
  1120.       #else
  1121.       DF_decode2(x2, { divide_0(); }, sign2=,exp2=,manthi2=,mantlo2=);
  1122.       DF_decode2(x1, { return x1; }, sign1=,exp1=,manthi1=,mantlo1=);
  1123.       #endif
  1124.       exp1 = exp1 - exp2; # Differenz der Exponenten
  1125.       sign1 = sign1 ^ sign2; # Ergebnis-Vorzeichen
  1126.       # Dividiere 2^54*mant1 durch mant2 oder (äquivalent)
  1127.       # 2^i*2^54*mant1 durch 2^i*mant2 für irgendein i mit 0 <= i <= 64-53 :
  1128.       # wähle i = 64-(DF_mant_len+1), also i+(DF_mant_len+2) = 65.
  1129.       #ifdef intQsize
  1130.       mant1 = mant1 << 1;
  1131.       mant2 = mant2 << (64-(DF_mant_len+1));
  1132.       manthi1 = (uint32)(mant1>>32); mantlo1 = (uint32)mant1;
  1133.       manthi2 = (uint32)(mant2>>32); mantlo2 = (uint32)mant2;
  1134.       #else
  1135.       manthi1 = (manthi1 << 1) | (mantlo1 >> 31); mantlo1 = mantlo1 << 1;
  1136.       manthi2 = (manthi2 << (64-(DF_mant_len+1))) | (mantlo2 >> ((DF_mant_len+1)-32)); mantlo2 = mantlo2 << (64-(DF_mant_len+1));
  1137.       #endif
  1138.      {var uintD mant1 [128/intDsize];
  1139.       var uintD mant2 [64/intDsize];
  1140.       #if (intDsize==32) || (intDsize==16) || (intDsize==8)
  1141.       set_32_Dptr(mant1,manthi1); set_32_Dptr(&mant1[32/intDsize],mantlo1);
  1142.         set_32_Dptr(&mant1[2*32/intDsize],0); set_32_Dptr(&mant1[3*32/intDsize],0);
  1143.       set_32_Dptr(mant2,manthi2); set_32_Dptr(&mant2[32/intDsize],mantlo2);
  1144.       #else
  1145.       {var reg1 uintD* ptr;
  1146.        ptr = &mant1[128/intDsize];
  1147.        doconsttimes(64/intDsize, { *--ptr = 0; } );
  1148.        doconsttimes(32/intDsize, { *--ptr = (uintD)mantlo1; mantlo1 = mantlo1>>intDsize; } );
  1149.        doconsttimes(32/intDsize, { *--ptr = (uintD)manthi1; manthi1 = manthi1>>intDsize; } );
  1150.       }
  1151.       {var reg1 uintD* ptr;
  1152.        ptr = &mant2[64/intDsize];
  1153.        doconsttimes(32/intDsize, { *--ptr = (uintD)mantlo2; mantlo2 = mantlo2>>intDsize; } );
  1154.        doconsttimes(32/intDsize, { *--ptr = (uintD)manthi2; manthi2 = manthi2>>intDsize; } );
  1155.       }
  1156.       #endif
  1157.       {var reg2 uintL mantlo;
  1158.        #ifdef intQsize
  1159.        var reg1 uint64 manthi;
  1160.        #else
  1161.        var reg1 uintL manthi;
  1162.        #endif
  1163.        {SAVE_NUM_STACK # num_stack retten
  1164.         var DS q;
  1165.         var DS r;
  1166.         UDS_divide(&mant1[0],128/intDsize,&mant1[128/intDsize],
  1167.                    &mant2[0],64/intDsize,&mant2[64/intDsize],
  1168.                    &q, &r
  1169.                   );
  1170.         # Es ist 2^53 <= q < 2^55, also q.len = ceiling(54/intDsize)=ceiling(55/intDsize),
  1171.         # und r=0 genau dann, wenn r.len=0.
  1172.         ASSERT(q.len==ceiling(54,intDsize))
  1173.         RESTORE_NUM_STACK # num_stack zurück
  1174.         {var reg3 uintD* ptr = q.MSDptr;
  1175.          manthi = get_max32_Dptr(23,ptr);
  1176.          mantlo = get_32_Dptr(&ptr[ceiling(23,intDsize)]);
  1177.         }
  1178.         # q = 2^32*manthi+mantlo.
  1179.         #ifdef intQsize
  1180.         manthi = (manthi<<32) | (uint64)mantlo;
  1181.         if (manthi >= bit(DF_mant_len+2))
  1182.           # Quotient >=2^54 -> 2 Bits wegrunden
  1183.           { var reg2 uint64 rounding_bits = manthi & (bit(2)-1);
  1184.             exp1 += 1; # Exponenten incrementieren
  1185.             manthi = manthi >> 2;
  1186.             if ( (rounding_bits < bit(1)) # 00,01 werden abgerundet
  1187.                  || ( (rounding_bits == bit(1)) # 10
  1188.                       && (r.len == 0) # und genau halbzahlig
  1189.                       && ((manthi & bit(0)) ==0) # -> round-to-even
  1190.                )    )
  1191.               # abrunden
  1192.               {}
  1193.               else
  1194.               # aufrunden
  1195.               { manthi += 1; }
  1196.           }
  1197.           else
  1198.           # Quotient <2^54 -> 1 Bit wegrunden
  1199.           { var reg2 uint64 rounding_bit = manthi & bit(0);
  1200.             manthi = manthi >> 1;
  1201.             if ( (rounding_bit == 0) # 0 wird abgerundet
  1202.                  || ( (r.len == 0) # genau halbzahlig
  1203.                       && ((manthi & bit(0)) ==0) # -> round-to-even
  1204.                )    )
  1205.               # abrunden
  1206.               {}
  1207.               else
  1208.               # aufrunden
  1209.               { manthi += 1;
  1210.                 if (manthi >= bit(DF_mant_len+1)) # rounding overflow?
  1211.                   { manthi = manthi>>1; exp1 = exp1+1; }
  1212.           }   }
  1213.         #else
  1214.         if (manthi >= bit(DF_mant_len-32+2))
  1215.           # Quotient >=2^54 -> 2 Bits wegrunden
  1216.           { var reg2 uintL rounding_bits = mantlo & (bit(2)-1);
  1217.             exp1 += 1; # Exponenten incrementieren
  1218.             mantlo = (mantlo >> 2) | (manthi << 30); manthi = manthi >> 2;
  1219.             if ( (rounding_bits < bit(1)) # 00,01 werden abgerundet
  1220.                  || ( (rounding_bits == bit(1)) # 10
  1221.                       && (r.len == 0) # und genau halbzahlig
  1222.                       && ((mantlo & bit(0)) ==0) # -> round-to-even
  1223.                )    )
  1224.               # abrunden
  1225.               {}
  1226.               else
  1227.               # aufrunden
  1228.               { mantlo += 1; if (mantlo==0) { manthi += 1; } }
  1229.           }
  1230.           else
  1231.           # Quotient <2^54 -> 1 Bit wegrunden
  1232.           { var reg2 uintL rounding_bit = mantlo & bit(0);
  1233.             mantlo = (mantlo >> 1) | (manthi << 31); manthi = manthi >> 1;
  1234.             if ( (rounding_bit == 0) # 0 wird abgerundet
  1235.                  || ( (r.len == 0) # genau halbzahlig
  1236.                       && ((mantlo & bit(0)) ==0) # -> round-to-even
  1237.                )    )
  1238.               # abrunden
  1239.               {}
  1240.               else
  1241.               # aufrunden
  1242.               { mantlo += 1;
  1243.                 if (mantlo==0)
  1244.                   { manthi += 1;
  1245.                     if (manthi >= bit(DF_mant_len-32+1)) # rounding overflow?
  1246.                       { manthi = manthi>>1; exp1 = exp1+1; }
  1247.           }   }   }
  1248.         #endif
  1249.        }
  1250.        #ifdef intQsize
  1251.        encode_DF(sign1,exp1,manthi, return);
  1252.        #else
  1253.        encode2_DF(sign1,exp1,manthi,mantlo, return);
  1254.        #endif
  1255.     }}}
  1256. #endif
  1257.  
  1258. # Liefert zu einem Double-Float x>=0 : (sqrt x), ein DF.
  1259. # DF_sqrt_DF(x)
  1260. # kann GC auslösen
  1261.   local object DF_sqrt_DF (object x);
  1262. # Methode:
  1263. # x = 0.0 -> Ergebnis 0.0
  1264. # Ergebnis-Vorzeichen := positiv,
  1265. # Ergebnis-Exponent := ceiling(e/2),
  1266. # Ergebnis-Mantisse:
  1267. #   Bilde aus [1,m51,...,m0,(55 Nullbits)] bei geradem e,
  1268. #         aus [0,1,m51,...,m0,(54 Nullbits)] bei ungeradem e
  1269. #   die Ganzzahl-Wurzel, eine 54-Bit-Zahl mit einer führenden 1.
  1270. #   Runde das letzte Bit weg:
  1271. #     Bit 0 = 0 -> abrunden,
  1272. #     Bit 0 = 1 und Wurzel exakt -> round-to-even,
  1273. #     Bit 0 = 1 und Rest >0 -> aufrunden.
  1274. #   Dabei um ein Bit nach rechts schieben.
  1275. #   Bei Aufrundung auf 2^53 (rounding overflow) Mantisse um 1 Bit nach rechts
  1276. #     schieben und Exponent incrementieren.
  1277. #ifdef intQsize # && (intDsize==32)
  1278.   local object DF_sqrt_DF(x)
  1279.     var reg5 object x;
  1280.     { # x entpacken:
  1281.       var reg4 sintWL exp;
  1282.       var reg1 uint64 mantx;
  1283.       DF_decode(x, { return x; },_EMA_,exp=,mantx=);
  1284.       # Um die 128-Bit-Ganzzahl-Wurzel ausnutzen zu können, fügen wir beim
  1285.       # Radikanden 74 bzw. 75 statt 54 bzw. 55 Nullbits an.
  1286.       if (exp & bit(0))
  1287.         # e ungerade
  1288.         { mantx = mantx << (63-(DF_mant_len+1)); exp = exp+1; }
  1289.         else
  1290.         # e gerade
  1291.         { mantx = mantx << (64-(DF_mant_len+1)); }
  1292.       exp = exp >> 1; # exp := exp/2
  1293.      {var uintD mant [128/intDsize];
  1294.       set_32_Dptr(mant,(uint32)(mantx>>32)); set_32_Dptr(&mant[32/intDsize],(uint32)mantx);
  1295.         set_32_Dptr(&mant[2*32/intDsize],0); set_32_Dptr(&mant[3*32/intDsize],0);
  1296.       {SAVE_NUM_STACK # num_stack retten
  1297.        var DS wurzel;
  1298.        var reg6 boolean exactp;
  1299.        UDS_sqrt(&mant[0],128/intDsize,&mant[128/intDsize], &wurzel, exactp=);
  1300.        # wurzel = isqrt(2^74_75 * mant), eine 64-Bit-Zahl.
  1301.        RESTORE_NUM_STACK # num_stack zurück
  1302.        {var reg3 uintD* ptr = wurzel.MSDptr;
  1303.         mantx = ((uint64)get_32_Dptr(ptr) << 32) | (uint64)get_32_Dptr(&ptr[32/intDsize]);
  1304.        }
  1305.        # Die hinteren 63-DF_mant_len Bits wegrunden:
  1306.        if ( ((mantx & bit(62-DF_mant_len)) ==0) # Bit 10 =0 -> abrunden
  1307.             || ( ((mantx & (bit(62-DF_mant_len)-1)) ==0) # Bit 10 =1 und Bits 9..0 >0 -> aufrunden
  1308.                  && exactp                   # Bit 10 =1 und Bits 9..0 =0, aber Rest -> aufrunden
  1309.                  # round-to-even, je nach Bit 11 :
  1310.                  && ((mantx & bit(63-DF_mant_len)) ==0)
  1311.           )    )
  1312.          # abrunden
  1313.          { mantx = mantx >> (63-DF_mant_len); }
  1314.          else
  1315.          # aufrunden
  1316.          { mantx = mantx >> (63-DF_mant_len);
  1317.            mantx += 1;
  1318.            if (mantx >= bit(DF_mant_len+1)) # rounding overflow?
  1319.              { mantx = mantx>>1; exp = exp+1; }
  1320.          }
  1321.       }}
  1322.       encode_DF(0,exp,mantx, return);
  1323.     }
  1324. #else
  1325.   local object DF_sqrt_DF(x)
  1326.     var reg5 object x;
  1327.     { # x entpacken:
  1328.       var reg4 sintWL exp;
  1329.       var reg1 uint32 manthi;
  1330.       var reg2 uint32 mantlo;
  1331.       DF_decode2(x, { return x; },_EMA_,exp=,manthi=,mantlo=);
  1332.       # Um die 128-Bit-Ganzzahl-Wurzel ausnutzen zu können, fügen wir beim
  1333.       # Radikanden 74 bzw. 75 statt 54 bzw. 55 Nullbits an.
  1334.       if (exp & bit(0))
  1335.         # e ungerade
  1336.         { manthi = (manthi << (63-(DF_mant_len+1))) | (mantlo >> ((DF_mant_len+1)-31));
  1337.           mantlo = mantlo << (63-(DF_mant_len+1));
  1338.           exp = exp+1;
  1339.         }
  1340.         else
  1341.         # e gerade
  1342.         { manthi = (manthi << (64-(DF_mant_len+1))) | (mantlo >> ((DF_mant_len+1)-32));
  1343.           mantlo = mantlo << (64-(DF_mant_len+1));
  1344.         }
  1345.       exp = exp >> 1; # exp := exp/2
  1346.      {var uintD mant [128/intDsize];
  1347.       #if (intDsize==32) || (intDsize==16) || (intDsize==8)
  1348.       set_32_Dptr(mant,manthi); set_32_Dptr(&mant[32/intDsize],mantlo);
  1349.         set_32_Dptr(&mant[2*32/intDsize],0); set_32_Dptr(&mant[3*32/intDsize],0);
  1350.       #else
  1351.       {var reg3 uintD* ptr;
  1352.        ptr = &mant[128/intDsize];
  1353.        doconsttimes(64/intDsize, { *--ptr = 0; } );
  1354.        doconsttimes(32/intDsize, { *--ptr = (uintD)mantlo; mantlo = mantlo>>intDsize; } );
  1355.        doconsttimes(32/intDsize, { *--ptr = (uintD)manthi; manthi = manthi>>intDsize; } );
  1356.       }
  1357.       #endif
  1358.       {SAVE_NUM_STACK # num_stack retten
  1359.        var DS wurzel;
  1360.        var reg6 boolean exactp;
  1361.        UDS_sqrt(&mant[0],128/intDsize,&mant[128/intDsize], &wurzel, exactp=);
  1362.        # wurzel = isqrt(2^74_75 * mant), eine 64-Bit-Zahl.
  1363.        RESTORE_NUM_STACK # num_stack zurück
  1364.        {var reg3 uintD* ptr = wurzel.MSDptr;
  1365.         manthi = get_32_Dptr(ptr); mantlo = get_32_Dptr(&ptr[32/intDsize]);
  1366.        }
  1367.        # Die hinteren 63-DF_mant_len Bits wegrunden:
  1368.        if ( ((mantlo & bit(62-DF_mant_len)) ==0) # Bit 10 =0 -> abrunden
  1369.             || ( ((mantlo & (bit(62-DF_mant_len)-1)) ==0) # Bit 10 =1 und Bits 9..0 >0 -> aufrunden
  1370.                  && exactp                   # Bit 10 =1 und Bits 9..0 =0, aber Rest -> aufrunden
  1371.                  # round-to-even, je nach Bit 11 :
  1372.                  && ((mantlo & bit(63-DF_mant_len)) ==0)
  1373.           )    )
  1374.          # abrunden
  1375.          { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
  1376.            manthi = manthi >> (63-DF_mant_len);
  1377.          }
  1378.          else
  1379.          # aufrunden
  1380.          { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
  1381.            manthi = manthi >> (63-DF_mant_len);
  1382.            mantlo += 1;
  1383.            if (mantlo==0)
  1384.              { manthi += 1;
  1385.                if (manthi >= bit(DF_mant_len-32+1)) # rounding overflow?
  1386.                  { manthi = manthi>>1; exp = exp+1; }
  1387.          }   }
  1388.       }}
  1389.       encode2_DF(0,exp,manthi,mantlo, return);
  1390.     }
  1391. #endif
  1392.  
  1393. # DF_to_I(x) wandelt ein Double-Float x, das eine ganze Zahl darstellt,
  1394. # in ein Integer um.
  1395. # kann GC auslösen
  1396.   local object DF_to_I (object x);
  1397. # Methode:
  1398. # Falls x=0.0, Ergebnis 0.
  1399. # Sonst (ASH Vorzeichen*Mantisse (e-53)).
  1400. #ifdef intQsize
  1401.   local object DF_to_I(x)
  1402.     var reg3 object x;
  1403.     { # x entpacken:
  1404.       var reg4 signean sign;
  1405.       var reg2 sintWL exp;
  1406.       var reg1 uint64 mant;
  1407.       DF_decode(x, { return Fixnum_0; }, sign=,exp=,mant=);
  1408.       exp = exp-(DF_mant_len+1);
  1409.       # mant mit Vorzeichen versehen:
  1410.       if (!(sign==0)) { mant = -mant; }
  1411.       # in ein Bignum umwandeln und shiften:
  1412.       return I_I_ash_I( Q_to_I(mant), L_to_FN(exp) );
  1413.     }
  1414. #else
  1415.   local object DF_to_I(x)
  1416.     var reg4 object x;
  1417.     { # x entpacken:
  1418.       var reg5 signean sign;
  1419.       var reg3 sintWL exp;
  1420.       var reg1 uint32 manthi;
  1421.       var reg2 uint32 mantlo;
  1422.       DF_decode2(x, { return Fixnum_0; }, sign=,exp=,manthi=,mantlo=);
  1423.       exp = exp-(DF_mant_len+1);
  1424.       # mant mit Vorzeichen versehen:
  1425.       if (!(sign==0))
  1426.         { manthi = -manthi; mantlo = -mantlo; if (!(mantlo==0)) { manthi -= 1; } }
  1427.       # in ein Bignum umwandeln und shiften:
  1428.       return I_I_ash_I( L2_to_I(manthi,mantlo), L_to_FN(exp) );
  1429.     }
  1430. #endif
  1431.  
  1432. # I_to_DF(x) wandelt ein Integer x in ein Double-Float um und rundet dabei.
  1433. # kann GC auslösen
  1434.   local object I_to_DF (object x);
  1435. # Methode:
  1436. # x=0 -> Ergebnis 0.0
  1437. # Merke Vorzeichen von x.
  1438. # x:=(abs x)
  1439. # Exponent:=(integer-length x)
  1440. #   Greife die 54 höchstwertigen Bits heraus (angeführt von einer 1).
  1441. #   Runde das letzte Bit weg:
  1442. #     Bit 0 = 0 -> abrunden,
  1443. #     Bit 0 = 1 und Rest =0 -> round-to-even,
  1444. #     Bit 0 = 1 und Rest >0 -> aufrunden.
  1445. #   Dabei um ein Bit nach rechts schieben.
  1446. #   Bei Aufrundung auf 2^53 (rounding overflow) Mantisse um 1 Bit nach rechts
  1447. #     schieben und Exponent incrementieren.
  1448.   local object I_to_DF(x)
  1449.     var reg7 object x;
  1450.     { if (eq(x,Fixnum_0)) { return DF_0; }
  1451.      {var reg8 signean sign = R_sign(x); # Vorzeichen
  1452.       if (!(sign==0)) { x = I_minus_I(x); } # bei x<0: x := (- x)
  1453.       {   var reg9 uintL exp = I_integer_length(x); # (integer-length x)
  1454.           # NDS zu x>0 bilden:
  1455.        {  var reg2 uintD* MSDptr;
  1456.           var reg5 uintC len;
  1457.           I_to_NDS_nocopy(x, MSDptr=,len=,_EMA_);
  1458.           # MSDptr/len/LSDptr ist die NDS zu x, len>0.
  1459.           # Führende Digits holen: Brauche DF_mant_len+1 Bits, dazu intDsize
  1460.           # Bits (die NDS kann mit bis zu intDsize Nullbits anfangen).
  1461.           # Dann werden diese Bits um (exp mod intDsize) nach rechts geschoben.
  1462.         { var reg4 uintD msd = *MSDptr++; # erstes Digit
  1463.           var reg1 uint32 msdd = 0; # weitere min(len-1,32/intDsize) Digits
  1464.           var reg1 uint32 msddf = 0; # weitere maximal 32/intDsize Digits
  1465.           #define NEXT_DIGIT(i)  \
  1466.             { if (--len == 0) goto ok;                            \
  1467.               msdd |= (uint32)(*MSDptr++) << (32-(i+1)*intDsize); \
  1468.             }
  1469.           DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
  1470.           #undef NEXT_DIGIT
  1471.           #define NEXT_DIGIT(i)  \
  1472.             { if (--len == 0) goto ok;                             \
  1473.               msddf |= (uint32)(*MSDptr++) << (32-(i+1)*intDsize); \
  1474.             }
  1475.           DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
  1476.           #undef NEXT_DIGIT
  1477.           --len; ok:
  1478.           # Die NDS besteht aus msd, msdd, msddf und len weiteren Digits.
  1479.           # Das höchste in 2^64*msd+2^32*msdd+msddf gesetzte Bit ist Bit Nummer
  1480.           # 63 + (exp mod intDsize).
  1481.          {var reg6 uintL shiftcount = exp % intDsize;
  1482.           #ifdef intQsize
  1483.           var reg3 uint64 mant = # führende 64 Bits
  1484.             (shiftcount==0
  1485.               ? (((uint64)msdd << 32) | (uint64)msddf)
  1486.               : (((uint64)msd << (64-shiftcount)) | ((uint64)msdd << (32-shiftcount)) | ((uint64)msddf >> shiftcount))
  1487.             );
  1488.           # Das höchste in mant gesetzte Bit ist Bit Nummer 63.
  1489.           if ( ((mant & bit(62-DF_mant_len)) ==0) # Bit 10 =0 -> abrunden
  1490.                || ( ((mant & (bit(62-DF_mant_len)-1)) ==0) # Bit 10 =1 und Bits 9..0 =0
  1491.                     && ((msddf & (bit(shiftcount)-1)) ==0) # und weitere Bits aus msddf =0
  1492.                     && (!test_loop_up(MSDptr,len)) # und alle weiteren Digits =0
  1493.                     # round-to-even, je nach Bit 11 :
  1494.                     && ((mant & bit(63-DF_mant_len)) ==0)
  1495.              )    )
  1496.             # abrunden
  1497.             { mant = mant >> (63-DF_mant_len); }
  1498.             else
  1499.             # aufrunden
  1500.             { mant = mant >> (63-DF_mant_len);
  1501.               mant += 1;
  1502.               if (mant >= bit(DF_mant_len+1)) # rounding overflow?
  1503.                 { mant = mant>>1; exp = exp+1; }
  1504.             }
  1505.           encode_DF(sign,(sintL)exp,mant, return);
  1506.           #else
  1507.           var reg3 uint32 manthi; # führende 32 Bits
  1508.           var reg3 uint32 mantlo; # nächste 32 Bits
  1509.           if (shiftcount==0)
  1510.             { manthi = msdd; mantlo = msddf; }
  1511.             else
  1512.             { manthi = ((uint32)msd << (32-shiftcount)) | (msdd >> shiftcount);
  1513.               mantlo = (msdd << (32-shiftcount)) | (msddf >> shiftcount);
  1514.             }
  1515.           # Das höchste in mant gesetzte Bit ist Bit Nummer 63.
  1516.           if ( ((mantlo & bit(62-DF_mant_len)) ==0) # Bit 10 =0 -> abrunden
  1517.                || ( ((mantlo & (bit(62-DF_mant_len)-1)) ==0) # Bit 10 =1 und Bits 9..0 =0
  1518.                     && ((msddf & (bit(shiftcount)-1)) ==0) # und weitere Bits aus msddf =0
  1519.                     && (!test_loop_up(MSDptr,len)) # und alle weiteren Digits =0
  1520.                     # round-to-even, je nach Bit 11 :
  1521.                     && ((mantlo & bit(63-DF_mant_len)) ==0)
  1522.              )    )
  1523.             # abrunden
  1524.             { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
  1525.               manthi = manthi >> (63-DF_mant_len);
  1526.             }
  1527.             else
  1528.             # aufrunden
  1529.             { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
  1530.               manthi = manthi >> (63-DF_mant_len);
  1531.               mantlo += 1;
  1532.               if (mantlo==0)
  1533.                 { manthi += 1;
  1534.                   if (manthi >= bit(DF_mant_len-32+1)) # rounding overflow?
  1535.                     { manthi = manthi>>1; exp = exp+1; }
  1536.             }   }
  1537.           encode2_DF(sign,(sintL)exp,manthi,mantlo, return);
  1538.           #endif
  1539.     }}}}}}
  1540.  
  1541. # RA_to_DF(x) wandelt eine rationale Zahl x in ein Double-Float um
  1542. # und rundet dabei.
  1543. # kann GC auslösen
  1544.   local object RA_to_DF (object x);
  1545. # Methode:
  1546. # x ganz -> klar.
  1547. # x = +/- a/b mit Integers a,b>0:
  1548. #   Seien n,m so gewählt, daß
  1549. #     2^(n-1) <= a < 2^n, 2^(m-1) <= b < 2^m.
  1550. #   Dann ist 2^(n-m-1) < a/b < 2^(n-m+1).
  1551. #   Berechne n=(integer-length a) und m=(integer-length b) und
  1552. #   floor(2^(-n+m+54)*a/b) :
  1553. #   Bei n-m>=54 dividiere a durch (ash b (n-m-54)),
  1554. #   bei n-m<54 dividiere (ash a (-n+m+54)) durch b.
  1555. #   Der erste Wert ist >=2^53, <2^55.
  1556. #   Falls er >=2^54 ist, runde 2 Bits weg,
  1557. #   falls er <2^54 ist, runde 1 Bit weg.
  1558.   local object RA_to_DF(x)
  1559.     var reg3 object x;
  1560.     { if (RA_integerp(x)) { return I_to_DF(x); }
  1561.       # x Ratio
  1562.       pushSTACK(TheRatio(x)->rt_den); # b
  1563.       x = TheRatio(x)->rt_num; # +/- a
  1564.      {var reg7 signean sign = R_sign(x); # Vorzeichen
  1565.       if (!(sign==0)) { x = I_minus_I(x); } # Betrag nehmen, liefert a
  1566.       pushSTACK(x);
  1567.       # Stackaufbau: b, a.
  1568.       {var reg4 sintL lendiff = I_integer_length(x) # (integer-length a)
  1569.                                 - I_integer_length(STACK_1); # (integer-length b)
  1570.        if (lendiff > DF_exp_high-DF_exp_mid) # Exponent >= n-m > Obergrenze ?
  1571.          { fehler_overflow(); } # -> Overflow
  1572.        if (lendiff < DF_exp_low-DF_exp_mid-2) # Exponent <= n-m+2 < Untergrenze ?
  1573.          { if (underflow_allowed())
  1574.              { fehler_underflow(); } # -> Underflow
  1575.              else
  1576.              { skipSTACK(2); return DF_0; }
  1577.          }
  1578.        { var reg5 object zaehler;
  1579.          var reg6 object nenner;
  1580.          if (lendiff >= DF_mant_len+2)
  1581.            # n-m-54>=0
  1582.            { nenner = I_I_ash_I(STACK_1,fixnum((uint32)(lendiff - (DF_mant_len+2)))); # (ash b n-m-54)
  1583.              zaehler = popSTACK(); # a
  1584.              skipSTACK(1);
  1585.            }
  1586.            else
  1587.            { zaehler = I_I_ash_I(popSTACK(),fixnum((uint32)((DF_mant_len+2) - lendiff))); # (ash a -n+m+54)
  1588.              nenner = popSTACK(); # b
  1589.            }
  1590.          # Division zaehler/nenner durchführen:
  1591.          I_I_divide_I_I(zaehler,nenner);
  1592.          # Stackaufbau: q, r.
  1593.          # 2^53 <= q < 2^55, also ist q Bignum mit ceiling(55/intDsize) Digits.
  1594.         {var reg3 uintD* ptr = &TheBignum(STACK_1)->data[0];
  1595.          #ifdef intQsize
  1596.          var reg1 uint64 mant =
  1597.            ((uint64)get_max32_Dptr(23,ptr) << 32)
  1598.            | (uint64)get_32_Dptr(&ptr[ceiling(23,intDsize)]);
  1599.          if (mant >= bit(DF_mant_len+2))
  1600.            # 2^54 <= q < 2^55, schiebe um 2 Bits nach rechts
  1601.            { var reg2 uint64 rounding_bits = mant & (bit(2)-1);
  1602.              lendiff = lendiff+1; # Exponent := n-m+1
  1603.              mant = mant >> 2;
  1604.              if ( (rounding_bits < bit(1)) # 00,01 werden abgerundet
  1605.                   || ( (rounding_bits == bit(1)) # 10
  1606.                        && (eq(STACK_0,Fixnum_0)) # und genau halbzahlig (r=0)
  1607.                        && ((mant & bit(0)) ==0) # -> round-to-even
  1608.                 )    )
  1609.                # abrunden
  1610.                goto ab;
  1611.                else
  1612.                # aufrunden
  1613.                goto auf;
  1614.            }
  1615.            else
  1616.            { var reg2 uint64 rounding_bit = mant & bit(0);
  1617.              mant = mant >> 1;
  1618.              if ( (rounding_bit == 0) # 0 wird abgerundet
  1619.                   || ( (eq(STACK_0,Fixnum_0)) # genau halbzahlig (r=0)
  1620.                        && ((mant & bit(0)) ==0) # -> round-to-even
  1621.                 )    )
  1622.                # abrunden
  1623.                goto ab;
  1624.                else
  1625.                # aufrunden
  1626.                goto auf;
  1627.            }
  1628.          auf:
  1629.          mant += 1;
  1630.          if (mant >= bit(DF_mant_len+1)) # rounding overflow?
  1631.            { mant = mant>>1; lendiff = lendiff+1; }
  1632.          ab:
  1633.          skipSTACK(2);
  1634.          # Fertig.
  1635.          encode_DF(sign,lendiff,mant, return);
  1636.          #else
  1637.          var reg1 uint32 manthi = get_max32_Dptr(23,ptr);
  1638.          var reg1 uint32 mantlo = get_32_Dptr(&ptr[ceiling(23,intDsize)]);
  1639.          if (manthi >= bit(DF_mant_len-32+2))
  1640.            # 2^54 <= q < 2^55, schiebe um 2 Bits nach rechts
  1641.            { var reg2 uintL rounding_bits = mantlo & (bit(2)-1);
  1642.              lendiff = lendiff+1; # Exponent := n-m+1
  1643.              mantlo = (mantlo >> 2) | (manthi << 30); manthi = manthi >> 2;
  1644.              if ( (rounding_bits < bit(1)) # 00,01 werden abgerundet
  1645.                   || ( (rounding_bits == bit(1)) # 10
  1646.                        && (eq(STACK_0,Fixnum_0)) # und genau halbzahlig (r=0)
  1647.                        && ((mantlo & bit(0)) ==0) # -> round-to-even
  1648.                 )    )
  1649.                # abrunden
  1650.                goto ab;
  1651.                else
  1652.                # aufrunden
  1653.                goto auf;
  1654.            }
  1655.            else
  1656.            { var reg2 uintL rounding_bit = mantlo & bit(0);
  1657.              mantlo = (mantlo >> 1) | (manthi << 31); manthi = manthi >> 1;
  1658.              if ( (rounding_bit == 0) # 0 wird abgerundet
  1659.                   || ( (eq(STACK_0,Fixnum_0)) # genau halbzahlig (r=0)
  1660.                        && ((mantlo & bit(0)) ==0) # -> round-to-even
  1661.                 )    )
  1662.                # abrunden
  1663.                goto ab;
  1664.                else
  1665.                # aufrunden
  1666.                goto auf;
  1667.            }
  1668.          auf:
  1669.          mantlo += 1;
  1670.          if (mantlo==0)
  1671.            { manthi += 1;
  1672.              if (manthi >= bit(DF_mant_len-32+1)) # rounding overflow?
  1673.                { manthi = manthi>>1; lendiff = lendiff+1; }
  1674.            }
  1675.          ab:
  1676.          skipSTACK(2);
  1677.          # Fertig.
  1678.          encode2_DF(sign,lendiff,manthi,mantlo, return);
  1679.          #endif
  1680.     }}}}}
  1681.  
  1682.