mytest01.c

大量数学算法,主要用于密码学中.实现十分完备,用起来十分方便.配有测试程序,大家可以自己试用,使用十分简单

#include "tommath.h"

//////////////////////////////////////////////////////////////////////////////////////////
//大整数fp_int的结构(二进)显示
//
/* 2-64进位制的数字  chars used in radix conversions */
const char *mp_s_rmap = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";

//大整数的定义：
//
//mp_digit被定义为无符号长整数，在32为机上，=32比特
typedef unsigned long      mp_digit;
//缺省分配给mp_digit的大小-决定了计算的最大精确度
#define MP_PREC            64     /* default digits of precision */
//每个mp_digit的32比特里仅用28比特进行计算
/* default case is 28-bit digits, defines MP_28BIT as a handy macro to test */
#define DIGIT_BIT          28 
#define MP_28BIT

//大整数mp_int结构
//
typedef struct  {
    int used, alloc, sign;
    mp_digit *dp;//每块(mp_digit,32bite)28bite
} mp_int;

//mp_int大整数结构的显示
//
void ddraw(mp_int *a)
{
  int x;
  //used-用了多少个mp_digit，alloc-分配多少mp_digit(缺省64)，sign-正负（0为非负，1为负）
  printf("%d, %d, %d ", a->used, a->alloc, a->sign);
  //显示使用的mp_digit
  for (x = a->used - 1; x >= 0; x--) { 
      //按16进显示
	  printf("%08lx ", a->dp[x]);
  }
  printf("\n");
}
///////////////////////////////////////////////////////////////////////
//大整数的输入输出
//
//不调用其它函数
// 
/* 为大整数结构a分配内存，缺省分配MP_PREC=64个mp_digit，
   因此，共有32*64比特用来表示大整数，而实际可处理28*64=1792比特 
*/
int mp_init (mp_int * a);

/* init an mp_init for a given size */
int mp_init_size (mp_int * a, int size);

/** 将大整数重新初始化，释放大整数所占内存 clear one (frees) **/
void mp_clear (mp_int * a);

/** 修剪未使用的mp_digit，调用本函数可剪除开头的那些被视为“使用-used”过的但为0的mp_digit */
void mp_clamp (mp_int * a);

/* 将a设置为0   set to zero */
void mp_zero (mp_int * a);
/* set to a digit (调用：mp_zero) */
void mp_set (mp_int * a, mp_digit b);

/* get the lower 32-bits of an mp_int */
unsigned long mp_get_int(mp_int * a) ;

/* 内存使用1：Reducing Memory Usage */
int mp_shrink (mp_int * a);
/* 内存使用2：Reducing Memory Usage
   在大整数a里增加数组dp的元素mp_digit:   Adding additional digits
   在大整数a里,used是表示a的实际mp_digit数，alloc是a可使用的mp_digit数，一般used<=alloc
   如果a需要更多的mp_digit数，可用本函数实现：
   1 如果size<=a->alloc，则本函数不作任何事
   2 如果size>a->alloc，则调用本函数后a的alloc变为新值，例如：
        mp_grow(&a, number.alloc + 20)
使a增加了20个mp_digit
*/
int mp_grow (mp_int * a, int size);
/* 右移位 shift right a certain amount of digits */
void mp_rshd (mp_int * a, int b);
/* returns the number of bits in an int */
int mp_count_bits (mp_int * a);
/* reverse an array, used for radix code */
void bn_reverse (unsigned char *s, int len);
/* swap the elements of two integers, for cases where you can't simply swap the mp_int pointers around  */
void mp_exch (mp_int * a, mp_int * b);
//
static int s_is_power_of_two(mp_digit b, int *p);

/////////////////////////////////////////////////////////////////
//调用上述函数之一
//
/* copy, b = a 本函数调用函数mp_grow */
int mp_copy (mp_int * a, mp_int * b);
/* 左移位  shift left a certain amount of digits   本函数调用:mp_grow */
int mp_lshd (mp_int * a, int b);
/* computes a = 2**b 调用 mp_grow */
int mp_2expt (mp_int * a, int b);
/* multiply by a digit 调用：mp_grow, mp_clamp */
int mp_mul_d (mp_int * a, mp_digit b, mp_int * c);
/* creates "a" then copies b into it    调用mp_init, mp_copy */
int mp_init_copy (mp_int * a, mp_int * b);
/* shift left by a certain bit count     调用：mp_copy, mp_grow, mp_lshd  */
int mp_mul_2d (mp_int * a, int b, mp_int * c);
/* divide by three (based on routine from MPI and the GMP\ manual) 
调用：mp_init_size,     mp_clamp,     mp_exch;  mp_clear;
*/
int mp_div_3 (mp_int * a, mp_int *c, mp_digit * d);
/* calc a value mod 2**b  调用：mp_zero, mp_copy,*/
int mp_mod_2d (mp_int * a, int b, mp_int * c);

///////////////////////////////////////////////////////////////////////
//调用上述两个层次
//
//mp_add_d与mp_sub_d互相调用
//
/* single digit addition  调用：mp_grow, mp_clamp  mp_sub_d */
int mp_add_d (mp_int * a, mp_digit b, mp_int * c);
/* single digit subtraction 调用：mp_grow, mp_clamp      mp_add_d */
int mp_sub_d (mp_int * a, mp_digit b, mp_int * c);
//
//////////////
/* set a 32-bit const   本函数调用函数:mp_zero, mp_clamp  mp_mul_2d */
int mp_set_int (mp_int * a, unsigned long b);
/* read a string [ASCII] in a given radix   本函数调用函数:mp_iszero  mp_mul_d, mp_add_d */
int mp_read_radix (mp_int * a, char *str, int radix);
/* shift right by a certain bit count (store quotient in c, optional remainder in d) 
  调用：mp_copy, mp_zero, mp_init, mp_clear, mp_rshd  mp_mod_2d */
int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d);
/* single digit division (based on routine from MPI) 
调用：mp_copy, mp_iszero, s_is_power_of_two mp_div_2d, mp_div_3 */
int mp_div_d (mp_int * a, mp_digit b, mp_int * c, mp_digit * d);
/* stores a bignum as a ASCII string in a given radix (2..64) 
  本函数调用函数:mp_iszero, mp_init_copy, mp_clear, bn_reverse mp_div_d */
int mp_toradix (mp_int * a, char *str, int radix);

///////////////////////////////////////////////////////////////////////////
//显示大整数
//调用mp_toradix(a, buf, k)将大整数a转换为十进制数字串放入buf,再用printf显示
//调用：mp_toradix
void n10draw(mp_int *a, char *name)
{
   char buf[16000];
   printf("%s ", name);
   mp_toradix(a, buf, 10);
   printf("%s\n", buf);
}
//16进制  调用：mp_toradix
void n16draw(mp_int *a, char *name)
{
   char buf[16000];
   printf("%s ", name);
   mp_toradix(a, buf, 16);
   printf("%s\n", buf);
}
//
static void draw10(mp_int *a)
{
   n10draw(a, "10进大整数a=");
}
static void draw16(mp_int *a)
{
   n16draw(a, "16进大整数a=");
}


///////////////////////////////////////////////
////////////////////////////////////////////////
//加、减法
/* high level addition (handles signs) */
int mp_add (mp_int * a, mp_int * b, mp_int * c);
/* high level subtraction (handles signs) */
int mp_sub (mp_int * a, mp_int * b, mp_int * c);

//用到：s_mp_sub,s_mp_sub,s_mp_add,
/* low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9 */
int s_mp_sub (mp_int * a, mp_int * b, mp_int * c);
/* low level addition, based on HAC pp.594, Algorithm 14.7 */
int s_mp_add (mp_int * a, mp_int * b, mp_int * c);
/* compare maginitude of two ints (unsigned) */
int mp_cmp_mag (mp_int * a, mp_int * b);

/////////////////////////////////////////////
//乘、除法
//
/* high level multiplication (handles sign) */
//调用：fast_s_mp_mul_digs, s_mp_mul (a, b, c)

/* 关键代码 Fast (comba) multiplier
 *
 * This is the fast column-array [comba] multiplier.  It is 
 * designed to compute the columns of the product first 
 * then handle the carries afterwards.  This has the effect 
 * of making the nested loops that compute the columns very
 * simple and schedulable on super-scalar processors.
 *
 * This has been modified to produce a variable number of 
 * digits of output so if say only a half-product is required 
 * you don't have to compute the upper half (a feature 
 * required for fast Barrett reduction).
 *
 * Based on Algorithm 14.12 on pp.595 of HAC.
 *
 */
/* size of comba arrays, should be at least 2 * 2**(BITS_PER_WORD - BITS_PER_DIGIT*2) */
#define MP_WARRAY               (1 << (sizeof(mp_word) * CHAR_BIT - 2 * DIGIT_BIT + 1))
int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs);
int mp_mul (mp_int * a, mp_int * b, mp_int * c);

/* integer signed division. 
 * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
 * HAC pp.598 Algorithm 14.20
 *
 * Note that the description in HAC is horribly 
 * incomplete.  For example, it doesn't consider 
 * the case where digits are removed from 'x' in 
 * the inner loop.  It also doesn't consider the 
 * case that y has fewer than three digits, etc..
 *
 * The overall algorithm is as described as 
 * 14.20 from HAC but fixed to treat these cases.
*/
/* compare two ints (signed)*/
int mp_cmp (mp_int * a, mp_int * b);
int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d);
//////////////////////
//模运算
//
/* c = a mod b, 0 <= c < b */
int mp_mod (mp_int * a, mp_int * b, mp_int * c);

///////////////////////
//概率素数
//


char ccc=' ';
///////////////////////////////////////////////////////////////////////////////////

int main(void)
{ 
  //2-64进大整数的赋值与二进、十进显示
  char *sss=malloc(sizeof(char)*1000000);
  //char *str=malloc(sizeof(char)*1000000); 

  // 测试 mp_get_int
  //大整数声明
  mp_int a, a1, a2, a3;
  unsigned long t;
  int i, n;
  int *x;
  unsigned long *u;

  printf("未初始化时整数*x与n,无符号长整数*u的地址\n");  
  printf("%08lx \n", &x);
  printf("%08lx \n", &u);
  printf("%08lx \n", &n);



  printf("未初始化时大整数a的结构(null)\n");  
  ddraw(&a);
  printf("\n");

  //大整数初始化
  mp_init(&a);   mp_init(&a1);   mp_init(&a2);   mp_init(&a3);
  printf("显示初始化后大整数a的结构：\n");  
  ddraw(&a);
  printf("\n");
    
  //
  //4294967296=2**32; 268435456-1=2**28-1
  printf("测试：使用mp_set给a赋mp_digit值2**28-1,2**28-2)\n");
  mp_set(&a,268435456-1);
  ddraw(&a);
  //mp_clear(&a);//
  //ddraw(&a); 
  mp_set(&a,17);
  ddraw(&a);
  printf("mp_set中268435456不可用！\n");
  mp_set(&a,268435456);//不可用！
  ddraw(&a);//0,64,0
  printf("\n");

  printf("测试：使用mp_zero给a赋0，向相当于重新初始化\n");
  mp_zero(&a);
  ddraw(&a);
  printf("mmmmm\n");


  printf("测试1 将随机产生32位无符号数t赋给大整数a：即a=t\n");
  //随机产生32位无符号数t
  srand(time(NULL));
  for(i=0;i<3;++i) {
    t = ((unsigned long)rand()*rand()+1)&0xFFFFFFFF;
    printf("t=%d\n",t);
	mp_set_int(&a,t);
    if (t!=mp_get_int(&a)) {//如果出错： 
      printf("(1)mp_get_int() bad result!\n");
      return 1;
    }
	ddraw(&a);draw10(&a);
  }
  printf("\n");
  
  //a=0
  mp_set_int(&a,0);
  if (mp_get_int(&a)!=0)
  { printf("(2)mp_get_int() bad result!\n");
    return 1;
  }
  printf("\n");
  
  //a=0xffffffff
  mp_set_int(&a,0xffffffff);
  if (mp_get_int(&a)!=0xffffffff)
  { printf("(3)mp_get_int() bad result!\n");
    return 1;
  }
  printf("\n");
  
  //注意：每块(mp_digit)28bite,16进中的7位
  t=268435456-1;//4294967296=2**32
  mp_set_int(&a,t);
  if (mp_get_int(&a)!=t)
  { printf("(4)mp_get_int() bad result!\n");
    return 1;
  }
  printf("268435456-1=2**28-1\n");
  ddraw(&a);draw10(&a);
  
  t=268435456;//268435456;//4294967296=2**32
  mp_set_int(&a,t);
  if (mp_get_int(&a)!=t)
  { printf("(4)mp_get_int() bad result!\n");
    return 1;
  }
  printf("268435456=2**28\n");
  ddraw(&a);draw10(&a);
   
  t=268435456+1;//268435456;//4294967296=2**32
  mp_set_int(&a,t);
  if (mp_get_int(&a)!=t)
  { printf("(4)mp_get_int() bad result!\n");
    return 1;
  }
  ddraw(&a);draw10(&a);

  t=4294967296-1;//=2**32-1
  mp_set_int(&a,t);
  if (mp_get_int(&a)!=t)
  { printf("(4)mp_get_int() bad result!\n");
    return 1;
  }
  ddraw(&a);draw10(&a);

   /////////////////////////
  //将radix进大整数字符串str赋给a
  mp_read_radix (&a, "yeuiewquiwyeuq68768126po000", 64);
  printf("yeuiewquiwyeuq68768126po000(64)=");
  draw10(&a);
  draw16(&a);

  
  
  //1  用fp_read_radix在程序中直接输入n进大整数sss,显示其二进结构
  printf("\n将11进数-a0a1234567890123456789ab01234567890赋给二进的a3显示:\n");
  n=11;
  mp_read_radix(&a3, "-a0a1234567890123456789ab01234567890", n); 
  ddraw(&a3);
  
  printf("\n---------------------------------------------------------------------\n");
  printf("二进a与n进sss的转换：\n\n");
   
  //2.1  用fp_read_radix键盘输入n进大整数sss并赋给a，并显示a
  printf("\n请输入进位制n=:\n");scanf("%d", &n);
  printf("\n请输入该进位制下的数=:\n");scanf("%s", &sss);
  mp_read_radix(&a, &sss, n);
  printf("\n显示此大整数fp_int a=:\n");
  ddraw(&a);//显示a的结构
  //2.2  将此二进的大整数a转换为十进字符串sss并显示  
  printf("\n将此二进的a转换为十进制字符串sss并显示:\n");
  mp_toradix(&a,&sss,10);
  printf(" a=%s\n",&sss);//显示字符串sss
  //2.3  将此二进的大整数a转换为23进字符串sss并显示  
  printf("\n将此二进的a转换为23进制字符串sss并显示:\n");
  mp_toradix(&a,&sss,23);
  printf(" a=%s\n",&sss);//显示字符串sss
 
  printf("\n---------------------------------------------------------------------\n");
  printf("计算2的幂:\n");
  printf("\n请输入幂n=:\n");scanf("%d", &n);
  //计算a2=2**n
  mp_2expt(&a, n);
  printf("\n将此二进的a2转换为十进制字符串sss并显示:\n");
  mp_toradix(&a,&sss,10);
  printf(" a=%s\n",&sss);//显示字符串sss,问题：n(6万)大时这个显示有问题，请研究！
  ddraw(&a);

  printf("\n---------------------------------------------------------------------\n");
  printf("加法c=a+b: \n");
  
  //1  用fp_read_radix在程序中直接输入n进大整数sss,显示其二进结构
  printf("\n将10进数111111112222222333333444444405555554赋给二进的a1:\n");
  n=10;
  mp_read_radix(&a1, "111111112222222333333444444405555554", n); 
  //ddraw(&a3);
  printf("\n将11进数111111112222222333333444444405555555赋给二进的a2:\n");
  n=10;
  mp_read_radix(&a2, "111111112222222333333444444405555555", n); 
  printf("\n计算a3=a1+a2, a=a1-a2\n");
  mp_add(&a1,&a2,&a3);
  mp_sub(&a1,&a2,&a);

  printf("\n计算a3=a1+a2并显示结果\n");
  printf("\n将此二进的a3转换为十进制字符串sss并显示:\n");
  mp_toradix(&a3,&sss,10);
  printf(" a3=%s\n",&sss);//显示字符串sss,问题：n(6万)大时这个显示有问题，请研究！
  ddraw(&a3);
  
  printf("\n计算a=a1-a2并显示结果\n");
  printf("\n将此二进的a转换为十进制字符串sss1并显示:\n");
  mp_toradix(&a,&sss,10);
  printf(" a=%s\n",&sss);//显示字符串sss,问题：n(6万)大时这个显示有问题，请研究！
  ddraw(&a);
  
  printf("\n\n贵州大学计算机软件与理论研究所：\n");
  printf("请击一键退出：\n"); scanf("%ccc", &ccc);       
  //研究：怎么释放sss,str?


}
/////////////////////////////////////////////////////////////////////
//函数定义
//

////////////////////////////////////////////////////////////////////
//不调用其它函数
//

///////////////////////////////////////////////////////////////////////
//不调用其它函数
// 
/* init a new mp_int 
   为大整数结构a分配内存，缺省分配MP_PREC=64个mp_digit，
   因此，共有32*64比特用来表示大整数，而实际可处理28*64=1792比特
   
	 本函数未调用其它函数
*/
//int mp_init (mp_int * a);
/* init a new mp_int */
int mp_init (mp_int * a)
{
  int i;

  /* allocate memory required and clear it */
  //OPT_CAST(x)定义为(x *)；XMALLOC定义为malloc-分配一段连接内存
  a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * MP_PREC);
  //即：  a->dp = (mp_digit *) malloc (sizeof (mp_digit) * MP_PREC);
  //如果内存分配失败
  if (a->dp == NULL) {
    return MP_MEM;//=-2 ：out of mem
  }
  //初始为0
  /* set the digits to zero */
  for (i = 0; i < MP_PREC; i++) {//MP_PREC=64
      a->dp[i] = 0;
  }

  /* set the used to zero, allocated digits to the default precision and sign to positive */
  a->used  = 0;      //一个mp_digit也没用过
  a->alloc = MP_PREC;//=64
  a->sign  = MP_ZPOS;//=0，非负整数

  return MP_OKAY;//=0: 分配成功 
}

/* init an mp_init for a given size 
   
	 本函数未调用其它函数
*/
int mp_init_size (mp_int * a, int size)
{
  int x;

  /* pad size so there are always extra digits */
  size += (MP_PREC * 2) - (size % MP_PREC);	
  
  /* alloc mem */
  a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * size);
  if (a->dp == NULL) {
    return MP_MEM;
  }

  /* set the members */
  a->used  = 0;
  a->alloc = size;
  a->sign  = MP_ZPOS;

  /* zero the digits */
  for (x = 0; x < size; x++) {
      a->dp[x] = 0;
  }

  return MP_OKAY;
}

/** 释放大整数所占内存 clear one (frees)
    
	  本函数未调用其它函数
**/
void mp_clear (mp_int * a)
{
  int i;

  /* only do anything if a hasn't been freed previously */
  if (a->dp != NULL) {
    /* first zero the digits */
    for (i = 0; i < a->used; i++) {
        a->dp[i] = 0;
    }
    /* free ram */
    XFREE(a->dp);//#define XFREE    free
    /* reset members to make debugging easier */
    a->dp    = NULL;
    a->alloc = a->used = 0;
    a->sign  = MP_ZPOS;
  }
}

/* 修剪未使用的mp_digit  trim unused digits 
 * 调用本函数可剪除开头的那些被视为“使用-used”过的但为0的mp_digit
 * This is used to ensure that leading zero digits are
 * trimed and the leading "used" digit will be non-zero
 * Typically very fast.  Also fixes the sign if there
 * are no more leading digits

  本函数未调用其它函数

 */
void mp_clamp (mp_int * a)
{
  /* decrease used while the most significant digit is zero. */
  while (a->used > 0 && a->dp[a->used - 1] == 0) {
    --(a->used);
  }

  /* reset the sign flag if used == 0 */
  if (a->used == 0) {
    a->sign = MP_ZPOS;
  }
}

/* 将a设置为0: a=0
  本函数未调用其它函数
*/
void mp_zero (mp_int * a)
{
  a->sign = MP_ZPOS;
  a->used = 0;
  //
  memset (a->dp, 0, sizeof (mp_digit) * a->alloc);
}

/* set to a digit 
  调用：mp_zero