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