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