2021SC@SDUSC 山东大学软件学院软件工程应用与实践

1、mrflsh2.c结构

mrflsh2.c的总体结构如下，，主要实现了expon()、 fexp()、flog()、fpowf()几个在miracl开源库中比较重要的函数，这一次的博客就是读一下这函数的功能。

 2、源代码 

static int expon(_MIPD_ big w,int n)
{  /* generator for C.F. of e */ 
    if (n==0) return 2;
    if (n%3==2) return 2*(n/3)+2;
    else return 1;
}

expon()方法的主要作用就是根据n的输入值获得不同的输出值。如果n等于0的话返回2，n除以3的余数为2的话返回2*(n/3)+2，否则返回1.

void fexp(_MIPD_ flash x,flash y)
{ /* calculates y=exp(x) */
    int i,n,nsq,m,sqrn,op[5];
    BOOL minus,rem;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;
    if (size(x)==0)
    {
        convert(_MIPP_ 1,y);
        return;
    }
    copy(x,y);

    MR_IN(54)

    minus=FALSE;
    if (size(y)<0)
    {
        minus=TRUE;
        negify(y,y);
    }
    ftrunc(_MIPP_ y,y,mr_mip->w9);
    n=size(y);
    if (n==MR_TOOBIG)
    {
        mr_berror(_MIPP_ MR_ERR_FLASH_OVERFLOW);
        MR_OUT
        return;
    }
    if (n==0) convert(_MIPP_ 1,y);
    else
    {
        build(_MIPP_ y,expon);
        if (minus)
        { /* underflow to zero - bit of a bodge */
            rem=mr_mip->ERCON;
            mr_mip->ERCON=TRUE;
            fpower(_MIPP_ y,n,y);
            mr_mip->ERCON=rem;
            if (mr_mip->ERNUM)
            {
                mr_mip->ERNUM=0;
                zero(y);
                MR_OUT
                return;
            }
        }
        else fpower(_MIPP_ y,n,y);
    }
    if (size(mr_mip->w9)==0)
    {
        if (minus) frecip(_MIPP_ y,y);
        MR_OUT
        return;
    }
    sqrn=isqrt(mr_mip->lg2b*mr_mip->workprec,mr_mip->lg2b);
    nsq=0;
    copy(mr_mip->w9,mr_mip->w8);
    frecip(_MIPP_ mr_mip->w9,mr_mip->w9);
    ftrunc(_MIPP_ mr_mip->w9,mr_mip->w9,mr_mip->w9);
    m=logb2(_MIPP_ mr_mip->w9);
    if (m<sqrn)
    { /* scale fraction down */
        nsq=sqrn-m;
        expb2(_MIPP_ nsq,mr_mip->w9);
        fdiv(_MIPP_ mr_mip->w8,mr_mip->w9,mr_mip->w8);
    }
    zero(mr_mip->w10);
    op[0]=0x4B;  /* set up for x/(C+y) */
    op[1]=1;
    op[2]=0;
    for (m=sqrn;m>0;m--)
    { /* Unwind C.F. expansion for exp(x)-1 */
        if (m%2==0) op[4]=2,op[3]=1;
        else        op[4]=m,op[3]=(-1);
        flop(_MIPP_ mr_mip->w8,mr_mip->w10,op,mr_mip->w10);
    }
    op[0]=0x2C;  /* set up for (x+2).y */
    op[1]=op[3]=1;
    op[2]=2;
    op[4]=0;
    for (i=0;i<nsq;i++)
    { /* now square it back up again */
        flop(_MIPP_ mr_mip->w10,mr_mip->w10,op,mr_mip->w10);
    }
    op[2]=1;
    flop(_MIPP_ mr_mip->w10,y,op,y);
    if (minus) frecip(_MIPP_ y,y);
    MR_OUT
}

fexp()方法的主要作用就是获得以输入数据flash数据类型x为幂，以e为底数的指数并且赋值给输入数据y。首先调用size方法将x转化为int数据类型，观察x是否为0，是的话调用convert方法把y赋值为1的flash数据类型，否则的话调用copy方法那x赋值给y。然后调用size方法看int数据类型的y是否小于0，是的话调用方法获得y的相反数，调用ftrunc方法把y切分为flsah类型的mr_mip->w9和big数据类型的y，调用size方法将y转化为int数据类型赋值给n。如果没有发生一个下溢到零位的错误的话，就等于fpower方法直接获得y的n次方赋值给y。之后该函数内部使用了多种其他函数来计算e的x次方，该方法的时间复杂度为O(n2.5)，可以说再相对较少的运行时间内就可以得出答案。

void flog(_MIPD_ flash x,flash y)
{ /* calculate y=log(x) to base e */
    BOOL hack;
    int op[5];
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    copy(x,y);
    if (mr_mip->ERNUM) return;
    if (size(y)==1)
    {
        zero(y);
        return;
    }

    MR_IN(55)

    if (size(y)<=0)
    {
        mr_berror(_MIPP_ MR_ERR_NEG_LOG);
        MR_OUT
        return;
    }
    hack=FALSE;
    if (mr_lent(y)<=2)
    { /* for 'simple' values of y */
        hack=TRUE;
        build(_MIPP_ mr_mip->w11,expon);
        fdiv(_MIPP_ y,mr_mip->w11,y);
    }
    op[0]=0x68;
    op[1]=op[3]=1;
    op[2]=(-1);
    op[4]=0;
    mr_mip->workprec=mr_mip->stprec;
    dconv(_MIPP_ log(fdsize(_MIPP_ y)),mr_mip->w11);
    while (mr_mip->workprec!=mr_mip->nib)
    { /* Newtons iteration w11=w11+(y-exp(w11))/exp(w11) */
        if (mr_mip->workprec<mr_mip->nib) mr_mip->workprec*=2;
        if (mr_mip->workprec>=mr_mip->nib) mr_mip->workprec=mr_mip->nib;
        else if (mr_mip->workprec*2>mr_mip->nib) mr_mip->workprec=(mr_mip->nib+1)/2;
        fexp(_MIPP_ mr_mip->w11,mr_mip->w12);
        flop(_MIPP_ y,mr_mip->w12,op,mr_mip->w12);
        fadd(_MIPP_ mr_mip->w12,mr_mip->w11,mr_mip->w11);
    }
    copy(mr_mip->w11,y);
    if (hack) fincr(_MIPP_ y,1,1,y);
    MR_OUT
}

flog()方法的主要作用就是获得以输入数据flash数据类型x为真数，以e为底数的对数并且赋值给输入数据y。首先调用mr_lent方法获得y的len属性，如果小于2的话，把hack赋值为true，调用build方法创建mr_mip->w11，调用fdiv方法获得y/mr_mip->w11并且赋值给y。然后给数组op赋予不同的初值。再把mr_mip->stprec赋值给mr_mip->workprec，调用fdsize方法把ycongflash形式转化为double形式，再调用log方法获得以y为真数，以e为底数的对数，然后调用dconv方法ba上述结果转化为flash形式再赋值给mr_mip->w11。接下来循环使用牛顿法计算w11+(y-exp(w11))/exp(w11)的近似值并且赋值给w11，最后调用copy方法把w11赋值给y，如果hack调用true的话，调用fincr方法使得y+1/1并且赋值给y。

void fpowf(_MIPD_ flash x,flash y,flash z)
{ /* raise flash number to flash power *
   *     z=x^y  -> z=exp(y.log(x))     */
    int n;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(56)

    n=size(y);
    if (mr_abs(n)<MR_TOOBIG)
    { /* y is small int */
        fpower(_MIPP_ x,n,z);
        MR_OUT
        return;
    }
    frecip(_MIPP_ y,mr_mip->w13);
    n=size(mr_mip->w13);
    if (mr_abs(n)<MR_TOOBIG)
    { /* 1/y is small int */
        froot(_MIPP_ x,n,z);
        MR_OUT
        return;
    }
    copy(x,z);
    flog(_MIPP_ z,z);
    fdiv(_MIPP_ z,mr_mip->w13,z);
    fexp(_MIPP_ z,z);
    MR_OUT
}

fpowf()方法的主要作用就是获得以输入数据x为底数，以y为幂的指数，并且把结果赋值给z。首先调用size方法把y转化为int类型赋值给n，然后调用mr_abs方法获得n的绝对值，如果n足够小的话，调用fpower方法直接获得x的n次方赋值给z并且退出。然后调用 frecip方法获得1/y并且赋值给mr_mip->w13，再调用size方法把mr_mip->w13转化为int类型赋值给n，如果n足够小的话，调用fpower方法直接获得x的n次方赋值给z并且退出。y或者1/y都不够小的话调用copy方法把x赋值给z，调用flog方法获得z的自然对数并且赋值给z，调用fdiv方法获得 z/mr_mip->w13并且赋值给z，再调用 fexp方法获得e的z次方并且赋值给z。