c_lip.c

密码大家Shoup写的数论算法c语言实现

   _lo = _b*_s; \
   _hi = (long) (((double) _b)*_ds); \
}

/* shift is -1..+3 */
#define zam_loop(a,t,nb) \
{ \
   long _a = (a), _t = (t), _nb = (nb); \
   long _vv; \
   double _yy; \
   _vv = _nb*_s; \
   _yy = ((double) _nb)*_ds; \
   _lo = _lo + _a + _t; \
   (t) = _hi + ((_lo - (_hi<<NTL_NBITS)) >> NTL_NBITS); \
   (a) = _lo & NTL_RADIXM; \
   _lo = _vv; \
   _hi = (long) _yy; \
}

/* shift is -1..+2 */
#define zsx_loop(a,t,nb) \
{ \
   long _t = (t), _nb = (nb); \
   long _vv; \
   double _yy; \
   _vv = _nb*_s; \
   _yy = ((double) _nb)*_ds; \
   _lo = _lo + _t; \
   (t) = _hi + ((_lo - (_hi<<NTL_NBITS)) >> NTL_NBITS); \
   (a) = _lo & NTL_RADIXM; \
   _lo = _vv; \
   _hi = (long) _yy; \
}

/* shift is -3..+1 */
#define zam_subloop(a,t,nb) \
{ \
   long _a = (a), _t = (t), _nb = (nb); \
   long _vv; \
   double _yy; \
   _vv = _nb*_s; \
   _yy = ((double) _nb)*_ds; \
   _lo = _a + _t - _lo; \
   (t) = ((_lo + (_hi<<NTL_NBITS)) >> NTL_NBITS) - _hi; \
   (a) = _lo & NTL_RADIXM; \
   _lo = _vv; \
   _hi = (long) _yy; \
}

/* shift is -1..+3 */
#define zam_finish(a,t) \
{ \
   long _a = (a), _t = (t); \
   _lo = _lo + _a + _t; \
   (t) = _hi + ((_lo - (_hi<<NTL_NBITS)) >> NTL_NBITS); \
   (a) = _lo & NTL_RADIXM; \
}

/* shift is -1..+2 */
#define zsx_finish(a,t) \
{ \
   long _t = (t); \
   _lo = _lo + _t; \
   (t) = _hi + ((_lo - (_hi<<NTL_NBITS)) >> NTL_NBITS); \
   (a) = _lo & NTL_RADIXM; \
}

/* shift is -3..+1 */
#define zam_subfinish(a,t) \
{ \
   long _a = (a), _t = (t); \
   _lo = _a + _t - _lo; \
   (t) = ((_lo + (_hi<<NTL_NBITS)) >> NTL_NBITS) - _hi; \
   (a) = _lo & NTL_RADIXM; \
}

#else

/* right shift is not arithmetic */

/* shift is 0..4 */
#define zaddmulp(a, b, d, t) \
{ \
   long _a = (a), _b = (b), _d = (d), _t = (t); \
   long _t1 =  _b*_d + _a + _t; \
   long _t2 = (long) ( ((double) _b)*(((double) _d)*NTL_FRADIX_INV) ) - 1; \
   (t) = _t2 + (((unsigned long)(_t1 - (_t2 << NTL_NBITS))) >> NTL_NBITS); \
   (a) = _t1 & NTL_RADIXM; \
}

#define zxmulp(a, b, d, t) \
{ \
   long _b = (b), _d = (d), _t = (t); \
   long _t1 =  _b*_d + _t; \
   long _t2 = (long) ( ((double) _b)*(((double) _d)*NTL_FRADIX_INV) ) - 1; \
   (t) = _t2 + (((unsigned long)(_t1 - (_t2 << NTL_NBITS))) >> NTL_NBITS); \
   (a) = _t1 & NTL_RADIXM; \
}

/* shift is 0..3 */
#define zaddmulpsq(a, b, t) \
{ \
   long _a = (a), _b = (b), _t = (t); \
   long _t1 = _b*_b + _a; \
   long _t2 = (long) ( ((double) _b)*(((double) _b)*NTL_FRADIX_INV) ) - 1; \
   (t) = _t2 + (((unsigned long)(_t1 - (_t2 << NTL_NBITS))) >> NTL_NBITS); \
   (a) = _t1 & NTL_RADIXM; \
}

#define zam_decl double _ds; long _hi, _lo, _s;

#define zam_init(b,s) \
{ \
   long _b = (b); \
   _s = (s); \
   _ds = _s*NTL_FRADIX_INV; \
   _lo = _b*_s; \
   _hi = (long) (((double) _b)*_ds); \
}

/* shift is 0..4 */
#define zam_loop(a,t,nb) \
{ \
   long _a = (a), _t = (t), _nb = (nb); \
   long _vv; \
   double _yy; \
   _vv = _nb*_s; \
   _yy = ((double) _nb)*_ds; \
   _hi--; \
   _lo = _lo + _a + _t; \
   (t) = _hi + (((unsigned long)(_lo - (_hi<<NTL_NBITS))) >> NTL_NBITS); \
   (a) = _lo & NTL_RADIXM; \
   _lo = _vv; \
   _hi = (long) _yy; \
}

/* shift is 0..3 */
#define zsx_loop(a,t,nb) \
{ \
   long _t = (t), _nb = (nb); \
   long _vv; \
   double _yy; \
   _vv = _nb*_s; \
   _yy = ((double) _nb)*_ds; \
   _hi--; \
   _lo = _lo + _t; \
   (t) = _hi + (((unsigned long)(_lo - (_hi<<NTL_NBITS))) >> NTL_NBITS); \
   (a) = _lo & NTL_RADIXM; \
   _lo = _vv; \
   _hi = (long) _yy; \
}

/* shift is 0..4 */
#define zam_subloop(a,t,nb) \
{ \
   long _a = (a), _t = (t), _nb = (nb); \
   long _vv; \
   double _yy; \
   _vv = _nb*_s; \
   _yy = ((double) _nb)*_ds; \
   _hi += 3; \
   _lo = _a + _t - _lo; \
   (t) = (((unsigned long)(_lo + (_hi<<NTL_NBITS))) >> NTL_NBITS) - _hi; \
   (a) = _lo & NTL_RADIXM; \
   _lo = _vv; \
   _hi = (long) _yy; \
}

/* shift is 0..4 */
#define zam_finish(a,t) \
{ \
   long _a = (a), _t = (t); \
   _lo = _lo + _a + _t; \
   _hi--; \
   (t) = _hi + (((unsigned long)(_lo - (_hi<<NTL_NBITS))) >> NTL_NBITS); \
   (a) = _lo & NTL_RADIXM; \
}

/* shift is 0..3 */
#define zsx_finish(a,t) \
{ \
   long _t = (t); \
   _lo = _lo + _t; \
   _hi--; \
   (t) = _hi + (((unsigned long)(_lo - (_hi<<NTL_NBITS))) >> NTL_NBITS); \
   (a) = _lo & NTL_RADIXM; \
}

/* shift is 0..4 */
#define zam_subfinish(a,t) \
{ \
   long _a = (a), _t = (t); \
   _hi += 3; \
   _lo = _a + _t - _lo; \
   (t) = (((unsigned long)(_lo + (_hi<<NTL_NBITS))) >> NTL_NBITS) - _hi; \
   (a) = _lo & NTL_RADIXM; \
}
#endif

#endif

#endif

#endif

static
void zaddmulone(long *lama, long *lamb)
{ 
   long lami; 
   long lams = 0; 
 
   lams = 0; 
   for (lami = (*lamb++); lami > 0; lami--) { 
      lams += (*lama + *lamb++); 
      *lama++ = lams & NTL_RADIXM; 
      lams >>= NTL_NBITS; 
   } 
   *lama += lams; 
}

#if (NTL_ARITH_RIGHT_SHIFT)
static
void zsubmulone(long *lama, long *lamb)
{ 
   long lami; 
   long lams = 0; 
 
   lams = 0; 
   for (lami = (*lamb++); lami > 0; lami--) { 
      lams += (*lama - *lamb++); 
      *lama++ = lams & NTL_RADIXM; 
      lams >>= NTL_NBITS; 
   } 
   *lama += lams; 
}
#else
static
void zsubmulone(long *lama, long *lamb)
{ 
   long lami; 
   long lams = 0; 
 
   lams = 0; 
   for (lami = (*lamb++); lami > 0; lami--) { 
      lams = *lama - *lamb++ - lams; 
      *lama++ = lams & NTL_RADIXM; 
      lams = (lams < 0);
   } 
   *lama -= lams; 
}

#endif

#if (defined(NTL_SINGLE_MUL))

static
void zaddmul(long ams, long *ama, long *amb) 
{ 
   long carry = 0; 
   long i = (*amb++);
   double dams = (double) ams;
   double xx;
   double yy;
   unsigned long lo_wd, lo;
   unsigned long hi_wd, hi;

   xx  =  ((double) (*amb++))*dams + DENORMALIZE;
   for (; i > 1; i--) { 
      yy =  ((double) (*amb++))*dams +DENORMALIZE;
      NTL_FetchHiLo(hi_wd, lo_wd, xx);
      lo = lo_wd & 0x3FFFFFF;
      hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
      lo = lo + (*ama) + carry;
      *ama = lo & 0x3FFFFFF;
      carry = hi + (lo >> 26);
      ama++; 
      xx = yy;
   } 

   NTL_FetchHiLo(hi_wd, lo_wd, xx);
   lo = lo_wd & 0x3FFFFFF;
   hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
   lo = lo + (*ama) + carry;
   *ama = lo & 0x3FFFFFF;
   carry = hi + (lo >> 26);
   ama++; 
   *ama += carry; 
}

static
void zsxmul(long ams, long *ama, long *amb) 
{ 
   long carry = 0; 
   long i = (*amb++);
   double dams = (double) ams;
   double xx;
   double yy;
   unsigned long lo_wd, lo;
   unsigned long hi_wd, hi;

   xx  =  ((double) (*amb++))*dams + DENORMALIZE;
   for (; i > 1; i--) { 
      yy =  ((double) (*amb++))*dams +DENORMALIZE;
      NTL_FetchHiLo(hi_wd, lo_wd, xx);
      lo = lo_wd & 0x3FFFFFF;
      hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
      lo = lo + carry;
      *ama = lo & 0x3FFFFFF;
      carry = hi + (lo >> 26);
      ama++; 
      xx = yy;
   } 

   NTL_FetchHiLo(hi_wd, lo_wd, xx);
   lo = lo_wd & 0x3FFFFFF;
   hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
   lo = lo + carry;
   *ama = lo & 0x3FFFFFF;
   carry = hi + (lo >> 26);
   ama++; 
   *ama = carry; 
}

static
void zaddmulsq(long ams, long *ama, long *amb) 
{ 
   long carry = 0; 
   long i = ams;
   double dams = (double) (*amb++);
   double xx;
   double yy;
   unsigned long lo, lo_wd;
   unsigned long hi, hi_wd;

   xx =  ((double) (*amb++))*dams + DENORMALIZE;
   for (; i > 1; i--) { 
      yy =  ((double) (*amb++))*dams + DENORMALIZE;
      NTL_FetchHiLo(hi_wd, lo_wd, xx);
      lo = lo_wd & 0x3FFFFFF;
      hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
      lo = lo + (*ama) + carry;
      *ama = lo & 0x3FFFFFF;
      carry = hi + (lo >> 26);
      ama++; 
      xx = yy;
   } 
   if (i==1) {
      NTL_FetchHiLo(hi_wd, lo_wd, xx);
      lo = lo_wd & 0x3FFFFFF;
      hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
      lo = lo + (*ama) + carry;
      *ama = lo & 0x3FFFFFF;
      carry = hi + (lo >> 26);
      ama++; 
   }
   *ama += carry; 
}

#else

#if (defined(NTL_AVOID_FLOAT) || defined(NTL_LONG_LONG))
static
void zaddmul(long lams, long *lama, long *lamb)
{
        long lami;
        long lamcarry = 0;

        for (lami = (*lamb++); lami > 0; lami--)
        {
                zaddmulp(*lama, *lamb, lams, lamcarry);
                lama++;
                lamb++;
        }
        *lama += lamcarry;
}

static
void zsxmul(long lams, long *lama, long *lamb)
{
        long lami;
        long lamcarry = 0;

        for (lami = (*lamb++); lami > 0; lami--)
        {
                zxmulp(*lama, *lamb, lams, lamcarry);
                lama++;
                lamb++;
        }
        *lama = lamcarry;
}

static
void zaddmulsq(long lsqi, long *lsqa, long *lsqb)
{
        long lsqs = *(lsqb);
        long lsqcarry = 0;

        lsqb++;
        for (; lsqi > 0; lsqi--)
        {
                zaddmulp(*lsqa, *lsqb, lsqs, lsqcarry);
                lsqa++;
                lsqb++;
        }
        *lsqa += lsqcarry;
}



#else

static
void zaddmul(long lams, long *lama, long *lamb)
{ 
   long lami = (*lamb++)-1; 
   long lamcarry = 0; 
   zam_decl;

   zam_init(*lamb, lams);
   lamb++;

 
   for (; lami > 0; lami--) { 
      zam_loop(*lama, lamcarry, *lamb);
      lama++; 
      lamb++; 
   } 
   zam_finish(*lama, lamcarry);
   lama++;
        *lama += lamcarry; 
}



static
void zsxmul(long lams, long *lama, long *lamb)
{ 
   long lami = (*lamb++)-1; 
   long lamcarry = 0; 
   zam_decl;

   zam_init(*lamb, lams);
   lamb++;

 
   for (; lami > 0; lami--) { 
      zsx_loop(*lama, lamcarry, *lamb);
      lama++; 
      lamb++; 
   } 
   zsx_finish(*lama, lamcarry);
   lama++;
   *lama = lamcarry; 
}



static
void zaddmulsq(long lsqi, long *lsqa, long *lsqb)
{ 
   long lsqs; 
   long lsqcarry; 
   zam_decl

   if (lsqi <= 0) return;

   lsqs = *lsqb;
   lsqcarry = 0;

   lsqb++; 
   zam_init(*lsqb, lsqs);
   lsqb++;
   lsqi--;
   for (; lsqi > 0; lsqi--) { 
      zam_loop(*lsqa, lsqcarry, *lsqb);
      lsqa++; 
      lsqb++; 
   } 
   zam_finish(*lsqa, lsqcarry);
   lsqa++;
   *lsqa += lsqcarry; 
}


#endif
#endif


#if (defined(NTL_SINGLE_MUL))

#if (NTL_ARITH_RIGHT_SHIFT)

static
void zsubmul(long ams, long *ama, long *amb) 
{ 
   long carry = 0; 
   long i = (*amb++);
   double dams = (double) ams;
   double xx;
   double yy;
   unsigned long lo_wd, lo;
   unsigned long hi_wd, hi;

   xx  =  ((double) (*amb++))*dams + DENORMALIZE;
   for (; i > 1; i--) { 
      yy =  ((double) (*amb++))*dams +DENORMALIZE;
      NTL_FetchHiLo(hi_wd, lo_wd, xx);
      lo = lo_wd & 0x3FFFFFF;
      hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
      lo = (*ama) + carry - lo;
      *ama = lo & 0x3FFFFFF;
      carry = (((long)lo) >> 26) - hi;
      ama++; 
      xx = yy;
   } 

   NTL_FetchHiLo(hi_wd, lo_wd, xx);
   lo = lo_wd & 0x3FFFFFF;
   hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
   lo = (*ama) + carry - lo;
   *ama = lo & 0x3FFFFFF;
   carry = (((long)lo) >> 26) - hi;
   ama++; 
   *ama += carry; 
}

#else

static
void zsubmul(long ams, long *ama, long *amb) 
{ 
   long carry = 0; 
   long i = (*amb++);
   double dams = (double) ams;
   double xx;
   double yy;
   unsigned long lo_wd, lo;
   unsigned long hi_wd, hi;

   xx  =  ((double) (*amb++))*dams + DENORMALIZE;
   for (; i > 1; i--) { 
      yy =  ((double) (*amb++))*dams +DENORMALIZE;
      NTL_FetchHiLo(hi_wd, lo_wd, xx);
      lo = lo_wd & 0x3FFFFFF;
      hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
      lo = (*ama) + carry - lo;
      *ama = lo & 0x3FFFFFF;
      carry = ((lo + (1L << 27)) >> 26) - hi - 2;
      ama++; 
      xx = yy;
   } 

   NTL_FetchHiLo(hi_wd, lo_wd, xx);
   lo = lo_wd & 0x3FFFFFF;
   hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
   lo = (*ama) + carry - lo;
   *ama = lo & 0x3FFFFFF;
   carry = ((lo + (1L << 27)) >> 26) - hi - 2;
   ama++; 
   *ama += carry; 
}

#endif

#else

#if (defined(NTL_LONG_LONG))

static
void zsubmul(long lams, long *lama, long *lamb)
{
        long lami;
        long lamcarry = 0;

        lams = -lams;

        for (lami = (*lamb++); lami > 0; lami--)
        {
                zaddmulp(*lama, *lamb, lams, lamcarry);
                lama++;
                lamb++;
        }
        *lama += lamcarry;
}

#elif (defined(NTL_AVOID_FLOAT))

static void
zsubmul(
        long r,
        _ntl_verylong a,
        _ntl_verylong b
        )
{
        long rd = NTL_RADIX - r;
        long i;
        long carry = NTL_RADIX;

        for (i = (*b++); i > 0; i--)
        {
                zaddmulp(*a, *b, rd, carry);
                a++;
                carry += NTL_RADIXM - (*b++);
        }
        *a += carry - NTL_RADIX; /* unnormalized */
}



#else

static
void zsubmul(long lams, long *lama, long *lamb)
{ 
   long lami = (*lamb++)-1; 
   long lamcarry = 0; 
   zam_decl;

   zam_init(*lamb, lams);
   lamb++;

 
   for (; lami > 0; lami--) { 
      zam_subloop(*lama, lamcarry, *lamb);
      lama++; 
      lamb++; 
   } 
   zam_subfinish(*lama, lamcarry);
   lama++;
   *lama += lamcarry; 
}

#endif

#endif




/*
	zdiv21 returns quot, numhigh so

	quot = (numhigh*NTL_RADIX + numlow)/denom;
	numhigh  = (numhigh*NTL_RADIX + numlow)%denom;

Assumes 0 <= numhigh < denom < NTL_RADIX and 0 <= numlow < NTL_RADIX.
*/



#define zdiv21(numhigh, numlow, denom, deninv, quot) \
{ \
   long lr21; \
   long lq21 = (long) (((NTL_FRADIX * (double) (numhigh)) \
          + (double) (numlow)) * (deninv)); \
   long lp21; \
   MulLo(lp21, lq21, denom); \
   lr21 = (numhigh << NTL_NBITS) + numlow - lp21; \
   if (lr21 < 0) { \
      lq21--; \
      lr21 += denom; \