home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / Misc / CLISP-1.LHA / CLISP960530-sr.lha / src / lfloat.d < prev    next >
Encoding:
Text File  |  1996-04-15  |  63.8 KB  |  1,435 lines
  1. # Grundfunktionen für Long-Floats
  2.  
  3. # Fehlermeldung bei zu langen Long-FLoats
  4.   nonreturning_function(local, fehler_LF_toolong, (void));
  5.   local void fehler_LF_toolong()
  6.     { 
  7.       //: DEUTSCH "Zu lange Long-Floats"
  8.       //: ENGLISH "long float too long"
  9.       //: FRANCAIS "LONG-FLOAT trop long"
  10.       fehler(error, GETTEXT("long float too long"));
  11.     }
  12.  
  13. # Entpacken eines Long-Float:
  14. # LF_decode(obj, zero_statement, sign=,exp=,mantMSDptr=,mantlen=,mantLSDptr=);
  15. # zerlegt ein Long-Float obj.
  16. # Ist obj=0.0, wird zero_statement ausgeführt.
  17. # Sonst: signean sign = Vorzeichen (0 = +, -1 = -),
  18. #        sintL exp = Exponent (vorzeichenbehaftet),
  19. #        UDS mantMSDptr/mantlen/mantLSDptr = Mantisse
  20. #          (>= 2^(intDsize*mantlen-1), < 2^(intDsize*mantlen)),
  21. #          mit mantlen>=LF_minlen.
  22.   #define LF_decode(obj, zero_statement, sign_zuweisung,exp_zuweisung,mantMSDptr_zuweisung,mantlen_zuweisung,mantLSDptr_zuweisung)  \
  23.     { var reg3 object _obj = (obj);                                               \
  24.       var reg1 Lfloat _x = TheLfloat(_obj);                                       \
  25.       var reg2 uintL uexp = _x->expo;                                             \
  26.       if (uexp==0)                                                                \
  27.         { mantlen_zuweisung _x->len; zero_statement } # e=0 -> Zahl 0.0           \
  28.         else                                                                      \
  29.         { exp_zuweisung (sintL)(uexp - LF_exp_mid);     # Exponent                \
  30.           sign_zuweisung R_sign(_obj);                  # Vorzeichen              \
  31.           unused (mantMSDptr_zuweisung &(_x->data[0])); # Mantissen-UDS           \
  32.           mantLSDptr_zuweisung &(_x->data[(uintP)( mantlen_zuweisung _x->len )]); \
  33.     }   }
  34.  
  35. # Einpacken eines Long-Float:
  36. # encode_LF0(len,erg_zuweisung) liefert ein Long-Float 0.0 mit len Digits.
  37. # > uintC len: Anzahl der Digits
  38. # < object erg: neues Long-Float 0.0 mit len Digits
  39. # kann GC auslösen
  40.   #define encode_LF0(len,erg_zuweisung)  \
  41.     { var reg2 uintC _len = (len);                                                 \
  42.       var reg1 object _erg = allocate_lfloat(_len,0,0); # Exponent 0, Vorzeichen + \
  43.       clear_loop_up(&TheLfloat(_erg)->data[0],_len); # Mantisse := 0               \
  44.       erg_zuweisung _erg;                                                          \
  45.     }
  46.  
  47. # Einpacken eines Long-Float:
  48. # encode_LF1s(sign,len,erg_zuweisung) liefert ein Long-Float +-1.0 mit len Digits.
  49. # > signean sign: Vorzeichen
  50. # > uintC len: Anzahl der Digits
  51. # < object erg: neues Long-Float +1.0 oder -1.0 mit len Digits
  52. # kann GC auslösen
  53.   #define encode_LF1s(sign,len,erg_zuweisung)  \
  54.     { var reg2 uintC _len = (len);                                                   \
  55.       var reg1 object _erg = allocate_lfloat(_len,LF_exp_mid+1,(sign)); # Exponent 1 \
  56.       TheLfloat(_erg)->data[0] = bit(intDsize-1); # Mantisse := 2^(intDsize*len-1)   \
  57.       clear_loop_up(&TheLfloat(_erg)->data[1],_len-1);                               \
  58.       erg_zuweisung _erg;                                                            \
  59.     }
  60.  
  61. # Einpacken eines Long-Float:
  62. # encode_LF1(len,erg_zuweisung) liefert ein Long-Float 1.0 mit len Digits.
  63. # > uintC len: Anzahl der Digits
  64. # < object erg: neues Long-Float 1.0 mit len Digits
  65. # kann GC auslösen
  66.   #define encode_LF1(len,erg_zuweisung)  encode_LF1s(0,len,erg_zuweisung)
  67.  
  68. # Einpacken eines Long-Float:
  69. # encode_LFu(sign,uexp,mantMSDptr,mantlen, erg_zuweisung) liefert ein Long-Float
  70. # > signean sign: Vorzeichen
  71. # > uintL exp: Exponent + LF_exp_mid
  72. # > uintD* mantMSDptr: Pointer auf eine NUDS mit gesetztem höchstem Bit
  73. # > uintC mantlen: Anzahl der Digits, >= LF_minlen
  74. # < object erg: neues Long-Float mit der UDS mantMSDptr/mantlen/.. als Mantisse
  75. # Der Exponent wird nicht auf Überlauf/Unterlauf getestet.
  76. # kann GC auslösen
  77.   #define encode_LFu(sign,uexp,mantMSDptr,mantlen,erg_zuweisung)  \
  78.     { var reg2 uintC _len = (mantlen);                                                 \
  79.       var reg1 object _erg = allocate_lfloat(_len,uexp,(sign)); # Exponent             \
  80.       copy_loop_up((mantMSDptr),&TheLfloat(_erg)->data[0],_len); # Mantisse übertragen \
  81.       erg_zuweisung _erg;                                                              \
  82.     }
  83.  
  84. # Einpacken eines Long-Float:
  85. # encode_LF(sign,exp,mantMSDptr,mantlen, erg_zuweisung) liefert ein Long-Float
  86. # > signean sign: Vorzeichen
  87. # > sintL exp: Exponent
  88. # > uintD* mantMSDptr: Pointer auf eine NUDS mit gesetztem höchstem Bit
  89. # > uintC mantlen: Anzahl der Digits, >= LF_minlen
  90. # < object erg: neues Long-Float mit der UDS mantMSDptr/mantlen/.. als Mantisse
  91. # Der Exponent wird nicht auf Überlauf/Unterlauf getestet.
  92. # kann GC auslösen
  93.   #define encode_LF(sign,exp,mantMSDptr,mantlen,erg_zuweisung)  \
  94.     encode_LFu(sign,LF_exp_mid+(uintL)(exp),mantMSDptr,mantlen,_EMA_ erg_zuweisung)
  95.  
  96. # Hash-Code eines Long-Float: Mischung aus Exponent, Länge, erste 32 Bit
  97.   global uint32 hashcode_lfloat (object obj);
  98.   global uint32 hashcode_lfloat(obj)
  99.     var reg1 object obj;
  100.     { return TheLfloat(obj)->expo + TheLfloat(obj)->len
  101.              + get_32_Dptr(&TheLfloat(obj)->data[0]);
  102.     }
  103.  
  104. # LF_zerop(x) stellt fest, ob ein Long-Float x = 0.0 ist.
  105.   #define LF_zerop(x)  (TheLfloat(x)->expo == 0)
  106.  
  107. # Liefert zu einem Long-Float x : (ftruncate x), ein LF.
  108. # LF_ftruncate_LF(x)
  109. # x wird zur 0 hin zur nächsten ganzen Zahl gerundet.
  110. # kann GC auslösen
  111.   local object LF_ftruncate_LF (object x);
  112. # Methode:
  113. # x = 0.0 oder e<=0 -> Ergebnis 0.0
  114. # 1<=e<=16n -> letzte (16n-e) Bits der Mantisse auf 0 setzen,
  115. #              Exponent und Vorzeichen beibehalten
  116. # e>=16n -> Ergebnis x
  117. #if 0
  118.   local object LF_ftruncate_LF(x)
  119.     var reg2 object x;
  120.     { var reg8 signean sign;
  121.       var reg1 sintL exp;
  122.       var reg9 uintD* mantMSDptr;
  123.       var reg4 uintC mantlen;
  124.       LF_decode(x, { return x; }, sign=,exp=,mantMSDptr=,mantlen=,_EMA_);
  125.       if (exp<=0) { encode_LF0(mantlen, return); } # e<=0 -> Ergebnis 0.0
  126.       if ((uintL)exp >= intDsize*(uintL)mantlen) # e>=16n -> x als Ergebnis
  127.         { return x; }
  128.         else
  129.         # 0 < e < 16n
  130.         # neue NUDS erzeugen mit e Bits aus mant und 16n-e Nullbits:
  131.         { SAVE_NUM_STACK # num_stack retten
  132.           var reg7 uintD* MSDptr;
  133.           num_stack_need(mantlen, MSDptr=,_EMA_);
  134.           { var reg5 uintC count = floor((uintL)exp,intDsize); # zu kopierende Digits, < mantlen
  135.             var reg6 uintC bitcount = ((uintL)exp) % intDsize; # zu kopierende Bits danach, >=0, <intDsize
  136.             var reg3 uintD* ptr =
  137.               copy_loop_up(mantMSDptr,MSDptr,count); # count ganze Digits kopieren
  138.             *ptr++ = mantMSDptr[count] & minus_bitm(intDsize-bitcount); # dann bitcount Bits kopieren
  139.             clear_loop_up(ptr,mantlen-count-1); # Rest mit Nullen füllen
  140.           }
  141.           RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  142.           encode_LF(sign,exp,MSDptr,mantlen, return);
  143.     }   }
  144. #else
  145.   local object LF_ftruncate_LF(x)
  146.     var reg4 object x;
  147.     { var reg2 uintC len = TheLfloat(x)->len;
  148.       var reg1 uintL uexp = TheLfloat(x)->expo;
  149.       if (uexp <= LF_exp_mid)
  150.         { if (uexp == 0) { return x; } # x=0.0 -> Ergebnis 0.0
  151.           encode_LF0(len, return); # e<=0 -> Ergebnis 0.0
  152.         }
  153.      {var reg3 uintL exp = uexp - LF_exp_mid;
  154.       if (exp >= intDsize*(uintL)len) # e>=16n -> x als Ergebnis
  155.         { return x; }
  156.       # 0 < e < 16n
  157.       pushSTACK(x);
  158.       {var reg5 object y = allocate_lfloat(len,uexp,R_sign(x)); # neues Long-Float
  159.        x = popSTACK();
  160.        # y_mant := NUDS mit e Bits aus x_mant und 16n-e Nullbits:
  161.        {var reg7 uintC count = floor(exp,intDsize); # zu kopierende Digits, < mantlen
  162.         var reg9 uintC bitcount = exp % intDsize; # zu kopierende Bits danach, >=0, <intDsize
  163.         var reg8 uintD* x_mantMSDptr = &TheLfloat(x)->data[0];
  164.         var reg6 uintD* ptr =
  165.           copy_loop_up(x_mantMSDptr,&TheLfloat(y)->data[0],count); # count ganze Digits kopieren
  166.         *ptr++ = x_mantMSDptr[count] & minus_bitm(intDsize-bitcount); # dann bitcount Bits kopieren
  167.         clear_loop_up(ptr,len-count-1); # Rest mit Nullen füllen
  168.        }
  169.        return y;
  170.     }}}
  171. #endif
  172.  
  173. # Liefert zu einem Long-Float x : (futruncate x), ein LF.
  174. # LF_futruncate_LF(x)
  175. # x wird von der 0 weg zur nächsten ganzen Zahl gerundet.
  176. # kann GC auslösen
  177.   local object LF_futruncate_LF (object x);
  178. # Methode:
  179. # x = 0.0 -> Ergebnis 0.0
  180. # e<=0 -> Ergebnis 1.0 oder -1.0, je nach Vorzeichen von x.
  181. # 1<=e<16n -> Greife die letzten (16n-e) Bits von x heraus.
  182. #             Sind sie alle =0 -> Ergebnis x.
  183. #             Sonst setze sie alle auf 0 und erhöhe dann die vorderen e Bits
  184. #             um 1.
  185. #             Kein Überlauf -> fertig.
  186. #             Sonst (Ergebnis eine Zweierpotenz): Mantisse := .1000...000,
  187. #               e:=e+1. (Test auf Überlauf wegen e<=16n überflüssig)
  188. # e>=16n -> Ergebnis x.
  189. #if 0
  190.   local object LF_futruncate_LF(x)
  191.     var reg1 object x;
  192.     { var reg9 signean sign;
  193.       var reg3 sintL exp;
  194.       var reg2 uintD* mantMSDptr;
  195.       var reg7 uintC mantlen;
  196.       LF_decode(x, { return x; }, sign=,exp=,mantMSDptr=,mantlen=,_EMA_);
  197.       if (exp<=0) { encode_LF1s(sign,mantlen, return); } # e<=0 -> Ergebnis +-1.0
  198.       if ((uintL)exp >= intDsize*(uintL)mantlen) # e>=16n -> x als Ergebnis
  199.         { return x; }
  200.         else
  201.         # 0 < e < 16n
  202.         { # Testen, ob alle hinteren 16n-e Bits =0 sind:
  203.           var reg6 uintC count = floor((uintL)exp,intDsize); # zu kopierende Digits, < mantlen
  204.           var reg5 uintC bitcount = ((uintL)exp) % intDsize; # zu kopierende Bits danach, >=0, <intDsize
  205.           var reg5 uintD mask = minus_bitm(intDsize-bitcount); # Maske mit bitcount Bits
  206.           var reg2 uintD* mantptr = &mantMSDptr[count];
  207.           if (   ((mantptr[0] & ~mask) ==0)
  208.               && !test_loop_up(&mantptr[1],mantlen-count-1)
  209.              )
  210.             { return x; }
  211.           # neue NUDS erzeugen mit e Bits aus mant mit Increment
  212.           # und 16n-e Nullbits:
  213.          {SAVE_NUM_STACK # num_stack retten
  214.           var reg8 uintD* MSDptr;
  215.           num_stack_need(mantlen, MSDptr=,_EMA_);
  216.           { var reg4 uintD* ptr =
  217.               copy_loop_up(mantMSDptr,MSDptr,count); # count ganze Digits kopieren
  218.             if ((ptr[0] = ((mantptr[0] & mask) - mask)) == 0) # dann bitcount Bits kopieren und incrementieren
  219.               { if (!( inc_loop_down(ptr,count) ==0)) # evtl. weiterincrementieren
  220.                   { MSDptr[0] = bit(intDsize-1); exp = exp+1; } # evtl. Exponenten erhöhen
  221.               }
  222.             clear_loop_up(&ptr[1],mantlen-count-1); # Rest mit Nullen füllen
  223.           }
  224.           RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  225.           encode_LF(sign,exp,MSDptr,mantlen, return);
  226.     }   }}
  227. #else
  228.   local object LF_futruncate_LF(x)
  229.     var reg1 object x;
  230.     { var reg3 uintC len = TheLfloat(x)->len;
  231.       var reg7 uintL uexp = TheLfloat(x)->expo;
  232.       if (uexp <= LF_exp_mid)
  233.         { if (uexp == 0) { return x; } # x=0.0 -> Ergebnis 0.0
  234.           encode_LF1s(R_sign(x),len, return); # e<=0 -> Ergebnis +-1.0
  235.         }
  236.      {var reg7 uintL exp = uexp - LF_exp_mid;
  237.       if (exp >= intDsize*(uintL)len) # e>=16n -> x als Ergebnis
  238.         { return x; }
  239.       # 0 < e < 16n
  240.       # Testen, ob alle hinteren 16n-e Bits =0 sind:
  241.       {var reg5 uintC count = floor(exp,intDsize); # zu kopierende Digits, < mantlen
  242.        var reg6 uintC bitcount = exp % intDsize; # zu kopierende Bits danach, >=0, <intDsize
  243.        var reg3 uintD mask = minus_bitm(intDsize-bitcount); # Maske mit bitcount Bits
  244.        {var reg2 uintD* mantptr = &TheLfloat(x)->data[count];
  245.         if (   ((mantptr[0] & ~mask) ==0)
  246.             && !test_loop_up(&mantptr[1],len-count-1)
  247.            )
  248.           { return x; }
  249.        }
  250.        # Nein -> neues Long-Float produzieren:
  251.        pushSTACK(x);
  252.        {var reg2 object y = allocate_lfloat(len,uexp,R_sign(x)); # neues Long-Float
  253.         x = popSTACK();
  254.         # y_mant := NUDS mit e Bits aus x_mant mit Increment und 16n-e Nullbits:
  255.         {var reg8 uintD* x_mantMSDptr = &TheLfloat(x)->data[0];
  256.          var reg9 uintD* y_mantMSDptr = &TheLfloat(y)->data[0];
  257.          var reg4 uintD* ptr =
  258.            copy_loop_up(x_mantMSDptr,y_mantMSDptr,count); # count ganze Digits kopieren
  259.          if ((ptr[0] = ((x_mantMSDptr[count] & mask) - mask)) == 0) # dann bitcount Bits kopieren und incrementieren
  260.            { if (!( inc_loop_down(ptr,count) ==0)) # evtl. weiterincrementieren
  261.                { y_mantMSDptr[0] = bit(intDsize-1); (TheLfloat(y)->expo)++; } # evtl. Exponenten erhöhen
  262.            }
  263.          clear_loop_up(&ptr[1],len-count-1); # Rest mit Nullen füllen
  264.         }
  265.         return y;
  266.     }}}}
  267. #endif
  268.  
  269. # Liefert zu einem Long-Float x : (fround x), ein LF.
  270. # LF_fround_LF(x)
  271. # x wird zur nächsten ganzen Zahl gerundet.
  272. # kann GC auslösen
  273.   local object LF_fround_LF (object x);
  274. # Methode:
  275. # x = 0.0 oder e<0 -> Ergebnis 0.0
  276. # 0<=e<16n -> letzte (16n-e) Bits der Mantisse wegrunden,
  277. #             Exponent und Vorzeichen beibehalten.
  278. # e>=16n -> Ergebnis x
  279. #if 0
  280.   local object LF_fround_LF(x)
  281.     var reg4 object x;
  282.     { var reg9 signean sign;
  283.       var reg5 sintL exp;
  284.       var reg1 uintD* mantMSDptr;
  285.       var reg7 uintC mantlen;
  286.       LF_decode(x, { return x; }, sign=,exp=,mantMSDptr=,mantlen=,_EMA_);
  287.       if (exp<0) { encode_LF0(mantlen, return); } # e<0 -> Ergebnis 0.0
  288.       if ((uintL)exp >= intDsize*(uintL)mantlen) # e>=16n -> x als Ergebnis
  289.         { return x; }
  290.         else
  291.         # 0 <= e < 16n
  292.         { # alle hinteren 16n-e Bits wegrunden:
  293.           var reg6 uintC count = floor((uintL)exp,intDsize); # zu kopierende Digits, < mantlen
  294.           var reg6 uintC bitcount = ((uintL)exp) % intDsize; # zu kopierende Bits danach, >=0, <intDsize
  295.           var reg3 uintD mask = minus_bit(intDsize-bitcount-1); # Maske mit bitcount+1 Bits
  296.           var reg1 uintD* mantptr = &mantMSDptr[count];
  297.           if ((mantptr[0] & -mask) ==0) goto ab; # Bit 16n-e-1 =0 -> abrunden
  298.           if (!((mantptr[0] & ~mask) ==0)) goto auf; # Bit 16n-e-1 =1 und Bits 16n-e-2..0 >0 -> aufrunden
  299.           if (test_loop_up(&mantptr[1],mantlen-count-1)) goto auf;
  300.           # round-to-even, je nach Bit 16n-e :
  301.           if (bitcount>0)
  302.             { if ((mantptr[0] & (-2*mask)) ==0) goto ab; else goto auf; }
  303.             elif (count>0)
  304.               { if ((mantptr[-1] & bit(0)) ==0) goto ab; else goto auf; }
  305.               else
  306.               # bitcount=0, count=0, also exp=0: Abrunden von +-0.5 zu 0.0
  307.               { encode_LF0(mantlen, return); }
  308.           ab: # abrunden
  309.           { SAVE_NUM_STACK # num_stack retten
  310.             var reg8 uintD* MSDptr;
  311.             num_stack_need(mantlen, MSDptr=,_EMA_);
  312.            {var reg2 uintD* ptr =
  313.               copy_loop_up(mantMSDptr,MSDptr,count); # count ganze Digits kopieren
  314.             *ptr++ = mantMSDptr[count] & mask; # dann bitcount Bits kopieren
  315.             clear_loop_up(ptr,mantlen-count-1); # Rest mit Nullen füllen
  316.             RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  317.             encode_LF(sign,exp,MSDptr,mantlen, return);
  318.           }}
  319.           auf: # aufrunden
  320.           { SAVE_NUM_STACK # num_stack retten
  321.             var reg8 uintD* MSDptr;
  322.             num_stack_need(mantlen, MSDptr=,_EMA_);
  323.            {var reg2 uintD* ptr =
  324.               copy_loop_up(mantMSDptr,MSDptr,count); # count ganze Digits kopieren
  325.             if ((ptr[0] = ((mantptr[0] & mask) - mask)) == 0) # dann bitcount Bits kopieren und incrementieren
  326.               { if (!( inc_loop_down(ptr,count) ==0)) # evtl. weiterincrementieren
  327.                   { MSDptr[0] = bit(intDsize-1); exp = exp+1; } # evtl. Exponenten erhöhen
  328.               }
  329.             clear_loop_up(&ptr[1],mantlen-count-1); # Rest mit Nullen füllen
  330.             RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  331.             encode_LF(sign,exp,MSDptr,mantlen, return);
  332.           }}
  333.     }   }
  334. #else
  335.  
  336.   local object LF_fround_LF(x)
  337.     var reg1 object x;
  338.     { var reg4 uintC len = TheLfloat(x)->len;
  339.       var reg5 uintL uexp = TheLfloat(x)->expo;
  340.       if (uexp < LF_exp_mid)
  341.         { if (uexp == 0) { return x; } # x=0.0 -> Ergebnis 0.0
  342.           encode_LF0(len, return); # e<0 -> Ergebnis 0.0
  343.         }
  344.      {var reg7 uintL exp = uexp - LF_exp_mid;
  345.       if (exp >= intDsize*(uintL)len) # e>=16n -> x als Ergebnis
  346.         { return x; }
  347.       # 0 <= e < 16n
  348.       # alle hinteren 16n-e Bits wegrunden:
  349.       {var reg6 uintC count = floor(exp,intDsize); # zu kopierende Digits, < mantlen
  350.        var reg10 uintC bitcount = exp % intDsize; # zu kopierende Bits danach, >=0, <intDsize
  351.        var reg3 uintD mask = minus_bit(intDsize-bitcount-1); # Maske mit bitcount+1 Bits
  352.        {var reg2 uintD* mantptr = &TheLfloat(x)->data[count];
  353. #if !((__GNUC__ == 2) && (__GNUC_MINOR__ == 7))
  354.         if ((mantptr[0] & -mask) ==0) goto ab; # Bit 16n-e-1 =0 -> abrunden
  355. #else
  356.         if ((mantptr[0] & ((~(mask))+1)) ==0) goto ab; # Bit 16n-e-1 =0 -> abrunden
  357. #endif
  358.         if (!((mantptr[0] & ~mask) ==0)) goto auf; # Bit 16n-e-1 =1 und Bits 16n-e-2..0 >0 -> aufrunden
  359.         if (test_loop_up(&mantptr[1],len-count-1)) goto auf;
  360.         # round-to-even, je nach Bit 16n-e :
  361.         if (bitcount>0)
  362.           { if ((mantptr[0] & (-2*mask)) ==0) goto ab; else goto auf; }
  363.           elif (count>0)
  364.             { if ((mantptr[-1] & bit(0)) ==0) goto ab; else goto auf; }
  365.             else
  366.             # bitcount=0, count=0, also exp=0: Abrunden von +-0.5 zu 0.0
  367.             { encode_LF0(len, return); }
  368.        }
  369.        ab: # abrunden
  370.          pushSTACK(x);
  371.          {var reg8 object y = allocate_lfloat(len,uexp,R_sign(x)); # neues Long-Float
  372.           x = popSTACK();
  373.           # y_mant := NUDS mit e Bits aus x_mant und 16n-e Nullbits:
  374.           {var reg9 uintD* x_mantMSDptr = &TheLfloat(x)->data[0];
  375.            var reg2 uintD* ptr =
  376.              copy_loop_up(x_mantMSDptr,&TheLfloat(y)->data[0],count); # count ganze Digits kopieren
  377.            *ptr++ = x_mantMSDptr[count] & mask; # dann bitcount Bits kopieren
  378.            clear_loop_up(ptr,len-count-1); # Rest mit Nullen füllen
  379.           }
  380.           return y;
  381.          }
  382.        auf: # aufrunden
  383.          pushSTACK(x);
  384.          {var reg8 object y = allocate_lfloat(len,uexp,R_sign(x)); # neues Long-Float
  385.           x = popSTACK();
  386.           # y_mant := NUDS mit e Bits aus x_mant mit Increment und 16n-e Nullbits:
  387.           {var reg9 uintD* x_mantMSDptr = &TheLfloat(x)->data[0];
  388.            var reg9 uintD* y_mantMSDptr = &TheLfloat(y)->data[0];
  389.            var reg2 uintD* ptr =
  390.              copy_loop_up(x_mantMSDptr,y_mantMSDptr,count); # count ganze Digits kopieren
  391.            if ((ptr[0] = ((x_mantMSDptr[count] & mask) - mask)) == 0) # dann bitcount Bits kopieren und incrementieren
  392.              { if (!( inc_loop_down(ptr,count) ==0)) # evtl. weiterincrementieren
  393.                  { y_mantMSDptr[0] = bit(intDsize-1); (TheLfloat(y)->expo)++; } # evtl. Exponenten erhöhen
  394.              }
  395.            clear_loop_up(&ptr[1],len-count-1); # Rest mit Nullen füllen
  396.           }
  397.           return y;
  398.          }
  399.     }}}
  400. #endif
  401. # Liefert zu einem Long-Float x : (- x), ein LF.
  402. # LF_minus_LF(x)
  403. # kann GC auslösen
  404.   local object LF_minus_LF (object x);
  405. # Methode:
  406. # Falls x=0.0, fertig. Sonst Vorzeichenbit umdrehen und Pointer beibehalten.
  407.   local object LF_minus_LF(x)
  408.     var reg1 object x;
  409.     { if (TheLfloat(x)->expo == 0)
  410.         { return x; }
  411.         else
  412.         #ifdef SPVW_MIXED
  413.         { return as_object(as_oint(x) ^ wbit(vorz_bit_o)); }
  414.         #else
  415.         { var reg3 uintC len = TheLfloat(x)->len;
  416.           pushSTACK(x);
  417.          {var reg2 object mx = allocate_lfloat(len,TheLfloat(x)->expo,~R_sign(x));
  418.           x = popSTACK();
  419.           copy_loop_up(&TheLfloat(x)->data[0],&TheLfloat(mx)->data[0],len);
  420.           return mx;
  421.         }}
  422.         #endif
  423.     }
  424.  
  425. # LF_LF_comp(x,y) vergleicht zwei Long-Floats x und y.
  426. # Ergebnis: 0 falls x=y, +1 falls x>y, -1 falls x<y.
  427.   local signean LF_LF_comp (object x, object y);
  428. # Methode:
  429. # x und y haben verschiedenes Vorzeichen ->
  430. #    x < 0 -> x < y
  431. #    x >= 0 -> x > y
  432. # x und y haben gleiches Vorzeichen ->
  433. #    x >=0 -> vergleiche x und y (die rechten 24 Bits)
  434. #    x <0 -> vergleiche y und x (die rechten 24 Bits)
  435.   local signean LF_LF_comp(x,y)
  436.     var reg1 object x;
  437.     var reg2 object y;
  438.     { if (!R_minusp(y))
  439.         # y>=0
  440.         { if (!R_minusp(x))
  441.             # y>=0, x>=0
  442.             { # Vergleiche Exponenten und Mantissen:
  443.               { var reg3 uintL x_uexp = TheLfloat(x)->expo;
  444.                 var reg4 uintL y_uexp = TheLfloat(y)->expo;
  445.                 if (x_uexp < y_uexp) return signean_minus; # x<y
  446.                 if (x_uexp > y_uexp) return signean_plus; # x>y
  447.               }
  448.               { var reg5 uintC x_len = TheLfloat(x)->len;
  449.                 var reg6 uintC y_len = TheLfloat(y)->len;
  450.                 var reg4 uintC len = (x_len<y_len ? x_len : y_len); # min(x_len,y_len)
  451.                 # len Digits vergleichen:
  452.                 var reg3 signean erg =
  453.                   compare_loop_up(&TheLfloat(x)->data[0],&TheLfloat(y)->data[0],len);
  454.                 if (!(erg==0)) { return erg; } # verschieden -> fertig
  455.                 # gemeinsames Teilstück war gleich
  456.                 if (x_len == y_len) { return signean_null; } # gleiche Länge -> fertig
  457.                 if (x_len > y_len)
  458.                   # x länger als y
  459.                   { if (test_loop_up(&TheLfloat(x)->data[y_len],x_len-y_len))
  460.                       { return signean_plus; } # x>y
  461.                       else
  462.                       { return signean_null; }
  463.                   }
  464.                   else
  465.                   # y länger als x
  466.                   { if (test_loop_up(&TheLfloat(y)->data[x_len],y_len-x_len))
  467.                       { return signean_minus; } # x<y
  468.                       else
  469.                       { return signean_null; }
  470.                   }
  471.             } }
  472.             else
  473.             # y>=0, x<0
  474.             { return signean_minus; } # x<y
  475.         }
  476.         else
  477.         { if (!R_minusp(x))
  478.             # y<0, x>=0
  479.             { return signean_plus; } # x>y
  480.             else
  481.             # y<0, x<0
  482.             { # Vergleiche Exponenten und Mantissen:
  483.               { var reg3 uintL x_uexp = TheLfloat(x)->expo;
  484.                 var reg4 uintL y_uexp = TheLfloat(y)->expo;
  485.                 if (x_uexp < y_uexp) return signean_plus; # |x|<|y| -> x>y
  486.                 if (x_uexp > y_uexp) return signean_minus; # |x|>|y| -> x<y
  487.               }
  488.               { var reg5 uintC x_len = TheLfloat(x)->len;
  489.                 var reg6 uintC y_len = TheLfloat(y)->len;
  490.                 var reg4 uintC len = (x_len<y_len ? x_len : y_len); # min(x_len,y_len)
  491.                 # len Digits vergleichen:
  492.                 var reg3 signean erg =
  493.                   compare_loop_up(&TheLfloat(y)->data[0],&TheLfloat(x)->data[0],len);
  494.                 if (!(erg==0)) { return erg; } # verschieden -> fertig
  495.                 # gemeinsames Teilstück war gleich
  496.                 if (x_len == y_len) { return signean_null; } # gleiche Länge -> fertig
  497.                 if (x_len > y_len)
  498.                   # x länger als y
  499.                   { if (test_loop_up(&TheLfloat(x)->data[y_len],x_len-y_len))
  500.                       { return signean_minus; } # |x|>|y| -> x<y
  501.                       else
  502.                       { return signean_null; }
  503.                   }
  504.                   else
  505.                   # y länger als x
  506.                   { if (test_loop_up(&TheLfloat(y)->data[x_len],y_len-x_len))
  507.                       { return signean_plus; } # |x|<|y| -> x>y
  508.                       else
  509.                       { return signean_null; }
  510.                   }
  511.             } }
  512.         }
  513.     }
  514.  
  515. # LF_shorten_LF(x,len) verkürzt ein Long-Float x auf gegebene Länge len
  516. # und rundet dabei.
  517. # > object x: ein Long-FLoat
  518. # > uintC len: gewünschte Länge (>= LF_minlen, < TheLfloat(x)->len)
  519. # < object ergebnis: verkürztes Long-Float
  520. # kann GC auslösen
  521.   local object LF_shorten_LF (object x, uintC len);
  522.   local object LF_shorten_LF(x,len)
  523.     var reg2 object x;
  524.     var reg4 uintC len;
  525.     { # x = 0.0 braucht nicht abgefangen zu werden, da bei Mantisse 0 dann
  526.       # sowieso abgerundet wird, die Mantisse also 0 bleibt.
  527.       pushSTACK(x);
  528.      {var reg3 object y = allocate_lfloat(len,TheLfloat(x)->expo,R_sign(x)); # neues LF
  529.       x = popSTACK();
  530.       { var reg5 uintC oldlen = TheLfloat(x)->len; # alte Länge, > len
  531.         # Mantisse von x nach y kopieren:
  532.         copy_loop_up(&TheLfloat(x)->data[0],&TheLfloat(y)->data[0],len);
  533.         # Entscheiden, ob auf- oder abrunden:
  534.        {var reg1 uintD* ptr = &TheLfloat(x)->data[len];
  535.         if ( ((sintD)ptr[0] >= 0) # nächstes Bit eine 0 -> abrunden
  536.              || ( ((ptr[0] & ((uintD)bit(intDsize-1)-1)) ==0) # eine 1 und alles weitere Nullen?
  537.                   && !test_loop_up(&ptr[1],oldlen-len-1)
  538.                   # round-to-even
  539.                   && ((ptr[-1] & bit(0)) ==0)
  540.            )    )
  541.           # abrunden
  542.           {}
  543.           else
  544.           # aufrunden
  545.           { if ( inc_loop_down(&TheLfloat(y)->data[len],len) )
  546.               # Übertrag durch Aufrunden
  547.               { TheLfloat(y)->data[0] = bit(intDsize-1); # Mantisse := 10...0
  548.                 # Exponent erhöhen:
  549.                 if (++(TheLfloat(y)->expo) == LF_exp_high+1) { fehler_overflow(); }
  550.           }   }
  551.       }}
  552.       return y;
  553.     }}
  554.  
  555. # LF_extend_LF(x,len) verlängert ein Long-Float x auf gegebene Länge len.
  556. # > object x: ein Long-FLoat
  557. # > uintC len: gewünschte Länge (> TheLfloat(x)->len)
  558. # < object ergebnis: verlängertes Long-Float
  559. # kann GC auslösen
  560.   local object LF_extend_LF (object x, uintC len);
  561.   local object LF_extend_LF(x,len)
  562.     var reg1 object x;
  563.     var reg4 uintC len;
  564.     { pushSTACK(x);
  565.      {var reg2 object y = allocate_lfloat(len,TheLfloat(x)->expo,R_sign(x)); # neues LF
  566.       x = popSTACK();
  567.       { var reg5 uintC oldlen = TheLfloat(x)->len; # alte Länge, < len
  568.         # Mantisse von x nach y kopieren:
  569.         var reg3 uintD* ptr =
  570.           copy_loop_up(&TheLfloat(x)->data[0],&TheLfloat(y)->data[0],oldlen);
  571.         # und mit Null-Digits ergänzen:
  572.         clear_loop_up(ptr,len-oldlen);
  573.       }
  574.       return y;
  575.     }}
  576.  
  577. # LF_to_LF(x,len) wandelt ein Long-Float x in ein Long-Float gegebener Länge
  578. # len um und rundet dabei nötigenfalls.
  579. # > object x: ein Long-FLoat
  580. # > uintC len: gewünschte Länge (>= LF_minlen)
  581. # < object ergebnis: Long-Float gegebener Länge
  582. # kann GC auslösen
  583.   local object LF_to_LF (object x, uintC len);
  584.   local object LF_to_LF(x,len)
  585.     var reg2 object x;
  586.     var reg3 uintC len;
  587.     { var reg1 uintC oldlen = TheLfloat(x)->len;
  588.       if (len < oldlen) { return LF_shorten_LF(x,len); }
  589.       if (len > oldlen) { return LF_extend_LF(x,len); }
  590.       # len = oldlen
  591.       return x;
  592.     }
  593.  
  594. # Liefert zu zwei gleichlangen Long-Float x und y : (+ x y), ein LF.
  595. # LF_LF_plus_LF(x,y)
  596. # kann GC auslösen
  597.   local object LF_LF_plus_LF (object x, object y);
  598. # Methode (nach [Knuth, II, Seminumerical Algorithms, Abschnitt 4.2.1., S.200]):
  599. # Falls e1<e2, vertausche x1 und x2.
  600. # Also e1 >= e2.
  601. # Falls e2=0, also x2=0.0, Ergebnis x1.
  602. # Falls e1 - e2 >= 16n+2, Ergebnis x1.
  603. # Erweitere die Mantissen rechts um 3 Bits (Bit -1 als Schutzbit, Bits -2,-3
  604. #   als Rundungsbits: 00 exakt, 01 1.Hälfte, 10 exakte Mitte, 11 2.Hälfte.)
  605. # Schiebe die Mantisse von x2 um e0-e1 Bits nach rechts. (Dabei die Rundung
  606. # ausführen: Bit -3 ist das logische Oder der Bits -3,-4,-5,...)
  607. # Falls x1,x2 selbes Vorzeichen haben: Addiere dieses zur Mantisse von x1.
  608. # Falls x1,x2 verschiedenes Vorzeichen haben: Subtrahiere dieses von der
  609. #   Mantisse von x1. <0 -> (Es war e1=e2) Vertausche die Vorzeichen, negiere.
  610. #                    =0 -> Ergebnis 0.0
  611. # Exponent ist e1.
  612. # Normalisiere, fertig.
  613.   local object LF_LF_plus_LF(x1,x2)
  614.     var reg9 object x1;
  615.     var reg9 object x2;
  616.     { var reg9 uintL uexp1 = TheLfloat(x1)->expo;
  617.       var reg9 uintL uexp2 = TheLfloat(x2)->expo;
  618.       if (uexp1 < uexp2)
  619.         # x1 und x2 vertauschen
  620.         { swap(reg1 object, x1,x2); swap(reg1 uintL, uexp1,uexp2); }
  621.       # uexp1 >= uexp2
  622.       if (uexp2==0) { return x1; } # x2=0.0 -> x1 als Ergebnis
  623.      {var reg9 uintC len = TheLfloat(x1)->len; # Länge n von x1 und x2
  624.       var reg9 uintL expdiff = uexp1-uexp2; # e1-e2
  625.       #if !defined(SPVW_MIXED)
  626.       if ((expdiff == 0) && !same_sign_p(x1,x2))
  627.         # verschiedene Vorzeichen, aber gleicher Exponent
  628.         { # Vorzeichen des Ergebnisses festlegen:
  629.           var reg1 signean erg = # Mantissen (je len Digits) vergleichen
  630.             compare_loop_up(&TheLfloat(x1)->data[0],&TheLfloat(x2)->data[0],len);
  631.           if (erg==0) # Mantissen gleich
  632.             { encode_LF0(len, return); } # Ergebnis 0.0
  633.           if (erg<0) # |x1| < |x2|
  634.             # x1 und x2 vertauschen, expdiff bleibt =0
  635.             { swap(reg1 object, x1,x2); swap(reg1 uintL, uexp1,uexp2); }
  636.         }
  637.       #endif
  638.       if (expdiff >= intDsize * (uintL)len + 2) # e1-e2 >= 16n+2 ?
  639.         { return x1; } # ja -> x1 als Ergebnis
  640.       # neues Long-Float allozieren:
  641.       pushSTACK(x1); pushSTACK(x2);
  642.       { var reg7 object y = allocate_lfloat(len,uexp1,R_sign(x1));
  643.         x2 = popSTACK(); x1 = popSTACK();
  644.        {var reg9 uintL i = floor(expdiff,intDsize); # e1-e2 div 16 (>=0, <=n)
  645.         var reg9 uintL j = expdiff % intDsize; # e1-e2 mod 16 (>=0, <16)
  646.         # Mantisse von x2 muß um intDsize*i+j Bits nach rechts geschoben werden.
  647.         var reg8 uintC x2_len = len - i; # n-i Digits von x2 gebraucht
  648.         # x2_len Digits um j Bits nach rechts schieben und dabei kopieren:
  649.         SAVE_NUM_STACK # num_stack retten
  650.         var reg9 uintD* x2_MSDptr;
  651.         var reg9 uintD* x2_LSDptr;
  652.         var reg2 uintD rounding_bits;
  653.         num_stack_need(x2_len, x2_MSDptr=,x2_LSDptr=); # x2_len Digits Platz
  654.         if (j==0)
  655.           { copy_loop_up(&TheLfloat(x2)->data[0],x2_MSDptr,x2_len); rounding_bits = 0; }
  656.           else
  657.           { rounding_bits = shiftrightcopy_loop_up(&TheLfloat(x2)->data[0],x2_MSDptr,x2_len,j,0); }
  658.         # x2_MSDptr/x2_len/x2_LSDptr sind die essentiellen Digits von x2.
  659.         # rounding_bits enthält die letzten j herausgeschobenen Bits.
  660.         # Aus rounding_bits und den nächsten i Digits die 3 Rundungsbits
  661.         # (als Bits intDsize-1..intDsize-3 von rounding_bits) aufbauen:
  662.         if (j>=2)
  663.           # j>=2 -> Bits -1,-2 sind OK, Bit -3 bestimmen:
  664.           { if ((rounding_bits & (bit(intDsize-3)-1)) ==0)
  665.               { if (test_loop_up(&TheLfloat(x2)->data[x2_len],i))
  666.                   { rounding_bits |= bit(intDsize-3); } # Rundungsbit -3 setzen
  667.               }
  668.               else
  669.               { rounding_bits |= bit(intDsize-3); # Rundungsbit -3 setzen
  670.                 rounding_bits &= bitm(intDsize)-bit(intDsize-3); # andere Bits löschen
  671.           }   }
  672.           else
  673.           # j<=3 -> Bits intDsize-4..0 von rounding_bits sind bereits Null.
  674.           # nächstes und weitere i-1 Digits heranziehen:
  675.           { if (i > 0) # i=0 -> Bits -1,-2,-3 sind OK.
  676.               { var reg1 uintD* ptr = &TheLfloat(x2)->data[x2_len];
  677.                 rounding_bits |= (ptr[0] >> j); # weitere relevante Bits des nächsten Digit dazu
  678.                 if ((rounding_bits & (bit(intDsize-3)-1)) ==0) # Alle Bits -3,-4,... =0 ?
  679.                   { if (   (!((ptr[0] & (bit(3)-1)) ==0)) # j (<=3) untere Bits von ptr[0] alle =0 ?
  680.                         || test_loop_up(&ptr[1],i-1)
  681.                        )
  682.                       { rounding_bits |= bit(intDsize-3); } # Rundungsbit -3 setzen
  683.                   }
  684.                   else
  685.                   { rounding_bits |= bit(intDsize-3); # Rundungsbit -3 setzen
  686.                     rounding_bits &= bitm(intDsize)-bit(intDsize-3); # andere Bits löschen
  687.           }   }   }
  688.         # x2 liegt in verschobener Form in der UDS x2_MSDptr/x2_len/x2_LSDptr
  689.         # vor, mit Rundungsbits in Bit intDsize-1..intDsize-3 von rounding_bits.
  690.         {var reg3 uintD* y_mantMSDptr = &TheLfloat(y)->data[0];
  691.          var reg4 uintD* y_mantLSDptr = &y_mantMSDptr[(uintP)len];
  692.          if (same_sign_p(x1,x2))
  693.            # gleiche Vorzeichen -> Mantissen addieren
  694.            { # erst rechten Mantissenteil (x2_len Digits) durch Addition:
  695.              var reg5 uintD carry =
  696.                add_loop_down(&TheLfloat(x1)->data[(uintP)len],x2_LSDptr,
  697.                              y_mantLSDptr, x2_len
  698.                             );
  699.              # dann linken Mantissenteil (i Digits) direkt kopieren:
  700.              var reg1 uintD* ptr =
  701.                copy_loop_up(&TheLfloat(x1)->data[0],y_mantMSDptr,i);
  702.              # dann Übertrag vom rechten zum linken Mantissenteil addieren:
  703.              if (!(carry==0))
  704.                { if ( inc_loop_down(ptr,i) )
  705.                    # Übertrag über das erste Digit hinaus
  706.                    { # Exponent von y incrementieren:
  707.                      if ( ++(TheLfloat(y)->expo) == LF_exp_high+1 ) { fehler_overflow(); }
  708.                      # normalisiere durch Schieben um 1 Bit nach rechts:
  709.                     {var reg6 uintD carry_rechts =
  710.                        shift1right_loop_up(y_mantMSDptr,len,-1);
  711.                      rounding_bits = rounding_bits>>1; # Rundungsbits mitschieben
  712.                      if (!(carry_rechts==0)) { rounding_bits |= bit(intDsize-1); }
  713.                }   }}
  714.            }
  715.            else
  716.            # verschiedene Vorzeichen -> Mantissen subtrahieren
  717.            { # erst rechten Mantissenteil (x2_len Digits) durch Subtraktion:
  718.              rounding_bits = -rounding_bits;
  719.              {var reg5 uintD carry =
  720.                 subx_loop_down(&TheLfloat(x1)->data[(uintP)len],x2_LSDptr,
  721.                                y_mantLSDptr, x2_len,
  722.                                (rounding_bits==0 ? 0 : -1L)
  723.                               );
  724.               # dann linken Mantissenteil (i Digits) direkt kopieren:
  725.               var reg1 uintD* ptr =
  726.                 copy_loop_up(&TheLfloat(x1)->data[0],y_mantMSDptr,i);
  727.               # dann Übertrag des rechten vom linken Mantissenteil subtrahieren:
  728.               if (!(carry==0))
  729.                 { if ( dec_loop_down(ptr,i) )
  730.                     # Übertrag über das erste Digit hinaus, also e1=e2
  731.                     #if !defined(SPVW_MIXED)
  732.                     { NOTREACHED } # diesen Fall haben wir schon behandelt
  733.                     #else
  734.                     { # Negieren:
  735.                       y = as_object(as_oint(y) ^ wbit(vorz_bit_o));
  736.                       rounding_bits = -rounding_bits;
  737.                       if (rounding_bits==0)
  738.                         # Negieren ohne Carry
  739.                         { neg_loop_down(y_mantLSDptr,len); }
  740.                         else
  741.                         # Negieren mit Carry von rechts
  742.                         { # not_loop_down(y_mantLSDptr,len); # oder
  743.                           not_loop_up(y_mantMSDptr,len);
  744.                         }
  745.                     }
  746.                     #endif
  747.                 }
  748.              }
  749.              # UDS y_mantMSDptr/len/y_mantLSDptr/rounding_bits normalisieren:
  750.              {var reg1 uintD* ptr = y_mantMSDptr;
  751.               var reg5 uintL k = 0;
  752.               var reg6 uintC count;
  753.               dotimesC(count,len,
  754.                 { if (!(ptr[0]==0)) goto nonzero_found;
  755.                   ptr++; k++;
  756.                 });
  757.               if (!(rounding_bits==0)) goto nonzero_found;
  758.               # Die UDS ist ganz Null. Also war e1=e2, keine Rundungsbits.
  759.               #if !defined(SPVW_MIXED)
  760.               { NOTREACHED } # diesen Fall haben wir schon behandelt
  761.               #else
  762.               TheLfloat(y)->expo = 0; # 0.0 als Ergebnis
  763.               return as_object(as_oint(y) & ~wbit(vorz_bit_o));
  764.               #endif
  765.               nonzero_found: # Digit /=0 gefunden
  766.               # UDS von ptr nach y_mantMSDptr um k Digits nach unten kopieren:
  767.               if (k>0)
  768.                 # mindestens ein führendes Nulldigit. Also war e1-e2 = 0 oder 1.
  769.                 { ptr = copy_loop_up(ptr,y_mantMSDptr,len-k); # len-k Digits verschieben
  770.                   *ptr++ = rounding_bits; # Rundungsbits als weiteres Digit
  771.                   clear_loop_up(ptr,k-1); # dann k-1 Nulldigits
  772.                   rounding_bits = 0; # und keine weiteren Rundungsbits
  773.                   # Exponenten um intDsize*k erniedrigen:
  774.                   k = intDsize*k;
  775.                  {var reg1 uintL uexp = TheLfloat(y)->expo;
  776.                   #if !(LF_exp_low==1)
  777.                   if (uexp < k+LF_exp_low)
  778.                   #else
  779.                   if (uexp <= k)
  780.                   #endif
  781.                     { if (underflow_allowed())
  782.                         { fehler_underflow(); }
  783.                         else
  784.                         { encode_LF0(len, return); } # Ergebnis 0.0
  785.                     }
  786.                   TheLfloat(y)->expo = uexp - k;
  787.                 }}
  788.              }
  789.              # NUDS y_mantMSDptr/len/y_mantLSDptr/rounding_bits normalisieren:
  790.              {var reg5 uintL s;
  791.               integerlengthD(y_mantMSDptr[0], s = intDsize - );
  792.               # s = Anzahl der führenden Nullbits im ersten Word (>=0, <intDsize)
  793.               if (s > 0)
  794.                 { # Muß die NUDS y_mantMSDptr/len/y_mantLSDptr/rounding_bits
  795.                   # um s Bits nach links schieben.
  796.                   # (Bei e1-e2>1 ist dabei zwangsläufig s=1.)
  797.                   if (s==1)
  798.                     { shift1left_loop_down(y_mantLSDptr,len);
  799.                       if (rounding_bits & bit(intDsize-1))
  800.                         { y_mantLSDptr[-1] |= bit(0); }
  801.                       rounding_bits = rounding_bits << 1;
  802.                     }
  803.                     else # s>1, also e1-e2 <= 1 <= s.
  804.                     { shiftleft_loop_down(y_mantLSDptr,len,s,rounding_bits>>(intDsize-s));
  805.                       rounding_bits = 0; # = rounding_bits << s;
  806.                     }
  807.                   # Exponenten um s erniedrigen:
  808.                  {var reg1 uintL uexp = TheLfloat(y)->expo;
  809.                   #if !(LF_exp_low==1)
  810.                   if (uexp < s+LF_exp_low)
  811.                   #else
  812.                   if (uexp <= s)
  813.                   #endif
  814.                     { if (underflow_allowed())
  815.                         { fehler_underflow(); }
  816.                         else
  817.                         { encode_LF0(len, return); } # Ergebnis 0.0
  818.                     }
  819.                   TheLfloat(y)->expo = uexp - s;
  820.                 }}
  821.            } }
  822.          # Hier enthält rounding_bits Bit -1 als Bit intDsize-1, Bit -2 als
  823.          # Bit intDsize-2, Bit -3 als Oder(Bits intDsize-3..0) !
  824.          # Runden. Dazu rounding_bits inspizieren:
  825.          if ((rounding_bits & bit(intDsize-1)) ==0) goto ab; # Bit -1 gelöscht -> abrunden
  826.          rounding_bits = rounding_bits<<1; # Bits -2,-3
  827.          if (!(rounding_bits==0)) goto auf; # Bit -2 oder Bit -3 gesetzt -> aufrunden
  828.          # round-to-even:
  829.          if ((y_mantLSDptr[-1] & bit(0)) ==0) goto ab;
  830.          auf: # aufrunden
  831.            if ( inc_loop_down(y_mantLSDptr,len) )
  832.              { # Übertrag durchs Aufrunden
  833.                y_mantMSDptr[0] = bit(intDsize-1); # Mantisse := 10...0
  834.                # Exponent erhöhen:
  835.                if (++(TheLfloat(y)->expo) == LF_exp_high+1) { fehler_overflow(); }
  836.              }
  837.          ab: # abrunden
  838.            ;
  839.         }
  840.         RESTORE_NUM_STACK # num_stack zurück
  841.         # y fertig.
  842.         return y;
  843.     }}}}
  844.  
  845. # Liefert zu zwei gleichlangen Long-Float x und y : (- x y), ein LF.
  846. # LF_LF_minus_LF(x,y)
  847. # kann GC auslösen
  848.   local object LF_LF_minus_LF (object x, object y);
  849. # Methode:
  850. # (- x1 x2) = (+ x1 (- x2))
  851.   local object LF_LF_minus_LF(x1,x2)
  852.     var reg2 object x1;
  853.     var reg1 object x2;
  854.     { if (TheLfloat(x2)->expo == 0)
  855.         { return x1; }
  856.         else
  857.         #ifdef SPVW_MIXED
  858.         { return LF_LF_plus_LF(x1, as_object(as_oint(x2) ^ wbit(vorz_bit_o)) ); }
  859.         #else
  860.         { var reg4 uintC len2 = TheLfloat(x2)->len;
  861.           pushSTACK(x1); pushSTACK(x2);
  862.          {var reg3 object mx2 = allocate_lfloat(len2,TheLfloat(x2)->expo,~R_sign(x2));
  863.           x2 = popSTACK();
  864.           copy_loop_up(&TheLfloat(x2)->data[0],&TheLfloat(mx2)->data[0],len2);
  865.           return LF_LF_plus_LF(popSTACK(),mx2);
  866.         }}
  867.         #endif
  868.     }
  869.  
  870. # Liefert zu zwei gleichlangen Long-Float x und y : (* x y), ein LF.
  871. # LF_LF_mal_LF(x,y)
  872. # kann GC auslösen
  873.   local object LF_LF_mal_LF (object x, object y);
  874. # Methode:
  875. # Falls x1=0.0 oder x2=0.0 -> Ergebnis 0.0
  876. # Sonst: Ergebnis-Vorzeichen = VZ von x1 xor VZ von x2.
  877. #        Ergebnis-Exponent = Summe der Exponenten von x1 und x2.
  878. #        Produkt der Mantissen bilden (2n Digits).
  879. #        Falls das führende Bit =0 ist: Mantissenprodukt um 1 Bit nach links
  880. #          schieben (die vorderen n+1 Digits genügen)
  881. #          und Exponent decrementieren.
  882. #        Runden auf n Digits liefert die Ergebnis-Mantisse.
  883.   local object LF_LF_mal_LF(x1,x2)
  884.     var reg9 object x1;
  885.     var reg9 object x2;
  886.     { var reg6 uintL uexp1 = TheLfloat(x1)->expo;
  887.       if (uexp1==0) { return x1; } # x1=0.0 -> Ergebnis 0.0
  888.      {var reg8 uintL uexp2 = TheLfloat(x2)->expo;
  889.       if (uexp2==0) { return x2; } # x2=0.0 -> Ergebnis 0.0
  890.       # Exponenten addieren:
  891.       # (uexp1-LF_exp_mid) + (uexp2-LF_exp_mid) = (uexp1+uexp2-LF_exp_mid)-LF_exp_mid
  892.       uexp1 = uexp1 + uexp2;
  893.       if (uexp1 >= uexp2)
  894.         # kein Carry
  895.         { if (uexp1 < LF_exp_mid+LF_exp_low)
  896.             { if (underflow_allowed())
  897.                 { fehler_underflow(); }
  898.                 else
  899.                 { encode_LF0(TheLfloat(x1)->len, return); } # Ergebnis 0.0
  900.         }   }
  901.         else
  902.         # Carry
  903.         { if (uexp1 > (uintL)(LF_exp_mid+LF_exp_high+1)) { fehler_overflow(); } }
  904.       uexp1 = uexp1 - LF_exp_mid;
  905.       # Nun ist LF_exp_low <= uexp1 <= LF_exp_high+1.
  906.       # neues Long-Float allozieren:
  907.       pushSTACK(x1); pushSTACK(x2);
  908.       {var reg2 uintC len = TheLfloat(x1)->len; # Länge n von x1 und x2
  909.        var reg1 object y = allocate_lfloat(len,uexp1,
  910.                                            R_sign(as_object(as_oint(x1) ^ as_oint(x2))) # Vorzeichen kombinieren
  911.                                           );
  912.        x2 = popSTACK(); x1 = popSTACK();
  913.        # Produkt bilden:
  914.        {var reg4 uintD* MSDptr;
  915.         {SAVE_NUM_STACK # num_stack retten
  916.          UDS_UDS_mal_UDS(len,&TheLfloat(x1)->data[(uintP)len],
  917.                          len,&TheLfloat(x2)->data[(uintP)len],
  918.                          MSDptr=,_EMA_,_EMA_);
  919.          RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  920.         }
  921.         {var reg3 uintD* midptr = &MSDptr[(uintP)len]; # Pointer in die Mitte der 2n Digits
  922.          if ((sintD)MSDptr[0] >= 0) # führendes Bit abtesten
  923.            { # erste n+1 Digits um 1 Bit nach links schieben:
  924.              shift1left_loop_down(&midptr[1],len+1);
  925.              # Exponenten decrementieren:
  926.              if ((TheLfloat(y)->expo)-- == LF_exp_low-1)
  927.                { if (underflow_allowed())
  928.                    { fehler_underflow(); }
  929.                    else
  930.                    { encode_LF0(len, return); } # Ergebnis 0.0
  931.                }
  932.            }
  933.          # erste Hälfte des Mantissenprodukts übertragen:
  934.          {var reg5 uintD* y_mantMSDptr = &TheLfloat(y)->data[0];
  935.           var reg7 uintD* y_mantLSDptr =
  936.             copy_loop_up(MSDptr,y_mantMSDptr,len);
  937.           # Runden:
  938.           if ( ((sintD)midptr[0] >= 0) # nächstes Bit =0 -> abrunden
  939.                || ( ((midptr[0] & ((uintD)bit(intDsize-1)-1)) ==0) # Bit =1, weitere Bits >0 -> aufrunden
  940.                     && !test_loop_up(&midptr[1],len-1)
  941.                     # round-to-even
  942.                     && ((midptr[-1] & bit(0)) ==0)
  943.              )    )
  944.             # abrunden
  945.             {}
  946.             else
  947.             # aufrunden
  948.             { if ( inc_loop_down(y_mantLSDptr,len) )
  949.                 { # Übertrag durchs Aufrunden (kann nur auftreten,
  950.                   # wenn vorhin um 1 Bit nach links geschoben wurde)
  951.                   y_mantMSDptr[0] = bit(intDsize-1); # Mantisse := 10...0
  952.                   (TheLfloat(y)->expo)++; # Exponent wieder zurück-erhöhen
  953.             }   }
  954.           # LF_exp_low <= exp <= LF_exp_high sicherstellen:
  955.           if (TheLfloat(y)->expo == LF_exp_high+1) { fehler_overflow(); }
  956.        }}}
  957.        return y;
  958.     }}}
  959.  
  960. # Liefert zu zwei gleichlangen Long-Float x und y : (/ x y), ein LF.
  961. # LF_LF_durch_LF(x,y)
  962. # kann GC auslösen
  963.   local object LF_LF_durch_LF (object x, object y);
  964. # Methode:
  965. # x2 = 0.0 -> Error
  966. # x1 = 0.0 -> Ergebnis 0.0
  967. # Sonst:
  968. # Ergebnis-Vorzeichen = xor der beiden Vorzeichen von x1 und x2
  969. # Ergebnis-Exponent = Differenz der beiden Exponenten von x1 und x2
  970. # Ergebnis-Mantisse = Mantisse mant1 / Mantisse mant2, gerundet.
  971. #   mant1/mant2 > 1/2, mant1/mant2 < 2;
  972. #   nach Rundung mant1/mant2 >=1/2, <=2*mant1<2.
  973. #   Bei mant1/mant2 >=1 brauche 16n-1 Nachkommabits,
  974. #   bei mant1/mant2 <1 brauche 16n Nachkommabits.
  975. #   Fürs Runden: brauche ein Rundungsbit (Rest gibt an, ob exakt).
  976. #   Brauche daher insgesamt 16n+1 Nachkommabits von mant1/mant2.
  977. #   Dividiere daher (als Unsigned Integers)
  978. #     2^16(n+1)*(2^16n*m0) durch (2^16n*m1).
  979. #   Falls der Quotient >=2^16(n+1) ist, schiebe ihn um 1 Bit nach rechts,
  980. #     erhöhe den Exponenten um 1 und runde das letzte Digit weg.
  981. #   Falls der Quotient <2^16(n+1) ist, runde das letzte Digit weg. Bei rounding
  982. #     overflow schiebe um 1 Bit nach rechts und erhöhe den Exponenten um 1.
  983.   # Workaround gcc-2.7.0 bug on i386.
  984.     #if defined(__GNUC__)
  985.       #if (__GNUC__ == 2)
  986.         #if (__GNUC_MINOR__ == 7)
  987.           #define workaround_gcc270_bug()  *&uexp1 = *&uexp1;
  988.         #endif
  989.       #endif
  990.     #endif
  991.     #ifndef workaround_gcc270_bug
  992.       #define workaround_gcc270_bug()
  993.     #endif
  994.   local object LF_LF_durch_LF(x1,x2)
  995.     var reg7 object x1;
  996.     var reg9 object x2;
  997.     { var reg8 uintL uexp2 = TheLfloat(x2)->expo;
  998.       if (uexp2==0) { divide_0(); } # x2=0.0 -> Error
  999.      {var reg6 uintL uexp1 = TheLfloat(x1)->expo;
  1000.       if (uexp1==0) { return x1; } # x1=0.0 -> Ergebnis 0.0
  1001.       # Exponenten subtrahieren:
  1002.       # (uexp1-LF_exp_mid) - (uexp2-LF_exp_mid) = (uexp1-uexp2+LF_exp_mid)-LF_exp_mid
  1003.       if (uexp1 >= uexp2)
  1004.         { uexp1 = uexp1 - uexp2; # kein Carry
  1005.           workaround_gcc270_bug();
  1006.           if (uexp1 > LF_exp_high-LF_exp_mid) { fehler_overflow(); }
  1007.           uexp1 = uexp1 + LF_exp_mid;
  1008.         }
  1009.         else
  1010.         { uexp1 = uexp1 - uexp2; # Carry
  1011.           workaround_gcc270_bug();
  1012.           if (uexp1 < (uintL)(LF_exp_low-1-LF_exp_mid))
  1013.             { if (underflow_allowed())
  1014.                 { fehler_underflow(); }
  1015.                 else
  1016.                 { encode_LF0(TheLfloat(x1)->len, return); } # Ergebnis 0.0
  1017.             }
  1018.           uexp1 = uexp1 + LF_exp_mid;
  1019.         }
  1020.       # Nun ist LF_exp_low-1 <= uexp1 <= LF_exp_high.
  1021.       # neues Long-Float allozieren:
  1022.       pushSTACK(x1); pushSTACK(x2);
  1023.       {var reg2 uintC len = TheLfloat(x1)->len; # Länge n von x1 und x2
  1024.        var reg1 object y = allocate_lfloat(len,uexp1,
  1025.                                            R_sign(as_object(as_oint(x1) ^ as_oint(x2))) # Vorzeichen kombinieren
  1026.                                           );
  1027.        x2 = popSTACK(); x1 = popSTACK();
  1028.        # Zähler bilden:
  1029.        {SAVE_NUM_STACK # num_stack retten
  1030.         var reg9 uintD* z_MSDptr;
  1031.         var reg9 uintL z_len;
  1032.         var reg9 uintD* z_LSDptr;
  1033.         z_len = 2*(uintL)len + 1;
  1034.         if ((intCsize < 32) && (z_len > (uintL)(bitc(intCsize)-1))) { fehler_LF_toolong(); }
  1035.         num_stack_need(z_len, z_MSDptr=,z_LSDptr=);
  1036.         {var reg3 uintD* ptr =
  1037.            copy_loop_up(&TheLfloat(x1)->data[0],z_MSDptr,len); # n Digits kopieren
  1038.          clear_loop_up(ptr,len+1); # und n+1 Null-Digits
  1039.         }
  1040.         # Quotienten bilden: 2n+1-Digit-Zahl durch n-Digit-Zahl dividieren
  1041.         {var DS q;
  1042.          var DS r;
  1043.          {var reg4 uintD* x2_mantMSDptr = &TheLfloat(x2)->data[0];
  1044.           UDS_divide(z_MSDptr,z_len,z_LSDptr,
  1045.                      x2_mantMSDptr,len,&x2_mantMSDptr[(uintP)len],
  1046.                      &q, &r
  1047.                     );
  1048.          }
  1049.          # q ist der Quotient mit n+1 oder n+2 Digits, r der Rest.
  1050.          RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  1051.          if (q.len > len+1)
  1052.            # Quotient hat n+2 Digits -> um 1 Bit nach rechts schieben:
  1053.            { var reg5 uintD* y_mantMSDptr = &TheLfloat(y)->data[0];
  1054.              var reg4 uintD carry_rechts =
  1055.                shiftrightcopy_loop_up(&q.MSDptr[1],y_mantMSDptr,len,1,
  1056.                                       /* carry links = q.MSDptr[0] = 1 */ 1 );
  1057.              # Exponenten incrementieren:
  1058.              if (++(TheLfloat(y)->expo) == LF_exp_high+1) { fehler_overflow(); }
  1059.              # Runden:
  1060.              if ( (carry_rechts == 0) # herausgeschobenes Bit =0 -> abrunden
  1061.                   || ( (q.LSDptr[-1]==0) # =1 und weitere Bits >0 oder Rest >0 -> aufrunden
  1062.                        && (r.len==0)
  1063.                        # round-to-even
  1064.                        && ((q.LSDptr[-2] & bit(1)) ==0)
  1065.                 )    )
  1066.                # abrunden
  1067.                {}
  1068.                else
  1069.                # aufrunden
  1070.                { inc_loop_down(&y_mantMSDptr[(uintP)len],len); }
  1071.            }
  1072.            else
  1073.            # Quotient hat n+1 Digits -> nur kopieren:
  1074.            { var reg4 uintD* y_mantMSDptr = &TheLfloat(y)->data[0];
  1075.              copy_loop_up(q.MSDptr,y_mantMSDptr,len);
  1076.              # Runden:
  1077.              if ( ((sintD)(q.LSDptr[-1]) >= 0) # nächstes Bit =0 -> abrunden
  1078.                   || ( ((q.LSDptr[-1] & ((uintD)bit(intDsize-1)-1)) ==0) # =1 und weitere Bits >0 oder Rest >0 -> aufrunden
  1079.                        && (r.len==0)
  1080.                        # round-to-even
  1081.                        && ((q.LSDptr[-2] & bit(0)) ==0)
  1082.                 )    )
  1083.                # abrunden
  1084.                {}
  1085.                else
  1086.                # aufrunden
  1087.                { if ( inc_loop_down(&y_mantMSDptr[(uintP)len],len) )
  1088.                    # Übertrag durchs Aufrunden
  1089.                    { y_mantMSDptr[0] = bit(intDsize-1); # Mantisse := 10...0
  1090.                      # Exponenten incrementieren:
  1091.                      if (++(TheLfloat(y)->expo) == LF_exp_high+1) { fehler_overflow(); }
  1092.                }   }
  1093.            }
  1094.         }
  1095.         # LF_exp_low <= exp <= LF_exp_high sicherstellen:
  1096.         if (TheLfloat(y)->expo == LF_exp_low-1)
  1097.           { if (underflow_allowed())
  1098.               { fehler_underflow(); }
  1099.               else
  1100.               { encode_LF0(len, return); } # Ergebnis 0.0
  1101.           }
  1102.        }
  1103.        return y;
  1104.     }}}
  1105.  
  1106. # Liefert zu einem Long-Float x>=0 : (sqrt x), ein LF.
  1107. # LF_sqrt_LF(x)
  1108. # kann GC auslösen
  1109.   local object LF_sqrt_LF (object x);
  1110. # Methode:
  1111. # x = 0.0 -> Ergebnis 0.0
  1112. # Ergebnis-Vorzeichen := positiv,
  1113. # Ergebnis-Exponent := ceiling(e/2),
  1114. # Ergebnis-Mantisse:
  1115. #   Erweitere die Mantisse (n Digits) um n+2 Nulldigits nach hinten.
  1116. #   Bei ungeradem e schiebe dies (oder nur die ersten n+1 Digits davon)
  1117. #     um 1 Bit nach rechts.
  1118. #   Bilde daraus die Ganzzahl-Wurzel, eine n+1-Digit-Zahl mit einer
  1119. #     führenden 1.
  1120. #   Runde das letzte Digit weg:
  1121. #     Bit 15 = 0 -> abrunden,
  1122. #     Bit 15 = 1, Rest =0 und Wurzel exakt -> round-to-even,
  1123. #     sonst aufrunden.
  1124. #   Bei rounding overflow Mantisse um 1 Bit nach rechts schieben
  1125. #     und Exponent incrementieren.
  1126.   local object LF_sqrt_LF(x)
  1127.     var reg4 object x;
  1128.     { var reg3 uintL uexp = TheLfloat(x)->expo;
  1129.       if (uexp==0) { return x; } # x=0.0 -> 0.0 als Ergebnis
  1130.      {var reg2 uintC len = TheLfloat(x)->len;
  1131.       # Radikanden bilden:
  1132.       SAVE_NUM_STACK # num_stack retten
  1133.       var reg7 uintD* r_MSDptr;
  1134.       var reg9 uintD* r_LSDptr;
  1135.       var reg6 uintL r_len = 2*(uintL)len+2; # Länge des Radikanden
  1136.       if ((intCsize < 32) && (r_len > (uintL)(bitc(intCsize)-1))) { fehler_LF_toolong(); }
  1137.       num_stack_need(r_len, r_MSDptr=,r_LSDptr=);
  1138.       uexp = uexp - LF_exp_mid + 1;
  1139.       if (uexp & bit(0))
  1140.         # Exponent gerade
  1141.         {var reg1 uintD* ptr =
  1142.            copy_loop_up(&TheLfloat(x)->data[0],r_MSDptr,len); # n Digits kopieren
  1143.          clear_loop_up(ptr,len+2); # n+2 Nulldigits anhängen
  1144.         }
  1145.         else
  1146.         # Exponent ungerade
  1147.         {var reg5 uintD carry_rechts = # n Digits kopieren und um 1 Bit rechts shiften
  1148.            shiftrightcopy_loop_up(&TheLfloat(x)->data[0],r_MSDptr,len,1,0);
  1149.          var reg1 uintD* ptr = &r_MSDptr[(uintP)len];
  1150.          *ptr++ = carry_rechts; # Übertrag und
  1151.          clear_loop_up(ptr,len+1); # n+1 Nulldigits anhängen
  1152.         }
  1153.       uexp = (sintL)((sintL)uexp >> 1); # Exponent halbieren
  1154.       uexp = uexp + LF_exp_mid;
  1155.       {# Ergebnis allozieren:
  1156.        var reg1 object y = allocate_lfloat(len,uexp,0);
  1157.        var reg5 uintD* y_mantMSDptr = &TheLfloat(y)->data[0];
  1158.        # Wurzel ziehen:
  1159.        var DS w;
  1160.        var reg8 boolean exactp;
  1161.        UDS_sqrt(r_MSDptr,r_len,r_LSDptr, &w, exactp=);
  1162.        # w ist die Ganzzahl-Wurzel, eine n+1-Digit-Zahl.
  1163.        RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  1164.        copy_loop_up(w.MSDptr,y_mantMSDptr,len); # NUDS nach y kopieren
  1165.        # Runden:
  1166.        if ( ((sintD)(w.LSDptr[-1]) >= 0) # nächstes Bit =0 -> abrunden
  1167.             || ( ((w.LSDptr[-1] & ((uintD)bit(intDsize-1)-1)) ==0) # =1 und weitere Bits >0 oder Rest >0 -> aufrunden
  1168.                  && exactp
  1169.                  # round-to-even
  1170.                  && ((w.LSDptr[-2] & bit(0)) ==0)
  1171.           )    )
  1172.          # abrunden
  1173.          {}
  1174.          else
  1175.          # aufrunden
  1176.          { if ( inc_loop_down(&y_mantMSDptr[(uintP)len],len) )
  1177.              # Übertrag durchs Aufrunden
  1178.              { y_mantMSDptr[0] = bit(intDsize-1); # Mantisse := 10...0
  1179.                (TheLfloat(y)->expo)++; # Exponenten incrementieren
  1180.          }   }
  1181.        return y;
  1182.     }}}
  1183.  
  1184. # LF_to_I(x) wandelt ein Long-Float x, das eine ganze Zahl darstellt,
  1185. # in ein Integer um.
  1186. # kann GC auslösen
  1187.   local object LF_to_I (object x);
  1188. # Methode:
  1189. # Falls x=0.0, Ergebnis 0.
  1190. # Sonst (ASH Vorzeichen*Mantisse (e-16n)).
  1191.   local object LF_to_I(x)
  1192.     var reg4 object x;
  1193.     { var reg5 uintL uexp = TheLfloat(x)->expo;
  1194.       if (uexp==0) { return Fixnum_0; } # x=0.0 -> Ergebnis 0
  1195.       # Mantisse zu einem Integer machen:
  1196.      {SAVE_NUM_STACK # num_stack retten
  1197.       var reg1 uintD* MSDptr;
  1198.       var reg8 uintD* LSDptr;
  1199.       var reg2 uintC len = TheLfloat(x)->len;
  1200.       var reg3 uintC len1 = len+1; # brauche 1 Digit mehr
  1201.       if (uintCoverflow(len1)) { fehler_LF_toolong(); }
  1202.       num_stack_need(len1, MSDptr=,LSDptr=);
  1203.       copy_loop_up(&TheLfloat(x)->data[0],&MSDptr[1],len); # Mantisse kopieren
  1204.       MSDptr[0] = 0; # und zusätzliches Nulldigit
  1205.       # Mantisse ist die UDS MSDptr/len1/LSDptr.
  1206.       if (R_minusp(x))
  1207.         # x<0 -> Mantisse negieren:
  1208.         { neg_loop_down(LSDptr,len1); }
  1209.       # Vorzeichen*Mantisse ist die DS MSDptr/len1/LSDptr.
  1210.       pushSTACK(DS_to_I(MSDptr,len1)); # Vorzeichen*Mantisse als Integer
  1211.       RESTORE_NUM_STACK # num_stack zurück
  1212.       # e-16n = uexp-LF_exp_mid-16n als Integer bilden:
  1213.       {var reg6 uintL sub = LF_exp_mid + intDsize*(uintL)len;
  1214.        var reg7 object shiftcount = UL_UL_minus_I(uexp,sub);
  1215.        # (ASH Vorzeichen*Mantisse (- e 16n)) durchführen:
  1216.        return I_I_ash_I(popSTACK(),shiftcount);
  1217.     }}}
  1218.  
  1219. # I_to_LF(x,len) wandelt ein Integer x in ein Long-Float mit len Digits um
  1220. # und rundet dabei.
  1221. # kann GC auslösen
  1222.   local object I_to_LF (object x, uintC len);
  1223. # Methode:
  1224. # x=0 -> Ergebnis 0.0
  1225. # Merke Vorzeichen von x.
  1226. # x:=(abs x)
  1227. # Exponent:=(integer-length x)
  1228. # Mantisse enthalte die höchstwertigen 16n Bits des Integers x (wobei die
  1229. #   führenden 16-(e mod 16) Nullbits zu streichen sind).
  1230. # Runde die weiteren Bits weg:
  1231. #   Kommen keine mehr -> abrunden,
  1232. #   nächstes Bit = 0 -> abrunden,
  1233. #   nächstes Bit = 1 und Rest =0 -> round-to-even,
  1234. #   nächstes Bit = 1 und Rest >0 -> aufrunden.
  1235. # Bei Aufrundung: rounding overflow -> Mantisse um 1 Bit nach rechts schieben
  1236. #   und Exponent incrementieren.
  1237.   local object I_to_LF(x,len)
  1238.     var reg7 object x;
  1239.     var reg3 uintC len;
  1240.     { if (eq(x,Fixnum_0)) { encode_LF0(len, return); } # x=0 -> Ergebnis 0.0
  1241.      {var reg9 signean sign = R_sign(x); # Vorzeichen von x
  1242.       if (!(sign==0)) { x = I_minus_I(x); } # Betrag von x nehmen
  1243.       {var reg9 uintL exp = I_integer_length(x); # (integer-length x) < intDsize*2^intCsize
  1244.        # Teste, ob exp <= LF_exp_high-LF_exp_mid :
  1245.        if (   (log2_intDsize+intCsize < 32)
  1246.            && ((uintL)(intDsize*bitc(intCsize)-1) <= (uintL)(LF_exp_high-LF_exp_mid))
  1247.           )
  1248.          {} # garantiert exp <= intDsize*2^intCsize-1 <= LF_exp_high-LF_exp_mid
  1249.          else
  1250.          { if (!(exp <= (uintL)(LF_exp_high-LF_exp_mid))) { fehler_overflow(); } }
  1251.        # Long-Float bauen:
  1252.        pushSTACK(x);
  1253.        {var reg4 object y = allocate_lfloat(len,exp+LF_exp_mid,sign);
  1254.         var reg1 uintD* y_mantMSDptr = &TheLfloat(y)->data[0];
  1255.         var reg2 uintD* x_MSDptr;
  1256.         var reg6 uintC x_len;
  1257.         I_to_NDS_nocopy(popSTACK(), x_MSDptr=,x_len=,_EMA_); # NDS zu x bilden, x_len>0
  1258.         # x_MSDptr/x_len/.. um (exp mod 16) Bits nach rechts shiften und in
  1259.         # y einfüllen (genauer: nur maximal len Digits davon):
  1260.         {var reg8 uintL shiftcount = exp % intDsize;
  1261.          # Die NDS fängt mit intDsize-shiftcount Nullbits an, dann kommt eine 1.
  1262.          if (x_len > len)
  1263.            { x_len -= 1+len;
  1264.              if (shiftcount>0)
  1265.                { var reg5 uintD carry_rechts =
  1266.                    shiftrightcopy_loop_up(&x_MSDptr[1],y_mantMSDptr,len,shiftcount,x_MSDptr[0]);
  1267.                  # Mantisse ist gefüllt. Runden:
  1268.                  if ( ((sintD)carry_rechts >= 0) # nächstes Bit =0 -> abrunden
  1269.                       || ( ((carry_rechts & ((uintD)bit(intDsize-1)-1)) ==0) # =1, Rest >0 -> aufrunden
  1270.                            && !test_loop_up(&x_MSDptr[1+(uintP)len],x_len)
  1271.                            # round-to-even
  1272.                            && ((y_mantMSDptr[(uintP)len-1] & bit(0)) ==0)
  1273.                     )    )
  1274.                    goto ab; # aufrunden
  1275.                    else
  1276.                    goto auf; # aufrunden
  1277.                }
  1278.                else
  1279.                { copy_loop_up(&x_MSDptr[1],y_mantMSDptr,len);
  1280.                  # Mantisse ist gefüllt. Runden:
  1281.                 {var reg5 uintD* ptr = &x_MSDptr[1+(uintP)len];
  1282.                  if ( (x_len==0) # keine Bits mehr -> abrunden
  1283.                       || ((sintD)ptr[0] >= 0) # nächstes Bit =0 -> abrunden
  1284.                       || ( ((ptr[0] & ((uintD)bit(intDsize-1)-1)) ==0) # =1, Rest >0 -> aufrunden
  1285.                            && !test_loop_up(&ptr[1],x_len-1)
  1286.                            # round-to-even
  1287.                            && ((ptr[-1] & bit(0)) ==0)
  1288.                     )    )
  1289.                    goto ab; # aufrunden
  1290.                    else
  1291.                    goto auf; # aufrunden
  1292.                }}
  1293.              auf: # aufrunden
  1294.                if ( inc_loop_down(&y_mantMSDptr[(uintP)len],len) )
  1295.                  # Übertrag durchs Aufrunden
  1296.                  { y_mantMSDptr[0] = bit(intDsize-1); # Mantisse := 10...0
  1297.                    # Exponenten incrementieren:
  1298.                    if (   (log2_intDsize+intCsize < 32)
  1299.                        && ((uintL)(intDsize*bitc(intCsize)-1) < (uintL)(LF_exp_high-LF_exp_mid))
  1300.                       )
  1301.                      # garantiert exp < intDsize*2^intCsize-1 <= LF_exp_high-LF_exp_mid
  1302.                      { (TheLfloat(y)->expo)++; } # jetzt exp <= LF_exp_high-LF_exp_mid
  1303.                      else
  1304.                      { if (++(TheLfloat(y)->expo) == LF_exp_high+1) { fehler_overflow(); } }
  1305.                  }
  1306.              ab: # abrunden
  1307.                ;
  1308.            }
  1309.            else # x_len <= len
  1310.            { var reg5 uintD carry_rechts;
  1311.              len -= x_len;
  1312.              x_len -= 1;
  1313.              if (shiftcount>0)
  1314.                { carry_rechts = shiftrightcopy_loop_up(&x_MSDptr[1],y_mantMSDptr,x_len,shiftcount,x_MSDptr[0]); }
  1315.                else
  1316.                { copy_loop_up(&x_MSDptr[1],y_mantMSDptr,x_len); carry_rechts = 0; }
  1317.             {var reg9 uintD* y_ptr = &y_mantMSDptr[x_len];
  1318.              *y_ptr++ = carry_rechts; # Carry als nächstes Digit
  1319.              clear_loop_up(y_ptr,len); # dann len-x_len Nulldigits
  1320.            }}
  1321.         }
  1322.         return y;
  1323.     }}}}
  1324.  
  1325. # RA_to_LF(x,len) wandelt eine rationale Zahl x in ein Long-Float
  1326. # mit len Digits um und rundet dabei.
  1327. # kann GC auslösen
  1328.   local object RA_to_LF (object x, uintC len);
  1329. # Methode:
  1330. # x ganz -> klar.
  1331. # x = +/- a/b mit Integers a,b>0:
  1332. #   Sei k,m so gewählt, daß
  1333. #     2^(k-1) <= a < 2^k, 2^(m-1) <= b < 2^m.
  1334. #   Dann ist 2^(k-m-1) < a/b < 2^(k-m+1).
  1335. #   Ergebnis-Vorzeichen := Vorzeichen von x.
  1336. #   Berechne k=(integer-length a) und m=(integer-length b).
  1337. #   Ergebnis-Exponent := k-m.
  1338. #   Ergebnis-Mantisse:
  1339. #     Berechne floor(2^(-k+m+16n+1)*a/b) :
  1340. #       Bei k-m>=16n+1 dividiere a durch (ash b (k-m-16n-1)),
  1341. #       bei k-m<16n+1 dividiere (ash a (-k+m+16n+1)) durch b.
  1342. #     Der erste Wert ist >=2^16n, <2^(16n+2).
  1343. #     Falls er >=2^(16n+1) ist, erhöhe Exponent um 1,
  1344. #       runde 2 Bits weg und schiebe dabei um 2 Bits nach rechts;
  1345. #     falls er <2^(16n+1) ist,
  1346. #       runde 1 Bit weg und schiebe dabei um 1 Bit nach rechts.
  1347.   local object RA_to_LF(x,len)
  1348.     var reg5 object x;
  1349.     var reg4 uintC len;
  1350.     { if (RA_integerp(x)) { return I_to_LF(x,len); }
  1351.       # x Ratio
  1352.       pushSTACK(TheRatio(x)->rt_den); # b
  1353.       x = TheRatio(x)->rt_num; # +/- a
  1354.      {var reg8 signean sign = R_sign(x); # Vorzeichen
  1355.       if (!(sign==0)) { x = I_minus_I(x); } # Betrag nehmen, liefert a
  1356.       pushSTACK(x);
  1357.       # Stackaufbau: b, a.
  1358.       {var reg3 sintL lendiff = I_integer_length(x) # (integer-length a)
  1359.                                 - I_integer_length(STACK_1); # (integer-length b)
  1360.        # |lendiff| < intDsize*2^intCsize. Da für LF-Exponenten ein sintL zur
  1361.        # Verfügung steht, braucht man keinen Test auf Overflow oder Underflow.
  1362.        {var reg6 uintL difflimit = intDsize*(uintL)len + 1; # 16n+1
  1363.         var reg1 object zaehler;
  1364.         var reg2 object nenner;
  1365.         if (lendiff > (sintL)difflimit)
  1366.           # 0 <= k-m-16n-1 < k < intDsize*2^intCsize
  1367.           { nenner = I_I_ash_I(STACK_1,
  1368.                                (log2_intDsize+intCsize<=oint_data_len # intDsize*2^intCsize <= 2^oint_data_len ?
  1369.                                  ? fixnum( (uintL)(lendiff - difflimit))
  1370.                                  : UL_to_I((uintL)(lendiff - difflimit))
  1371.                               ));
  1372.             zaehler = popSTACK(); # a
  1373.             skipSTACK(1);
  1374.           }
  1375.           else
  1376.           # 0 < -k+m+16n+1 <= m+1 + 16n < intDsize*2^intCsize + intDsize*2^intCsize
  1377.           { var reg7 object shiftcount = # -k+m+16n+1
  1378.               (log2_intDsize+intCsize+1<=oint_data_len # 2*intDsize*2^intCsize <= 2^oint_data_len ?
  1379.                 ? fixnum( (uintL)(difflimit - lendiff))
  1380.                 : UL_to_I((uintL)(difflimit - lendiff))
  1381.               );
  1382.             zaehler = I_I_ash_I(popSTACK(),shiftcount); # (ash a -k+m+16n+1)
  1383.             nenner = popSTACK(); # b
  1384.           }
  1385.         # Division zaehler/nenner durchführen:
  1386.         I_I_divide_I_I(zaehler,nenner);
  1387.        }
  1388.        # Stackaufbau: q, r.
  1389.        # 2^16n <= q < 2^(16n+2), also ist q Bignum mit n+1 Digits.
  1390.        {var reg7 object y = allocate_lfloat(len,lendiff+LF_exp_mid,sign); # neues Long-Float
  1391.         var reg2 uintD* y_mantMSDptr = &TheLfloat(y)->data[0];
  1392.         {var reg1 uintD* q_MSDptr = &TheBignum(STACK_1)->data[0];
  1393.          if (q_MSDptr[0] == 1) # erstes Digit =1 oder =2,3 ?
  1394.            # 2^16n <= q < 2^(16n+1), also 2^(k-m-1) < a/b < 2^(k-m).
  1395.            { # Mantisse mit einer Schiebeschleife um 1 Bit nach rechts füllen:
  1396.              var reg6 uintD rounding_bit =
  1397.                shiftrightcopy_loop_up(&q_MSDptr[1],y_mantMSDptr,len,1,1);
  1398.              if ( (rounding_bit == 0) # herausgeschobenes Bit =0 -> abrunden
  1399.                   || ( eq(STACK_0,Fixnum_0) # =1 und Rest r > 0 -> aufrunden
  1400.                        # round-to-even
  1401.                        && ((y_mantMSDptr[(uintP)len-1] & bit(0)) ==0)
  1402.                 )    )
  1403.                goto ab; # abrunden
  1404.                else
  1405.                goto auf; # aufrunden
  1406.            }
  1407.            else
  1408.            # 2^(16n+1) <= q < 2^(16n+2), also 2^(k-m) < a/b < 2^(k-m+1).
  1409.            { # Mantisse mit einer Schiebeschleife um 2 Bit nach rechts füllen:
  1410.              var reg6 uintD rounding_bits =
  1411.                shiftrightcopy_loop_up(&q_MSDptr[1],y_mantMSDptr,len,2,q_MSDptr[0]);
  1412.              (TheLfloat(y)->expo)++; # Exponenten incrementieren auf k-m+1
  1413.              if ( ((sintD)rounding_bits >= 0) # herausgeschobenes Bit =0 -> abrunden
  1414.                   || ( ((rounding_bits & bit(intDsize-2)) ==0) # =1 und nächstes Bit =1 oder Rest r > 0 -> aufrunden
  1415.                        && eq(STACK_0,Fixnum_0)
  1416.                        # round-to-even
  1417.                        && ((y_mantMSDptr[(uintP)len-1] & bit(0)) ==0)
  1418.                 )    )
  1419.                goto ab; # abrunden
  1420.                else
  1421.                goto auf; # aufrunden
  1422.            }
  1423.         }
  1424.         auf: # aufrunden
  1425.           { if ( inc_loop_down(&y_mantMSDptr[(uintP)len],len) )
  1426.               # Übertrag durchs Aufrunden
  1427.               { y_mantMSDptr[0] = bit(intDsize-1); # Mantisse := 10...0
  1428.                 (TheLfloat(y)->expo)++; # Exponenten incrementieren
  1429.           }   }
  1430.         ab: # abrunden
  1431.         skipSTACK(2);
  1432.         return y;
  1433.     }}}}
  1434.  
  1435.