src/base/cl_low.h

   1 // Low-level arithmetic: operations on 16-bit and 32-bit words
   2
   3 #ifndef _CL_LOW_H
   4 #define _CL_LOW_H
   5
   6 namespace cln {
   7
   8 // Determines the sign of a 16-bit number.
   9 // sign_of(wert)
  10 // > wert: eine 16-Bit-Zahl
  11 // < sint16 ergebnis: 0 falls wert>=0, -1 falls wert<0.
  12 inline sint16 sign_of (sint16 wert)
  13 {
  14 #if defined(__sparc64__)
  15         return (sint64)wert >> 63;
  16 #elif defined(__sparc__) || defined(__arm__)
  17         return (sint32)wert >> 31;
  18 #else
  19         return (wert >= 0 ? 0 : -1);
  20 #endif
  21 }
  22
  23 // Determines the sign of a 32-bit number.
  24 // sign_of(wert)
  25 // > wert: eine 32-Bit-Zahl
  26 // < sint32 ergebnis: 0 falls wert>=0, -1 falls wert<0.
  27 inline sint32 sign_of (sint32 wert)
  28 {
  29 #if defined(__sparc64__)
  30         return (sint64)wert >> 63;
  31 #elif defined(__sparc__) || defined(__arm__)
  32         return wert >> 31;
  33 #else
  34         return (wert >= 0 ? 0 : -1);
  35 #endif
  36 }
  37
  38 #ifdef HAVE_FAST_LONGLONG
  39
  40 // Determines the sign of a 64-bit number.
  41 // sign_of(wert)
  42 // > wert: eine 64-Bit-Zahl
  43 // < sint64 ergebnis: 0 falls wert>=0, -1 falls wert<0.
  44 inline sint64 sign_of (sint64 wert)
  45 {
  46         return wert >> 63;
  47 }
  48
  49 #endif /* HAVE_FAST_LONGLONG */
  50
  51
  52 // High-Word einer 32-Bit-Zahl bestimmen
  53 // high16(wert)
  54 inline uint16 high16 (uint32 wert)
  55 {
  56         return wert >> 16;
  57 }
  58
  59 // Low-Word einer 32-Bit-Zahl bestimmen
  60 // low16(wert)
  61 inline uint16 low16 (uint32 wert)
  62 {
  63         return (uint16)wert;
  64 }
  65
  66 // Eine 32-Bit-Zahl aus ihrem High-Word und ihrem Low-Word bestimmen:
  67 // highlow32(uint16 high, uint16 low)
  68 inline uint32 highlow32 (uint16 high, uint16 low)
  69 {
  70         return ((uint32)high << 16) | (uint32)low;
  71 }
  72
  73 // Eine 32-Bit-Zahl aus ihrem High-Word und ihrem Low-Word 0 bestimmen:
  74 // highlow32_0(uint16 high)
  75 inline uint32 highlow32_0 (uint16 high)
  76 {
  77         return (uint32)high << 16;
  78 }
  79
  80 // High-Word einer 64-Bit-Zahl bestimmen
  81 // high32(wert)
  82 inline uint32 high32 (uint64 wert)
  83 {
  84         return wert >> 32;
  85 }
  86
  87 // Low-Word einer 64-Bit-Zahl bestimmen
  88 // low32(wert)
  89 inline uint32 low32 (uint64 wert)
  90 {
  91         return (uint32)wert;
  92 }
  93
  94 // Eine 64-Bit-Zahl aus ihrem High-Word und ihrem Low-Word bestimmen:
  95 // highlow64(uint32 high, uint32 low)
  96 inline uint64 highlow64 (uint32 high, uint32 low)
  97 {
  98         return ((uint64)high << 32) | (uint64)low;
  99 }
 100
 101 // Eine 64-Bit-Zahl aus ihrem High-Word und ihrem Low-Word 0 bestimmen:
 102 // highlow64_0(uint32 high)
 103 inline uint64 highlow64_0 (uint32 high)
 104 {
 105         return (uint64)high << 32;
 106 }
 107
 108
 109 // Multipliziert zwei 16-Bit-Zahlen miteinander und liefert eine 32-Bit-Zahl:
 110 // mulu16(arg1,arg2)
 111 // > arg1, arg2 : zwei 16-Bit-Zahlen
 112 // < ergebnis: eine 32-Bit-Zahl
 113 #if defined(__GNUC__) && defined(__sparc__) && !defined(__sparc64__) && defined(FAST_DOUBLE)
 114 // Ist das schneller als mulu16_ ??
 115 inline uint32 mulu16 (uint16 arg1, uint16 arg2)
 116 {
 117         union { double f; uint32 i[2]; } __fi;
 118         __fi.f = (double)(sint32)arg1 * (double)(sint32)arg2
 119                  + (double)(4503599627370496.0L); // + 2^52, zum Normalisieren
 120         return __fi.i[1]; // untere 32 Bit herausholen (benutzt CL_CPU_BIG_ENDIAN_P !)
 121 }
 122 #elif defined(__GNUC__) && defined(__sparc64__) && !defined(NO_ASM)
 123 inline uint32 mulu16 (uint16 arg1, uint16 arg2)
 124 {
 125         uint64 _prod;
 126         __asm__("umul %1,%2,%0"
 127                 : "=r" (_prod)
 128                 : "r" (arg1), "r" (arg2)
 129                );
 130         return _prod;
 131 }
 132 #elif defined(__GNUC__) && (defined(__i386__) || defined(__x86_64__)) && !defined(NO_ASM)
 133 inline uint32 mulu16 (uint16 arg1, uint16 arg2)
 134 {
 135         uint16 _hi;
 136         uint16 _lo;
 137         __asm__("mulw %2"
 138                 : "=d" /* %dx */ (_hi), "=a" /* %ax */ (_lo)
 139                 : "rm" (arg1), "1" /* %eax */ (arg2)
 140                );
 141         return highlow32(_hi,_lo);
 142 }
 143 #elif (defined(__sparc__) || defined(__sparc64__)) && !defined(NO_ASM)
 144   extern "C" uint32 mulu16_ (uint16 arg1, uint16 arg2);
 145   #define mulu16  mulu16_  // extern in Assembler
 146 #else
 147 inline uint32 mulu16 (uint16 arg1, uint16 arg2)
 148 {
 149         return arg1 * arg2;
 150 }
 151 #endif
 152
 153 // Multipliziert zwei 24-Bit-Zahlen zusammen und liefert eine 48-Bit-Zahl.
 154 // mulu24(arg1,arg2,hi=,lo=);
 155 // > arg1, arg2 : zwei 24-Bit-Zahlen
 156 // < 2^32*hi+lo : eine 48-Bit-Zahl
 157 #if defined(__sparc__) && !defined(__sparc64__) && defined(FAST_DOUBLE)
 158   #define mulu24(x,y,hi_zuweisung,lo_zuweisung)  \
 159     { var uint32 _x = (x);                                      \
 160       var uint32 _y = (y);                                      \
 161       var union { double f; uint32 i[2]; uint16 s[4]; } __fi;   \
 162       __fi.f = (double)(sint32)(_x)*(double)(sint32)(_y)        \
 163                + (double)(4503599627370496.0L); /* + 2^52, zum Normalisieren */\
 164       cl_unused (hi_zuweisung __fi.s[1]); /* mittlere 16 Bit herausholen, (benutzt CL_CPU_BIG_ENDIAN_P !) */\
 165       lo_zuweisung __fi.i[1]; /* untere 32 Bit herausholen (benutzt CL_CPU_BIG_ENDIAN_P !)    */\
 166     }
 167 #else
 168   #define mulu24  mulu32
 169 #endif
 170
 171 // Multipliziert zwei 32-Bit-Zahlen miteinander und liefert eine 32-Bit-Zahl:
 172 // mulu32_unchecked(arg1,arg2)
 173 // > arg1, arg2 : zwei 32-Bit-Zahlen
 174 // < ergebnis : eine 32-Bit-Zahl
 175 // Es wird vorausgesetzt, daß arg1*arg2 < 2^32.
 176 #if defined(__GNUC__) && defined(__sparc64__) && !defined(NO_ASM)
 177 inline uint32 mulu32_unchecked (uint32 arg1, uint32 arg2)
 178 {
 179         uint64 _prod;
 180         __asm__("umul %1,%2,%0"
 181                 : "=r" (_prod)
 182                 : "r" (arg1), "r" (arg2)
 183                );
 184         return _prod;
 185 }
 186 #elif defined(__sparc__) && !defined(NO_ASM)
 187   extern "C" uint32 mulu32_unchecked (uint32 x, uint32 y); // extern in Assembler
 188 #else
 189   // Wir können dafür auch die Bibliotheksroutine des C-Compilers nehmen:
 190   inline uint32 mulu32_unchecked (uint32 arg1, uint32 arg2)
 191   {
 192         return arg1 * arg2;
 193   }
 194 #endif
 195
 196 #ifdef __arm__
 197 // Returns the r0 part of a 64-bit (uint64) return value returned in r0,r1.
 198 // retval64_r0(value)
 199 // Returns the r1 part of a 64-bit (uint64) return value returned in r0,r1.
 200 // retval64_r1(value)
 201   #ifdef __ARMEL__
 202     /* little-endian ARM: uint64 retval = 2^32*r1+r0 */
 203     #define retval64_r0(value)  ((uint32)(uint64)(value))
 204     #define retval64_r1(value)  ((uint32)((uint64)(value)>>32))
 205   #else
 206     /* big-endian ARM: uint64 retval = 2^32*r0+r1 */
 207     #define retval64_r0(value)  ((uint32)((uint64)(value)>>32))
 208     #define retval64_r1(value)  ((uint32)(uint64)(value))
 209   #endif
 210 #endif
 211
 212 // Multipliziert zwei 32-Bit-Zahlen miteinander und liefert eine 64-Bit-Zahl:
 213 // mulu32(arg1,arg2,hi=,lo=);
 214 // > arg1, arg2 : zwei 32-Bit-Zahlen
 215 // < 2^32*hi+lo : eine 64-Bit-Zahl
 216 #if defined(__GNUC__) && defined(__arm__) && !defined(NO_ASM)
 217   extern "C" uint64 mulu32_ (uint32 arg1, uint32 arg2);
 218 #else
 219   extern "C" uint32 mulu32_ (uint32 arg1, uint32 arg2); // -> Low-Teil
 220 #endif
 221 #ifdef _MSC_VER
 222   // Workaround MSVC compiler bug: extern "C" results in wrong symbols, when
 223   // declared inside a namespace!
 224 } extern "C" uint32 mulu32_high; namespace cln {        // -> High-Teil
 225 #else
 226   extern "C" uint32 mulu32_high;                        // -> High-Teil
 227 #endif
 228 #if defined(__GNUC__) && defined(__m68k__) && !defined(NO_ASM)
 229   #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
 230     ({ var uint32 _x = (x);       \
 231        var uint32 _y = (y);       \
 232        var uint32 _hi;            \
 233        var uint32 _lo;            \
 234        __asm__("mulul %3,%0:%1" : "=d" (_hi), "=d"(_lo) : "1" (_x), "dm" (_y) ); \
 235        cl_unused (hi_zuweisung _hi); \
 236        lo_zuweisung _lo;          \
 237      })
 238 #elif defined(__GNUC__) && defined(__m68k__)
 239   #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
 240     ({ var uint32 _x = (x);                                             \
 241        var uint32 _y = (y);                                             \
 242        var uint16 _x1 = high16(_x);                                     \
 243        var uint16 _x0 = low16(_x);                                      \
 244        var uint16 _y1 = high16(_y);                                     \
 245        var uint16 _y0 = low16(_y);                                      \
 246        var uint32 _hi = mulu16(_x1,_y1); /* obere Portion */            \
 247        var uint32 _lo = mulu16(_x0,_y0); /* untere Portion */           \
 248        {var uint32 _mid = mulu16(_x0,_y1); /* 1. mittlere Portion */    \
 249         _hi += high16(_mid); _mid = highlow32_0(low16(_mid));           \
 250         _lo += _mid; if (_lo < _mid) { _hi += 1; } /* 64-Bit-Addition */\
 251        }                                                                \
 252        {var uint32 _mid = mulu16(_x1,_y0); /* 2. mittlere Portion */    \
 253         _hi += high16(_mid); _mid = highlow32_0(low16(_mid));           \
 254         _lo += _mid; if (_lo < _mid) { _hi += 1; } /* 64-Bit-Addition */\
 255        }                                                                \
 256        cl_unused (hi_zuweisung _hi);                                    \
 257        lo_zuweisung _lo;                                                \
 258      })
 259 #elif defined(__GNUC__) && defined(__sparc64__) && !defined(NO_ASM)
 260   #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
 261     ({ var uint64 _prod;                                \
 262        __asm__("umul %1,%2,%0"                                  \
 263                : "=r" (_prod)                                   \
 264                : "r" ((uint32)(x)), "r" ((uint32)(y))           \
 265               );                                                \
 266        cl_unused (hi_zuweisung (uint32)(_prod>>32));            \
 267        lo_zuweisung (uint32)(_prod);                            \
 268      })
 269 #elif defined(__GNUC__) && defined(__sparc__) && !defined(NO_ASM)
 270   #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
 271     ({ lo_zuweisung mulu32_(x,y); /* extern in Assembler */     \
 272       {var uint32 _hi __asm__("%g1");                   \
 273        cl_unused (hi_zuweisung _hi);                            \
 274      }})
 275 #elif defined(__GNUC__) && defined(__arm__) && !defined(NO_ASM)
 276   #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
 277     ({ var register uint64 _prod = mulu32_(x,y);        \
 278        hi_zuweisung retval64_r1(_prod);         \
 279        lo_zuweisung retval64_r0(_prod);         \
 280      })
 281 #elif defined(__GNUC__) && (defined(__i386__) || defined(__x86_64__)) && !defined(NO_ASM)
 282   #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
 283     ({ var uint32 _hi;                                  \
 284        var uint32 _lo;                                  \
 285        __asm__("mull %2"                                         \
 286                : "=d" /* %edx */ (_hi), "=a" /* %eax */ (_lo)    \
 287                : "g" ((uint32)(x)), "1" /* %eax */ ((uint32)(y)) \
 288               );                                                 \
 289        cl_unused (hi_zuweisung _hi); lo_zuweisung _lo;              \
 290      })
 291 #elif defined(__GNUC__) && defined(__mips__) && !defined(NO_ASM)
 292   #if __mips_isa_rev >= 6
 293     #define MULTU_HI_LO "mulu %1,%3,%2 ; muhu %0,%3,%2"
 294   #else
 295     #define MULTU_HI_LO "multu %3,%2 ; mfhi %0 ; mflo %1"
 296   #endif
 297   #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
 298     ({ var uint32 _hi;                       \
 299        var uint32 _lo;                       \
 300        __asm__(MULTU_HI_LO                            \
 301                : "=r" (_hi), "=r" (_lo)               \
 302                : "r" ((uint32)(x)), "r" ((uint32)(y)) \
 303               );                                      \
 304        cl_unused (hi_zuweisung _hi); lo_zuweisung _lo;   \
 305      })
 306 #elif defined(__GNUC__) && !defined(__arm__)
 307   #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
 308     ({ var uint64 _prod = (uint64)(uint32)(x) * (uint64)(uint32)(y); \
 309        cl_unused (hi_zuweisung (uint32)(_prod>>32));                             \
 310        lo_zuweisung (uint32)(_prod);                                          \
 311      })
 312 #else
 313   #define mulu32(x,y,hi_zuweisung,lo_zuweisung)  \
 314     { lo_zuweisung mulu32_(x,y); cl_unused (hi_zuweisung mulu32_high); }
 315   #if (defined(__m68k__) || defined(__sparc__) || defined(__sparc64__) || defined(__arm__) || (defined(__i386__) && !defined(MICROSOFT)) || defined(__x86_64__) || defined(__mips__) || defined(__hppa__)) && !defined(NO_ASM)
 316     // mulu32_ extern in Assembler
 317     #if defined(__sparc__) || defined(__sparc64__)
 318       extern "C" uint32 _get_g1 (void);
 319       #define mulu32_high  (_get_g1()) // Rückgabe im Register %g1
 320     #elif !defined(__hppa__)
 321       #define NEED_VAR_mulu32_high
 322     #endif
 323   #else
 324     #define NEED_FUNCTION_mulu32_
 325   #endif
 326 #endif
 327
 328 #ifdef HAVE_FAST_LONGLONG
 329
 330 // Multipliziert zwei 32-Bit-Zahlen miteinander und liefert eine 64-Bit-Zahl:
 331 // mulu32_w(arg1,arg2)
 332 // > arg1, arg2 : zwei 32-Bit-Zahlen
 333 // < result : eine 64-Bit-Zahl
 334 #if defined(__GNUC__) && defined(__sparc64__) && !defined(NO_ASM)
 335   // Prefer the umul instruction over the mulx instruction (overkill).
 336   #define mulu32_w(x,y)  \
 337     ({ var uint64 _prod;                                \
 338        __asm__("umul %1,%2,%0"                                  \
 339                : "=r" (_prod)                                   \
 340                : "r" ((uint32)(x)), "r" ((uint32)(y))           \
 341               );                                                \
 342        _prod;                                                   \
 343      })
 344 #elif defined(__GNUC__)
 345   #define mulu32_w(x,y)  ((uint64)(uint32)(x) * (uint64)(uint32)(y))
 346 #else
 347   extern "C" uint64 mulu32_w (uint32 arg1, uint32 arg2);
 348   #define NEED_FUNCTION_mulu32_w
 349 #endif
 350
 351 // Multipliziert zwei 64-Bit-Zahlen miteinander und liefert eine 128-Bit-Zahl:
 352 // mulu64(arg1,arg2,hi=,lo=);
 353 // > arg1, arg2 : zwei 64-Bit-Zahlen
 354 // < 2^64*hi+lo : eine 128-Bit-Zahl
 355   extern "C" uint64 mulu64_ (uint64 arg1, uint64 arg2); // -> Low-Teil
 356 #ifdef _MSC_VER
 357   // Workaround MSVC compiler bug.
 358 } extern "C" uint64 mulu64_high; namespace cln {        // -> High-Teil
 359 #else
 360   extern "C" uint64 mulu64_high;                        // -> High-Teil
 361 #endif
 362 #if defined(__GNUC__) && defined(__alpha__) && !defined(NO_ASM)
 363   #define mulu64(x,y,hi_zuweisung,lo_zuweisung)  \
 364     ({ var uint64 _x = (x);     \
 365        var uint64 _y = (y);     \
 366        var uint64 _hi;          \
 367        var uint64 _lo;          \
 368        __asm__("mulq %1,%2,%0"          \
 369                : "=r" (_lo)             \
 370                : "r" (_x), "r" (_y)     \
 371               );                        \
 372        __asm__("umulh %1,%2,%0"         \
 373                : "=r" (_hi)             \
 374                : "r" (_x), "r" (_y)     \
 375               );                        \
 376        cl_unused (hi_zuweisung _hi);    \
 377        lo_zuweisung _lo;                \
 378      })
 379 #elif defined(__GNUC__) && defined(__sparc64__) && !defined(NO_ASM)
 380   #define mulu64(x,y,hi_zuweisung,lo_zuweisung)  \
 381     ({ lo_zuweisung mulu64_(x,y); /* extern in Assembler */     \
 382       {var uint64 _hi __asm__("%g2");                   \
 383        cl_unused (hi_zuweisung _hi);                            \
 384      }})
 385 #elif defined(__GNUC__) && defined(__x86_64__) && !defined(NO_ASM)
 386   #define mulu64(x,y,hi_zuweisung,lo_zuweisung)  \
 387     ({ var uint64 _hi;                                   \
 388        var uint64 _lo;                                   \
 389        __asm__("mulq %2"                                          \
 390                : "=d" /* %rdx */ (_hi), "=a" /* %rax */ (_lo)     \
 391                : "rm" ((uint64)(x)), "1" /* %rax */ ((uint64)(y)) \
 392               );                                                  \
 393        cl_unused (hi_zuweisung _hi); lo_zuweisung _lo;               \
 394      })
 395 #elif defined(__GNUC__) && defined(__ia64__) && !defined(NO_ASM)
 396   #define mulu64(x,y,hi_zuweisung,lo_zuweisung)  \
 397     ({ var uint64 _x = (x);                                         \
 398        var uint64 _y = (y);                                         \
 399        var uint64 _hi;                                              \
 400        __asm__("xma.hu %0 = %1, %2, f0"                                     \
 401                : "=f" (_hi)                                                 \
 402                : "f" ((uint64)(_x)), "f" ((uint64)(_y))                     \
 403               );                                                            \
 404        cl_unused (hi_zuweisung _hi); lo_zuweisung ((uint64)(_x)*(uint64)(_y)); \
 405      })
 406 #else
 407   #define mulu64(x,y,hi_zuweisung,lo_zuweisung)  \
 408     { lo_zuweisung mulu64_(x,y); cl_unused (hi_zuweisung mulu64_high); }
 409   #if defined(__sparc64__) && !defined(NO_ASM)
 410     // mulu64_ extern in Assembler
 411     extern "C" uint64 _get_g2 (void);
 412     #define mulu64_high  (_get_g2()) // Rückgabe im Register %g2
 413   #else
 414     #define NEED_FUNCTION_mulu64_
 415   #endif
 416 #endif
 417
 418 #endif /* HAVE_FAST_LONGLONG */
 419
 420
 421 // Dividiert eine 16-Bit-Zahl durch eine 16-Bit-Zahl und
 422 // liefert einen 16-Bit-Quotienten und einen 16-Bit-Rest.
 423 // divu_1616_1616(x,y,q=,r=);
 424 // > uint16 x: Zähler
 425 // > uint16 y: Nenner
 426 // < uint16 q: floor(x/y)
 427 // < uint16 r: x mod y
 428 // < x = q*y+r
 429   #define divu_1616_1616(x,y,q_zuweisung,r_zuweisung)  \
 430     { var uint16 __x = (x);                                     \
 431       var uint16 __y = (y);                                     \
 432       q_zuweisung floor(__x,__y);                               \
 433       r_zuweisung (__x % __y);                                  \
 434     }
 435
 436 // Dividiert eine 32-Bit-Zahl durch eine 16-Bit-Zahl und
 437 // liefert einen 16-Bit-Quotienten und einen 16-Bit-Rest.
 438 // divu_3216_1616(x,y,q=,r=);
 439 // > uint32 x: Zähler
 440 // > uint16 y: Nenner
 441 // > Es sei bekannt, daß 0 <= x < 2^16*y .
 442 // < uint16 q: floor(x/y)
 443 // < uint16 r: x mod y
 444 // < x = q*y+r
 445 #if defined(__sparc__)
 446   extern "C" uint32 divu_3216_1616_ (uint32 x, uint16 y); // -> Quotient q, Rest r
 447 #elif defined(__GNUC__) && defined(__arm__) && !defined(NO_ASM)
 448   extern "C" uint64 divu_3216_1616_ (uint32 x, uint16 y); // -> Quotient q, Rest r
 449 #else
 450   extern "C" uint16 divu_3216_1616_ (uint32 x, uint16 y); // -> Quotient q
 451 #ifdef _MSC_VER
 452   // Workaround MSVC compiler bug.
 453 } extern "C" uint16 divu_16_rest; namespace cln {         // -> Rest r
 454 #else
 455   extern "C" uint16 divu_16_rest;                         // -> Rest r
 456 #endif
 457 #endif
 458 #if defined(__GNUC__) && defined(__sparc64__) && !defined(NO_ASM)
 459   #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
 460     ({var uint32 __x = (x);             \
 461       var uint16 __y = (y);             \
 462       var uint64 __q;                   \
 463       var uint64 __r;                   \
 464       __asm__ __volatile__ (            \
 465         "wr %%g0,%%g0,%%y\n\t"          \
 466         "udiv %2,%3,%0\n\t"             \
 467         "umul %0,%3,%1\n\t"             \
 468         "sub %2,%1,%1"                  \
 469         : "=&r" (__q), "=&r" (__r)      \
 470         : "r" (__x), "r" (__y));        \
 471       cl_unused (q_zuweisung (uint16)__q); \
 472       r_zuweisung (uint16)__r;          \
 473      })
 474 #elif defined(__GNUC__) && (defined(__sparc__) || defined(__sparc64__)) && !defined(NO_ASM)
 475   #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
 476     ({ var uint32 __qr = divu_3216_1616_(x,y); /* extern in Assembler */\
 477        cl_unused (q_zuweisung low16(__qr));                             \
 478        r_zuweisung high16(__qr);                                        \
 479      })
 480 #elif defined(__GNUC__) && defined(__m68k__) && !defined(NO_ASM)
 481   #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
 482     ({var uint32 __x = (x);                                             \
 483       var uint16 __y = (y);                                             \
 484       var uint32 __qr;                                                  \
 485       __asm__ __volatile__ ("                                           \
 486         divu %2,%0                                                      \
 487         " : "=d" (__qr) : "0" (__x), "dm" (__y));                       \
 488       cl_unused (q_zuweisung low16(__qr));                                      \
 489       r_zuweisung high16(__qr);                                         \
 490      })
 491 #elif defined(__GNUC__) && (defined(__i386__) || defined(__x86_64__)) && !defined(NO_ASM)
 492   #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
 493     ({var uint32 __x = (x);                                             \
 494       var uint16 __y = (y);                                             \
 495       var uint16 __q;                                                   \
 496       var uint16 __r;                                                   \
 497       __asm__("divw %4"                                                 \
 498               : "=a" /* %ax */ (__q), "=d" /* %dx */ (__r)              \
 499               : "1" /* %dx */ ((uint16)(high16(__x))), "0" /* %ax */ ((uint16)(low16(__x))), "rm" (__y) \
 500              );                                                         \
 501       cl_unused (q_zuweisung __q);                                              \
 502       r_zuweisung __r;                                                  \
 503      })
 504 #elif defined(__GNUC__) && defined(__arm__) && !defined(NO_ASM)
 505   #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
 506     { var uint64 __q = divu_3216_1616_(x,y); /* extern in Assembler */  \
 507       q_zuweisung retval64_r0(__q);     \
 508       r_zuweisung retval64_r1(__q);     \
 509     }
 510 #elif defined(__GNUC__) && !defined(__arm__)
 511   #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
 512     ({var uint32 __x = (x);                                             \
 513       var uint16 __y = (y);                                             \
 514       var uint16 __q = floor(__x,__y);                                  \
 515       cl_unused (q_zuweisung __q);                                              \
 516       r_zuweisung (__x - __q * __y);                                    \
 517      })
 518 #elif (defined(__sparc__) || defined(__sparc64__)) && !defined(NO_ASM)
 519   #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
 520     { var uint32 __qr = divu_3216_1616_(x,y); /* extern in Assembler */ \
 521       cl_unused (q_zuweisung low16(__qr));                                      \
 522       r_zuweisung high16(__qr);                                         \
 523     }
 524 #else
 525   #define divu_3216_1616(x,y,q_zuweisung,r_zuweisung)  \
 526     { cl_unused (q_zuweisung divu_3216_1616_(x,y)); r_zuweisung divu_16_rest; }
 527   #define NEED_FUNCTION_divu_3216_1616_
 528 #endif
 529
 530 // Dividiert eine 32-Bit-Zahl durch eine 16-Bit-Zahl und
 531 // liefert einen 32-Bit-Quotienten und einen 16-Bit-Rest.
 532 // divu_3216_3216(x,y,q=,r=);
 533 // > uint32 x: Zähler
 534 // > uint16 y: Nenner
 535 // Es sei bekannt, daß y>0.
 536 // < uint32 q: floor(x/y)
 537 // < uint16 r: x mod y
 538 // < x = q*y+r
 539   extern "C" uint32 divu_3216_3216_ (uint32 x, uint16 y); // -> Quotient q
 540 #ifdef _MSC_VER
 541   // Workaround MSVC compiler bug.
 542 } extern "C" uint16 divu_16_rest; namespace cln {         // -> Rest r
 543 #else
 544   extern "C" uint16 divu_16_rest;                         // -> Rest r
 545 #endif
 546 #if defined(__GNUC__) && defined(__sparc64__) && !defined(NO_ASM)
 547   #define divu_3216_3216(x,y,q_zuweisung,r_zuweisung)  \
 548     ({var uint32 __x = (x);        \
 549       var uint16 __y = (y);        \
 550       var uint64 __q;              \
 551       var uint64 __r;              \
 552       __asm__ __volatile__ (       \
 553         "wr %%g0,%%g0,%%y\n\t"     \
 554         "udiv %2,%3,%0\n\t"        \
 555         "umul %0,%3,%1\n\t"        \
 556         "sub %2,%1,%1"             \
 557         : "=&r" (__q), "=&r" (__r) \
 558         : "r" (__x), "r" (__y));   \
 559       q_zuweisung (uint32)__q;     \
 560       r_zuweisung (uint16)__r;     \
 561      })
 562 #elif defined(__sparc__) || defined(__sparc64__) || defined(__i386__) || defined(__x86_64__)
 563   #define divu_3216_3216  divu_3232_3232
 564 #else
 565   // Methode: (beta = 2^16)
 566   // x = x1*beta+x0 schreiben.
 567   // Division mit Rest: x1 = q1*y + r1, wobei 0 <= x1 < beta <= beta*y.
 568   // Also 0 <= q1 < beta, 0 <= r1 < y.
 569   // Division mit Rest: (r1*beta+x0) = q0*y + r0, wobei 0 <= r1*beta+x0 < beta*y.
 570   // Also 0 <= q0 < beta, 0 <= r0 < y
 571   // und x = x1*beta+x0 = (q1*beta+q0)*y + r0.
 572   // Setze q := q1*beta+q0 und r := r0.
 573   #define divu_3216_3216(x,y,q_zuweisung,r_zuweisung)  \
 574     { var uint32 _x = (x);                                              \
 575       var uint16 _y = (y);                                              \
 576       var uint16 _q1;                                                   \
 577       var uint16 _q0;                                                   \
 578       var uint16 _r1;                                                   \
 579       divu_3216_1616(high16(_x),_y, _q1 = , _r1 = );                    \
 580       divu_3216_1616(highlow32(_r1,low16(_x)),_y, _q0 = , r_zuweisung); \
 581       q_zuweisung highlow32(_q1,_q0);                                   \
 582     }
 583 #endif
 584
 585 // Dividiert eine 32-Bit-Zahl durch eine 32-Bit-Zahl und
 586 // liefert einen 32-Bit-Quotienten und einen 32-Bit-Rest.
 587 // divu_3232_3232(x,y,q=,r=);
 588 // > uint32 x: Zähler
 589 // > uint32 y: Nenner
 590 // Es sei bekannt, daß y>0.
 591 // < uint32 q: floor(x/y)
 592 // < uint32 r: x mod y
 593 // < x = q*y+r
 594   extern "C" uint32 divu_3232_3232_ (uint32 x, uint32 y); // -> Quotient q
 595 #ifdef _MSC_VER
 596   // Workaround MSVC compiler bug.
 597 } extern "C" uint32 divu_32_rest; namespace cln {         // -> Rest r
 598 #else
 599   extern "C" uint32 divu_32_rest;                         // -> Rest r
 600 #endif
 601 #if defined(__GNUC__) && defined(__sparc64__) && !defined(NO_ASM)
 602   #define divu_3232_3232(x,y,q_zuweisung,r_zuweisung)  \
 603     ({var uint32 __x = (x);        \
 604       var uint32 __y = (y);        \
 605       var uint64 __q;              \
 606       var uint64 __r;              \
 607       __asm__ __volatile__ (       \
 608         "wr %%g0,%%g0,%%y\n\t"     \
 609         "udiv %2,%3,%0\n\t"        \
 610         "umul %0,%3,%1\n\t"        \
 611         "sub %2,%1,%1"             \
 612         : "=&r" (__q), "=&r" (__r) \
 613         : "r" (__x), "r" (__y));   \
 614       q_zuweisung (uint32)__q;     \
 615       r_zuweisung (uint32)__r;     \
 616      })
 617   #define divu_3232_3232_(x,y) divu_6432_3232_(0,x,y)
 618 #elif defined(__sparc__) || defined(__sparc64__) || defined(__i386__) || defined(__x86_64__)
 619   #define divu_3232_3232(x,y,q_zuweisung,r_zuweisung)  \
 620     divu_6432_3232(0,x,y,q_zuweisung,r_zuweisung)
 621   #define divu_3232_3232_(x,y) divu_6432_3232_(0,x,y)
 622 #else
 623   // Methode: (beta = 2^n = 2^16, n = 16)
 624   // Falls y < beta, handelt es sich um eine 32-durch-16-Bit-Division.
 625   // Falls y >= beta:
 626   // Quotient  q = floor(x/y) < beta  (da 0 <= x < beta^2, y >= beta).
 627   // y habe genau n+k Bits (1 <= k <= n), d.h. 2^(n+k-1) <= y < 2^(n+k).
 628   // Schreibe  x = 2^k*x1 + x0  mit  x1 := floor(x/2^k)
 629   // und       y = 2^k*y1 + y0  mit  y1 := floor(y/2^k)
 630   // und bilde den Näherungs-Quotienten floor(x1/y1)
 631   // oder (noch besser) floor(x1/(y1+1)).
 632   // Wegen 0 <= x1 < 2^(2n) und 0 < 2^(n-1) <= y1 < 2^n
 633   // und  x1/(y1+1) <= x/y < x1/(y1+1) + 2
 634   // (denn x1/(y1+1) = (x1*2^k)/((y1+1)*2^k) <= (x1*2^k)/y <= x/y
 635   // und x/y - x1/(y1+1) = (x+x*y1-x1*y)/(y*(y1+1))
 636   // = (x+x0*y1-x1*y0)/(y*(y1+1)) <= (x+x0*y1)/(y*(y1+1))
 637   // <= x/(y*(y1+1)) + x0/y
 638   // <= 2^(2n)/(2^(n+k-1)*(2^(n-1)+1)) + 2^k/2^(n+k-1)
 639   // = 2^(n-k+1)/(2^(n-1)+1) + 2^(1-n) <= 2^n/(2^(n-1)+1) + 2^(1-n) < 2 )
 640   // gilt  floor(x1/(y1+1)) <= floor(x/y) <= floor(x1/(y1+1)) + 2  .
 641   // Man bildet also  q:=floor(x1/(y1+1))  (ein Shift um n Bit oder
 642   // eine (2n)-durch-n-Bit-Division, mit Ergebnis q <= floor(x/y) < beta)
 643   // und x-q*y und muß hiervon noch höchstens 2 mal y abziehen und q
 644   // incrementieren, um den Quotienten  q = floor(x/y)  und den Rest
 645   // x-floor(x/y)*y  der Division zu bekommen.
 646   #define divu_3232_3232(x,y,q_zuweisung,r_zuweisung)  \
 647     { var uint32 _x = (x);                                              \
 648       var uint32 _y = (y);                                              \
 649       if (_y <= (uint32)(bit(16)-1))                                    \
 650         { var uint16 _q1;                                               \
 651           var uint16 _q0;                                               \
 652           var uint16 _r1;                                               \
 653           divu_3216_1616(high16(_x),_y, _q1 = , _r1 = );                \
 654           divu_3216_1616(highlow32(_r1,low16(_x)),_y, _q0 = , r_zuweisung); \
 655           q_zuweisung highlow32(_q1,_q0);                               \
 656         }                                                               \
 657         else                                                            \
 658         { var uint32 _x1 = _x; /* x1 := x */                            \
 659           var uint32 _y1 = _y; /* y1 := y */                            \
 660           var uint16 _q;                                                \
 661           do { _x1 = floor(_x1,2); _y1 = floor(_y1,2); } /* k erhöhen */\
 662              until (_y1 <= (uint32)(bit(16)-1)); /* bis y1 < beta */    \
 663           { var uint16 _y2 = low16(_y1)+1; /* y1+1 bilden */            \
 664             if (_y2==0)                                                 \
 665               { _q = high16(_x1); } /* y1+1=beta -> ein Shift */        \
 666               else                                                      \
 667               { divu_3216_1616(_x1,_y2,_q=,); } /* Division von x1 durch y1+1 */\
 668           }                                                             \
 669           /* _q = q = floor(x1/(y1+1)) */                               \
 670           /* x-q*y bilden (eine 16-mal-32-Bit-Multiplikation ohne Überlauf): */\
 671           _x -= highlow32_0(mulu16(_q,high16(_y))); /* q * high16(y) * beta */\
 672           /* gefahrlos, da q*high16(y) <= q*y/beta <= x/beta < beta */  \
 673           _x -= mulu16(_q,low16(_y)); /* q * low16(y) */                \
 674           /* gefahrlos, da q*high16(y)*beta + q*low16(y) = q*y <= x */  \
 675           /* Noch höchstens 2 mal y abziehen: */                       \
 676           if (_x >= _y)                                                 \
 677             { _q += 1; _x -= _y;                                        \
 678               if (_x >= _y)                                             \
 679                 { _q += 1; _x -= _y; }                                  \
 680             }                                                           \
 681           r_zuweisung _x;                                               \
 682           q_zuweisung (uint32)(_q);                                     \
 683     }   }
 684   #define NEED_FUNCTION_divu_3232_3232_
 685 #endif
 686
 687 // Dividiert eine 64-Bit-Zahl durch eine 32-Bit-Zahl und
 688 // liefert einen 32-Bit-Quotienten und einen 32-Bit-Rest.
 689 // divu_6432_3232(xhi,xlo,y,q=,r=);
 690 // > uint32 xhi,xlo: x = 2^32*xhi+xlo = Zähler
 691 // > uint32 y: Nenner
 692 // > Es sei bekannt, daß 0 <= x < 2^32*y .
 693 // < uint32 q: floor(x/y)
 694 // < uint32 r: x mod y
 695 // < x = q*y+r
 696 extern "C" uint32 divu_6432_3232_ (uint32 xhi, uint32 xlo, uint32 y); // -> Quotient q
 697 #ifdef _MSC_VER
 698   // Workaround MSVC compiler bug.
 699 } extern "C" uint32 divu_32_rest; namespace cln {                       // -> Rest r
 700 #else
 701   extern "C" uint32 divu_32_rest;                                       // -> Rest r
 702 #endif
 703 #if defined(__GNUC__) && defined(__m68k__) && !defined(NO_ASM)
 704   #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
 705     ({var uint32 __xhi = (xhi);                                         \
 706       var uint32 __xlo = (xlo);                                         \
 707       var uint32 __y = (y);                                             \
 708       var uint32 __q;                                                   \
 709       var uint32 __r;                                                   \
 710       __asm__ __volatile__ ("                                           \
 711         divul %4,%1:%0                                                  \
 712         " : "=d" (__q), "=d" (__r) : "1" (__xhi), "0" (__xlo), "dm" (__y)); \
 713       cl_unused (q_zuweisung __q);                                              \
 714       r_zuweisung __r;                                                  \
 715      })
 716   #define divu_6432_3232_(xhi,xlo,y) \
 717     ({var uint32 ___q; divu_6432_3232(xhi,xlo,y,___q=,); ___q; })
 718 #elif defined(__GNUC__) && defined(__sparc64__) && !defined(NO_ASM)
 719   #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
 720     ({var uint32 __xhi = (xhi);    \
 721       var uint32 __xlo = (xlo);    \
 722       var uint32 __y = (y);        \
 723       var uint64 __q;              \
 724       var uint64 __r;              \
 725       __asm__ __volatile__ (       \
 726         "wr %2,%%g0,%%y\n\t"       \
 727         "udiv %3,%4,%0\n\t"        \
 728         "umul %0,%4,%1\n\t"        \
 729         "sub %3,%1,%1"             \
 730         : "=&r" (__q), "=&r" (__r) \
 731         : "r" (__xhi), "r" (__xlo), "r" (__y)); \
 732       cl_unused (q_zuweisung (uint32)__q); \
 733       r_zuweisung (uint32)__r;     \
 734      })
 735 #elif defined(__GNUC__) && (defined(__sparc__) || defined(__sparc64__)) && !defined(NO_ASM)
 736   #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
 737     ({ var uint32 _q = divu_6432_3232_(xhi,xlo,y); /* extern in Assembler */\
 738        var uint32 _r __asm__("%g1");                                \
 739        cl_unused (q_zuweisung _q); r_zuweisung _r;                                  \
 740      })
 741 #elif defined(__GNUC__) && (defined(__i386__) || defined(__x86_64__)) && !defined(NO_ASM)
 742   #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
 743     ({var uint32 __xhi = (xhi);                                         \
 744       var uint32 __xlo = (xlo);                                         \
 745       var uint32 __y = (y);                                             \
 746       var uint32 __q;                                                   \
 747       var uint32 __r;                                                   \
 748       __asm__ __volatile__ (                                            \
 749          "divl %4"                                                      \
 750          : "=a" /* %eax */ (__q), "=d" /* %edx */ (__r)                 \
 751          : "1" /* %edx */ (__xhi), "0" /* %eax */ (__xlo), "rm" (__y)   \
 752          );                                                             \
 753       cl_unused (q_zuweisung __q);                                              \
 754       r_zuweisung __r;                                                  \
 755      })
 756   #define divu_6432_3232_(xhi,xlo,y) \
 757     ({var uint32 ___q; divu_6432_3232(xhi,xlo,y,___q=,); ___q; })
 758 #elif defined(__GNUC__) && !defined(__arm__)
 759   #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung) \
 760     ({var uint32 __xhi = (xhi);                                         \
 761       var uint32 __xlo = (xlo);                                         \
 762       var uint64 __x = ((uint64)__xhi << 32) | (uint64)__xlo;           \
 763       var uint32 __y = (y);                                             \
 764       var uint32 __q = floor(__x,(uint64)__y);                          \
 765       cl_unused (q_zuweisung __q); r_zuweisung __xlo - __q * __y;               \
 766      })
 767   #define divu_6432_3232_(xhi,xlo,y) \
 768     ({var uint32 ___q; divu_6432_3232(xhi,xlo,y,___q=,); ___q; })
 769 #else
 770   #define divu_6432_3232(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
 771     { cl_unused (q_zuweisung divu_6432_3232_(xhi,xlo,y)); r_zuweisung divu_32_rest; }
 772   #if (defined(__m68k__) || defined(__sparc__) || defined(__sparc64__) || (defined(__i386__) && !defined(MICROSOFT)) || defined(__x86_64__) || defined(__hppa__)) && !defined(NO_ASM)
 773     // divu_6432_3232_ extern in Assembler
 774     #if defined(__sparc__) || defined(__sparc64__)
 775       extern "C" uint32 _get_g1 (void);
 776       #define divu_32_rest  (_get_g1()) // Rückgabe im Register %g1
 777     #else
 778       #define NEED_VAR_divu_32_rest
 779     #endif
 780   #else
 781     #define NEED_FUNCTION_divu_6432_3232_
 782   #endif
 783 #endif
 784
 785 #ifdef HAVE_FAST_LONGLONG
 786
 787 // Dividiert eine 64-Bit-Zahl durch eine 32-Bit-Zahl und
 788 // liefert einen 32-Bit-Quotienten und einen 32-Bit-Rest.
 789 // divu_6432_3232_w(x,y,q=,r=);
 790 // > uint64 x: Zähler
 791 // > uint32 y: Nenner
 792 // > Es sei bekannt, daß 0 <= x < 2^32*y .
 793 // < uint32 q: floor(x/y)
 794 // < uint32 r: x mod y
 795 // < x = q*y+r
 796 #if defined(__GNUC__) && defined(__sparc64__) && !defined(NO_ASM)
 797   // Prefer the udiv and umul instructions over the udivx and mulx instructions
 798   // (overkill).
 799   #define divu_6432_3232_w(x,y,q_zuweisung,r_zuweisung)  \
 800     ({var uint64 __x = (x);           \
 801       var uint32 __xhi = high32(__x); \
 802       var uint32 __xlo = low32(__x);  \
 803       var uint32 __y = (y);           \
 804       var uint64 __q;                 \
 805       var uint64 __r;                 \
 806       __asm__ __volatile__ (          \
 807         "wr %2,%%g0,%%y\n\t"          \
 808         "udiv %3,%4,%0\n\t"           \
 809         "umul %0,%4,%1\n\t"           \
 810         "sub %3,%1,%1"                \
 811         : "=&r" (__q), "=&r" (__r)    \
 812         : "r" (__xhi), "r" (__xlo), "r" (__y)); \
 813       q_zuweisung (uint32)__q;        \
 814       r_zuweisung (uint32)__r;        \
 815      })
 816 #elif defined(__GNUC__) && (defined(__alpha__) || defined(__ia64__) || defined(__mips64__) || defined(__sparc64__))
 817   // On __alpha__, computing the remainder by multiplication is just two
 818   // instructions, compared to the __remqu (libc) function call for the %
 819   // operator.
 820   // On __ia64__, computing the remainder by multiplication is just four
 821   // instructions, compared to the __umoddi3 (libgcc) function call for the %
 822   // operator.
 823   // On __mips64__, computing the remainder by multiplication is just two
 824   // instructions, compared to the __umoddi3 (libgcc) function call for the %
 825   // operator.
 826   // On __sparc64__, computing the remainder by multiplication uses a 32-bit
 827   // multiplication instruction, compared to a 64-bit multiplication when the %
 828   // operator is used.
 829   #define divu_6432_3232_w(x,y,q_zuweisung,r_zuweisung)  \
 830     ({var uint64 __x = (x);                                             \
 831       var uint32 __y = (y);                                             \
 832       var uint32 __q = floor(__x,(uint64)__y);                          \
 833       q_zuweisung __q; r_zuweisung (uint32)__x - __q * __y;             \
 834      })
 835 #elif defined(__GNUC__) && defined(__x86_64__)
 836   // On __x86_64__, gcc 4.0 performs both quotient and remainder computation
 837   // in a single instruction.
 838   #define divu_6432_3232_w(x,y,q_zuweisung,r_zuweisung)  \
 839     ({var uint64 __x = (x);                                             \
 840       var uint32 __y = (y);                                             \
 841       var uint32 __q = floor(__x,(uint64)__y);                          \
 842       q_zuweisung __q; r_zuweisung __x % (uint64)__y;                   \
 843      })
 844 #else
 845   #define divu_6432_3232_w(x,y,q_zuweisung,r_zuweisung)  \
 846     { var uint64 __x = (x);                                               \
 847       divu_6432_3232(high32(__x),low32(__x),(y),q_zuweisung,r_zuweisung); \
 848     }
 849 #endif
 850
 851 // Dividiert eine 64-Bit-Zahl durch eine 32-Bit-Zahl und
 852 // liefert einen 64-Bit-Quotienten und einen 32-Bit-Rest.
 853 // divu_6432_6432(x,y,q=,r=);
 854 // > uint64 x: Zähler
 855 // > uint32 y: Nenner
 856 // > Es sei bekannt, daß y>0.
 857 // < uint64 q: floor(x/y)
 858 // < uint32 r: x mod y
 859 // < x = q*y+r
 860 #if defined(__GNUC__) && (defined(__alpha__) || defined(__ia64__) || defined(__mips64__) || defined(__sparc64__))
 861   // On __alpha__, computing the remainder by multiplication is just two
 862   // instructions, compared to the __remqu (libc) function call for the %
 863   // operator.
 864   // On __ia64__, computing the remainder by multiplication is just four
 865   // instructions, compared to the __umoddi3 (libgcc) function call for the %
 866   // operator.
 867   // On __mips64__, computing the remainder by multiplication is just two
 868   // instructions, compared to the __umoddi3 (libgcc) function call for the %
 869   // operator.
 870   // On __sparc64__, computing the remainder by multiplication uses a 32-bit
 871   // multiplication instruction, compared to a 64-bit multiplication when the %
 872   // operator is used.
 873   #define divu_6432_6432(x,y,q_zuweisung,r_zuweisung)  \
 874     ({var uint64 _x = (x);                    \
 875       var uint32 _y = (y);                    \
 876       var uint64 _q;                          \
 877       q_zuweisung _q = floor(_x,(uint64)_y);  \
 878       r_zuweisung low32(_x) - low32(_q) * _y; \
 879      })
 880 #elif defined(__GNUC__) && defined(__x86_64__)
 881   // On __x86_64__, gcc 4.0 performs both quotient and remainder computation
 882   // in a single instruction.
 883   #define divu_6432_6432(x,y,q_zuweisung,r_zuweisung)  \
 884     ({var uint64 _x = (x);               \
 885       var uint32 _y = (y);               \
 886       q_zuweisung floor(_x,(uint64)_y);  \
 887       r_zuweisung _x % (uint64)_y;       \
 888      })
 889 #else
 890   // Methode: (beta = 2^32)
 891   // x = x1*beta+x0 schreiben.
 892   // Division mit Rest: x1 = q1*y + r1, wobei 0 <= x1 < beta <= beta*y.
 893   // Also 0 <= q1 < beta, 0 <= r1 < y.
 894   // Division mit Rest: (r1*beta+x0) = q0*y + r0, wobei 0 <= r1*beta+x0 < beta*y.
 895   // Also 0 <= q0 < beta, 0 <= r0 < y
 896   // und x = x1*beta+x0 = (q1*beta+q0)*y + r0.
 897   // Setze q := q1*beta+q0 und r := r0.
 898   #if defined(__GNUC__)
 899     #define divu_6432_6432(x,y,q_zuweisung,r_zuweisung)  \
 900       ({var uint64 _x = (x);            \
 901         var uint32 _y = (y);            \
 902         var uint32 _q1;                 \
 903         var uint32 _q0;                 \
 904         var uint32 _r1;                 \
 905         divu_6432_3232(0,high32(_x),_y, _q1 = , _r1 = ); \
 906         divu_6432_3232(_r1,low32(_x),_y, _q0 = , r_zuweisung); \
 907         q_zuweisung highlow64(_q1,_q0); \
 908        })
 909   #else
 910     #define divu_6432_6432(x,y,q_zuweisung,r_zuweisung)  \
 911       {var uint64 _x = (x);            \
 912        var uint32 _y = (y);            \
 913        var uint32 _q1;                 \
 914        var uint32 _q0;                 \
 915        var uint32 _r1;                 \
 916        divu_6432_3232(0,high32(_x),_y, _q1 = , _r1 = ); \
 917        divu_6432_3232(_r1,low32(_x),_y, _q0 = , r_zuweisung); \
 918        q_zuweisung highlow64(_q1,_q0); \
 919       }
 920   #endif
 921 #endif
 922
 923 // Dividiert eine 64-Bit-Zahl durch eine 64-Bit-Zahl und
 924 // liefert einen 64-Bit-Quotienten und einen 64-Bit-Rest.
 925 // divu_6464_6464(x,y,q=,r=);
 926 // > uint64 x: Zähler
 927 // > uint64 y: Nenner
 928 // > Es sei bekannt, daß y>0.
 929 // < uint64 q: floor(x/y)
 930 // < uint64 r: x mod y
 931 // < x = q*y+r
 932 #if defined(__GNUC__) && (defined(__alpha__) || defined(__ia64__) || defined(__mips64__) || defined(__sparc64__))
 933   // On __alpha__, computing the remainder by multiplication is just two
 934   // instructions, compared to the __remqu (libc) function call for the %
 935   // operator.
 936   // On __ia64__, computing the remainder by multiplication is just four
 937   // instructions, compared to the __umoddi3 (libgcc) function call for the %
 938   // operator.
 939   // On __mips64__, computing the remainder by multiplication is just two
 940   // instructions, compared to the __umoddi3 (libgcc) function call for the %
 941   // operator.
 942   // On __sparc64__, it doesn't matter.
 943   #define divu_6464_6464(x,y,q_zuweisung,r_zuweisung)  \
 944     ({var uint64 _x = (x);           \
 945       var uint64 _y = (y);           \
 946       var uint64 _q;                 \
 947       q_zuweisung _q = floor(_x,_y); \
 948       r_zuweisung _x - _q * _y;      \
 949      })
 950 #elif defined(__GNUC__) && (defined(__sparc64__) || defined(__x86_64__))
 951   // On __sparc64__, it doesn't matter.
 952   // On __x86_64__, gcc 4.0 performs both quotient and remainder computation
 953   // in a single instruction.
 954   #define divu_6464_6464(x,y,q_zuweisung,r_zuweisung)  \
 955     ({var uint64 _x = (x);      \
 956       var uint64 _y = (y);      \
 957       q_zuweisung floor(_x,_y); \
 958       r_zuweisung _x % _y;      \
 959      })
 960 #else
 961   // For unknown CPUs, we don't know whether gcc's __udivdi3 function plus a
 962   // multiplication is slower or faster than our own divu_6464_6464_ routine.
 963   // Anyway, call our own routine.
 964   extern "C" uint64 divu_6464_6464_ (uint64 x, uint64 y); // -> Quotient q
 965 #ifdef _MSC_VER
 966   // Workaround MSVC compiler bug.
 967 } extern "C" uint64 divu_64_rest; namespace cln {         // -> Rest r
 968 #else
 969   extern "C" uint64 divu_64_rest;                         // -> Rest r
 970 #endif
 971   #define divu_6464_6464(x,y,q_zuweisung,r_zuweisung)  \
 972     { q_zuweisung divu_6464_6464_(x,y); r_zuweisung divu_64_rest; }
 973   #define NEED_VAR_divu_64_rest
 974   #define NEED_FUNCTION_divu_6464_6464_
 975 #endif
 976
 977 // Dividiert eine 128-Bit-Zahl durch eine 64-Bit-Zahl und
 978 // liefert einen 64-Bit-Quotienten und einen 64-Bit-Rest.
 979 // divu_12864_6464(xhi,xlo,y,q=,r=);
 980 // > uint64 xhi,xlo: x = 2^64*xhi+xlo = Zähler
 981 // > uint64 y: Nenner
 982 // > Es sei bekannt, daß 0 <= x < 2^64*y .
 983 // < uint64 q: floor(x/y)
 984 // < uint64 r: x mod y
 985 // < x = q*y+r
 986   extern "C" uint64 divu_12864_6464_ (uint64 xhi, uint64 xlo, uint64 y); // -> Quotient q
 987 #ifdef _MSC_VER
 988   // Workaround MSVC compiler bug.
 989 } extern "C" uint64 divu_64_rest; namespace cln {                        // -> Rest r
 990 #else
 991   extern "C" uint64 divu_64_rest;                                        // -> Rest r
 992 #endif
 993 #if defined(__GNUC__) && defined(__x86_64__) && !defined(NO_ASM)
 994   #define divu_12864_6464(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
 995     ({var uint64 __xhi = (xhi);                                         \
 996       var uint64 __xlo = (xlo);                                         \
 997       var uint64 __y = (y);                                             \
 998       var uint64 __q;                                                   \
 999       var uint64 __r;                                                   \
1000       __asm__ __volatile__ (                                            \
1001          "divq %4"                                                      \
1002          : "=a" /* %rax */ (__q), "=d" /* %rdx */ (__r)                 \
1003          : "1" /* %rdx */ (__xhi), "0" /* %rax */ (__xlo), "rm" (__y)   \
1004          );                                                             \
1005       q_zuweisung __q;                                                  \
1006       r_zuweisung __r;                                                  \
1007      })
1008   #define divu_12864_64364_(xhi,xlo,y) \
1009     ({var uint64 ___q; divu_12864_6464(xhi,xlo,y,___q=,); ___q; })
1010 #else
1011   #define divu_12864_6464(xhi,xlo,y,q_zuweisung,r_zuweisung)  \
1012     { q_zuweisung divu_12864_6464_(xhi,xlo,y); r_zuweisung divu_64_rest; }
1013   #define NEED_VAR_divu_64_rest
1014   #define NEED_FUNCTION_divu_12864_6464_
1015 #endif
1016
1017 #endif /* HAVE_FAST_LONGLONG */
1018
1019
1020 // Zieht die Ganzzahl-Wurzel aus einer 32-Bit-Zahl und
1021 // liefert eine 16-Bit-Wurzel und einen Rest.
1022 // isqrt_32_16(x,y=,sqrtp=);
1023 // > uint32 x: Radikand, >= 2^30, < 2^32
1024 // < uint16 y: floor(sqrt(x)), >= 2^15, < 2^16
1025 // < boolean sqrtp: /=0, falls x=y^2
1026   // Methode:
1027   // y := 2^16 als Anfangswert,
1028   // y := floor((y + floor(x/y))/2) als nächster Wert,
1029   // solange z := floor(x/y) < y, setze y := floor((y+z)/2).
1030   // y ist fertig; x=y^2 genau dann, wenn z=y und die letzte Division aufging.
1031   // (Beweis:
1032   //  1. Die Folge der y ist streng monoton fallend.
1033   //  2. Stets gilt y >= floor(sqrt(x)) (denn für alle y>0 ist
1034   //     y + x/y >= 2*sqrt(x) und daher  floor((y + floor(x/y))/2) =
1035   //     floor(y/2 + x/(2*y)) >= floor(sqrt(x)) ).
1036   //  3. Am Schluß gilt x >= y^2.
1037   // )
1038   #define isqrt_32_16(x,y_zuweisung,sqrtp_zuweisung)  \
1039     { var uint32 _x = (x);                                              \
1040       var uint16 _x1 = high16(_x);                                      \
1041       var uint16 _y = floor(_x1,2) | bit(16-1);                         \
1042       loop                                                              \
1043         { var uint16 _z;                                                \
1044           var uint16 _r;                                                \
1045           if (_x1 >= _y) /* Division _x/_y ergäbe Überlauf -> _z > _y */\
1046             { cl_unused (sqrtp_zuweisung FALSE); break; }                       \
1047           divu_3216_1616(_x,_y, _z=,_r=); /* Dividiere _x/_y */         \
1048           if (_z >= _y)                                                 \
1049             { cl_unused (sqrtp_zuweisung (_z == _y) && (_r == 0)); break; } \
1050           _y = floor((uint16)(_z+_y),2) | bit(16-1); /* _y muß >= 2^15 bleiben */\
1051         }                                                               \
1052       y_zuweisung _y;                                                   \
1053     }
1054
1055 // Zieht die Ganzzahl-Wurzel aus einer 64-Bit-Zahl und
1056 // liefert eine 32-Bit-Wurzel und einen Rest.
1057 // isqrt_64_32(xhi,xlo,y=,sqrtp=);
1058 // > uint32 xhi,xlo: Radikand x = 2^32*xhi+xlo, >= 2^62, < 2^64
1059 // < uint32 y: floor(sqrt(x)), >= 2^31, < 2^32
1060 // < boolean sqrtp: /=0, falls x=y^2
1061 #if defined(__sparc__) || defined(__sparc64__) || defined(__m68k__) || defined(__hppa__)
1062   // Methode:
1063   // y := 2^32 als Anfangswert,
1064   // y := floor((y + floor(x/y))/2) als nächster Wert,
1065   // solange z := floor(x/y) < y, setze y := floor((y+z)/2).
1066   // y ist fertig; x=y^2 genau dann, wenn z=y und die letzte Division aufging.
1067   // (Beweis:
1068   //  1. Die Folge der y ist streng monoton fallend.
1069   //  2. Stets gilt y >= floor(sqrt(x)) (denn für alle y>0 ist
1070   //     y + x/y >= 2*sqrt(x) und daher  floor((y + floor(x/y))/2) =
1071   //     floor(y/2 + x/(2*y)) >= floor(sqrt(x)) ).
1072   //  3. Am Schluß gilt x >= y^2.
1073   // )
1074   #define isqrt_64_32(xhi,xlo,y_zuweisung,sqrtp_zuweisung)  \
1075     { var uint32 _xhi = (xhi);                                          \
1076       var uint32 _xlo = (xlo);                                          \
1077       var uint32 _y = floor(_xhi,2) | bit(32-1);                        \
1078       loop                                                              \
1079         { var uint32 _z;                                                \
1080           var uint32 _rest;                                             \
1081           if (_xhi >= _y) /* Division _x/_y ergäbe Überlauf -> _z > _y */\
1082             { sqrtp_zuweisung FALSE; break; }                           \
1083           divu_6432_3232(_xhi,_xlo,_y, _z=,_rest=); /* Dividiere _x/_y */\
1084           if (_z >= _y)                                                 \
1085             { sqrtp_zuweisung (_z == _y) && (_rest == 0); break; }      \
1086           _y = floor(_z+_y,2) | bit(32-1); /* _y muß >= 2^31 bleiben */        \
1087         }                                                               \
1088       y_zuweisung _y;                                                   \
1089     }
1090 #else
1091   // Methode:
1092   // Wie bei UDS_sqrt mit n=2.
1093   // y = 2^16*yhi + ylo ansetzen.
1094   // Dann muß
1095   //   yhi = floor(y/2^16) = floor(floor(sqrt(x))/2^16)
1096   //       = floor(sqrt(x)/2^16) = floor(sqrt(x/2^32)) = isqrt(xhi)
1097   // sein. Es folgt yhi >= 2^15.
1098   // Danach sucht man das größte ylo >=0 mit
1099   // x - 2^32*yhi^2 >= 2*2^16*yhi*ylo + ylo^2.
1100   // Dazu setzen wir  xhi*2^32+xlo := x - 2^32*yhi^2
1101   // (also xhi := xhi - yhi^2, das ist >=0, <=2*yhi).
1102   // Die Schätzung für die zweite Ziffer
1103   //     ylo' := min(2^16-1,floor((xhi*2^32+xlo)/(2*2^16*yhi)))
1104   // erfüllt ylo'-1 <= ylo <= ylo', ist also um höchstens 1 zu groß.
1105   // (Beweis: Rechte Ungleichung klar, da  ylo < 2^16  und
1106   //   xhi*2^32+xlo >= 2*2^16*yhi*ylo + ylo^2 >= 2*2^16*yhi*ylo
1107   //   ==> (xhi*2^32+xlo)/(2*2^16*yhi) >= ylo  gelten muß.
1108   //   Linke Ungleichung: Falls floor(...)>=2^16, ist
1109   //   xhi*2^32+xlo >= 2*2^16*2^16*yhi >= 2*2^16*yhi*(2^16-1) + 2^32
1110   //                >= 2*2^16*yhi*(2^16-1) + (2^16-1)^2
1111   //   und xhi*2^32+xlo < 2*2^16*2^16*yhi + (2^16)^2, also
1112   //   ylo = 2^16-1 = ylo'.
1113   //   Sonst ist ylo' = floor((xhi*2^32+xlo)/(2*2^16*yhi)), also
1114   //   xhi*2^32+xlo >= 2*2^16*yhi*ylo' >= 2*2^16*yhi*(ylo'-1) + 2^32
1115   //                >= 2*2^16*yhi*(ylo'-1) + (ylo'-1)^2,
1116   //   also ylo >= ylo'-1 nach Definition von ylo.)
1117   #define isqrt_64_32(xhi,xlo,y_zuweisung,sqrtp_zuweisung)  \
1118     { var uint32 _xhi = (xhi);                                          \
1119       var uint32 _xlo = (xlo);                                          \
1120       var uint16 _yhi;                                                  \
1121       var uint16 _ylo;                                                  \
1122       /* erste Ziffer berechnen: */                                     \
1123       isqrt_32_16(_xhi,_yhi=,); /* yhi := isqrt(xhi) */                 \
1124       _xhi -= mulu16(_yhi,_yhi); /* jetzt 0 <= xhi <= 2*yhi */          \
1125       /* x = 2^32*yhi^2 + 2^32*xhi + xlo */                             \
1126       /* Schätzung für die zweite Ziffer berechnen: */                        \
1127       /* ylo := min(2^16-1,floor((xhi*2^32+xlo)/(2*2^16*yhi))) bilden: */\
1128      {var uint32 _z = (_xhi << 15) | (_xlo >> 17); /* < 2^15*(2*yhi+1) */\
1129       var uint32 _r = highlow32_0(_yhi);                                \
1130       if (_z >= _r)                                                     \
1131         { _ylo = bit(16)-1; _r = _z - _r + (uint32)_yhi; }              \
1132         else                                                            \
1133         { divu_3216_1616(_z,_yhi, _ylo=,_r=); }                         \
1134       /* x = 2^32*yhi^2 + 2*2^16*yhi*ylo + 2^17*r + (xlo mod 2^17), */  \
1135       /* 0 <= r < yhi + 2^15 */                                         \
1136       _xlo = (_r << 17) | (_xlo & (bit(17)-1));                         \
1137       /* x = 2^32*yhi^2 + 2*2^16*yhi*ylo + 2^32*floor(r/2^15) + xlo */  \
1138       _z = mulu16(_ylo,_ylo); /* z = ylo^2 */                           \
1139       /* Versuche vom Rest 2^32*floor(r/2^15) + xlo  z zu subtrahieren. */\
1140       /* Falls Rest >= z (d.h. r>=2^15 oder xlo>=z), ist ylo fertig, */ \
1141       /* und es gilt x=y^2 genau dann, wenn r<2^15 und xlo=z. */        \
1142       /* Sonst (d.h. r<2^15 und xlo<z), muß man ylo erniedrigen. Dazu */\
1143       /* setzt man  ylo := ylo-1, z := z-(2*ylo+1), */                  \
1144       /* Rest := Rest + 2^17*yhi = xlo + 2^17*yhi >= 2^32 > z, also x>y^2. */\
1145       if (_r < bit(15))                                                 \
1146         { if (_xlo < _z)                                                \
1147             { _ylo -= 1; cl_unused (sqrtp_zuweisung FALSE); }           \
1148             else                                                        \
1149             { cl_unused (sqrtp_zuweisung (_xlo == _z)); }                       \
1150         }                                                               \
1151         else                                                            \
1152         { cl_unused (sqrtp_zuweisung FALSE); }                          \
1153       y_zuweisung highlow32(_yhi,_ylo);                                 \
1154     }}
1155 #endif
1156
1157 #ifdef HAVE_FAST_LONGLONG
1158
1159 // Zieht die Ganzzahl-Wurzel aus einer 128-Bit-Zahl und
1160 // liefert eine 64-Bit-Wurzel und einen Rest.
1161 // isqrt_128_64(xhi,xlo,y=,sqrtp=);
1162 // > uint64 xhi,xlo: Radikand x = 2^64*xhi+xlo, >= 2^126, < 2^128
1163 // < uint64 y: floor(sqrt(x)), >= 2^63, < 2^64
1164 // < boolean sqrtp: /=0, falls x=y^2
1165   // Methode:
1166   // Wie bei UDS_sqrt mit n=2.
1167   // y = 2^32*yhi + ylo ansetzen.
1168   // Dann muß
1169   //   yhi = floor(y/2^32) = floor(floor(sqrt(x))/2^32)
1170   //       = floor(sqrt(x)/2^32) = floor(sqrt(x/2^64)) = isqrt(xhi)
1171   // sein. Es folgt yhi >= 2^31.
1172   // Danach sucht man das größte ylo >=0 mit
1173   // x - 2^64*yhi^2 >= 2*2^32*yhi*ylo + ylo^2.
1174   // Dazu setzen wir  xhi*2^64+xlo := x - 2^64*yhi^2
1175   // (also xhi := xhi - yhi^2, das ist >=0, <=2*yhi).
1176   // Die Schätzung für die zweite Ziffer
1177   //     ylo' := min(2^32-1,floor((xhi*2^64+xlo)/(2*2^32*yhi)))
1178   // erfüllt ylo'-1 <= ylo <= ylo', ist also um höchstens 1 zu groß.
1179   // (Beweis: Rechte Ungleichung klar, da  ylo < 2^32  und
1180   //   xhi*2^64+xlo >= 2*2^32*yhi*ylo + ylo^2 >= 2*2^32*yhi*ylo
1181   //   ==> (xhi*2^64+xlo)/(2*2^32*yhi) >= ylo  gelten muß.
1182   //   Linke Ungleichung: Falls floor(...)>=2^32, ist
1183   //   xhi*2^64+xlo >= 2*2^32*2^32*yhi >= 2*2^32*yhi*(2^32-1) + 2^64
1184   //                >= 2*2^32*yhi*(2^32-1) + (2^32-1)^2
1185   //   und xhi*2^64+xlo < 2*2^32*2^32*yhi + (2^32)^2, also
1186   //   ylo = 2^32-1 = ylo'.
1187   //   Sonst ist ylo' = floor((xhi*2^64+xlo)/(2*2^32*yhi)), also
1188   //   xhi*2^64+xlo >= 2*2^32*yhi*ylo' >= 2*2^32*yhi*(ylo'-1) + 2^64
1189   //                >= 2*2^32*yhi*(ylo'-1) + (ylo'-1)^2,
1190   //   also ylo >= ylo'-1 nach Definition von ylo.)
1191   #define isqrt_128_64(x_hi,x_lo,y_zuweisung,sqrtp_zuweisung)  \
1192     { var uint64 xhi = (x_hi);                                          \
1193       var uint64 xlo = (x_lo);                                          \
1194       var uint32 yhi;                                                   \
1195       var uint32 ylo;                                                   \
1196       /* erste Ziffer berechnen: */                                     \
1197       isqrt_64_32(high32(xhi),low32(xhi),yhi=,); /* yhi := isqrt(xhi) */\
1198       xhi -= mulu32_w(yhi,yhi); /* jetzt 0 <= xhi <= 2*yhi */           \
1199       /* x = 2^64*yhi^2 + 2^64*xhi + xlo */                             \
1200       /* Schätzung für die zweite Ziffer berechnen: */                        \
1201       /* ylo := min(2^32-1,floor((xhi*2^64+xlo)/(2*2^32*yhi))) bilden: */\
1202      {var uint64 z = (xhi << 31) | (xlo >> 33); /* < 2^31*(2*yhi+1) */  \
1203       var uint64 r = highlow64_0(yhi);                                  \
1204       if (z >= r)                                                       \
1205         { ylo = bit(32)-1; r = z - r + (uint64)yhi; }                   \
1206         else                                                            \
1207         { divu_6432_3232_w(z,yhi, ylo=,r=); }                           \
1208       /* x = 2^64*yhi^2 + 2*2^32*yhi*ylo + 2^33*r + (xlo mod 2^33), */  \
1209       /* 0 <= r < yhi + 2^31 */                                         \
1210       xlo = (r << 33) | (xlo & (bit(33)-1));                            \
1211       /* x = 2^64*yhi^2 + 2*2^32*yhi*ylo + 2^64*floor(r/2^31) + xlo */  \
1212       z = mulu32_w(ylo,ylo); /* z = ylo^2 */                            \
1213       /* Versuche vom Rest 2^64*floor(r/2^31) + xlo  z zu subtrahieren. */\
1214       /* Falls Rest >= z (d.h. r>=2^31 oder xlo>=z), ist ylo fertig, */ \
1215       /* und es gilt x=y^2 genau dann, wenn r<2^31 und xlo=z. */        \
1216       /* Sonst (d.h. r<2^31 und xlo<z), muß man ylo erniedrigen. Dazu */\
1217       /* setzt man  ylo := ylo-1, z := z-(2*ylo+1), */                  \
1218       /* Rest := Rest + 2^33*yhi = xlo + 2^33*yhi >= 2^64 > z, also x>y^2. */\
1219       if (r < bit(31))                                                  \
1220         { if (xlo < z)                                                  \
1221             { ylo -= 1; sqrtp_zuweisung FALSE; }                        \
1222             else                                                        \
1223             { sqrtp_zuweisung (xlo == z); }                             \
1224         }                                                               \
1225         else                                                            \
1226         { sqrtp_zuweisung FALSE; }                                      \
1227       y_zuweisung highlow64(yhi,ylo);                                   \
1228     }}
1229
1230 #endif /* HAVE_FAST_LONGLONG */
1231
1232 // Zieht die Ganzzahl-Wurzel aus einer 32-Bit-Zahl und
1233 // liefert eine 16-Bit-Wurzel.
1234 // isqrt(x)
1235 // > uintL x : Radikand, >=0, <2^32
1236 // < uintL ergebnis : Wurzel, >=0, <2^16
1237   extern uintL isqrt (uintL x);
1238
1239 // Extracts integer root of a 64-bit number and returns a 32-bit number.
1240 // isqrt(x)
1241 // > uintQ x : radicand, >=0, <2^64
1242 // < uintL result : square root, >=0, <2^32
1243   extern uintL isqrt (uintQ x);
1244
1245 // Sorry for this. We need an isqrt function taking uintC arguments but we
1246 // cannot use overloading since this would lead to ambiguities with any of the
1247 // two signatures above.
1248   inline uintL isqrtC (uintC x)
1249   {
1250 #if (intCsize==32)
1251       return isqrt((uintL)x);
1252 #else
1253       return isqrt((uintQ)x);
1254 #endif
1255   }
1256
1257
1258 // Zieht die Ganzzahl-Wurzel aus einer 64-Bit-Zahl und
1259 // liefert eine 32-Bit-Wurzel.
1260 // isqrt(x1,x0)
1261 // > uintL2 x = x1*2^32+x0 : Radikand, >=0, <2^64
1262 // < uintL ergebnis : Wurzel, >=0, <2^32
1263   extern uintL isqrt (uintL x1, uintL x0);
1264
1265
1266 // Bits einer 8-Bit-Zahl zählen:
1267 // integerlength8(digit,size=);
1268 // setzt size auf die höchste in digit vorkommende Bitnummer.
1269 // > digit: ein uint8 >0
1270 // < size: >0, <=8, mit 2^(size-1) <= digit < 2^size
1271 #if defined(__GNUC__) && defined(__m68k__) && !defined(NO_ASM)
1272   #define integerlength8(digit,size_zuweisung)  \
1273     { var uintL _zero_counter; /* zählt die führenden Nullbits in digit            */\
1274       __asm__("bfffo %1{#0:#8},%0" : "=d" (_zero_counter) : "dm" ((uint8)(digit)) ); \
1275       size_zuweisung (8-_zero_counter);                                              \
1276     }
1277 #elif defined(__sparc__) && !defined(__sparc64__)
1278   #define integerlength8(digit,size_zuweisung)  \
1279     integerlength32((uint32)(digit),size_zuweisung) // siehe unten
1280 #elif defined(__GNUC__) && defined(__i386__) && !defined(NO_ASM)
1281   #define integerlength8(digit,size_zuweisung)  \
1282     integerlength16((uint16)(digit),size_zuweisung)
1283 #else
1284   #define integerlength8(digit,size_zuweisung)  \
1285     { var uintC _bitsize = 1;                                   \
1286       var uintL _x8 = (uint8)(digit);                           \
1287       /* _x8 hat höchstens 8 Bits.                             */\
1288       if (_x8 >= bit(4)) { _x8 = _x8>>4; _bitsize += 4; }               \
1289       /* _x8 hat höchstens 4 Bits.                             */\
1290       if (_x8 >= bit(2)) { _x8 = _x8>>2; _bitsize += 2; }               \
1291       /* _x8 hat höchstens 2 Bits.                             */\
1292       if (_x8 >= bit(1)) { /* _x8 = _x8>>1; */ _bitsize += 1; } \
1293       /* _x8 hat höchstens 1 Bit. Dieses Bit muß gesetzt sein. */\
1294       size_zuweisung _bitsize;                                  \
1295     }
1296 #endif
1297
1298 // Bits einer 16-Bit-Zahl zählen:
1299 // integerlength16(digit,size=);
1300 // setzt size auf die höchste in digit vorkommende Bitnummer.
1301 // > digit: ein uint16 >0
1302 // < size: >0, <=16, mit 2^(size-1) <= digit < 2^size
1303 #if defined(__GNUC__) && defined(__m68k__) && !defined(NO_ASM)
1304   #define integerlength16(digit,size_zuweisung)  \
1305     { var uintL _zero_counter; /* zählt die führenden Nullbits in digit              */\
1306       __asm__("bfffo %1{#0:#16},%0" : "=d" (_zero_counter) : "dm" ((uint16)(digit)) ); \
1307       size_zuweisung (16-_zero_counter);                                               \
1308     }
1309 #elif defined(__sparc__) && !defined(__sparc64__)
1310   #define integerlength16(digit,size_zuweisung)  \
1311     integerlength32((uint32)(digit),size_zuweisung) // siehe unten
1312 #elif defined(__GNUC__) && defined(__i386__) && !defined(NO_ASM)
1313   #define integerlength16(digit,size_zuweisung)  \
1314     { var uintW _one_position; /* Position der führenden 1                 */\
1315       __asm__("bsrw %1,%0" : "=r" (_one_position) : "r" ((uint16)(digit)) ); \
1316       size_zuweisung (1+_one_position);                                      \
1317     }
1318 // Die weiteren kommen von gcc/longlong.h :
1319 #elif defined(__GNUC__) && defined(__ibm032__) && !defined(NO_ASM) // RT/ROMP
1320   #define integerlength16(digit,size_zuweisung)  \
1321     { var uintL _zero_counter; /* zählt die führenden Nullbits in digit   */\
1322       __asm__("clz %0,%1" : "=r" (_zero_counter) : "r" ((uint32)(digit)) ); \
1323       size_zuweisung (16-_zero_counter);                                    \
1324     }
1325 #else
1326   #define integerlength16(digit,size_zuweisung)  \
1327     { var uintC _bitsize = 1;                                           \
1328       var uintWL _x16 = (uint16)(digit);                                        \
1329       /* _x16 hat höchstens 16 Bits.                                   */\
1330       if (_x16 >= bit(8)) { _x16 = _x16>>8; _bitsize += 8; }            \
1331       /* _x16 hat höchstens 8 Bits.                                    */\
1332       if (_x16 >= bit(4)) { _x16 = _x16>>4; _bitsize += 4; }            \
1333       /* _x16 hat höchstens 4 Bits.                                    */\
1334       if (_x16 >= bit(2)) { _x16 = _x16>>2; _bitsize += 2; }            \
1335       /* _x16 hat höchstens 2 Bits.                                    */\
1336       if (_x16 >= bit(1)) { /* _x16 = _x16>>1; */ _bitsize += 1; }              \
1337       /* _x16 hat höchstens 1 Bit. Dieses Bit muß gesetzt sein.        */\
1338       size_zuweisung _bitsize;                                          \
1339     }
1340 #endif
1341
1342 // Bits einer 32-Bit-Zahl zählen:
1343 // integerlength32(digit,size=);
1344 // setzt size auf die höchste in digit vorkommende Bitnummer.
1345 // > digit: ein uint32 >0
1346 // < size: >0, <=32, mit 2^(size-1) <= digit < 2^size
1347 #if defined(__GNUC__) && defined(__m68k__) && !defined(NO_ASM)
1348   #define integerlength32(digit,size_zuweisung)  \
1349     { var uintL _zero_counter; /* zählt die führenden Nullbits in digit              */\
1350       __asm__("bfffo %1{#0:#32},%0" : "=d" (_zero_counter) : "dm" ((uint32)(digit)) ); \
1351       size_zuweisung (32-_zero_counter);                                               \
1352     }
1353 #elif defined(__sparc__) && !defined(__sparc64__) && defined(FAST_DOUBLE)
1354   #define integerlength32(digit,size_zuweisung)  \
1355     {var union { double f; uint32 i[2]; } __fi;                         \
1356      const int df_mant_len = 52;  /* mantissa bits (excl. hidden bit) */\
1357      const int df_exp_mid = 1022; /* exponent bias */                   \
1358      /* Bilde 2^52 + digit:                                           */\
1359      __fi.i[0] = (uint32)(df_mant_len+1+df_exp_mid) << (df_mant_len-32); /* Vorzeichen 0, Exponent 53 */\
1360      __fi.i[1] = (digit); /* untere 32 Bits setzen (benutzt CL_CPU_BIG_ENDIAN_P !) */\
1361      /* subtrahiere 2^52:                                             */\
1362      __fi.f = __fi.f - (double)(4503599627370496.0L);                   \
1363      /* Hole davon den Exponenten:                                    */\
1364      size_zuweisung ((__fi.i[0] >> (df_mant_len-32)) - df_exp_mid);     \
1365     }
1366 #elif defined(__GNUC__) && defined(__i386__) && !defined(NO_ASM)
1367   #define integerlength32(digit,size_zuweisung)  \
1368     { var uintL _one_position; /* Position der führenden 1                  */\
1369       __asm__("bsrl %1,%0" : "=r" (_one_position) : "rm" ((uint32)(digit)) ); \
1370       size_zuweisung (1+_one_position);                                       \
1371     }
1372 #elif defined(__hppa__) && !defined(NO_ASM)
1373   #define integerlength32(digit,size_zuweisung)  \
1374     size_zuweisung length32(digit);
1375   extern "C" uintL length32 (uintL digit); // extern in Assembler
1376 // Die weiteren kommen von gcc/longlong.h :
1377 #elif defined(__GNUC__) && (defined(__a29k__) || defined(___AM29K__)) && !defined(NO_ASM)
1378   #define integerlength32(digit,size_zuweisung)  \
1379     { var uintL _zero_counter; /* zählt die führenden Nullbits in digit   */\
1380       __asm__("clz %0,%1" : "=r" (_zero_counter) : "r" ((uint32)(digit)) ); \
1381       size_zuweisung (32-_zero_counter);                                    \
1382     }
1383 #elif defined(__GNUC__) && defined(__gmicro__) && !defined(NO_ASM)
1384   #define integerlength32(digit,size_zuweisung)  \
1385     { var uintL _zero_counter; /* zählt die führenden Nullbits in digit      */\
1386       __asm__("bsch/1 %1,%0" : "=g" (_zero_counter) : "g" ((uint32)(digit)) ); \
1387       size_zuweisung (32-_zero_counter);                                       \
1388     }
1389 #elif defined(__GNUC__) && defined(__rs6000__) && !defined(NO_ASM)
1390  #ifdef _AIX
1391   // old assembler syntax
1392   #define integerlength32(digit,size_zuweisung)  \
1393     { var uintL _zero_counter; /* zählt die führenden Nullbits in digit     */\
1394       __asm__("cntlz %0,%1" : "=r" (_zero_counter) : "r" ((uint32)(digit)) ); \
1395       size_zuweisung (32-_zero_counter);                                      \
1396     }
1397  #else
1398   // new assembler syntax
1399   #define integerlength32(digit,size_zuweisung)  \
1400     { var uintL _zero_counter; /* zählt die führenden Nullbits in digit      */\
1401       __asm__("cntlzw %0,%1" : "=r" (_zero_counter) : "r" ((uint32)(digit)) ); \
1402       size_zuweisung (32-_zero_counter);                                       \
1403     }
1404  #endif
1405 #elif defined(__GNUC__) && defined(__m88k__) && !defined(NO_ASM)
1406   #define integerlength32(digit,size_zuweisung)  \
1407     { var uintL _one_position; /* Position der führenden 1                */\
1408       __asm__("ff1 %0,%1" : "=r" (_one_position) : "r" ((uint32)(digit)) ); \
1409       size_zuweisung (1+_one_position);                                     \
1410     }
1411 #elif defined(__GNUC__) && defined(__ibm032__) && !defined(NO_ASM) // RT/ROMP
1412   #define integerlength32(digit,size_zuweisung)  \
1413     { var uintL _x32 = (uint32)(digit);                         \
1414       if (_x32 >= bit(16))                                      \
1415         { integerlength16(_x32>>16,size_zuweisung 16 + ); }     \
1416         else                                                    \
1417         { integerlength16(_x32,size_zuweisung); }               \
1418     }
1419 #else
1420   #define integerlength32(digit,size_zuweisung)  \
1421     { var uintC _bitsize = 1;                                           \
1422       var uintL _x32 = (uint32)(digit);                                 \
1423       /* _x32 hat höchstens 32 Bits.                                   */\
1424       if (_x32 >= bit(16)) { _x32 = _x32>>16; _bitsize += 16; }         \
1425       /* _x32 hat höchstens 16 Bits.                                   */\
1426       if (_x32 >= bit(8)) { _x32 = _x32>>8; _bitsize += 8; }            \
1427       /* _x32 hat höchstens 8 Bits.                                    */\
1428       if (_x32 >= bit(4)) { _x32 = _x32>>4; _bitsize += 4; }            \
1429       /* _x32 hat höchstens 4 Bits.                                    */\
1430       if (_x32 >= bit(2)) { _x32 = _x32>>2; _bitsize += 2; }            \
1431       /* _x32 hat höchstens 2 Bits.                                    */\
1432       if (_x32 >= bit(1)) { /* _x32 = _x32>>1; */ _bitsize += 1; }      \
1433       /* _x32 hat höchstens 1 Bit. Dieses Bit muß gesetzt sein.        */\
1434       size_zuweisung _bitsize;                                          \
1435     }
1436   #define GENERIC_INTEGERLENGTH32
1437 #endif
1438
1439 // Bits einer 64-Bit-Zahl zählen:
1440 // integerlength64(digit,size=);
1441 // setzt size auf die höchste in digit vorkommende Bitnummer.
1442 // > digit: ein uint64 >0
1443 // < size: >0, <=64, mit 2^(size-1) <= digit < 2^size
1444 #ifdef GENERIC_INTEGERLENGTH32
1445   #define integerlength64(digit,size_zuweisung)  \
1446     { var uintC _bitsize = 1;                                           \
1447       var uint64 _x64 = (uint64)(digit);                                \
1448       /* _x64 hat höchstens 64 Bits.                                   */\
1449       if (_x64 >= bit(32)) { _x64 = _x64>>32; _bitsize += 32; }         \
1450       /* _x64 hat höchstens 32 Bits.                                   */\
1451       if (_x64 >= bit(16)) { _x64 = _x64>>16; _bitsize += 16; }         \
1452       /* _x64 hat höchstens 16 Bits.                                   */\
1453       if (_x64 >= bit(8)) { _x64 = _x64>>8; _bitsize += 8; }            \
1454       /* _x64 hat höchstens 8 Bits.                                    */\
1455       if (_x64 >= bit(4)) { _x64 = _x64>>4; _bitsize += 4; }            \
1456       /* _x64 hat höchstens 4 Bits.                                    */\
1457       if (_x64 >= bit(2)) { _x64 = _x64>>2; _bitsize += 2; }            \
1458       /* _x64 hat höchstens 2 Bits.                                    */\
1459       if (_x64 >= bit(1)) { /* _x64 = _x64>>1; */ _bitsize += 1; }      \
1460       /* _x64 hat höchstens 1 Bit. Dieses Bit muß gesetzt sein.        */\
1461       size_zuweisung _bitsize;                                          \
1462     }
1463 #else
1464   #define integerlength64(digit,size_zuweisung)  \
1465     { var uint64 _x64 = (digit);                                        \
1466       var uintC _bitsize64 = 0;                                         \
1467       var uint32 _x32_from_integerlength64;                             \
1468       if (_x64 >= (1ULL << 32)) {                                       \
1469         _x32_from_integerlength64 = _x64>>32; _bitsize64 += 32;         \
1470       } else {                                                          \
1471         _x32_from_integerlength64 = _x64;                               \
1472       }                                                                 \
1473       integerlength32(_x32_from_integerlength64, size_zuweisung _bitsize64 + ); \
1474     }
1475 #endif
1476
1477 // Bits einer uintC-Zahl zählen:
1478 // integerlengthC(digit,size=);
1479 // setzt size auf die höchste in digit vorkommende Bitnummer.
1480 // > digit: ein uintC >0
1481 // < size: >0, <=intCsize, mit 2^(size-1) <= digit < 2^size
1482   #if (intCsize==32)
1483     #define integerlengthC  integerlength32
1484   #endif
1485   #if (intCsize==64)
1486     #define integerlengthC  integerlength64
1487   #endif
1488
1489 // Hintere Nullbits eines 32-Bit-Wortes zählen:
1490 // ord2_32(digit,count=);
1491 // setzt size auf die kleinste in digit vorkommende Bitnummer.
1492 // > digit: ein uint32 >0
1493 // < count: >=0, <32, mit 2^count | digit, digit/2^count ungerade
1494   #if defined(__GNUC__) && defined(__i386__) && !defined(NO_ASM)
1495     #define ord2_32(digit,count_zuweisung)  \
1496       { var uintL _one_position; /* Position der letzten 1                    */\
1497         __asm__("bsfl %1,%0" : "=r" (_one_position) : "rm" ((uint32)(digit)) ); \
1498         count_zuweisung _one_position;                                          \
1499       }
1500     #define FAST_ORD2
1501   #elif defined(__sparc__) && !defined(__sparc64__)
1502     #define ord2_32(digit,count_zuweisung)  \
1503     { var uint32 n = (digit);                                             \
1504       n = n | -n;                                                         \
1505       n = (n<<4) + n;                                                     \
1506       n = (n<<6) + n;                                                     \
1507       n = n - (n<<16); /* or  n = n ^ (n<<16);  or  n = n &~ (n<<16);  */ \
1508       /* static const char ord2_tab [64] = {-1,0,1,12,2,6,-1,13,3,-1,7,-1,-1,-1,-1,14,10,4,-1,-1,8,-1,-1,25,-1,-1,-1,-1,-1,21,27,15,31,11,5,-1,-1,-1,-1,-1,9,-1,-1,24,-1,-1,20,26,30,-1,-1,-1,-1,23,-1,19,29,-1,22,18,28,17,16,-1}; */ \
1509       /* count_zuweisung ord2_tab[n>>26];                              */ \
1510       count_zuweisung "\377\000\001\014\002\006\377\015\003\377\007\377\377\377\377\016\012\004\377\377\010\377\377\031\377\377\377\377\377\025\033\017\037\013\005\377\377\377\377\377\011\377\377\030\377\377\024\032\036\377\377\377\377\027\377\023\035\377\026\022\034\021\020"[n>>26]; \
1511     }
1512     #define FAST_ORD2
1513   #else
1514     // Sei n = ord2(x). Dann ist logxor(x,x-1) = 2^n + (2^n-1) = 2^(n+1)-1.
1515     // Also  (ord2 x) = (1- (integer-length (logxor x (1- x)))) .
1516     #define ord2_32(digit,count_zuweisung)  \
1517       { var uint32 _digit = (digit) ^ ((digit) - 1);    \
1518         integerlength32(_digit,count_zuweisung -1 + )   \
1519       }
1520   #endif
1521
1522 // Hintere Nullbits eines 64-Bit-Wortes zählen:
1523 // ord2_64(digit,count=);
1524 // setzt size auf die kleinste in digit vorkommende Bitnummer.
1525 // > digit: ein uint64 >0
1526 // < count: >=0, <64, mit 2^count | digit, digit/2^count ungerade
1527   // Sei n = ord2(x). Dann ist logxor(x,x-1) = 2^n + (2^n-1) = 2^(n+1)-1.
1528   // Also  (ord2 x) = (1- (integer-length (logxor x (1- x)))) .
1529   #define ord2_64(digit,count_zuweisung)  \
1530     { var uint64 _digit = (digit) ^ ((digit) - 1);      \
1531       integerlength64(_digit,count_zuweisung -1 + )     \
1532     }
1533
1534
1535 // Bits eines Wortes zählen.
1536 // logcount_NN();
1537 // > xNN: ein uintNN
1538 // < xNN: Anzahl der darin gesetzten Bits
1539   // Bits von x8 zählen: (Input x8, Output x8)
1540   #define logcount_8()  \
1541     ( /* x8 besteht aus 8 1-Bit-Zählern (0,1).        */\
1542       x8 = (x8 & 0x55U) + ((x8 & 0xAAU) >> 1),          \
1543       /* x8 besteht aus 4 2-Bit-Zählern (0,1,2).      */\
1544       x8 = (x8 & 0x33U) + ((x8 & 0xCCU) >> 2),          \
1545       /* x8 besteht aus 2 4-Bit-Zählern (0,1,2,3,4).  */\
1546       x8 = (x8 & 0x0FU) + (x8 >> 4)                     \
1547       /* x8 besteht aus 1 8-Bit-Zähler (0,...,8).     */\
1548     )
1549   // Bits von x16 zählen: (Input x16, Output x16)
1550   #define logcount_16()  \
1551     ( /* x16 besteht aus 16 1-Bit-Zählern (0,1).      */\
1552       x16 = (x16 & 0x5555U) + ((x16 & 0xAAAAU) >> 1),   \
1553       /* x16 besteht aus 8 2-Bit-Zählern (0,1,2).     */\
1554       x16 = (x16 & 0x3333U) + ((x16 & 0xCCCCU) >> 2),   \
1555       /* x16 besteht aus 4 4-Bit-Zählern (0,1,2,3,4). */\
1556       x16 = (x16 & 0x0F0FU) + ((x16 & 0xF0F0U) >> 4),   \
1557       /* x16 besteht aus 2 8-Bit-Zählern (0,...,8).   */\
1558       x16 = (x16 & 0x00FFU) + (x16 >> 8)                \
1559       /* x16 besteht aus 1 16-Bit-Zähler (0,...,16).  */\
1560     )
1561   // Bits von x32 zählen: (Input x32, Output x32)
1562   #define logcount_32()  \
1563     ( /* x32 besteht aus 32 1-Bit-Zählern (0,1).              */\
1564       x32 = (x32 & 0x55555555UL) + ((x32 & 0xAAAAAAAAUL) >> 1), \
1565       /* x32 besteht aus 16 2-Bit-Zählern (0,1,2).            */\
1566       x32 = (x32 & 0x33333333UL) + ((x32 & 0xCCCCCCCCUL) >> 2), \
1567       /* x32 besteht aus 8 4-Bit-Zählern (0,1,2,3,4).         */\
1568       x32 = high16(x32)+low16(x32),                             \
1569       /* x32 besteht aus 4 4-Bit-Zählern (0,...,8).           */\
1570       x32 = (x32 & 0x0F0FU) + ((x32 & 0xF0F0U) >> 4),           \
1571       /* x32 besteht aus 2 8-Bit-Zählern (0,...,16).          */\
1572       x32 = (x32 & 0x00FFU) + (x32 >> 8)                        \
1573       /* x32 besteht aus 1 16-Bit-Zähler (0,...,32).          */\
1574     )
1575   // Bits von x64 zählen: (Input x64, Output x64)
1576   #define logcount_64()  \
1577     ( /* x64 besteht aus 64 1-Bit-Zählern (0,1).                             */\
1578       x64 = (x64 & 0x5555555555555555ULL) + ((x64 & 0xAAAAAAAAAAAAAAAAULL) >> 1),\
1579       /* x64 besteht aus 32 2-Bit-Zählern (0,1,2).                           */\
1580       x64 = (x64 & 0x3333333333333333ULL) + ((x64 & 0xCCCCCCCCCCCCCCCCULL) >> 2),\
1581       /* x64 besteht aus 16 4-Bit-Zählern (0,1,2,3,4).                       */\
1582       x64 = (uint32)(x64 + (x64 >> 32)),                                       \
1583       /* x64 besteht aus 8 4-Bit-Zählern (0,...,8).                          */\
1584       x64 = (x64 & 0x0F0F0F0FUL) + ((x64 & 0xF0F0F0F0UL) >> 4),                \
1585       /* x64 besteht aus 4 8-Bit-Zählern (0,...,16).                         */\
1586       x64 = (x64 & 0x00FF00FFU) + ((x64 & 0xFF00FF00U) >> 8),                  \
1587       /* x64 besteht aus 2 16-Bit-Zählern (0,...,32).                        */\
1588       x64 = (x64 & 0x0000FFFFU) + (x64 >> 16)                                  \
1589       /* x64 besteht aus 1 16-Bit-Zähler (0,...,64).                         */\
1590     )
1591
1592 }  // namespace cln
1593
1594 #endif /* _CL_LOW_H */