home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / Misc / CLISP-1.LHA / CLISP960530-sr.lha / src / arilev0.d < prev    next >
Encoding:
Text File  |  1996-04-15  |  47.5 KB  |  1,040 lines
  1. # Arithmetik, Level 0
  2. # operiert auf einzelnen 16-Bit-Wörtern und 32-Bit-Wörtern (unsigned).
  3.  
  4. # Vorzeichen einer 32-Bit-Zahl bestimmen
  5. # sign_of_sint32(wert)
  6. # > wert: eine 32-Bit-Zahl
  7. # < sint16 ergebnis: 0 falls wert>=0, -1 falls wert<0.
  8.   extern sint16 sign_of_sint32 (sint32 wert);
  9. # im 68000-Assembler (Input D0.L, Output D0.W):
  10. #   SWAP D0  # Vorzeichen nach Bit 15 schieben
  11. #   EXT.L D0 # nach Bit 31..16 kopieren
  12. #   SWAP D0  # nach Bit 15..0 schieben
  13.   #if defined(GNU) && defined(MC680X0) && !defined(NO_ASM)
  14.     #define sign_of_sint32(wert)  \
  15.       ({var sint32 __wert = (wert);  \
  16.         var sint16 __ergebnis;        \
  17.         __asm__ ("\
  18.           swap %0;  \
  19.           extl %0;  \
  20.           swap %0   \
  21.           " : "=d" (__ergebnis) : "0" (__wert));  \
  22.         __ergebnis;                               \
  23.        })
  24.   #elif defined(SPARC) || defined(ARM)
  25.     #define sign_of_sint32(wert)  (((sint32)(wert)) >> 31)
  26.   #else
  27.     #define sign_of_sint32(wert)  ((sint32)(wert) >= 0 ? 0 : -1)
  28.   #endif
  29.  
  30. # Vorzeichen einer 16-Bit-Zahl bestimmen
  31. # sign_of_sint16(wert)
  32. # > wert: eine 16-Bit-Zahl
  33. # < sint16 ergebnis: 0 falls wert>=0, -1 falls wert<0.
  34.   extern sint16 sign_of_sint16 (sint16 wert);
  35. # im 68000-Assembler (Input D0.W, Output D0.W):
  36. #   EXT.L D0 # Vorzeichen nach Bit 31..16 kopieren
  37. #   SWAP D0  # nach Bit 15..0 schieben
  38.   #if defined(GNU) && defined(MC680X0) && !defined(NO_ASM)
  39.     #define sign_of_sint16(wert)  \
  40.       ({var sint16 __wert = (wert);  \
  41.         var sint16 __ergebnis;        \
  42.         __asm__ ("\
  43.           extl %0;  \
  44.           swap %0   \
  45.           " : "=d" (__ergebnis) : "0" (__wert));  \
  46.         __ergebnis;                               \
  47.        })
  48.   #elif defined(SPARC)
  49.     #define sign_of_sint16(wert)  (((sint32)(sint16)(wert)) >> 31)
  50.   #else
  51.     #define sign_of_sint16(wert)  ((sint16)(wert) >= 0 ? 0 : -1)
  52.   #endif
  53.  
  54. # High-Word einer 32-Bit-Zahl bestimmen
  55. # high16(wert)
  56.   extern uint16 high16 (uint32 wert);
  57. # im 68000-Assembler (Input D0.L, Output D0.W):
  58. #   SWAP D0
  59.   #ifdef GNU
  60.     #if defined(MC680X0) && !defined(NO_ASM)
  61.       #define high16(wert)  \
  62.         ({var uint32 __wert = (wert);  \
  63.           var uint16 __ergebnis;        \
  64.           __asm__ ("\
  65.             swap %0   \
  66.             " : "=d" (__ergebnis) : "0" (__wert));  \
  67.           __ergebnis;                               \
  68.          })
  69.     #endif
  70.   #endif
  71.   #ifndef high16
  72.     #define high16(wert)  ((uint16)((uint32)(wert)>>16))
  73.   #endif
  74.  
  75. # Low-Word einer 32-Bit-Zahl bestimmen
  76. # low16(wert)
  77.   extern uint16 low16 (uint32 wert);
  78.   #define low16(wert)  ((uint16)(uint32)(wert))
  79.  
  80. # Eine 32-Bit-Zahl aus ihrem High-Word und ihrem Low-Word bestimmen:
  81. # highlow32(uint16 high, uint16 low)
  82.   extern uint32 highlow32 (uint16 high, uint16 low);
  83. # im 68000-Assembler (Input D0.W,D1.W, Output D0.L):
  84. #   SWAP D0
  85. #   MOVE.W D1,D0
  86.   #ifdef GNU
  87.     #if defined(MC680X0) && !defined(NO_ASM)
  88.       #define highlow32(high,low)  \
  89.         ({var uint16 __high = (high);  \
  90.           var uint16 __low = (low);    \
  91.           var uint32 __ergebnis;       \
  92.           __asm__ __volatile__ ("\
  93.             swap %0;     \
  94.             movew %2,%0  \
  95.             " : "=&d" (__ergebnis) : "0" (__high), "g" (__low));  \
  96.           __ergebnis;                                             \
  97.          })
  98.     #endif
  99.   #endif
  100.   #ifndef highlow32
  101.     #define highlow32(high,low)  \
  102.       (((uint32)(uint16)(high) << 16) | (uint32)(uint16)(low))
  103.   #endif
  104.  
  105. # Eine 32-Bit-Zahl aus ihrem High-Word und ihrem Low-Word 0 bestimmen:
  106. # highlow32_0(uint16 high)
  107.   extern uint32 highlow32_0 (uint16 high);
  108.   # define highlow32_0(high)  highlow32(high,0)
  109. # im 68000-Assembler (Input D0.W,D1.W, Output D0.L):
  110. #   SWAP D0
  111. #   CLR.W D0
  112.   #ifdef GNU
  113.     #if defined(MC680X0) && !defined(NO_ASM)
  114.       #define highlow32_0(high)  \
  115.         ({var uint16 __high = (high);  \
  116.           var uint32 __ergebnis;       \
  117.           __asm__ __volatile__ ("\
  118.             swap %0;     \
  119.             clrw %0      \
  120.             " : "=d" (__ergebnis) : "0" (__high));  \
  121.           __ergebnis;                               \
  122.          })
  123.     #endif
  124.   #endif
  125.   #ifndef highlow32_0
  126.     #define highlow32_0(high)  ((uint32)(uint16)(high) << 16)
  127.   #endif
  128.  
  129. # Multipliziert zwei 16-Bit-Zahlen miteinander und liefert eine 32-Bit-Zahl:
  130. # mulu16(arg1,arg2)
  131. # > arg1, arg2 : zwei 16-Bit-Zahlen
  132. # < ergebnis: eine 32-Bit-Zahl
  133.   extern uint32 mulu16 (uint16 arg1, uint16 arg2);
  134. # in 68000-Assembler (Input D0.W, D1.W, Output D0.L):
  135. #   MULU D1,D0
  136.   #ifdef GNU
  137.     #if defined(SPARC) && defined(FAST_DOUBLE) # Ist das schneller als _mulu16 ??
  138.       #define mulu16(arg1,arg2)  \
  139.         ({var union { double f; uint32 i[2]; } __fi; \
  140.           __fi.f = (double)(sint32)(uint16)(arg1)*(double)(sint32)(uint16)(arg2) \
  141.                    + (double)(4503599627370496.0L); # + 2^52, zum Normalisieren  \
  142.           __fi.i[1]; # untere 32 Bit herausholen (benutzt BIG_ENDIAN_P !)        \
  143.          })
  144.     #elif defined(I80Z86) && !defined(NO_ASM)
  145.       #define mulu16(arg1,arg2)  \
  146.         ({ var register uint16 _hi;                                         \
  147.            var register uint16 _lo;                                         \
  148.            __asm__("mulw %2"                                                \
  149.                    : "=d" /* %dx */ (_hi), "=a" /* %ax */ (_lo)             \
  150.                    : "rm" ((uint16)(arg1)), "1" /* %eax */ ((uint16)(arg2)) \
  151.                   );                                                        \
  152.            highlow32(_hi,_lo);                                              \
  153.          })
  154.     #endif
  155.   #else
  156.     #if defined(SPARC)
  157.       #define mulu16  mulu16_  # extern in Assembler
  158.     #endif
  159.   #endif
  160.   #ifndef mulu16
  161.     #define mulu16(arg1,arg2)  ((uint32)(uint16)(arg1)*(uint32)(uint16)(arg2))
  162.   #endif
  163.  
  164. # Multipliziert zwei 24-Bit-Zahlen zusammen und liefert eine 48-Bit-Zahl.
  165. # mulu24(arg1,arg2,hi=,lo=);
  166. # > arg1, arg2 : zwei 24-Bit-Zahlen
  167. # < 2^32*hi+lo : eine 48-Bit-Zahl
  168.   #if defined(SPARC) && defined(FAST_DOUBLE)
  169.     #define mulu24(x,y,hi_zuweisung,lo_zuweisung)  \
  170.       { var reg1 uint32 _x = (x);                                                       \
  171.         var reg2 uint32 _y = (y);                                                       \
  172.         var union { double f; uint32 i[2]; uint16 s[4]; } __fi;                         \
  173.         __fi.f = (double)(sint32)(_x)*(double)(sint32)(_y)                              \
  174.                  + (double)(4503599627370496.0L); # + 2^52, zum Normalisieren           \
  175.         hi_zuweisung __fi.s[1]; # mittlere 16 Bit herausholen, (benutzt BIG_ENDIAN_P !) \
  176.         lo_zuweisung __fi.i[1]; # untere 32 Bit herausholen (benutzt BIG_ENDIAN_P !)    \
  177.       }
  178.   #elif defined(MC680X0) && !defined(MC680Y0)
  179.     # Methode:
  180.     # Sei x = x1*2^16+x0, y = y1*2^16+y0  mit
  181.     # 0 <= x0,y0 < 2^16, 0 <= x1,y1 < 2^8 . Dann ist das Produkt x*y
  182.     # >=0 und < 2^48. Es belegt also 3 16-Bit-Worte.
  183.     # In  x * y = x1*y1*2^32 + (x1*y0+x0*y1)*2^16 + x0*y0
  184.     # bestimmen die ersten beiden Summanden und das High-Word des dritten
  185.     # Summanden die beiden höherwertigen Words des Ergebnisses (und
  186.     # hierbei gibt es keinen Überlauf, da das Produkt in 3 Words paßt!),
  187.     # das Low-Word des dritten Summanden ist auch das des Ergebnisses.
  188.     #define mulu24(x,y,hi_zuweisung,lo_zuweisung)  \
  189.       { var reg1 uint32 _x = (x);                                  \
  190.         var reg2 uint32 _y = (y);                                  \
  191.         var reg4 uint32 _erg21; # Teilresultat für Words 2 und 1 des Ergebnisses \
  192.         var reg3 uint32 _erg10; # Teilresultat für Words 1 und 0 des Ergebnisses \
  193.         { var reg3 uint16 _x1 = high16(_x);                        \
  194.           var reg4 uint16 _y1 = high16(_y);                        \
  195.           _erg21 = highlow32_0(mulu16(_x1,_y1))                    \
  196.                   + mulu16(_x1,low16(_y)) + mulu16(low16(_x),_y1); \
  197.         }                                                          \
  198.         _erg10 = mulu16(low16(_x),low16(_y));                      \
  199.         # Teilresultate kombinieren:                               \
  200.         _erg21 += high16(_erg10);                                  \
  201.         hi_zuweisung high16(_erg21);                               \
  202.         lo_zuweisung highlow32(low16(_erg21),low16(_erg10));       \
  203.       }
  204.   #else
  205.     #define mulu24  mulu32
  206.   #endif
  207.  
  208. # Multipliziert zwei 32-Bit-Zahlen miteinander und liefert eine 64-Bit-Zahl:
  209. # mulu32(arg1,arg2,hi=,lo=);
  210. # > arg1, arg2 : zwei 32-Bit-Zahlen
  211. # < 2^32*hi+lo : eine 64-Bit-Zahl
  212.   extern uint32 mulu32_ (uint32 arg1, uint32 arg2); # -> Low-Teil
  213.   extern uint32 mulu32_high;                        # -> High-Teil
  214. # in 68000-Assembler (Input D0.L,D1.L, Output D0.L,D1.L, verändert D2-D4):
  215. #   ; D0.L = 2^16*a+b, D1.L = 2^16*c+d -> Produkt
  216. #   ; (2^16*a+b)*(2^16*c+d) = 2^32*a*c + 2^16*(a*d+b*c) + b*d
  217. #   MOVE.L D0,D2 ! SWAP D2
  218. #   MOVE.L D1,D3 ! SWAP D1
  219. #   MOVE.L D1,D4 ! MULU D2,D1 ; a*c
  220. #   MULU D3,D2 ; a*d
  221. #   MULU D0,D4 ; b*c
  222. #   MULU D3,D0 ; b*d
  223. #   CLR.L D3 ; Hilfsregister für Zero-Extend
  224. #   SWAP D2 ! MOVE.W D2,D3 ! ADD.L D3,D1 ; high16(a*d) zu D1.L addieren
  225. #   SWAP D4 ! MOVE.W D4,D3 ! ADD.L D3,D1 ; high16(b*c) zu D1.L addieren
  226. #   CLR.W D2 ! ADD.L D2,D0 ! BCC.S \1 ! ADDQ.L #1,D1 ! \1: ; 2^16*low16(a*d) zu D0.L addieren
  227. #   CLR.W D4 ! ADD.L D4,D0 ! BCC.S \2 ! ADDQ.L #1,D1 ! \2: ; 2^16*low16(b*c) zu D0.L addieren
  228. #   ; D0.L = lo, D1.L = hi fertig.
  229. # in 68020-Assembler (Input D0.L,D1.L, Output D0.L,D1.L):
  230. #   MULU.L D1,D1:D0
  231.   #ifdef GNU
  232.     #ifdef MC680X0
  233.       #if !(defined(MC680Y0) && !defined(NO_ASM))
  234.         #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
  235.           ({ var reg4 uint32 _x = (x);                                      \
  236.              var reg4 uint32 _y = (y);                                      \
  237.              var reg5 uint16 _x1 = high16(_x);                              \
  238.              var reg5 uint16 _x0 = low16(_x);                               \
  239.              var reg5 uint16 _y1 = high16(_y);                              \
  240.              var reg5 uint16 _y0 = low16(_y);                               \
  241.              var reg3 uint32 _hi = mulu16(_x1,_y1); # obere Portion         \
  242.              var reg2 uint32 _lo = mulu16(_x0,_y0); # untere Portion        \
  243.              {var reg1 uint32 _mid = mulu16(_x0,_y1); # 1. mittlere Portion \
  244.               _hi += high16(_mid); _mid = highlow32_0(low16(_mid));         \
  245.               _lo += _mid; if (_lo < _mid) { _hi += 1; } # 64-Bit-Addition  \
  246.              }                                                              \
  247.              {var reg1 uint32 _mid = mulu16(_x1,_y0); # 2. mittlere Portion \
  248.               _hi += high16(_mid); _mid = highlow32_0(low16(_mid));         \
  249.               _lo += _mid; if (_lo < _mid) { _hi += 1; } # 64-Bit-Addition  \
  250.              }                                                              \
  251.              hi_zuweisung _hi;                                              \
  252.              lo_zuweisung _lo;                                              \
  253.            })
  254.       #else
  255.         #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
  256.           ({ var reg3 uint32 _x = (x); \
  257.              var reg4 uint32 _y = (y); \
  258.              var reg2 uint32 _hi;      \
  259.              var reg1 uint32 _lo;      \
  260.              __asm__("mulul %3,%0:%1" : "=d" (_hi), "=d"(_lo) : "1" (_x), "dm" (_y) ); \
  261.              hi_zuweisung _hi;         \
  262.              lo_zuweisung _lo;         \
  263.            })
  264.       #endif
  265.     #elif defined(SPARC)
  266.       #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
  267.         ({ lo_zuweisung mulu32_(x,y); # extern in Assembler \
  268.           {var register uint32 _hi __asm__("%g1");          \
  269.            hi_zuweisung _hi;                                \
  270.          }})
  271.     #elif defined(ARM)
  272.       #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
  273.         ({ lo_zuweisung mulu32_(x,y); # extern in Assembler \
  274.           {var register uint32 _hi __asm__("%r1"/*"%a2"*/); \
  275.            hi_zuweisung _hi;                                \
  276.          }})
  277.     #elif defined(I80Z86) && !defined(NO_ASM)
  278.       #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
  279.         ({ var register uint32 _hi;                                  \
  280.            var register uint32 _lo;                                  \
  281.            __asm__("mull %2"                                         \
  282.                    : "=d" /* %edx */ (_hi), "=a" /* %eax */ (_lo)    \
  283.                    : "g" ((uint32)(x)), "1" /* %eax */ ((uint32)(y)) \
  284.                   );                                                 \
  285.            hi_zuweisung _hi; lo_zuweisung _lo;                       \
  286.          })
  287.     #elif defined(MIPS) && !defined(NO_ASM)
  288.       #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
  289.         ({ var register uint32 _hi;                       \
  290.            var register uint32 _lo;                       \
  291.            __asm__("multu %3,%2 ; mfhi %0 ; mflo %1"      \
  292.                    : "=r" (_hi), "=r" (_lo)               \
  293.                    : "r" ((uint32)(x)), "r" ((uint32)(y)) \
  294.                   );                                      \
  295.            hi_zuweisung _hi; lo_zuweisung _lo;            \
  296.          })
  297.     #elif defined(HAVE_LONGLONG)
  298.       #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
  299.         ({ var register uint64 _prod = (uint64)(x) * (uint64)(y); \
  300.            hi_zuweisung (uint32)(_prod>>32);                      \
  301.            lo_zuweisung (uint32)(_prod);                          \
  302.          })
  303.     #endif
  304.   #endif
  305.   #if defined(WATCOM) && defined(I80Z86) && !defined(NO_ASM)
  306.     #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
  307.       { var register uint32 _hi;                  \
  308.         var register uint32 _lo;                  \
  309.         _lo = mulu32_(x,y), _hi = mulu32_high_(); \
  310.         hi_zuweisung _hi; lo_zuweisung _lo;       \
  311.       }
  312.     extern uint32 mulu32_high_ (void);
  313.     #pragma aux mulu32_ = 0xF7 0xE2 /* mull %edx */ parm [eax] [edx] value [eax] modify [eax edx];
  314.     #pragma aux mulu32_high_ = /* */ value [edx] modify [];
  315.   #endif
  316.   #ifndef mulu32
  317.     #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
  318.       { lo_zuweisung mulu32_(x,y); hi_zuweisung mulu32_high; }
  319.     #if defined(MC680X0) || defined(SPARC) || defined(ARM) || (defined(I80Z86) && !defined(WATCOM) && !defined(MICROSOFT)) || defined(MIPS) || defined(HPPA) || defined(VAX)
  320.       # mulu32_ extern in Assembler
  321.       #if defined(SPARC)
  322.         #define mulu32_high  (uint32)(_get_g1()) # Rückgabe im Register %g1
  323.       #elif defined(LISPARIT) && !defined(HPPA) # In arihppa.d ist mulu32_high bereits definiert.
  324.         global uint32 mulu32_high;
  325.       #endif
  326.     #else
  327.       #ifdef LISPARIT
  328.       global uint32 mulu32_high;
  329.       global uint32 mulu32_(x,y)
  330.         var reg4 uint32 x;
  331.         var reg4 uint32 y;
  332.         { var reg5 uint16 x1 = high16(x);
  333.           var reg5 uint16 x0 = low16(x);
  334.           var reg5 uint16 y1 = high16(y);
  335.           var reg5 uint16 y0 = low16(y);
  336.           var reg3 uint32 hi = mulu16(x1,y1); # obere Portion
  337.           var reg2 uint32 lo = mulu16(x0,y0); # untere Portion
  338.           {var reg1 uint32 mid = mulu16(x0,y1); # 1. mittlere Portion
  339.            hi += high16(mid); mid = highlow32_0(low16(mid));
  340.            lo += mid; if (lo < mid) { hi += 1; } # 64-Bit-Addition
  341.           }
  342.           {var reg1 uint32 mid = mulu16(x1,y0); # 2. mittlere Portion
  343.            hi += high16(mid); mid = highlow32_0(low16(mid));
  344.            lo += mid; if (lo < mid) { hi += 1; } # 64-Bit-Addition
  345.           }
  346.           mulu32_high = hi; return lo;
  347.         }
  348.       #endif
  349.     #endif
  350.   #endif
  351.  
  352. # Multipliziert zwei 32-Bit-Zahlen miteinander und liefert eine 32-Bit-Zahl:
  353. # mulu32_unchecked(arg1,arg2)
  354. # > arg1, arg2 : zwei 32-Bit-Zahlen
  355. # < ergebnis : eine 32-Bit-Zahl
  356. # Es wird vorausgesetzt, daß arg1*arg2 < 2^32.
  357.   #if (defined(GNU) && defined(MC680X0) && !defined(MC680Y0))
  358.     extern uint32 mulu32_unchecked (uint32 x, uint32 y);
  359.     #ifdef LISPARIT
  360.     global uint32 mulu32_unchecked(x,y)
  361.       var reg2 uint32 x;
  362.       var reg3 uint32 y;
  363.       { # Methode:
  364.         # Falls x>=2^16 und y>=2^16 wäre, wäre das Produkt zu groß.
  365.         # Falls x<2^16 : y = y1*2^16+y0 schreiben, (x*y1)*2^16 + (x*y0) bilden.
  366.         # Falls y<2^16 : x = x1*2^16+x0 schreiben, (x1*y)*2^16 + (x0*y) bilden.
  367.         # Falls sogar x<2^16 und y<2^16: nur x*y bilden.
  368.         var reg4 uint16 x1 = high16(x);
  369.         var reg5 uint16 y1 = high16(y);
  370.         if (x1==0)
  371.           if (y1==0)
  372.             return mulu16((uint16)(x),(uint16)(y));
  373.             else
  374.             return highlow32_0(mulu16((uint16)(x),y1))
  375.                    + mulu16((uint16)(x),low16(y));
  376.           else
  377.           return highlow32_0(mulu16(x1,(uint16)(y)))
  378.                  + mulu16(low16(x),(uint16)(y));
  379.       }
  380.     #endif
  381.   #elif defined(SPARC)
  382.     extern uint32 mulu32_unchecked (uint32 x, uint32 y); # extern in Assembler
  383.   #else
  384.     # Wir können dafür auch die Bibliotheksroutine des C-Compilers nehmen:
  385.     #define mulu32_unchecked(x,y)  ((uint32)((uint32)(x)*(uint32)(y)))
  386.   #endif
  387.  
  388. # Dividiert eine 16-Bit-Zahl durch eine 16-Bit-Zahl und
  389. # liefert einen 16-Bit-Quotienten und einen 16-Bit-Rest.
  390. # divu_1616_1616(x,y,q=,r=);
  391. # > uint16 x: Zähler
  392. # > uint16 y: Nenner
  393. # < uint16 q: floor(x/y)
  394. # < uint16 r: x mod y
  395. # < x = q*y+r
  396.   #define divu_1616_1616(x,y,q_zuweisung,r_zuweisung)  \
  397.     { var reg1 uint16 __x = (x);   \
  398.       var reg2 uint16 __y = (y);   \
  399.       q_zuweisung floor(__x,__y);  \
  400.       r_zuweisung (__x % __y);     \
  401.     }
  402.  
  403. # Dividiert eine 32-Bit-Zahl durch eine 16-Bit-Zahl und
  404. # liefert einen 16-Bit-Quotienten und einen 16-Bit-Rest.
  405. # divu_3216_1616(x,y,q=,r=);
  406. # > uint32 x: Zähler
  407. # > uint16 y: Nenner
  408. # > Es sei bekannt, daß 0 <= x < 2^16*y .
  409. # < uint16 q: floor(x/y)
  410. # < uint16 r: x mod y
  411. # < x = q*y+r
  412.   extern uint16 divu_3216_1616_ (uint32 x, uint16 y); # -> Quotient q
  413.   extern uint16 divu_16_rest;                         # -> Rest r
  414. # im 68000-Assembler (Input D0.L,D1.W, Output D0.W,D1.W):
  415. #   DIVU D1,D0    ; D0.L=x / D1.W=y -> q=D0.W, r=D0.H.W
  416. #   MOVE.L D0,D1
  417. #   SWAP D1
  418.   #ifdef GNU
  419.     #if defined(SPARC)
  420.       #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
  421.         ({ var reg1 uint32 __qr = divu_3216_1616_(x,y); # extern in Assembler \
  422.            q_zuweisung low16(__qr);  \
  423.            r_zuweisung high16(__qr); \
  424.          })
  425.     #elif defined(MC680X0) && !defined(NO_ASM)
  426.       #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
  427.         ({var uint32 __x = (x);      \
  428.           var uint16 __y = (y);      \
  429.           var uint32 __qr;           \
  430.           __asm__ __volatile__ ("\
  431.             divu %2,%0   \
  432.             " : "=d" (__qr) : "0" (__x), "dm" (__y));  \
  433.           q_zuweisung low16(__qr);   \
  434.           r_zuweisung high16(__qr);  \
  435.          })
  436.     #elif defined(I80Z86) && !defined(NO_ASM)
  437.       #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
  438.         ({var uint32 __x = (x);  \
  439.           var uint16 __y = (y);  \
  440.           var uint16 __q;        \
  441.           var uint16 __r;        \
  442.           __asm__("divw %4"      \
  443.                   : "=a" /* %ax */ (__q), "=d" /* %dx */ (__r) \
  444.                   : "1" /* %dx */ ((uint16)(high16(__x))), "0" /* %ax */ ((uint16)(low16(__x))), "rm" (__y) \
  445.                  );              \
  446.           q_zuweisung __q;       \
  447.           r_zuweisung __r;       \
  448.          })
  449.     #else
  450.       #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
  451.         ({var uint32 __x = (x);            \
  452.           var uint16 __y = (y);            \
  453.           var uint16 __q = floor(__x,__y); \
  454.           q_zuweisung __q;                 \
  455.           r_zuweisung (__x - __q * __y);   \
  456.          })
  457.     #endif
  458.   #else
  459.     #if defined(SPARC)
  460.       #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
  461.         { var reg1 uint32 __qr = divu_3216_1616_(x,y); # extern in Assembler \
  462.           q_zuweisung low16(__qr);  \
  463.           r_zuweisung high16(__qr); \
  464.         }
  465.     #elif defined(ARM)
  466.       #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
  467.         { q_zuweisung divu_3216_1616_(x,y); # extern in Assembler \
  468.           r_zuweisung divu_16_rest;                               \
  469.         }
  470.       #ifdef LISPARIT
  471.       global uint16 divu_16_rest;
  472.       #endif
  473.     #else
  474.       #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
  475.         { q_zuweisung divu_3216_1616_(x,y); r_zuweisung divu_16_rest; }
  476.       #ifdef LISPARIT
  477.       global uint16 divu_16_rest;
  478.       global uint16 divu_3216_1616_(x,y)
  479.         var reg3 uint32 x;
  480.         var reg2 uint16 y;
  481.         { var reg1 uint16 q = floor(x,(uint32)y);
  482.           divu_16_rest = x - (uint32)q * (uint32)y;
  483.           return q;
  484.         }
  485.       #endif
  486.     #endif
  487.   #endif
  488.  
  489. # Dividiert eine 32-Bit-Zahl durch eine 16-Bit-Zahl und
  490. # liefert einen 32-Bit-Quotienten und einen 16-Bit-Rest.
  491. # divu_3216_3216(x,y,q=,r=);
  492. # > uint32 x: Zähler
  493. # > uint16 y: Nenner
  494. # Es sei bekannt, daß y>0.
  495. # < uint32 q: floor(x/y)
  496. # < uint16 r: x mod y
  497. # < x = q*y+r
  498.   extern uint32 divu_3216_3216_ (uint32 x, uint16 y); # -> Quotient q
  499.   extern uint16 divu_16_rest;                         # -> Rest r
  500. # im 68000-Assembler (Input D0.L,D1.W, Output D0.L,D1.W, verändert D2,D3):
  501. #   MOVE.L D0,D2
  502. #   CLR.W D2
  503. #   SWAP D2        ; D2.L = D2.W := D0.H.W = high16(x)
  504. #   DIVU D1,D2     ; durch y dividieren
  505. #   MOVE.W D2,D3   ; Quotient nach D3.W
  506. #   MOVE.W D0,D2   ; Rest (in D2.H.W) mit D0.W = low16(x) kombinieren
  507. #   DIVU D1,D2     ; und wieder durch y dividieren
  508. #   MOVE.W D3,D0   ; ersten Quotiententeil
  509. #   SWAP D0        ; mal 2^16
  510. #   MOVE.W D2,D0   ; plus zweiten Quotiententeil, liefert q
  511. #   SWAP D2
  512. #   MOVE.W D2,D1   ; r = Rest der zweiten Division
  513. # oder (Input D0.L,D1.W, Output D0.L,D1.W, verändert D2,D3):
  514. #   MOVE.L D0,D2   ; x retten
  515. #   CLR.W D0
  516. #   SWAP D0        ; D0.L = D0.W := high16(x)
  517. #   DIVU D1,D0     ; durch y dividieren
  518. #   MOVE.W D0,D3   ; Quotient nach D3.W
  519. #   MOVE.W D2,D0   ; Rest (in D0.H.W) mit D2.W = low16(x) kombinieren
  520. #   DIVU D1,D0     ; und wieder durch y dividieren
  521. #   SWAP D0
  522. #   MOVE.W D0,D1   ; r = Rest der zweiten Division
  523. #   MOVE.W D3,D0
  524. #   SWAP D0        ; beide Quotienten kombinieren, liefert q
  525.   #if defined(SPARC) || defined(I80Z86)
  526.     #define divu_3216_3216  divu_3232_3232
  527.   #elif 1
  528.     # Methode: (beta = 2^16)
  529.     # x = x1*beta+x0 schreiben.
  530.     # Division mit Rest: x1 = q1*y + r1, wobei 0 <= x1 < beta <= beta*y.
  531.     # Also 0 <= q1 < beta, 0 <= r1 < y.
  532.     # Division mit Rest: (r1*beta+x0) = q0*y + r0, wobei 0 <= r1*beta+x0 < beta*y.
  533.     # Also 0 <= q0 < beta, 0 <= r0 < y
  534.     # und x = x1*beta+x0 = (q1*beta+q0)*y + r0.
  535.     # Setze q := q1*beta+q0 und r := r0.
  536.     #ifdef GNU
  537.       #define divu_3216_3216(x,y,q_zuweisung,r_zuweisung)  \
  538.         ({var uint32 _x = (x);     \
  539.           var uint16 _y = (y);     \
  540.           var uint16 _q1;          \
  541.           var uint16 _q0;          \
  542.           var uint16 _r1;          \
  543.           divu_3216_1616(high16(_x),_y, _q1 = , _r1 = ); \
  544.           divu_3216_1616(highlow32(_r1,low16(_x)),_y, _q0 = , _EMA_ r_zuweisung); \
  545.           q_zuweisung highlow32(_q1,_q0);  \
  546.          })
  547.     #else
  548.       #define divu_3216_3216(x,y,q_zuweisung,r_zuweisung)  \
  549.         {var reg1 uint32 _x = (x);     \
  550.          var reg2 uint16 _y = (y);     \
  551.          var reg3 uint16 _q1;          \
  552.          var reg4 uint16 _q0;          \
  553.          var reg5 uint16 _r1;          \
  554.          divu_3216_1616(high16(_x),_y, _q1 = , _r1 = ); \
  555.          divu_3216_1616(highlow32(_r1,low16(_x)),_y, _q0 = , _EMA_ r_zuweisung); \
  556.          q_zuweisung highlow32(_q1,_q0);  \
  557.         }
  558.     #endif
  559.   #else
  560.     #define divu_3216_3216(x,y,q_zuweisung,r_zuweisung)  \
  561.       { q_zuweisung divu_3216_3216_(x,y); r_zuweisung divu_16_rest; }
  562.     #if 0
  563.       # divu_3216_3216_ extern in Assembler
  564.     #else
  565.       #ifdef LISPARIT
  566.       global uint32 divu_3216_3216_(x,y)
  567.         var reg1 uint32 x;
  568.         var reg2 uint16 y;
  569.         { var reg4 uint16 q1;
  570.           var reg5 uint16 q0;
  571.           var reg3 uint16 r1;
  572.           divu_3216_1616(high16(x),y, q1 = , r1 = );
  573.           divu_3216_1616(highlow32(r1,low16(x)),y, q0 = , divu_16_rest =);
  574.           return highlow32(q1,q0);
  575.         }
  576.       #endif
  577.     #endif
  578.   #endif
  579.  
  580. # Dividiert eine 32-Bit-Zahl durch eine 32-Bit-Zahl und
  581. # liefert einen 32-Bit-Quotienten und einen 32-Bit-Rest.
  582. # divu_3232_3232(x,y,q=,r=);
  583. # > uint32 x: Zähler
  584. # > uint32 y: Nenner
  585. # Es sei bekannt, daß y>0.
  586. # < uint32 q: floor(x/y)
  587. # < uint32 r: x mod y
  588. # < x = q*y+r
  589.   extern uint32 divu_3232_3232_ (uint32 x, uint32 y); # -> Quotient q
  590.   extern uint32 divu_32_rest;                         # -> Rest r
  591.   #if defined(SPARC) || defined(I80Z86) || defined(HPPA_DIV_WORKS)
  592.     #define divu_3232_3232(x,y,q_zuweisung,r_zuweisung)  \
  593.       divu_6432_3232(0,x,y,_EMA_ q_zuweisung,_EMA_ r_zuweisung)
  594.     #define divu_3232_3232_(x,y) divu_6432_3232_(0,x,y)
  595.   #elif 1
  596.     # Methode: (beta = 2^n = 2^16, n = 16)
  597.     # Falls y < beta, handelt es sich um eine 32-durch-16-Bit-Division.
  598.     # Falls y >= beta:
  599.     # Quotient  q = floor(x/y) < beta  (da 0 <= x < beta^2, y >= beta).
  600.     # y habe genau n+k Bits (1 <= k <= n), d.h. 2^(n+k-1) <= y < 2^(n+k).
  601.     # Schreibe  x = 2^k*x1 + x0  mit  x1 := floor(x/2^k)
  602.     # und       y = 2^k*y1 + y0  mit  y1 := floor(y/2^k)
  603.     # und bilde den Näherungs-Quotienten floor(x1/y1)
  604.     # oder (noch besser) floor(x1/(y1+1)).
  605.     # Wegen 0 <= x1 < 2^(2n) und 0 < 2^(n-1) <= y1 < 2^n
  606.     # und  x1/(y1+1) <= x/y < x1/(y1+1) + 2
  607.     # (denn x1/(y1+1) = (x1*2^k)/((y1+1)*2^k) <= (x1*2^k)/y <= x/y
  608.     # und x/y - x1/(y1+1) = (x+x*y1-x1*y)/(y*(y1+1))
  609.     # = (x+x0*y1-x1*y0)/(y*(y1+1)) <= (x+x0*y1)/(y*(y1+1))
  610.     # <= x/(y*(y1+1)) + x0/y
  611.     # <= 2^(2n)/(2^(n+k-1)*(2^(n-1)+1)) + 2^k/2^(n+k-1)
  612.     # = 2^(n-k+1)/(2^(n-1)+1) + 2^(1-n) <= 2^n/(2^(n-1)+1) + 2^(1-n) < 2 )
  613.     # gilt  floor(x1/(y1+1)) <= floor(x/y) <= floor(x1/(y1+1)) + 2  .
  614.     # Man bildet also  q:=floor(x1/(y1+1))  (ein Shift um n Bit oder
  615.     # eine (2n)-durch-n-Bit-Division, mit Ergebnis q <= floor(x/y) < beta)
  616.     # und x-q*y und muß hiervon noch höchstens 2 mal y abziehen und q
  617.     # incrementieren, um den Quotienten  q = floor(x/y)  und den Rest
  618.     # x-floor(x/y)*y  der Division zu bekommen.
  619.     #define divu_3232_3232(x,y,q_zuweisung,r_zuweisung)  \
  620.       { var uint32 _x = (x);                                                    \
  621.         var uint32 _y = (y);                                                    \
  622.         if (_y <= (uint32)(bit(16)-1))                                          \
  623.           { var uint16 _q1;                                                     \
  624.             var uint16 _q0;                                                     \
  625.             var uint16 _r1;                                                     \
  626.             divu_3216_1616(high16(_x),_y, _q1 = , _r1 = );                      \
  627.             divu_3216_1616(highlow32(_r1,low16(_x)),_y, _q0 = , _EMA_ r_zuweisung); \
  628.             q_zuweisung highlow32(_q1,_q0);                                     \
  629.           }                                                                     \
  630.           else                                                                  \
  631.           { var uint32 _x1 = _x; # x1 := x                                      \
  632.             var uint32 _y1 = _y; # y1 := y                                      \
  633.             var uint16 _q;                                                      \
  634.             do { _x1 = floor(_x1,2); _y1 = floor(_y1,2); } # k erhöhen          \
  635.                until (_y1 <= (uint32)(bit(16)-1)); # bis y1 < beta              \
  636.             { var uint16 _y2 = low16(_y1)+1; # y1+1 bilden                      \
  637.               if (_y2==0)                                                       \
  638.                 { _q = high16(_x1); } # y1+1=beta -> ein Shift                  \
  639.                 else                                                            \
  640.                 { divu_3216_1616(_x1,_y2,_q=,_EMA_); } # Division von x1 durch y1+1  \
  641.             }                                                                   \
  642.             # _q = q = floor(x1/(y1+1))                                         \
  643.             # x-q*y bilden (eine 16-mal-32-Bit-Multiplikation ohne Überlauf):   \
  644.             _x -= highlow32_0(mulu16(_q,high16(_y))); # q * high16(y) * beta    \
  645.             # gefahrlos, da q*high16(y) <= q*y/beta <= x/beta < beta            \
  646.             _x -= mulu16(_q,low16(_y)); # q * low16(y)                          \
  647.             # gefahrlos, da q*high16(y)*beta + q*low16(y) = q*y <= x            \
  648.             # Noch höchstens 2 mal y abziehen:                                  \
  649.             if (_x >= _y)                                                       \
  650.               { _q += 1; _x -= _y;                                              \
  651.                 if (_x >= _y)                                                   \
  652.                   { _q += 1; _x -= _y;                                          \
  653.               }   }                                                             \
  654.             r_zuweisung _x;                                                     \
  655.             q_zuweisung (uint32)(_q);                                           \
  656.       }   }
  657.     #ifdef LISPARIT
  658.     # Dies dient nur noch als Hilfsfunktion für arilev1.d.
  659.     # Die Rückgabe des Restes in divu_32_rest ist also hier nicht nötig.
  660.     global uint32 divu_3232_3232_(x,y)
  661.       var reg2 uint32 x;
  662.       var reg1 uint32 y;
  663.       { var reg3 uint32 q;
  664.         divu_3232_3232(x,y,q=,_EMA_);
  665.         return q;
  666.       }
  667.     #endif
  668.   #else
  669.     #define divu_3232_3232(x,y,q_zuweisung,r_zuweisung)  \
  670.       { q_zuweisung divu_3232_3232_(x,y); r_zuweisung divu_32_rest; }
  671.     #if 0
  672.       # divu_3232_3232_ extern in Assembler
  673.     #else
  674.       #ifdef LISPARIT
  675.       global uint32 divu_3232_3232_(x,y)
  676.         var reg2 uint32 x;
  677.         var reg1 uint32 y;
  678.         { var reg3 uint32 q = floor(x,y);
  679.           divu_32_rest = x - q*y;
  680.           return q;
  681.         }
  682.       #endif
  683.     #endif
  684.   #endif
  685.  
  686. # Dividiert eine 64-Bit-Zahl durch eine 32-Bit-Zahl und
  687. # liefert einen 32-Bit-Quotienten und einen 32-Bit-Rest.
  688. # divu_6432_3232(xhi,xlo,y,q=,r=);
  689. # > uint32 xhi,xlo: x = 2^32*xhi+xlo = Zähler
  690. # > uint32 y: Nenner
  691. # > Es sei bekannt, daß 0 <= x < 2^32*y .
  692. # < uint32 q: floor(x/y)
  693. # < uint32 r: x mod y
  694. # < x = q*y+r
  695.   extern uint32 divu_6432_3232_ (uint32 xhi, uint32 xlo, uint32 y); # -> Quotient q
  696.   extern uint32 divu_32_rest;                                       # -> Rest r
  697.   #ifdef GNU
  698.     #if defined(MC680Y0) && !defined(NO_ASM)
  699.       #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
  700.         ({var uint32 __xhi = (xhi);  \
  701.           var uint32 __xlo = (xlo);  \
  702.           var uint32 __y = (y);      \
  703.           var uint32 __q;            \
  704.           var uint32 __r;            \
  705.           __asm__ __volatile__ ("\
  706.             divul %4,%1:%0   \
  707.             " : "=d" (__q), "=d" (__r) : "1" (__xhi), "0" (__xlo), "dm" (__y));  \
  708.           q_zuweisung __q;           \
  709.           r_zuweisung __r;           \
  710.          })
  711.       #define divu_6432_3232_(xhi,xlo,y) \
  712.         ({var reg1 uint32 ___q; divu_6432_3232(xhi,xlo,y,___q=,_EMA_); ___q; })
  713.     #elif defined(SPARC)
  714.       #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
  715.         ({ var reg1 uint32 _q = divu_6432_3232_(xhi,xlo,y); # extern in Assembler \
  716.            var register uint32 _r __asm__("%g1");                                 \
  717.            q_zuweisung _q; r_zuweisung _r;                                        \
  718.          })
  719.     #elif defined(ARM)
  720.       #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
  721.         ({ var reg1 uint32 _q = divu_6432_3232_(xhi,xlo,y); # extern in Assembler \
  722.            var register uint32 _r __asm__("%r1"/*"%a2"*/);                        \
  723.            q_zuweisung _q; r_zuweisung _r;                                        \
  724.          })
  725.     #elif defined(I80Z86) && !defined(NO_ASM)
  726.       #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
  727.         ({var uint32 __xhi = (xhi);  \
  728.           var uint32 __xlo = (xlo);  \
  729.           var uint32 __y = (y);      \
  730.           var uint32 __q;            \
  731.           var uint32 __r;            \
  732.           __asm__ __volatile__ (     \
  733.              "divl %4"               \
  734.              : "=a" /* %eax */ (__q), "=d" /* %edx */ (__r)               \
  735.              : "1" /* %edx */ (__xhi), "0" /* %eax */ (__xlo), "rm" (__y) \
  736.              );                      \
  737.           q_zuweisung __q;           \
  738.           r_zuweisung __r;           \
  739.          })
  740.       #define divu_6432_3232_(xhi,xlo,y) \
  741.         ({var reg1 uint32 ___q; divu_6432_3232(xhi,xlo,y,___q=,_EMA_); ___q; })
  742.     #elif defined(HAVE_LONGLONG) && !defined(HPPA_DIV_WORKS)
  743.       #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung) \
  744.         ({var reg1 uint32 __xhi = (xhi);                           \
  745.           var reg1 uint32 __xlo = (xlo);                           \
  746.           var reg1 uint64 __x = (uint64)__xhi<<32 | (uint64)__xlo; \
  747.           var reg1 uint32 __y = (y);                               \
  748.           var reg1 uint32 __q = floor(__x,(uint64)__y);            \
  749.           q_zuweisung __q; r_zuweisung __xlo - __q * __y;          \
  750.          })
  751.       #define divu_6432_3232_(xhi,xlo,y) \
  752.         ({var reg1 uint32 ___q; divu_6432_3232(xhi,xlo,y,___q=,_EMA_); ___q; })
  753.     #endif
  754.   #endif
  755.   #if defined(WATCOM) && defined(I80Z86) && !defined(NO_ASM)
  756.     #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
  757.       { var uint32 __xhi = (xhi);  \
  758.         var uint32 __xlo = (xlo);  \
  759.         var uint32 __y = (y);      \
  760.         var uint32 __q;            \
  761.         var uint32 __r;            \
  762.         __q = divu_6432_3232_(__xhi,__xlo,__y); __r = divu_6432_3232_rest(); \
  763.         q_zuweisung __q;           \
  764.         r_zuweisung __r;           \
  765.       }
  766.     extern uint32 divu_6432_3232_rest (void);
  767.     #pragma aux divu_6432_3232_ = 0xF7 0xF1 /* divl %ecx */ parm [edx] [eax] [ecx] value [eax] modify [eax edx];
  768.     #pragma aux divu_6432_3232_rest = /* */ value [edx] modify [];
  769.   #endif
  770.   #ifndef divu_6432_3232
  771.     #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
  772.       { q_zuweisung divu_6432_3232_(xhi,xlo,y); r_zuweisung divu_32_rest; }
  773.     #if defined(MC680Y0) || defined(SPARC) || defined(ARM) || (defined(I80Z86) && !defined(WATCOM) && !defined(MICROSOFT)) || defined(HPPA)
  774.       # divu_6432_3232_ extern in Assembler
  775.       #if defined(SPARC)
  776.         #define divu_32_rest  (uint32)(_get_g1()) # Rückgabe im Register %g1
  777.       #elif defined(LISPARIT)
  778.         global uint32 divu_32_rest;
  779.       #endif
  780.     #else
  781.       #ifdef LISPARIT
  782.       # Methode:
  783.       # Wie UDS_divide mit intDsize=16, a_len=4, b_len=2.
  784.       global uint32 divu_32_rest;
  785.       global uint32 divu_6432_3232_(xhi,xlo,y)
  786.         var reg9 uint32 xhi;
  787.         var reg6 uint32 xlo;
  788.         var reg1 uint32 y;
  789.         { if (y <= (uint32)(bit(16)-1))
  790.             # 48-durch-16-Bit-Division,
  791.             # aufgebaut aus zwei 32-durch-16-Bit-Divisionen:
  792.             { var reg4 uint16 q1;
  793.               var reg3 uint16 q0;
  794.               var reg2 uint16 r1;
  795.               divu_3216_1616(highlow32(low16(xhi),high16(xlo)),y, q1=,r1=);
  796.               divu_3216_1616(highlow32(r1,low16(xlo)),y, q0=, divu_32_rest=(uint32) );
  797.               return highlow32(q1,q0);
  798.             }
  799.           # y>=2^16
  800.          {# y shiften:
  801.           var reg10 uintL s = 0;
  802.           while ((sint32)y >= 0) { y = y<<1; s++; }
  803.           # x entsprechend shiften:
  804.           if (!(s==0))
  805.             { xhi = (xhi << s) | (xlo >> (32-s)); xlo = xlo << s; }
  806.           # 64-durch-32-Bit-Division,
  807.           # aufgebaut aus zwei 48-durch-32-Bit-Divisionen.
  808.           # Methode für eine 48-durch-32-Bit-Division x/y mit 0 <= x < 2^16*y :
  809.           # (beta = 2^n = 2^16, n = 16)
  810.           # Wir wissen beta^2/2 <= y < beta^2, Quotient  q = floor(x/y) < beta.
  811.           # Schreibe  x = beta*x1 + x0  mit  x1 := floor(x/beta)
  812.           # und       y = beta*y1 + y0  mit  y1 := floor(y/beta)
  813.           # und bilde den Näherungs-Quotienten floor(x1/y1)
  814.           # oder (noch besser) floor(x1/(y1+1)).
  815.           # Wegen 0 <= x1 < 2^(2n) und 0 < 2^(n-1) <= y1 < 2^n
  816.           # und  x1/(y1+1) <= x/y < x1/(y1+1) + 2
  817.           # (denn x1/(y1+1) = (x1*beta)/((y1+1)*beta) <= (x1*beta)/y <= x/y
  818.           # und x/y - x1/(y1+1) = (x+x*y1-x1*y)/(y*(y1+1))
  819.           # = (x+x0*y1-x1*y0)/(y*(y1+1)) <= (x+x0*y1)/(y*(y1+1))
  820.           # <= x/(y*(y1+1)) + x0/y = (x/y)/(y1+1) + x0/y
  821.           # <= 2^n/(2^(n-1)+1) + 2^n/2^(2n-1) = 2^n/(2^(n-1)+1) + 2^(1-n) < 2 )
  822.           # gilt  floor(x1/(y1+1)) <= floor(x/y) <= floor(x1/(y1+1)) + 2  .
  823.           # Man bildet also  q:=floor(x1/(y1+1))  (ein Shift um n Bit oder
  824.           # eine (2n)-durch-n-Bit-Division, mit Ergebnis q <= floor(x/y) < beta)
  825.           # und x-q*y und muß hiervon noch höchstens 2 mal y abziehen und q
  826.           # incrementieren, um den Quotienten  q = floor(x/y)  und den Rest
  827.           # x-floor(x/y)*y  der Division zu bekommen.
  828.           { var reg2 uint16 y1_1 = high16(y)+1; # y1+1
  829.             var reg7 uint16 q1;
  830.             var reg8 uint16 q0;
  831.             var reg3 uint32 r;
  832.             # 2^16*xhi+high16(xlo) durch y dividieren:
  833.            {var reg5 uint16 r16;
  834.             var reg4 uint32 r2;
  835.             if (y1_1==0)
  836.               { q1 = high16(xhi); r16 = low16(xhi); }
  837.               else
  838.               { divu_3216_1616(xhi,y1_1, q1=,r16=); }
  839.             # q1 = floor(xhi/(y1+1)), r16 = xhi - (y1+1)*q1 (>=0, <=y1)
  840.             # Bilde r := (2^16*xhi+high16(xlo)) - y*q1
  841.             #          = 2^16*(xhi-y1*q1) + high16(xlo) - y0*q1
  842.             #          = 2^16*r16 + 2^16*q1 + high16(xlo) - y0*q1 (>=0)
  843.             # Dies ist < 2^16*y1 + 2^32 <= y + 2^32 <= 3*y, kann überlaufen!
  844.             r = highlow32(r16,high16(xlo)); # 2^16*r16 + high16(xlo) < 2^32
  845.             r2 = highlow32_0(q1) - mulu16(low16(y),q1); # 2^16*q1 - y0*q1 < 2^32
  846.             # 0 <= r+r2 < 3*y. Bei der Addition auf Carry testen!
  847.             # Carry -> jedenfalls y <= r+r2 < y + 2^32 <= 3*y.
  848.             # kein Carry -> jedenfalls 0 <= r+r2 < 2^32 <= 2*y.
  849.             if ((r += r2) < r2) # addieren, r >= 2^32 ?
  850.               { q1 += 1; r -= y; }
  851.             # jetzt noch 0 <= r < 2^32 <= 2*y
  852.             if (r >= y)
  853.               { q1 += 1; r -= y; }
  854.            }# Quotient q1, Rest r fertig.
  855.             # 2^16*r+low16(xlo) durch y dividieren:
  856.            {var reg5 uint16 r16;
  857.             var reg4 uint32 r2;
  858.             if (y1_1==0)
  859.               { q0 = high16(r); r16 = low16(r); }
  860.               else
  861.               { divu_3216_1616(r,y1_1, q0=,r16=); }
  862.             # q0 = floor(r/(y1+1)), r16 = r - (y1+1)*q0 (>=0, <=y1)
  863.             # Bilde r := (2^16*r+low16(xlo)) - y*q0
  864.             #          = 2^16*(r-y1*q0) + low16(xlo) - y0*q0
  865.             #          = 2^16*r16 + 2^16*q0 + low16(xlo) - y0*q0 (>=0)
  866.             # Dies ist < 2^16*y1 + 2^32 <= y + 2^32 <= 3*y, kann überlaufen!
  867.             r = highlow32(r16,low16(xlo)); # 2^16*r16 + low16(xlo) < 2^32
  868.             r2 = highlow32_0(q0) - mulu16(low16(y),q0); # 2^16*q0 - y0*q0 < 2^32
  869.             # 0 <= r+r2 < 3*y. Bei der Addition auf Carry testen!
  870.             # Carry -> jedenfalls y <= r+r2 < y + 2^32 <= 3*y.
  871.             # kein Carry -> jedenfalls 0 <= r+r2 < 2^32 <= 2*y.
  872.             if ((r += r2) < r2) # addieren, r >= 2^32 ?
  873.               { q0 += 1; r -= y; }
  874.             # jetzt noch 0 <= r < 2^32 <= 2*y
  875.             if (r >= y)
  876.               { q0 += 1; r -= y; }
  877.            }# Quotient q0, Rest r fertig.
  878.             divu_32_rest = r >> s; # Rest
  879.             return highlow32(q1,q0); # Quotient
  880.         }}}
  881.       #endif
  882.     #endif
  883.   #endif
  884.  
  885. # Zieht die Ganzzahl-Wurzel aus einer 32-Bit-Zahl und
  886. # liefert eine 16-Bit-Wurzel und einen Rest.
  887. # isqrt_32_16(x,y=,sqrtp=);
  888. # > uint32 x: Radikand, >= 2^30, < 2^32
  889. # < uint16 y: floor(sqrt(x)), >= 2^15, < 2^16
  890. # < boolean sqrtp: /=0, falls x=y^2
  891.   # Methode:
  892.   # y := 2^16 als Anfangswert,
  893.   # y := floor((y + floor(x/y))/2) als nächster Wert,
  894.   # solange z := floor(x/y) < y, setze y := floor((y+z)/2).
  895.   # y ist fertig; x=y^2 genau dann, wenn z=y und die letzte Division aufging.
  896.   # (Beweis:
  897.   #  1. Die Folge der y ist streng monoton fallend.
  898.   #  2. Stets gilt y >= floor(sqrt(x)) (denn für alle y>0 ist
  899.   #     y + x/y >= 2*sqrt(x) und daher  floor((y + floor(x/y))/2) =
  900.   #     floor(y/2 + x/(2*y)) >= floor(sqrt(x)) ).
  901.   #  3. Am Schluß gilt x >= y^2.
  902.   # )
  903.   #define isqrt_32_16(x,y_zuweisung,sqrtp_zuweisung)  \
  904.     { var reg4 uint32 _x = (x);                                          \
  905.       var reg3 uint16 _x1 = high16(_x);                                  \
  906.       var reg1 uint16 _y = floor(_x1,2) | bit(16-1);                     \
  907.       loop                                                               \
  908.         { var reg2 uint16 _z;                                            \
  909.           var reg5 uint16 _r;                                            \
  910.           if (_x1 >= _y) # Division _x/_y ergäbe Überlauf -> _z > _y     \
  911.             { unused (sqrtp_zuweisung FALSE); break; }                   \
  912.           divu_3216_1616(_x,_y, _z=,_r=); # Dividiere _x/_y              \
  913.           if (_z >= _y)                                                  \
  914.             { unused (sqrtp_zuweisung (_z == _y) && (_r == 0)); break; } \
  915.           _y = floor((uint16)(_z+_y),2) | bit(16-1); # _y muß >= 2^15 bleiben  \
  916.         }                                                                \
  917.       y_zuweisung _y;                                                    \
  918.     }
  919.  
  920. # Zieht die Ganzzahl-Wurzel aus einer 64-Bit-Zahl und
  921. # liefert eine 32-Bit-Wurzel und einen Rest.
  922. # isqrt_64_32(xhi,xlo,y=,sqrtp=);
  923. # > uint32 xhi,xlo: Radikand x = 2^32*xhi+xlo, >= 2^62, < 2^64
  924. # < uint32 y: floor(sqrt(x)), >= 2^31, < 2^32
  925. # < boolean sqrtp: /=0, falls x=y^2
  926.   #if (defined(SPARC) || defined(MC680Y0) || defined(HPPA))
  927.     # Methode:
  928.     # y := 2^32 als Anfangswert,
  929.     # y := floor((y + floor(x/y))/2) als nächster Wert,
  930.     # solange z := floor(x/y) < y, setze y := floor((y+z)/2).
  931.     # y ist fertig; x=y^2 genau dann, wenn z=y und die letzte Division aufging.
  932.     # (Beweis:
  933.     #  1. Die Folge der y ist streng monoton fallend.
  934.     #  2. Stets gilt y >= floor(sqrt(x)) (denn für alle y>0 ist
  935.     #     y + x/y >= 2*sqrt(x) und daher  floor((y + floor(x/y))/2) =
  936.     #     floor(y/2 + x/(2*y)) >= floor(sqrt(x)) ).
  937.     #  3. Am Schluß gilt x >= y^2.
  938.     # )
  939.     #define isqrt_64_32(xhi,xlo,y_zuweisung,sqrtp_zuweisung)  \
  940.       { var reg3 uint32 _xhi = (xhi);                                   \
  941.         var reg4 uint32 _xlo = (xlo);                                   \
  942.         var reg1 uint32 _y = floor(_xhi,2) | bit(32-1);                 \
  943.         loop                                                            \
  944.           { var reg2 uint32 _z;                                         \
  945.             var reg5 uint32 _rest;                                      \
  946.             if (_xhi >= _y) # Division _x/_y ergäbe Überlauf -> _z > _y \
  947.               { sqrtp_zuweisung FALSE; break; }                         \
  948.             divu_6432_3232(_xhi,_xlo,_y, _z=,_rest=); # Dividiere _x/_y \
  949.             if (_z >= _y)                                               \
  950.               { sqrtp_zuweisung (_z == _y) && (_rest == 0); break; }    \
  951.             _y = floor(_z+_y,2) | bit(32-1); # _y muß >= 2^31 bleiben   \
  952.           }                                                             \
  953.         y_zuweisung _y;                                                 \
  954.       }
  955.   #else
  956.     # Methode:
  957.     # Wie bei UDS_sqrt mit n=2.
  958.     # y = 2^16*yhi + ylo ansetzen.
  959.     # Dann muß
  960.     #   yhi = floor(y/2^16) = floor(floor(sqrt(x))/2^16)
  961.     #       = floor(sqrt(x)/2^16) = floor(sqrt(x/2^32)) = isqrt(xhi)
  962.     # sein. Es folgt yhi >= 2^15.
  963.     # Danach sucht man das größte ylo >=0 mit
  964.     # x - 2^32*yhi^2 >= 2*2^16*yhi*ylo + ylo^2.
  965.     # Dazu setzen wir  xhi*2^32+xlo := x - 2^32*yhi^2
  966.     # (also xhi := xhi - yhi^2, das ist >=0, <=2*yhi).
  967.     # Die Schätzung für die zweite Ziffer
  968.     #     ylo' := min(2^16-1,floor((xhi*2^32+xlo)/(2*2^16*yhi)))
  969.     # erfüllt ylo'-1 <= ylo <= ylo', ist also um höchstens 1 zu groß.
  970.     # (Beweis: Rechte Ungleichung klar, da  ylo < 2^16  und
  971.     #   xhi*2^32+xlo >= 2*2^16*yhi*ylo + ylo^2 >= 2*2^16*yhi*ylo
  972.     #   ==> (xhi*2^32+xlo)/(2*2^16*yhi) >= ylo  gelten muß.
  973.     #   Linke Ungleichung: Falls floor(...)>=2^16, ist
  974.     #   xhi*2^32+xlo >= 2*2^16*2^16*yhi >= 2*2^16*yhi*(2^16-1) + 2^32
  975.     #                >= 2*2^16*yhi*(2^16-1) + (2^16-1)^2
  976.     #   und xhi*2^32+xlo < 2*2^16*2^16*yhi + (2^16)^2, also
  977.     #   ylo = 2^16-1 = ylo'.
  978.     #   Sonst ist ylo' = floor((xhi*2^32+xlo)/(2*2^16*yhi)), also
  979.     #   xhi*2^32+xlo >= 2*2^16*yhi*ylo' >= 2*2^16*yhi*(ylo'-1) + 2^32
  980.     #                >= 2*2^16*yhi*(ylo'-1) + (ylo'-1)^2,
  981.     #   also ylo >= ylo'-1 nach Definition von ylo.)
  982.     #define isqrt_64_32(xhi,xlo,y_zuweisung,sqrtp_zuweisung)  \
  983.       { var reg4 uint32 _xhi = (xhi);                                        \
  984.         var reg3 uint32 _xlo = (xlo);                                        \
  985.         var reg6 uint16 _yhi;                                                \
  986.         var reg5 uint16 _ylo;                                                \
  987.         # erste Ziffer berechnen:                                            \
  988.         isqrt_32_16(_xhi,_yhi=,_EMA_); # yhi := isqrt(xhi)                        \
  989.         _xhi -= mulu16(_yhi,_yhi); # jetzt 0 <= xhi <= 2*yhi                 \
  990.         # x = 2^32*yhi^2 + 2^32*xhi + xlo                                    \
  991.         # Schätzung für die zweite Ziffer berechnen:                         \
  992.         # ylo := min(2^16-1,floor((xhi*2^32+xlo)/(2*2^16*yhi))) bilden:      \
  993.        {var reg1 uint32 _z = (_xhi << 15) | (_xlo >> 17); # < 2^15*(2*yhi+1) \
  994.         var reg2 uint32 _r = highlow32_0(_yhi);                              \
  995.         if (_z >= _r)                                                        \
  996.           { _ylo = bit(16)-1; _r = _z - _r + (uint32)_yhi; }                 \
  997.           else                                                               \
  998.           { divu_3216_1616(_z,_yhi, _ylo=,_r=); }                            \
  999.         # x = 2^32*yhi^2 + 2*2^16*yhi*ylo + 2^17*r + (xlo mod 2^17),         \
  1000.         # 0 <= r < yhi + 2^15                                                \
  1001.         _xlo = (_r << 17) | (_xlo & (bit(17)-1));                            \
  1002.         # x = 2^32*yhi^2 + 2*2^16*yhi*ylo + 2^32*floor(r/2^15) + xlo         \
  1003.         _z = mulu16(_ylo,_ylo); # z = ylo^2                                  \
  1004.         # Versuche vom Rest 2^32*floor(r/2^15) + xlo  z zu subtrahieren.     \
  1005.         # Falls Rest >= z (d.h. r>=2^15 oder xlo>=z), ist ylo fertig,        \
  1006.         # und es gilt x=y^2 genau dann, wenn r<2^15 und xlo=z.               \
  1007.         # Sonst (d.h. r<2^15 und xlo<z), muß man ylo erniedrigen. Dazu       \
  1008.         # setzt man  ylo := ylo-1, z := z-(2*ylo+1),                         \
  1009.         # Rest := Rest + 2^17*yhi = xlo + 2^17*yhi >= 2^32 > z, also x>y^2.  \
  1010.         if (_r < bit(15))                                                    \
  1011.           { if (_xlo < _z)                                                   \
  1012.               { _ylo -= 1; sqrtp_zuweisung FALSE; }                          \
  1013.               else                                                           \
  1014.               { sqrtp_zuweisung (_xlo == _z); }                              \
  1015.           }                                                                  \
  1016.           else                                                               \
  1017.           { sqrtp_zuweisung FALSE; }                                         \
  1018.         y_zuweisung highlow32(_yhi,_ylo);                                    \
  1019.       }}
  1020.   #endif
  1021.  
  1022. # Eine 32-Bit-Zahl aus zwei aufeinanderfolgenden 16-Bit-Digits einer UDS
  1023. # zusammensetzen: highlow32_at(ptr)
  1024.   #if BIG_ENDIAN_P && defined(MC680X0)
  1025.     # ptr als 32-Bit-Pointer auffassen und darauf zugreifen
  1026.     #define highlow32_at(ptr)  (*(uint32*)(ptr))
  1027.   #else
  1028.     #define highlow32_at(ptr)  highlow32(((uint16*)(ptr))[0],((uint16*)(ptr))[1])
  1029.   #endif
  1030.  
  1031. # Eine 32-Bit-Zahl in zwei aufeinanderfolgende 16-Bit-Digits einer UDS abspeichern:
  1032. # set_highlow32_at(ptr,value32); wobei ptr und value32 Variablen.
  1033.   #if BIG_ENDIAN_P && defined(MC680X0)
  1034.     # ptr als 32-Bit-Pointer auffassen und darauf zugreifen
  1035.     #define set_highlow32_at(ptr,value32)  (*(uint32*)(ptr)=(value32))
  1036.   #else
  1037.     #define set_highlow32_at(ptr,value32)  (((uint16*)(ptr))[0]=high16(value32),((uint16*)(ptr))[1]=low16(value32))
  1038.   #endif
  1039.  
  1040.