digit.h

A Library of Efficient Data Types and Algorithms,封装了常用的ADT及其相关算法的软件包

/*******************************************************************************
+
+  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_NAMESPACE

typedef 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_asm
inline digit digMul(digit a, digit b, digit* high);
#else
digit digMul(digit a, digit b, digit* high);
#endif


inline digit digDivTwoByOne(digit high, digit low, digit D, digit* q);
#ifdef WORD_LENGTH_32
inline digit digDivTwoByOne32(digit high, digit low, digit D, digit* q);
#endif
inline 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_asm
inline 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_32
inline 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;
}
#endif


void 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