mrcore.c

miracl-大数运算库,大家使用有什么问题请多多提意见

/*
 *   MIRACL Core module - contains initialisation code and general purpose 
 *   utilities
 *   mrcore.c
 *
 *   Copyright (c) 1988-2005 Shamus Software Ltd.
 */

#include "miracl.h"
#include <stdlib.h>
#include <string.h>

/*** Multi-Threaded Support ***/

#ifndef MR_GENERIC_MT

  #ifdef MR_WINDOWS_MT
    #include <windows.h>
    DWORD mr_key;   

    miracl *get_mip()
    {
        return (miracl *)TlsGetValue(mr_key); 
    }

    void mr_init_threading()
    {
        mr_key=TlsAlloc();
    }

    void mr_end_threading()
    {
        TlsFree(mr_key);
    }

  #endif

  #ifdef MR_UNIX_MT
    #include <pthread.h>
    pthread_key_t mr_key;

    miracl *get_mip()
    {
        return (miracl *)pthread_getspecific(mr_key); 
    }

    void mr_init_threading()
    {
        pthread_key_create(&mr_key,(void(*)(void *))NULL);
    }

    void mr_end_threading()
    {
        pthread_key_delete(mr_key);
    }
  #endif

  #ifndef MR_WINDOWS_MT
    #ifndef MR_UNIX_MT
      #ifdef MR_STATIC
          miracl mip;
          miracl *mr_mip=&mip;
      #else
          miracl *mr_mip=NULL;     /* MIRACL's one and only global variable */
      #endif
      miracl *get_mip()
      {
          return (miracl *)mr_mip; 
      }
    #endif
  #endif

#endif

/* See Advanced Windows by Jeffrey Richter, Chapter 12 for methods for
   creating different instances of this global for each executing thread 
   when using Windows '95/NT
*/

#ifdef MR_STATIC

static const int mr_small_primes[]=
{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,
107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,
223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,
337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,
457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,
593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,
719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,
857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,
997,0};

#endif

#ifndef MR_STRIPPED_DOWN
#ifndef MR_NO_STANDARD_IO

static char *names[] =
{"your program","innum","otnum","jack","normalise",
"multiply","divide","incr","decr","premult",
"subdiv","fdsize","egcd","cbase",
"cinnum","cotnum","nroot","power",
"powmod","bigdig","bigrand","nxprime","isprime",
"mirvar","mad","multi_inverse","putdig",
"add","subtract","mirsys","xgcd",
"fpack","dconv","mr_shift","mround","fmul",
"fdiv","fadd","fsub","fcomp","fconv",
"frecip","fpmul","fincr","","ftrunc",
"frand","sftbit","build","logb2","expint",
"fpower","froot","fpi","fexp","flog","fpowf",
"ftan","fatan","fsin","fasin","fcos","facos",
"ftanh","fatanh","fsinh","fasinh","fcosh",
"facosh","flop","gprime","powltr","fft_mult",
"crt_init","crt","otstr","instr","cotstr","cinstr","powmod2",
"prepare_monty","nres","redc","nres_modmult","nres_powmod",
"nres_moddiv","nres_powltr","divisible","remain",
"fmodulo","nres_modadd","nres_modsub","nres_negate",
"ecurve_init","ecurve_add","ecurve_mult",
"epoint_init","epoint_set","epoint_get","nres_powmod2",
"nres_sqroot","sqroot","nres_premult","ecurve_mult2",
"ecurve_sub","trial_division","nxsafeprime","nres_lucas","lucas",
"brick_init","pow_brick","set_user_function",
"nres_powmodn","powmodn","ecurve_multn",
"ebrick_init","mul_brick","epoint_norm","nres_multi_inverse","",
"nres_dotprod","epoint_negate","ecurve_multi_add",
"ecurve2_init","","epoint2_set","epoint2_norm","epoint2_get",
"epoint2_comp","ecurve2_add","epoint2_negate","ecurve2_sub",
"ecurve2_multi_add","ecurve2_mult","ecurve2_multn","ecurve2_mult2",
"ebrick2_init","mul2_brick","prepare_basis","strong_bigrand",
"bytes_to_big","big_to_bytes","set_io_buffer_size",
"epoint_getxyz","epoint_double_add","nres_double_inverse",
"double_inverse","epoint_x","hamming","expb2","bigbits"};
/* 0 - 150 (151 in all) */

#endif
#endif


#ifdef MR_NOASM

/* C only versions of muldiv/muldvd/muldvd2/muldvm */
/* Note that mr_large should be twice the size of mr_small */

mr_small muldiv(mr_small a,mr_small b,mr_small c,mr_small m,mr_small *rp)
{
    mr_small q;
    mr_large ldres,p=(mr_large)a*b+c;
    q=(mr_small)(MR_LROUND(p/m));
    *rp=(mr_small)(p-(mr_large)q*m);
    return q;
}

#ifdef MR_FP_ROUNDING

mr_small imuldiv(mr_small a,mr_small b,mr_small c,mr_small m,mr_large im,mr_small *rp)
{
    mr_small q;
    mr_large ldres,p=(mr_large)a*b+c;
    q=(mr_small)MR_LROUND(p*im);
    *rp=(mr_small)(p-(mr_large)q*m);
    return q;
}

#endif

#ifndef MR_NOFULLWIDTH

mr_small muldvm(mr_small a,mr_small c,mr_small m,mr_small *rp)
{
    mr_small q;
    union doubleword dble;
    dble.h[MR_BOT]=c;
    dble.h[MR_TOP]=a;

    q=(mr_small)(dble.d/m);
    *rp=(mr_small)(dble.d-(mr_large)q*m);
    return q;
}

mr_small muldvd(mr_small a,mr_small b,mr_small c,mr_small *rp)
{
    union doubleword dble;
    dble.d=(mr_large)a*b+c;

    *rp=dble.h[MR_BOT];
    return dble.h[MR_TOP];
}

void muldvd2(mr_small a,mr_small b,mr_small *c,mr_small *rp)
{
    union doubleword dble;
    dble.d=(mr_large)a*b+*c+*rp;
    *rp=dble.h[MR_BOT];
    *c=dble.h[MR_TOP];
}

#endif
#endif

#ifdef MR_NOFULLWIDTH

/* no FULLWIDTH working, so supply dummies */

/*

mr_small muldvd(mr_small a,mr_small b,mr_small c,mr_small *rp)
{
    return (mr_small)0;
}

mr_small muldvm(mr_small a,mr_small c,mr_small m,mr_small *rp)
{
    return (mr_small)0;
}

void muldvd2(mr_small a,mr_small b,mr_small *c,mr_small *rp)
{
}

*/

#endif

#ifndef MR_NO_STANDARD_IO

static void mputs(char *s)
{ /* output a string */
    int i=0;
    while (s[i]!=0) fputc((int)s[i++],stdout);
}

#endif

void mr_berror(_MIPD_ int nerr)
{  /*  Big number error routine  */
#ifndef MR_STRIPPED_DOWN
int i;
#endif

#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif

if (mr_mip->ERCON)
{
    mr_mip->ERNUM=nerr;
    return;
}
#ifndef MR_NO_STANDARD_IO

#ifndef MR_STRIPPED_DOWN
mputs("\nMIRACL error from routine ");
if (mr_mip->depth<MR_MAXDEPTH) mputs(names[mr_mip->trace[mr_mip->depth]]);
else                           mputs("???");
fputc('\n',stdout);

for (i=mr_mip->depth-1;i>=0;i--)
{
    mputs("              called from ");
    if (i<MR_MAXDEPTH) mputs(names[mr_mip->trace[i]]);
    else               mputs("???");
    fputc('\n',stdout);
}

switch (nerr)
{
case 1 :
mputs("Number base too big for representation\n");
break;
case 2 :
mputs("Division by zero attempted\n");
break;
case 3 : 
mputs("Overflow - Number too big\n");
break;
case 4 :
mputs("Internal result is negative\n");
break;
case 5 : 
mputs("Input format error\n");
break;
case 6 :
mputs("Illegal number base\n");
break;
case 7 : 
mputs("Illegal parameter usage\n");
break;
case 8 :
mputs("Out of space\n");
break;
case 9 :
mputs("Even root of a negative number\n");
break;
case 10:
mputs("Raising integer to negative power\n");
break;
case 11:
mputs("Attempt to take illegal root\n");
break;
case 12:
mputs("Integer operation attempted on Flash number\n");
break;
case 13:
mputs("Flash overflow\n");
break;
case 14:
mputs("Numbers too big\n");
break;
case 15:
mputs("Log of a non-positive number\n");
break;
case 16:
mputs("Flash to double conversion failure\n");
break;
case 17:
mputs("I/O buffer overflow\n");
break;
case 18:
mputs("MIRACL not initialised - no call to mirsys()\n");
break;
case 19:
mputs("Illegal modulus \n");
break;
case 20:
mputs("No modulus defined\n");
break;
case 21:
mputs("Exponent too big\n");
break;
case 22:
mputs("Unsupported Feature - check mirdef.h\n");
break;
case 23:
mputs("Specified double length type isn't double length\n");
break;
case 24:
mputs("Specified basis is NOT irreducible\n");
break;
case 25:
mputs("Unable to control Floating-point rounding\n");
break;
case 26:
mputs("Base must be binary (MR_ALWAYS_BINARY defined in mirdef.h ?)\n");
break;
case 27:
mputs("No irreducible basis defined\n");
break;
case 28:
mputs("Composite modulus\n");
break;
default:
mputs("Undefined error\n");
break;
}
exit(0);
#else
mputs("MIRACL error\n");
exit(0);
#endif

#endif
}

#ifndef MR_STRIPPED_DOWN

void mr_track(_MIPDO_ )
{ /* track course of program execution *
   * through the MIRACL routines       */

#ifndef MR_NO_STANDARD_IO

    int i;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    for (i=0;i<mr_mip->depth;i++) fputc('-',stdout);
    fputc('>',stdout);
    mputs(names[mr_mip->trace[mr_mip->depth]]);
    fputc('\n',stdout);
#endif
}

#endif

mr_small brand(_MIPDO_ )
{ /* Marsaglia & Zaman random number generator */
    int i,k;
    mr_unsign32 pdiff,t;
    mr_small r;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->lg2b>32)
    { /* underlying type is > 32 bits. Assume <= 64 bits */
        mr_mip->rndptr+=2;
        if (mr_mip->rndptr<NK-1)
        {
            r=(mr_small)mr_mip->ira[mr_mip->rndptr];
            r=mr_shiftbits(r,mr_mip->lg2b-32);
            r+=(mr_small)mr_mip->ira[mr_mip->rndptr+1];
            return r;
        }
    }
    else
    {
        mr_mip->rndptr++;
        if (mr_mip->rndptr<NK) return (mr_small)mr_mip->ira[mr_mip->rndptr];
    }
    mr_mip->rndptr=0;
    for (i=0,k=NK-NJ;i<NK;i++,k++)
    { /* calculate next NK values */
        if (k==NK) k=0;
        t=mr_mip->ira[k];
        pdiff=t - mr_mip->ira[i] - mr_mip->borrow;
        if (pdiff<t) mr_mip->borrow=0;
        if (pdiff>t) mr_mip->borrow=1;
        mr_mip->ira[i]=pdiff; 
    }
    if (mr_mip->lg2b>32)
    { /* double up */
        r=(mr_small)mr_mip->ira[0];
        r=mr_shiftbits(r,mr_mip->lg2b-32);
        r+=(mr_small)mr_mip->ira[1];
        return r;
    }
    else return (mr_small)(mr_mip->ira[0]);
}

void irand(_MIPD_ mr_unsign32 seed)
{ /* initialise random number system */
    int i,in;
    mr_unsign32 t,m=1L;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    mr_mip->borrow=0L;
    mr_mip->rndptr=0;
    mr_mip->ira[0]^=seed;
    for (i=1;i<NK;i++)
    { /* fill initialisation vector */
        in=(NV*i)%NK;
        mr_mip->ira[in]=m; 
        t=m;
        m=seed-m;
        seed=t;
    }
    for (i=0;i<1000;i++) brand(_MIPPO_ ); /* "warm-up" & stir the generator */
}

mr_small mr_shiftbits(mr_small x,int n)
{
#ifdef MR_FP
    int i;
    mr_small dres;
    if (n==0) return x;
    if (n>0)
    {
        for (i=0;i<n;i++) x=x+x;
        return x;
    }
    n=-n;
    for (i=0;i<n;i++) x=MR_DIV(x,2.0);
    return x;
#else
    if (n==0) return x;
    if (n>0) x<<=n;
    else x>>=n;
    return x;
#endif

}

mr_small mr_setbase(_MIPD_ mr_small nb)
{  /* set base. Pack as many digits as  *
    * possible into each computer word  */
    mr_small temp;
#ifdef MR_FP
    mr_small dres;
#endif
#ifndef MR_NOFULLWIDTH
    BOOL fits;
    int bits;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    fits=FALSE;
    bits=MIRACL;
    while (bits>1) 
    {
        bits/=2;
        temp=((mr_small)1<<bits);
        if (temp==nb)
        {
            fits=TRUE;
            break;
        }
        if (temp<nb || (bits%2)!=0) break;
    }
    if (fits)
    {
        mr_mip->apbase=nb;
        mr_mip->pack=MIRACL/bits;
        mr_mip->base=0;
        return 0;
    }
#endif
    mr_mip->apbase=nb;
    mr_mip->pack=1;
    mr_mip->base=nb;
    if (mr_mip->base==0) return 0;
    temp=MR_DIV(MAXBASE,nb);
    while (temp>=nb)
    {
        temp=MR_DIV(temp,nb);
        mr_mip->base*=nb;
        mr_mip->pack++;
    }
#ifdef MR_FP_ROUNDING
    mr_mip->inverse_base=mr_invert(mr_mip->base);
    return mr_mip->inverse_base;
#else
    return 0;
#endif
}

#ifdef MR_FLASH

BOOL fit(big x,big y,int f)
{ /* returns TRUE if x/y would fit flash format of length f */
    int n,d;
    n=(int)(x->len&(MR_OBITS));
    d=(int)(y->len&(MR_OBITS));
    if (n==1 && x->w[0]==1) n=0;
    if (d==1 && y->w[0]==1) d=0;
    if (n+d<=f) return TRUE;
    return FALSE;
}

#endif

int mr_lent(flash x)
{ /* return length of big or flash in words */
    mr_unsign32 lx;
    lx=(x->len&(MR_OBITS));
#ifdef MR_FLASH
    return (int)((lx&(MR_MSK))+((lx>>(MR_BTS))&(MR_MSK)));
#else
    return (int)lx;
#endif
}

void zero(flash x)
{ /* set big/flash number to zero */
    int i,n;
    mr_small *g;
    if (x==NULL) return;
#ifdef MR_FLASH
    n=mr_lent(x);
#else
    n=(x->len&MR_OBITS);
#endif
    g=x->w;

    for (i=0;i<n;i++)
        g[i]=0;

    x->len=0;

}

void convert(_MIPD_ int n ,big x)
{  /*  convert integer n to big number format  */
    int m;
    mr_unsign32 s;
#ifdef MR_FP
    mr_small dres;
#endif
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    zero(x);
    if (n==0) return;
    s=0;
    if (n<0)
    {
        s=MR_MSBIT;
        n=(-n);
    }
    m=0;
    if (mr_mip->base==0)
    {
#ifndef MR_NOFULLWIDTH
#if MR_IBITS > MIRACL
        while (n>0)
        {
            x->w[m++]=(mr_small)(n%((mr_small)1<<(MIRACL)));
            n/=((mr_small)1<<(MIRACL));
        }
#else
        x->w[m++]=(mr_small)n;
#endif
#endif
    }
    else while (n>0)
    {
        x->w[m++]=MR_REMAIN((mr_small)n,mr_mip->base);
        n/=mr_mip->base;
    }
    x->len=(m|s);
}

void uconvert(_MIPD_ unsigned int n ,big x)
{  /*  convert integer n to big number format  */
    int m;
#ifdef MR_FP
    mr_small dres;
#endif
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    zero(x);
    if (n==0) return;