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