/*******************************************************************************++  LEDA 4.5  +++  digit.h+++  Copyright (c) 1995-2004+  by Algorithmic Solutions Software GmbH+  All rights reserved.+ *******************************************************************************/// $Revision: 1.3 $  $Date: 2004/02/06 11:19:32 $#ifndef DIGIT_H#define DIGIT_H#include <LEDA/std/math.h>#ifdef DIGIT_CHECKS_ON#include <LEDA/std/assert.h>#endif#if !defined(WORD_LENGTH_32) && !defined(WORD_LENGTH_64)#error WORD_LENGTH must be either 32 or 64#endif#include <LEDA/numbers/digit_asm.h>#include <LEDA/numbers/numbers_common.h>LEDA_BEGIN_NAMESPACEtypedef unsigned long digit;typedef unsigned int  digit_sz;extern __exportF digit digPower(digit x, digit_sz n);inline digit_sz digLeadingZeros(digit x);inline digit_sz digTrailingZeros(digit x);inline digit_sz digNextPower(digit x);inline digit digModReduce(digit x, digit p) { return x%p; } static const digit_sz  DIGIT_SIZE    = sizeof(digit);static const digit_sz  DIGIT_LENGTH  = CHAR_BIT*sizeof(digit);static const digit_sz  HALF_LENGTH   = DIGIT_LENGTH/2;static const digit MAX_DIGIT     = ~0UL;static const digit MIN_NDIGIT    = 1UL<<(DIGIT_LENGTH-1);static const digit HALF_BASE     = 1UL<<HALF_LENGTH;static const digit MAX_HALFDIGIT = HALF_BASE-1;static const digit LEFT_MASK     = MAX_DIGIT-MAX_HALFDIGIT;static const digit RIGHT_MASK    = MAX_HALFDIGIT;#define LOW_WORD(w)  (w & RIGHT_MASK)#define HIGH_WORD(w) (w >> HALF_LENGTH)static const double digit_base = ldexp(1.0,DIGIT_LENGTH);static const double digit_pow32 = ldexp(1.0,32);static const double digit_pow64 = ldexp(1.0,64);static const double digit_roundup_factor   = 1+ldexp(8.0,-53);static const double digit_rounddown_factor = 1-ldexp(8.0,-53);static const digit_sz digit_decimal_base_length10 = (digit_sz) log10(digit_base);static const digit_sz digit_decimal_base_length2  = digit_sz(digit_decimal_base_length10/log10(2.0));static const digit digit_decimal_base = digPower(10,digit_decimal_base_length10);static const digit_sz  digit_shift_decimal_to_normal=DIGIT_LENGTH-1-digit_decimal_base_length2;static const digit digit_shifted_decimal_base=digit_decimal_base<<digit_shift_decimal_to_normal;#ifndef digMul_asminline digit digMul(digit a, digit b, digit* high);#elsedigit digMul(digit a, digit b, digit* high);#endifinline digit digDivTwoByOne(digit high, digit low, digit D, digit* q);#ifdef WORD_LENGTH_32inline digit digDivTwoByOne32(digit high, digit low, digit D, digit* q);#endifinline digit digDivTwoByOneNormal(digit high, digit low, digit D, digit* q);inline digit digDivThreeByTwoNormal             (digit A1, digit A2, digit A3, digit B1, digit B2);char* digToDecimal(digit x, char* cp);char* digFromDecimal(digit& x, char* cp);inline digit_sz digLeadingZeros(digit x){  digit_sz rshift=0;  digit y = x>>HALF_LENGTH;  if (y) {    rshift=HALF_LENGTH;    x=y;  }    while(x) {    y=x;    x>>=2;    rshift+=2;  }  if (y==1)     rshift--;       return DIGIT_LENGTH-rshift;}inline digit_sz digTrailingZeros(digit x){  digit_sz lshift=0;  digit y = x<<HALF_LENGTH;  if (y) {    lshift=HALF_LENGTH;    x=y;  }    while(x) {    y=x;    x<<=2;    lshift+=2;  }  if (y==(digit(1)<<(DIGIT_LENGTH-1)))     lshift--;       return DIGIT_LENGTH-lshift;}inline digit_sz digNextPower(digit s){ digit p = 1;  while (p < s) p <<= 1;  return digit_sz(p);}#ifndef digMul_asminline digit digMul(digit a, digit b, digit* high)  // computes a*b = (high,low) and returns low{ digit al = LOW_WORD(a);   // low part of a  digit ah = HIGH_WORD(a);  // high part of a  digit bl = LOW_WORD(b);   // low part of b  digit bh = HIGH_WORD(b);  // low part of b  digit c,L,H;  c = bl*al;             L = LOW_WORD(c);  c = HIGH_WORD(c) + bl*ah;  H = HIGH_WORD(c);  c = LOW_WORD(c)  + bh*al;  L |= (LOW_WORD(c) << HALF_LENGTH);  *high =  H + HIGH_WORD(c) + bh*ah;  return L;}#endif#ifdef WORD_LENGTH_32inline digit digDivTwoByOne32(digit high, digit low, digit D, digit* q){   // precondition: digit size = 32 bit, high < D  // *q = [high,low] / D, returns remainder  // we can compute it with double precision floating point division  double H = high*digit_base;  double L = low;  digit Q = digit(((H+L)/D)*digit_rounddown_factor);  double Q_high = (double) (Q & LEFT_MASK);  double Q_low  = (double) (Q & RIGHT_MASK);  double diff = ((H - D*Q_high) - D*Q_low) + L;  if (diff >= D)     {      Q++;      diff-=D;    }  *q = Q;  return (digit) diff;}#endifvoid digDivThreeHalvesByTwoSpecial(digit A12, digit a3, digit B, digit b1, digit b2, digit& q, digit& R);inline void digDivThreeHalvesByTwoNormal(digit A12, digit a3, digit B, digit b1, digit b2, digit& q, digit& R){  // small letters = half-digits  // capital letters = digits#ifdef DIGIT_CHECKS_ON  assert(b1>>(HALF_LENGTH-1)!=0);  assert(A12<B);#endif  q = A12/b1;  if (HIGH_WORD(q))    digDivThreeHalvesByTwoSpecial(A12,a3,B,b1,b2,q,R);  else {    digit r = A12-q*b1;    digit A = (r<<HALF_LENGTH)+a3;    digit D = q*b2;    R = A-D;    if (A<D) // remainder negative => q too large at one by least        {          q--;          R+=B; // R = B + R; -B<=R<B holds                    if (R>B) // q still too large            {              q--;              R+=B; // this is the correct R            }        }    }#ifdef DIGIT_CHECKS_ON#ifdef WORD_LENGTH_32   double xd = double(A12)*HALF_BASE+a3;   double bd = B;   assert (q==digit(xd/bd));   assert (R==digit(xd-bd*q));#else   digit Ph,Pl;   Pl  = digMul(B,q,&Ph);   Pl += R;   if (Pl < R)     Ph++;   digit _A12 = (Ph<<HALF_LENGTH)+HIGH_WORD(Pl);   digit _a3 = LOW_WORD(Pl);   assert(_A12==A12);   assert(_a3==a3);#endif#endif}inline digit digDivTwoByOneNormal(digit A12, digit A34, digit B12, digit* Q){  digit b1 = HIGH_WORD(B12);  digit b2 =  LOW_WORD(B12);  digit a3 = HIGH_WORD(A34);  digit a4 =  LOW_WORD(A34);  digit q1, q2, R1;  digit R;  digDivThreeHalvesByTwoNormal(A12,a3,B12,b1,b2,q1,R1);  digDivThreeHalvesByTwoNormal(R1 ,a4,B12,b1,b2,q2,R );  *Q = (q1<<HALF_LENGTH)+q2;#ifdef DIGIT_CHECKS_ON  digit _A12, _A34;  _A34=digMul(*Q,B12,&_A12);  _A34+=R;  if (_A34<R)    _A12++;  assert(_A12==A12);  assert(_A34==A34);#endif  return R;}inline digit digDivThreeByTwoNormal             (digit A1, digit A2, digit A3, digit B1, digit B2){  digit Q,P1,P2,R1,R2;  R1 = digDivTwoByOneNormal(A1,A2,B1,&Q);   P2 =  digMul(Q,B2,&P1);     R2 = A3;  if (P1>R1||(P1==R1&&P2>R2))    {      Q--;      // [R1,R2]+=[B1,B2]      R2+=B2;      if (R2<B2)         R1++;      R1+=B1;      if (R1>=B1) // no overflow in R1        if (P1>R1||(P1==R1&&P2>R2))          Q--;    }  return Q;}inline digit digDivTwoByOne(digit A12, digit A34, digit B12, digit* Q){  digit_sz lshift = digLeadingZeros(B12);  if (lshift) {    B12<<=lshift;    A12<<=lshift;    A12+=(A34>>(DIGIT_LENGTH-lshift));    A34<<=lshift;  }  return digDivTwoByOneNormal(A12,A34,B12,Q)>>lshift;}bool digIsPrime(digit p);bool digIsPseudoPrime(digit p, digit K);digit digGeneratePrimeWithBound(digit a);digit digGeneratePseudoPrimeWithBound(digit a, digit K);int digFromDouble(double x, digit& hi, digit& lo);LEDA_END_NAMESPACE#endif