15701338753

Committed 17 Jun 2025 07:49AM UTC coverage: 50.382% (-0.004%) from 50.386%

Build # 15701338753

Build Type

Pull #662

github

Committed by

web-flow

Commit Message

Merge f1f68c050 into 207386593

Pull Request Pull Request #662: Cleanup: change VOID into void

Run Details

178 of 245 new or added lines in 34 files covered. (72.65%)

2 existing lines in 1 file now uncovered.

41784 of 82935 relevant lines covered (50.38%)

2640008.85 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

29.1

/sources/evaluate.c

/*
          #[ includes :
*/

#include "form3.h"
#include <stdarg.h>
#include <gmp.h>
#include <mpfr.h>

#define EXTRAPRECISION 4

#define RND MPFR_RNDN

#define auxr1 (((mpfr_t *)(AT.auxr_))[0])
#define auxr2 (((mpfr_t *)(AT.auxr_))[1])
#define auxr3 (((mpfr_t *)(AT.auxr_))[2])
#define auxr4 (((mpfr_t *)(AT.auxr_))[3])
#define auxr5 (((mpfr_t *)(AT.auxr_))[4])

mpf_t ln2;

int PackFloat(WORD *,mpf_t);
int UnpackFloat(mpf_t,WORD *);
void RatToFloat(mpf_t, UWORD *, int);
void FormtoZ(mpz_t,UWORD *,WORD);

/*
          #] includes : 
          #[ mpfr_ :
                 #[ IntegerToFloatr :

                Converts a Form long integer to a mpfr_ float of default size.
                We assume that sizeof(unsigned long int) = 2*sizeof(UWORD).
*/

void IntegerToFloatr(mpfr_t result, UWORD *formlong, int longsize)
{
        mpz_t z;
        mpz_init(z);
        FormtoZ(z,formlong,longsize);
        mpfr_set_z(result,z,RND);
        mpz_clear(z);
}

/*
                 #] IntegerToFloatr : 
                 #[ RatToFloatr :

                Converts a Form rational to a gmp float of default size.
*/

void RatToFloatr(mpfr_t result, UWORD *formrat, int ratsize)
{
        GETIDENTITY
        int nnum, nden;
        UWORD *num, *den;
        int sgn = 0;
        if ( ratsize < 0 ) { ratsize = -ratsize; sgn = 1; }
        nnum = nden = (ratsize-1)/2;
        num = formrat; den = formrat+nnum;
        while ( nnum > 1 && num[nnum-1] == 0 ) { nnum--; }
        while ( nden > 1 && den[nden-1] == 0 ) { nden--; }
        IntegerToFloatr(auxr4,num,nnum);
        IntegerToFloatr(auxr5,den,nden);
        mpfr_div(result,auxr4,auxr5,RND);
        if ( sgn > 0 ) mpfr_neg(result,result,RND);
}

/*
                 #] RatToFloatr : 
                 #[ SetfFloatPrecision :

        The AC.DefaultPrecision is in bits.
        prec is in limbs.
        fprec is NOT in bits for mpfr_ which is slightly less than in mpf_
        We also set the auxr1,auxr2,auxr3,auxr4,auxr5 variables.
*/

void SetfFloatPrecision(LONG prec)
{
/*
        mpfr_prec_t fprec = prec * mp_bits_per_limb - EXTRAPRECISION;
*/
        mpfr_prec_t fprec = prec + 1;
        mpfr_set_default_prec(fprec);
#ifdef WITHPTHREADS
        int totnum = AM.totalnumberofthreads, id;
        mpfr_t *a;
  #ifdef WITHSORTBOTS
        totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads);
  #endif
        if ( AB[0]->T.auxr_ ) {
            for ( id = 0; id < totnum; id++ ) {
                        a = (mpfr_t *)AB[id]->T.auxr_;
                        mpfr_clears(a[0],a[1],a[2],a[3],a[4],(mpfr_ptr)0);
                        M_free(AB[id]->T.auxr_,"AB[id]->T.auxr_");
                        AB[id]->T.auxr_ = 0;
                }
        }
    for ( id = 0; id < totnum; id++ ) {
                AB[id]->T.auxr_ = (void *)Malloc1(sizeof(mpfr_t)*5,"AB[id]->T.auxr_");
                a = (mpfr_t *)AB[id]->T.auxr_;
/*
                We work here with a[0] etc because the aux1 etc contain B which
                in the current routine would be AB[0] only
*/
                mpfr_inits2(fprec,a[0],a[1],a[2],a[3],a[4],(mpfr_ptr)0);
        }
#else
        if ( AT.auxr_ ) {
                mpfr_clears(auxr1,auxr2,auxr3,auxr4,auxr5,(mpfr_ptr)0);
                M_free(AT.auxr_,"AT.auxr_");
                AT.auxr_ = 0;
        }
        AT.auxr_ = (void *)Malloc1(sizeof(mpfr_t)*5,"AT.auxr_");
        mpfr_inits2(fprec,auxr1,auxr2,auxr3,auxr4,auxr5,(mpfr_ptr)0);
#endif
}

/*
                 #] SetfFloatPrecision : 
                 #[ ClearfFloat :
*/

void ClearfFloat(void)
{
#ifdef WITHPTHREADS
        int totnum = AM.totalnumberofthreads, id;
        mpfr_t *a;
  #ifdef WITHSORTBOTS
        totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads);
  #endif
        if ( AB[0]->T.auxr_ ) {
            for ( id = 0; id < totnum; id++ ) {
                        a = (mpfr_t *)AB[id]->T.auxr_;
                        mpfr_clears(a[0],a[1],a[2],a[3],a[4],(mpfr_ptr)0);
                        M_free(AB[id]->T.auxr_,"AB[id]->T.auxr_");
                        AB[id]->T.auxr_ = 0;
                }
/*
                M_free(AB[0]->T.auxr_,"AB[0]->T.auxr_");
                AB[0]->T.auxr_ = 0;
*/
        }
#else
        if ( AT.auxr_ ) {
                mpfr_clears(auxr1,auxr2,auxr3,auxr4,auxr5,(mpfr_ptr)0);
                M_free(AT.auxr_,"AT.auxr_");
                AT.auxr_ = 0;
        }
#endif
/*
        if ( AS.delta_1 ) {
                mpf_clear(ln2);
                AS.delta_1 = 0;
        }
*/
}

/*
                 #] ClearfFloat : 
                 #[ GetFloatArgument :

                Convert an argument to an mpfr float if possible.
                Return value: 0: was converted successfully.
                              -1: could not convert.
                Note: arguments with more than one term should be evaluated first
                inside an Argument environment. If there is still more than one
                term remaining we get the code: "could not convert".
                par tells which argument we need. If it is negative it should also
                be the last argument.
*/

int GetFloatArgument(PHEAD mpfr_t f_out,WORD *fun,int par)
{
        WORD *term, *tn, *t, *tstop, first, ncoef, *arg, *argp, abspar;
        arg = fun+FUNHEAD;
        abspar = ABS(par);
        while ( arg < fun+fun[1] && abspar > 1 ) { abspar--; NEXTARG(arg) }
        if ( arg >= fun+fun[1] || abspar!= 1 ) return(-1);
        if ( par < 0 ) {
                argp = arg; NEXTARG(argp); if ( argp != fun+fun[1] ) return(-1);
        }
        if ( *arg < 0 ) {
                if ( *arg == -SNUMBER ) {
                        mpfr_set_si(f_out,(LONG)(arg[1]),RND);
                }
                else if ( *arg == -SYMBOL && ( arg[1] == PISYMBOL ) ) {
                        mpfr_const_pi(f_out,RND);
                }
                else if ( *arg == -INDEX && arg[1] >= 0 && arg[1] < AM.OffsetIndex ) {
                        mpfr_set_ui(f_out,(ULONG)(arg[1]),RND);
                }
                else { return(-1); }
                return(0);
        }
        else if ( arg[0] != arg[ARGHEAD]+ARGHEAD ) {        /* more than one term */
                return(-1);
        }
        term = arg+ARGHEAD;
        tn = term+*term;
        tstop = tn-ABS(tn[-1]);
        t = term+1;
        first = 1;
        while ( t <= tstop ) {
                if ( t == tstop ) { /* Fraction */
                        if ( first ) RatToFloatr(f_out,(UWORD *)tstop,tn[-1]);
                        else {
                                RatToFloatr(auxr4,(UWORD *)tstop,tn[-1]);
                                mpfr_mul(f_out,f_out,auxr4,RND);
                        }
                        return(0);
                }
                else if ( *t == FLOATFUN ) {
                        if ( t+t[1] != tstop ) return(-1);
                        if ( UnpackFloat(aux5,t) < 0 ) return(-1);
                        if ( first ) {
                                mpfr_set_f(f_out,aux5,RND);
                                first = 0;
                        }
                        else {
                                mpfr_set_f(auxr4,aux5,RND);
                                mpfr_mul(f_out,f_out,auxr4,RND);
                        }
                        first = 0;
                        if ( tn[-1] < 0 ) { /* change sign */
                                ncoef = -tn[-1];
                                mpfr_neg(f_out,f_out,RND);
                        }
                        else ncoef = tn[-1];
                        if ( ncoef == 3 && tn[-2] == 1 && tn[-3] == 1 ) return(0);
/*
                        If the argument was properly normalized we are not supposed to come here.
*/
                        MLOCK(ErrorMessageLock);
                        MesPrint("Unnormalized argument in GetFloatArgument: %a",*term,term);
                        MUNLOCK(ErrorMessageLock);
                        Terminate(-1);
                }
                else if ( t[1] == SYMBOL && t[2] == 4 && t[3] == PISYMBOL && t[4] == 1 ) {
                        if ( first ) {
                                mpfr_const_pi(f_out,RND);
                                first = 0;
                        }
                        else {
                                mpfr_const_pi(auxr5,RND);
                                mpfr_mul(f_out,f_out,auxr5,RND);
                        }
                }
                else {        /* This we do not / cannot do */
                        return(-1);
                }
                t += t[1];
        }
        return(0);
}

/*
                 #] GetFloatArgument : 
                 #[ GetPiArgument :

                Tests for sin,cos,tan whether the argument is a simple
                multiple of pi or can be reduced to such.
                Return value: -1: the answer is no.
                0-23: we have 0, 1/12*pi_,...,23/12*pi_
                With success most values can be worked out by simple means.
*/

int GetPiArgument(PHEAD WORD *arg)
{
        UWORD *coef, *co, *tco, twelve, *numer, *denom, *c, *d;
        WORD i, ii, i2, iflip, nc, nd, *t, *tn, *tstop;
        int rem;
/*
        One: determine whether there is a rational coefficient and a pi_:
*/
        if ( *arg == -SNUMBER && arg[1] == 0 ) return(0);
        if ( *arg == -SYMBOL && arg[1] == PISYMBOL ) return(12);
        if ( *arg < 0 ) return(-1);
        if ( arg[ARGHEAD] != *arg-ARGHEAD ) return(-1);
        t = arg+ARGHEAD+1;
        tn = arg+*arg; tstop = tn - ABS(tn[-1]);
        if ( *t != SYMBOL || t[1] != 4 || t[2] != PISYMBOL || t[3] != 1
                || t+t[1] != tstop ) return(-1);
/*
        The denominator must be a divisor of 12
        1: copy the coefficient
*/
        co = coef = NumberMalloc("GetPiArgument");
        tco = (UWORD *)tstop;
        i = tn[-1]; if ( i < 0 ) { i = -i; iflip = 1; } else iflip = 0;
        ii = (i-1)/2;
        twelve = 12;
        NCOPY(co,tco,i);
        co = coef;
        Mully(BHEAD co,&ii,&twelve,1);
/*
        Now the denominator should be 1
        i2 = i = (ii-1)/2;
*/
        i2 = i = ii;
        numer = co; denom = numer + i;
        if ( i > 1 ) {
                while ( i2 > 0 ) {
                        if ( denom[i2--] != 0 ) {
                                NumberFree(co,"GetPiArgument");
                                return(-1);
                        }
                }
        }
        if ( *denom != 1 ) {
                NumberFree(co,"GetPiArgument");
                return(-1);
        }
/*
        Now we need the numerator modulus 24.
*/
        if ( i == 1 ) {
                rem = *numer % 24;
        }
        else {
                c = NumberMalloc("GetPiArgument");
                d = NumberMalloc("GetPiArgument");
                twelve *= 2;
                DivLong(numer,i,&twelve,1,c,&nc,d,&nd);
                rem = *d % 24;
                NumberFree(d,"GetPiArgument");
                NumberFree(c,"GetPiArgument");
        }
        NumberFree(co,"GetPiArgument");
        if ( iflip && rem != 0 ) rem = 24-rem;
        return(rem);
}

/*
                 #] GetPiArgument : 
                 #[ EvaluateFun :

                What we need to do is:
                1: look for a function to be treated.
                2: make sure its argument is treatable.
                3: call the proper mpfr_ function.
                4: accumulate the result.
                For some functions we have to insert 'smart' shortcuts as is
                the case with sin_(pi_) or sin_(pi_/6) of sqrt(4/9) etc.
                Otherwise we may have to insert a value for pi_ first.
                There are several types of arguments:
                a: (short) integers.
                b: rationals.
                c: floats.
                We accumulate the result(s) in auxr2. The argument comes in aux5
                and a result in auxr3 which then gets multiplied into auxr2.
                In the end we have to combine auxr2 with whatever coefficient
                existed already.

                When the float system is started we need for aux only 3 variables
                per thread. auxr1-auxr3. This should be done separately.
                The main problem is the conversion of mpfr_t to float_ and/or mpf_t
                and back.
*/


int EvaluateFun(PHEAD WORD *term, WORD level, WORD *pars)
{
        WORD *t, *tstop, *tt, *tnext, *newterm, i;
        WORD *oldworkpointer = AT.WorkPointer, nsize, nsgn;
        int retval = 0, first = 1, pimul;

        tstop = term + *term; tstop -= ABS(tstop[-1]);
        if ( AT.WorkPointer < term+*term ) AT.WorkPointer = term + *term;
        t = term+1;
        mpfr_set_ui(auxr2,1L,RND);
        while ( t < tstop ) {
                if ( pars[2] == *t ) {        /* have to do this one if possible */
TestArgument:
/*
                        There must be a single argument, except for the AGM function
*/
                        tnext = t+t[1]; tt = t+FUNHEAD; NEXTARG(tt);
                        if( *t == SYMBOL) {
                                for ( WORD* ti = t+2; ti < t+t[1]; ti+=2 ) {
                                        if( ( *ti == PISYMBOL || *ti == EESYMBOL  || *ti == EMSYMBOL )
                                         && ( pars[2] == ALLFUNCTIONS || pars[3] == *ti ) ) {
                                                if ( *ti == PISYMBOL )
                                                        mpfr_const_pi(auxr3,RND);
                                                else if ( *ti == EESYMBOL ) {
                                                        mpfr_set_ui(auxr3,1,RND);
                                                        mpfr_exp(auxr3,auxr3,RND);
                                                }
                                                else if ( *ti == EMSYMBOL )
                                                        mpfr_const_euler(auxr3,RND);
                                                if ( ti[1] != 1 )
                                                        mpfr_pow_si(auxr3,auxr3,ti[1],RND);
                                                mpfr_mul(auxr2,auxr2,auxr3,RND);
                                                ti[1] = 0;
                                        }
                                }
                                first = 0;
                                goto nextfun;
                        }
                        if ( tt != tnext && *t != AGMFUNCTION ) goto nextfun;
                        if ( *t == SINFUNCTION ) {
                                pimul = GetPiArgument(BHEAD t+FUNHEAD);
                                if ( pimul >= 0 && pimul < 24 ) {
                                        if ( pimul == 0 || pimul == 12 ) goto getout;
                                        if ( pimul > 12 ) { pimul = 24-pimul; nsgn = -1; }
                                        else nsgn = 1;
                                        if ( pimul > 6 ) pimul = 12-pimul;
                                        switch ( pimul ) {
                                                case 1:
label1:
                                                        mpfr_sqrt_ui(auxr3,3L,RND);
                                                        mpfr_sub_ui(auxr3,auxr3,1L,RND);
                                                        mpfr_sqrt_ui(auxr1,2L,RND);
                                                        mpfr_div_ui(auxr1,auxr1,4L,RND);
                                                        mpfr_mul(auxr3,auxr3,auxr1,RND);
                                                        if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
                                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                                        break;
                                                case 2:
label2:
                                                        mpfr_div_ui(auxr2,auxr2,2L,RND);
                                                        if ( nsgn < 0 ) mpfr_neg(auxr2,auxr2,RND);
                                                        break;
                                                case 3:
label3:
                                                        mpfr_sqrt_ui(auxr3,2L,RND);
                                                        mpfr_div_ui(auxr3,auxr3,2L,RND);
                                                        if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
                                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                                        break;
                                                case 4:
label4:
                                                        mpfr_sqrt_ui(auxr3,3L,RND);
                                                        mpfr_div_ui(auxr3,auxr3,2L,RND);
                                                        if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
                                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                                        break;
                                                case 5:
label5:
                                                        mpfr_sqrt_ui(auxr3,3L,RND);
                                                        mpfr_add_ui(auxr3,auxr3,1L,RND);
                                                        mpfr_sqrt_ui(auxr1,2L,RND);
                                                        mpfr_div_ui(auxr1,auxr1,4L,RND);
                                                        mpfr_mul(auxr3,auxr3,auxr1,RND);
                                                        if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
                                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                                        break;
                                                case 6:
label6:
                                                        if ( nsgn < 0 ) mpfr_neg(auxr2,auxr2,RND);
                                                        break;
                                        }
                                        *t = 0; first = 0;
                                        goto nextfun;
                                }
                        }
                        else if ( *t == COSFUNCTION ) {
                                pimul = GetPiArgument(BHEAD t+FUNHEAD);
                                if ( pimul >= 0 && pimul < 24 ) {
                                        if ( pimul > 12 ) pimul = 24-12;
                                        if ( pimul > 6 ) { pimul = 12-pimul; nsgn = -1; }
                                        else nsgn = 1;
                                        if ( pimul == 6 ) goto getout;
                                        switch ( pimul ) {
                                                case 0: goto label6;
                                                case 1: goto label5;
                                                case 2: goto label4;
                                                case 3: goto label3;
                                                case 4: goto label2;
                                                case 5: goto label1;
                                        }
                                }
                        }
                        else if ( *t == TANFUNCTION ) {
                                pimul = GetPiArgument(BHEAD t+FUNHEAD);
                                if ( pimul >= 0 && pimul < 24 ) {
                                        if ( pimul == 6 || pimul == 18 ) goto nextfun;
                                        if ( pimul == 0 || pimul == 12 ) goto getout;
                                        if ( pimul > 12 ) { pimul = 24-pimul; nsgn = -1; }
                                        else nsgn = 1;
                                        if ( pimul > 6 ) { pimul = 12-pimul; nsgn = -nsgn; }
                                        switch ( pimul ) {
                                                case 1:
                                                        mpfr_sqrt_ui(auxr3,3L,RND);
                                                        mpfr_sub_ui(auxr3,auxr3,2L,RND);
                                                        if ( nsgn > 0 ) mpfr_neg(auxr3,auxr3,RND);
                                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                                        break;
                                                case 2:
                                                        mpfr_sqrt_ui(auxr3,3L,RND);
                                                        mpfr_div_ui(auxr3,auxr3,3L,RND);
                                                        if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
                                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                                        break;
                                                case 3:
                                                        break;
                                                case 4:
                                                        mpfr_sqrt_ui(auxr3,3L,RND);
                                                        if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
                                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                                        break;
                                                case 5:
                                                        mpfr_sqrt_ui(auxr3,3L,RND);
                                                        mpfr_add_ui(auxr3,auxr3,2L,RND);
                                                        if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);
                                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                                        break;
                                        }
                                        *t = 0;
                                        goto nextfun;
                                }
                        }

                        if ( *t == AGMFUNCTION ) {
                                if ( GetFloatArgument(BHEAD auxr1,t,1) < 0 ) goto nextfun;
                                if ( GetFloatArgument(BHEAD auxr3,t,-2) < 0 ) goto nextfun;
                        }
                        else if ( GetFloatArgument(BHEAD auxr1,t,-1) < 0 ) goto nextfun;
                        nsgn = mpfr_sgn(auxr1);

                        switch ( *t ) {
                                case SQRTFUNCTION:
                                        if ( nsgn < 0 ) goto nextfun;
                                        if ( nsgn == 0 ) goto getout;
                                        else mpfr_sqrt(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case LNFUNCTION:
                                        if ( nsgn <= 0 ) goto nextfun;
                                        else mpfr_log(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case LI2FUNCTION: /* should be between -1 and +1 */
                                        mpfr_abs(auxr3,auxr1,RND);
                                        if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;
                                        mpfr_li2(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case GAMMAFUN:
/*
                                        We cannot do this when the argument is a non-positive integer
*/
                                        if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] <= 0 ) goto nextfun;
                                        if ( t[FUNHEAD] == t[1]-FUNHEAD && ABS(t[t[1]-1]) ==
                                                t[FUNHEAD]-ARGHEAD-1 && t[t[1]-1] < 0 ) {
                                                nsize = (-t[t[1]-1]-1)/2;
                                                if ( t[t[1]-1-nsize] == 1 ) {
                                                        for ( i = 1; i < nsize; i++ ) {
                                                                if ( t[t[1]-1-nsize+i] != 0 ) break;
                                                        }
                                                        if ( i >= nsize ) goto nextfun;
                                                }
                                        }
                                        mpfr_gamma(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case AGMFUNCTION:
                                        mpfr_agm(auxr3,auxr1,auxr3,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case SINHFUNCTION:
                                        mpfr_sinh(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case COSHFUNCTION:
                                        mpfr_cosh(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case TANHFUNCTION:
                                        mpfr_tanh(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case ASINHFUNCTION:
                                        mpfr_asinh(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case ACOSHFUNCTION:
                                        mpfr_acosh(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case ATANHFUNCTION:
                                        if ( nsgn == 0 ) goto getout;
                                        mpfr_abs(auxr3,auxr1,RND);
                                        if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;
                                        mpfr_atanh(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case ASINFUNCTION:
                                        if ( nsgn == 0 ) goto getout;
                                        mpfr_abs(auxr3,auxr1,RND);
                                        if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;
                                        mpfr_asin(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case ACOSFUNCTION:
                                        if ( nsgn == 0 ) goto getout;
                                        mpfr_abs(auxr3,auxr1,RND);
                                        if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;
                                        mpfr_acos(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case ATANFUNCTION:
                                        mpfr_atan(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case SINFUNCTION:
                                        mpfr_sin(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case COSFUNCTION:
                                        mpfr_cos(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                case TANFUNCTION:
                                        mpfr_tan(auxr3,auxr1,RND);
                                        mpfr_mul(auxr2,auxr2,auxr3,RND);
                                break;
                                default:
                                        goto nextfun;
                                break;
                        }
                        first = 0;
                        *t = 0;
                        goto nextfun;
                }
                else if ( pars[2] == ALLFUNCTIONS ) {
/*
                        Now we have to test whether this is one of our functions
*/
                        if ( t[1] == FUNHEAD ) goto nextfun;
                        switch ( *t ) {
                                case SQRTFUNCTION:
                                case LNFUNCTION:
                                case LI2FUNCTION:
                                case LINFUNCTION:
                                case ASINFUNCTION:
                                case ACOSFUNCTION:
                                case ATANFUNCTION:
                                case SINHFUNCTION:
                                case COSHFUNCTION:
                                case TANHFUNCTION:
                                case ASINHFUNCTION:
                                case ACOSHFUNCTION:
                                case ATANHFUNCTION:
                                case SINFUNCTION:
                                case COSFUNCTION:
                                case TANFUNCTION:
                                case AGMFUNCTION:
                                case SYMBOL:
                                        goto TestArgument;
                                case MZV:
                                case EULER:
                                case MZVHALF:
                                        goto nextfun;
                                default:
                                        MLOCK(ErrorMessageLock);
                                        MesPrint("Function in evaluate statement not yet implemented.");
                                        MUNLOCK(ErrorMessageLock);
                                        break;
                        }
                }
                else goto nextfun;
nextfun:
                t += t[1];
        }
        if ( first == 1 ) return(Generator(BHEAD term,level));
        mpfr_get_f(aux4,auxr2,RND);
/*
        Step 3:
        Now the regular coefficient, if it is not 1/1.
        We have two cases: size +- 3, or bigger.
*/
        
        nsize = term[*term-1];
        if ( nsize < 0 ) { nsize = -nsize; nsgn = -1; }
        else nsgn = 1;
        if ( aux4->_mp_size < 0 ) {
                aux4->_mp_size = -aux4->_mp_size;
                nsgn = -nsgn;
        }

        if ( nsize == 3 ) {
                if ( tstop[0] != 1 ) {
                        mpf_mul_ui(aux4,aux4,(ULONG)((UWORD)tstop[0]));
                }
                if ( tstop[1] != 1 ) {
                        mpf_div_ui(aux4,aux4,(ULONG)((UWORD)tstop[1]));
                }
        }
        else {
                RatToFloat(aux5,(UWORD *)tstop,nsize);
                mpf_mul(aux4,aux4,aux5);
        }
/*
        Now we have to locate possible other float_ functions.
        Note possible incompatibilities between the mpf and mpfr formats.
*/
        t = term+1;
        while ( t < tstop ) {
                if ( *t == FLOATFUN ) {
                        UnpackFloat(aux5,t);
                        mpf_mul(aux4,aux4,aux5);
                }
                t += t[1];
        }
/*
        Now we should compose the new term in the WorkSpace.
*/
        t = term+1;
        newterm = AT.WorkPointer;
        tt = newterm+1;
        while ( t < tstop ) {
                if ( *t == 0 || *t == FLOATFUN ) t += t[1];
                else {
                        i = t[1]; NCOPY(tt,t,i);
                }
        }
        PackFloat(tt,aux4);
        tt += tt[1];
        *tt++ = 1; *tt++ = 1; *tt++ = 3*nsgn;
        *newterm = tt-newterm;
        AT.WorkPointer = tt;
        retval = Generator(BHEAD newterm,level);
getout:
        AT.WorkPointer = oldworkpointer;
        return(retval);
}

/*
                 #] EvaluateFun : 
          #] mpfr_ : 
*/

1	/*
2	#[ includes :
3	*/
4
5	#include "form3.h"
6	#include <stdarg.h>
7	#include <gmp.h>
8	#include <mpfr.h>
9
10	#define EXTRAPRECISION 4
11
12	#define RND MPFR_RNDN
13
14	#define auxr1 (((mpfr_t *)(AT.auxr_))[0])
15	#define auxr2 (((mpfr_t *)(AT.auxr_))[1])
16	#define auxr3 (((mpfr_t *)(AT.auxr_))[2])
17	#define auxr4 (((mpfr_t *)(AT.auxr_))[3])
18	#define auxr5 (((mpfr_t *)(AT.auxr_))[4])
19
20	mpf_t ln2;
21
22	int PackFloat(WORD *,mpf_t);
23	int UnpackFloat(mpf_t,WORD *);
24	void RatToFloat(mpf_t, UWORD *, int);
25	void FormtoZ(mpz_t,UWORD *,WORD);
26
27	/*
28	#] includes :
29	#[ mpfr_ :
30	#[ IntegerToFloatr :
31
32	Converts a Form long integer to a mpfr_ float of default size.
33	We assume that sizeof(unsigned long int) = 2*sizeof(UWORD).
34	*/
35
36	void IntegerToFloatr(mpfr_t result, UWORD *formlong, int longsize)	×
37	{
38	mpz_t z;	×
39	mpz_init(z);	×
40	FormtoZ(z,formlong,longsize);	×
41	mpfr_set_z(result,z,RND);	×
42	mpz_clear(z);	×
43	}	×
44
45	/*
46	#] IntegerToFloatr :
47	#[ RatToFloatr :
48
49	Converts a Form rational to a gmp float of default size.
50	*/
51
52	void RatToFloatr(mpfr_t result, UWORD *formrat, int ratsize)	×
53	{
54	GETIDENTITY
55	int nnum, nden;	×
56	UWORD num, den;	×
57	int sgn = 0;	×
58	if ( ratsize < 0 ) { ratsize = -ratsize; sgn = 1; }	×
59	nnum = nden = (ratsize-1)/2;	×
60	num = formrat; den = formrat+nnum;	×
61	while ( nnum > 1 && num[nnum-1] == 0 ) { nnum--; }	×
62	while ( nden > 1 && den[nden-1] == 0 ) { nden--; }	×
63	IntegerToFloatr(auxr4,num,nnum);	×
64	IntegerToFloatr(auxr5,den,nden);	×
65	mpfr_div(result,auxr4,auxr5,RND);	×
66	if ( sgn > 0 ) mpfr_neg(result,result,RND);	×
67	}	×
68
69	/*
70	#] RatToFloatr :
71	#[ SetfFloatPrecision :
72
73	The AC.DefaultPrecision is in bits.
74	prec is in limbs.
75	fprec is NOT in bits for mpfr_ which is slightly less than in mpf_
76	We also set the auxr1,auxr2,auxr3,auxr4,auxr5 variables.
77	*/
78
79	void SetfFloatPrecision(LONG prec)	24✔
80	{
81	/*
82	mpfr_prec_t fprec = prec * mp_bits_per_limb - EXTRAPRECISION;
83	*/
84	mpfr_prec_t fprec = prec + 1;	24✔
85	mpfr_set_default_prec(fprec);	24✔
86	#ifdef WITHPTHREADS
87	int totnum = AM.totalnumberofthreads, id;	16✔
88	mpfr_t *a;	16✔
89	#ifdef WITHSORTBOTS
90	totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads);	16✔
91	#endif
92	if ( AB[0]->T.auxr_ ) {	16✔
93	for ( id = 0; id < totnum; id++ ) {
94	a = (mpfr_t *)AB[id]->T.auxr_;
95	mpfr_clears(a[0],a[1],a[2],a[3],a[4],(mpfr_ptr)0);
96	M_free(AB[id]->T.auxr_,"AB[id]->T.auxr_");
97	AB[id]->T.auxr_ = 0;
98	}
99	}
100	for ( id = 0; id < totnum; id++ ) {	96✔
101	AB[id]->T.auxr_ = (void )Malloc1(sizeof(mpfr_t)5,"AB[id]->T.auxr_");	80✔
102	a = (mpfr_t *)AB[id]->T.auxr_;	80✔
103	/*
104	We work here with a[0] etc because the aux1 etc contain B which
105	in the current routine would be AB[0] only
106	*/
107	mpfr_inits2(fprec,a[0],a[1],a[2],a[3],a[4],(mpfr_ptr)0);	80✔
108	}
109	#else
110	if ( AT.auxr_ ) {	8✔
111	mpfr_clears(auxr1,auxr2,auxr3,auxr4,auxr5,(mpfr_ptr)0);
112	M_free(AT.auxr_,"AT.auxr_");
113	AT.auxr_ = 0;
114	}
115	AT.auxr_ = (void )Malloc1(sizeof(mpfr_t)5,"AT.auxr_");	8✔
116	mpfr_inits2(fprec,auxr1,auxr2,auxr3,auxr4,auxr5,(mpfr_ptr)0);	8✔
117	#endif
118	}	24✔
119
120	/*
121	#] SetfFloatPrecision :
122	#[ ClearfFloat :
123	*/
124
NEW 125	void ClearfFloat(void)	×
126	{
127	#ifdef WITHPTHREADS
128	int totnum = AM.totalnumberofthreads, id;
129	mpfr_t *a;
130	#ifdef WITHSORTBOTS
131	totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads);
132	#endif
133	if ( AB[0]->T.auxr_ ) {
134	for ( id = 0; id < totnum; id++ ) {
135	a = (mpfr_t *)AB[id]->T.auxr_;
136	mpfr_clears(a[0],a[1],a[2],a[3],a[4],(mpfr_ptr)0);
137	M_free(AB[id]->T.auxr_,"AB[id]->T.auxr_");
138	AB[id]->T.auxr_ = 0;
139	}
140	/*
141	M_free(AB[0]->T.auxr_,"AB[0]->T.auxr_");
142	AB[0]->T.auxr_ = 0;
143	*/
144	}
145	#else
146	if ( AT.auxr_ ) {
147	mpfr_clears(auxr1,auxr2,auxr3,auxr4,auxr5,(mpfr_ptr)0);
148	M_free(AT.auxr_,"AT.auxr_");
149	AT.auxr_ = 0;
150	}
151	#endif
152	/*
153	if ( AS.delta_1 ) {
154	mpf_clear(ln2);
155	AS.delta_1 = 0;
156	}
157	*/
158	}	×
159
160	/*
161	#] ClearfFloat :
162	#[ GetFloatArgument :
163
164	Convert an argument to an mpfr float if possible.
165	Return value: 0: was converted successfully.
166	-1: could not convert.
167	Note: arguments with more than one term should be evaluated first
168	inside an Argument environment. If there is still more than one
169	term remaining we get the code: "could not convert".
170	par tells which argument we need. If it is negative it should also
171	be the last argument.
172	*/
173
174	int GetFloatArgument(PHEAD mpfr_t f_out,WORD *fun,int par)	12✔
175	{
176	WORD term, tn, t, tstop, first, ncoef, arg, argp, abspar;	12✔
177	arg = fun+FUNHEAD;	12✔
178	abspar = ABS(par);	12✔
179	while ( arg < fun+fun[1] && abspar > 1 ) { abspar--; NEXTARG(arg) }	12✔
180	if ( arg >= fun+fun[1] \|\| abspar!= 1 ) return(-1);	12✔
181	if ( par < 0 ) {	12✔
182	argp = arg; NEXTARG(argp); if ( argp != fun+fun[1] ) return(-1);	12✔
183	}
184	if ( *arg < 0 ) {	12✔
185	if ( *arg == -SNUMBER ) {	12✔
186	mpfr_set_si(f_out,(LONG)(arg[1]),RND);	6✔
187	}
188	else if ( *arg == -SYMBOL && ( arg[1] == PISYMBOL ) ) {	6✔
189	mpfr_const_pi(f_out,RND);	6✔
190	}
191	else if ( *arg == -INDEX && arg[1] >= 0 && arg[1] < AM.OffsetIndex ) {	×
192	mpfr_set_ui(f_out,(ULONG)(arg[1]),RND);	×
193	}
194	else { return(-1); }
195	return(0);	12✔
196	}
197	else if ( arg[0] != arg[ARGHEAD]+ARGHEAD ) { /* more than one term */	×
198	return(-1);
199	}
200	term = arg+ARGHEAD;	×
201	tn = term+*term;	×
202	tstop = tn-ABS(tn[-1]);	×
203	t = term+1;	×
204	first = 1;	×
205	while ( t <= tstop ) {	×
206	if ( t == tstop ) { /* Fraction */	×
207	if ( first ) RatToFloatr(f_out,(UWORD *)tstop,tn[-1]);	×
208	else {
209	RatToFloatr(auxr4,(UWORD *)tstop,tn[-1]);	×
210	mpfr_mul(f_out,f_out,auxr4,RND);	×
211	}
212	return(0);	×
213	}
214	else if ( *t == FLOATFUN ) {	×
215	if ( t+t[1] != tstop ) return(-1);	×
216	if ( UnpackFloat(aux5,t) < 0 ) return(-1);	×
217	if ( first ) {	×
218	mpfr_set_f(f_out,aux5,RND);	×
219	first = 0;	×
220	}
221	else {
222	mpfr_set_f(auxr4,aux5,RND);	×
223	mpfr_mul(f_out,f_out,auxr4,RND);	×
224	}
225	first = 0;	×
226	if ( tn[-1] < 0 ) { /* change sign */	×
227	ncoef = -tn[-1];	×
228	mpfr_neg(f_out,f_out,RND);	×
229	}
230	else ncoef = tn[-1];
231	if ( ncoef == 3 && tn[-2] == 1 && tn[-3] == 1 ) return(0);	×
232	/*
233	If the argument was properly normalized we are not supposed to come here.
234	*/
235	MLOCK(ErrorMessageLock);	×
236	MesPrint("Unnormalized argument in GetFloatArgument: %a",*term,term);	×
237	MUNLOCK(ErrorMessageLock);	×
238	Terminate(-1);	×
239	}
240	else if ( t[1] == SYMBOL && t[2] == 4 && t[3] == PISYMBOL && t[4] == 1 ) {	×
241	if ( first ) {	×
242	mpfr_const_pi(f_out,RND);	×
243	first = 0;	×
244	}
245	else {
246	mpfr_const_pi(auxr5,RND);	×
247	mpfr_mul(f_out,f_out,auxr5,RND);	×
248	}
249	}
250	else { /* This we do not / cannot do */
251	return(-1);
252	}
253	t += t[1];	×
254	}
255	return(0);
256	}
257
258	/*
259	#] GetFloatArgument :
260	#[ GetPiArgument :
261
262	Tests for sin,cos,tan whether the argument is a simple
263	multiple of pi or can be reduced to such.
264	Return value: -1: the answer is no.
265	0-23: we have 0, 1/12pi_,...,23/12pi_
266	With success most values can be worked out by simple means.
267	*/
268
269	int GetPiArgument(PHEAD WORD *arg)	×
270	{
271	UWORD coef, co, tco, twelve, numer, denom, c, *d;	×
272	WORD i, ii, i2, iflip, nc, nd, t, tn, *tstop;	×
273	int rem;	×
274	/*
275	One: determine whether there is a rational coefficient and a pi_:
276	*/
277	if ( *arg == -SNUMBER && arg[1] == 0 ) return(0);	×
278	if ( *arg == -SYMBOL && arg[1] == PISYMBOL ) return(12);	×
279	if ( *arg < 0 ) return(-1);	×
280	if ( arg[ARGHEAD] != *arg-ARGHEAD ) return(-1);	×
281	t = arg+ARGHEAD+1;	×
282	tn = arg+*arg; tstop = tn - ABS(tn[-1]);	×
283	if ( *t != SYMBOL \|\| t[1] != 4 \|\| t[2] != PISYMBOL \|\| t[3] != 1	×
284	\|\| t+t[1] != tstop ) return(-1);	×
285	/*
286	The denominator must be a divisor of 12
287	1: copy the coefficient
288	*/
289	co = coef = NumberMalloc("GetPiArgument");	×
290	tco = (UWORD *)tstop;	×
291	i = tn[-1]; if ( i < 0 ) { i = -i; iflip = 1; } else iflip = 0;	×
292	ii = (i-1)/2;	×
293	twelve = 12;	×
294	NCOPY(co,tco,i);	×
295	co = coef;	×
296	Mully(BHEAD co,&ii,&twelve,1);	×
297	/*
298	Now the denominator should be 1
299	i2 = i = (ii-1)/2;
300	*/
301	i2 = i = ii;	×
302	numer = co; denom = numer + i;	×
303	if ( i > 1 ) {	×
304	while ( i2 > 0 ) {	×
305	if ( denom[i2--] != 0 ) {	×
306	NumberFree(co,"GetPiArgument");	×
307	return(-1);	×
308	}
309	}
310	}
311	if ( *denom != 1 ) {	×
312	NumberFree(co,"GetPiArgument");	×
313	return(-1);	×
314	}
315	/*
316	Now we need the numerator modulus 24.
317	*/
318	if ( i == 1 ) {	×
319	rem = *numer % 24;	×
320	}
321	else {
322	c = NumberMalloc("GetPiArgument");	×
323	d = NumberMalloc("GetPiArgument");	×
324	twelve *= 2;	×
325	DivLong(numer,i,&twelve,1,c,&nc,d,&nd);	×
326	rem = *d % 24;	×
327	NumberFree(d,"GetPiArgument");	×
328	NumberFree(c,"GetPiArgument");	×
329	}
330	NumberFree(co,"GetPiArgument");	×
331	if ( iflip && rem != 0 ) rem = 24-rem;	×
332	return(rem);
333	}
334
335	/*
336	#] GetPiArgument :
337	#[ EvaluateFun :
338
339	What we need to do is:
340	1: look for a function to be treated.
341	2: make sure its argument is treatable.
342	3: call the proper mpfr_ function.
343	4: accumulate the result.
344	For some functions we have to insert 'smart' shortcuts as is
345	the case with sin_(pi_) or sin_(pi_/6) of sqrt(4/9) etc.
346	Otherwise we may have to insert a value for pi_ first.
347	There are several types of arguments:
348	a: (short) integers.
349	b: rationals.
350	c: floats.
351	We accumulate the result(s) in auxr2. The argument comes in aux5
352	and a result in auxr3 which then gets multiplied into auxr2.
353	In the end we have to combine auxr2 with whatever coefficient
354	existed already.
355
356	When the float system is started we need for aux only 3 variables
357	per thread. auxr1-auxr3. This should be done separately.
358	The main problem is the conversion of mpfr_t to float_ and/or mpf_t
359	and back.
360	*/
361
362
363	int EvaluateFun(PHEAD WORD term, WORD level, WORD pars)	114✔
364	{
365	WORD t, tstop, tt, tnext, *newterm, i;	114✔
366	WORD *oldworkpointer = AT.WorkPointer, nsize, nsgn;	114✔
367	int retval = 0, first = 1, pimul;	114✔
368
369	tstop = term + *term; tstop -= ABS(tstop[-1]);	114✔
370	if ( AT.WorkPointer < term+term ) AT.WorkPointer = term + term;	114✔
371	t = term+1;	114✔
372	mpfr_set_ui(auxr2,1L,RND);	114✔
373	while ( t < tstop ) {	348✔
374	if ( pars[2] == t ) { / have to do this one if possible */	234✔
375	TestArgument:	54✔
376	/*
377	There must be a single argument, except for the AGM function
378	*/
379	tnext = t+t[1]; tt = t+FUNHEAD; NEXTARG(tt);	120✔
380	if( *t == SYMBOL) {	120✔
381	for ( WORD* ti = t+2; ti < t+t[1]; ti+=2 ) {	240✔
382	if( ( ti == PISYMBOL \|\| ti == EESYMBOL \|\| *ti == EMSYMBOL )	132✔
383	&& ( pars[2] == ALLFUNCTIONS \|\| pars[3] == *ti ) ) {	108✔
384	if ( *ti == PISYMBOL )	72✔
385	mpfr_const_pi(auxr3,RND);	42✔
386	else if ( *ti == EESYMBOL ) {	30✔
387	mpfr_set_ui(auxr3,1,RND);	18✔
388	mpfr_exp(auxr3,auxr3,RND);	18✔
389	}
390	else if ( *ti == EMSYMBOL )	12✔
391	mpfr_const_euler(auxr3,RND);	12✔
392	if ( ti[1] != 1 )	72✔
393	mpfr_pow_si(auxr3,auxr3,ti[1],RND);	6✔
394	mpfr_mul(auxr2,auxr2,auxr3,RND);	72✔
395	ti[1] = 0;	72✔
396	}
397	}
398	first = 0;	108✔
399	goto nextfun;	108✔
400	}
401	if ( tt != tnext && *t != AGMFUNCTION ) goto nextfun;	12✔
402	if ( *t == SINFUNCTION ) {	12✔
403	pimul = GetPiArgument(BHEAD t+FUNHEAD);	×
404	if ( pimul >= 0 && pimul < 24 ) {	×
405	if ( pimul == 0 \|\| pimul == 12 ) goto getout;	×
406	if ( pimul > 12 ) { pimul = 24-pimul; nsgn = -1; }	×
407	else nsgn = 1;
408	if ( pimul > 6 ) pimul = 12-pimul;	×
409	switch ( pimul ) {	×
410	case 1:
411	label1:	×
412	mpfr_sqrt_ui(auxr3,3L,RND);	×
413	mpfr_sub_ui(auxr3,auxr3,1L,RND);	×
414	mpfr_sqrt_ui(auxr1,2L,RND);	×
415	mpfr_div_ui(auxr1,auxr1,4L,RND);	×
416	mpfr_mul(auxr3,auxr3,auxr1,RND);	×
417	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
418	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
419	break;	×
420	case 2:
421	label2:	×
422	mpfr_div_ui(auxr2,auxr2,2L,RND);	×
423	if ( nsgn < 0 ) mpfr_neg(auxr2,auxr2,RND);	×
424	break;
425	case 3:
426	label3:	×
427	mpfr_sqrt_ui(auxr3,2L,RND);	×
428	mpfr_div_ui(auxr3,auxr3,2L,RND);	×
429	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
430	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
431	break;	×
432	case 4:
433	label4:	×
434	mpfr_sqrt_ui(auxr3,3L,RND);	×
435	mpfr_div_ui(auxr3,auxr3,2L,RND);	×
436	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
437	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
438	break;	×
439	case 5:
440	label5:	×
441	mpfr_sqrt_ui(auxr3,3L,RND);	×
442	mpfr_add_ui(auxr3,auxr3,1L,RND);	×
443	mpfr_sqrt_ui(auxr1,2L,RND);	×
444	mpfr_div_ui(auxr1,auxr1,4L,RND);	×
445	mpfr_mul(auxr3,auxr3,auxr1,RND);	×
446	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
447	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
448	break;	×
449	case 6:
450	label6:	×
451	if ( nsgn < 0 ) mpfr_neg(auxr2,auxr2,RND);	×
452	break;
453	}
454	*t = 0; first = 0;	×
455	goto nextfun;	×
456	}
457	}
458	else if ( *t == COSFUNCTION ) {	12✔
459	pimul = GetPiArgument(BHEAD t+FUNHEAD);	×
460	if ( pimul >= 0 && pimul < 24 ) {	×
461	if ( pimul > 12 ) pimul = 24-12;	×
462	if ( pimul > 6 ) { pimul = 12-pimul; nsgn = -1; }	×
463	else nsgn = 1;
464	if ( pimul == 6 ) goto getout;	×
465	switch ( pimul ) {	×
466	case 0: goto label6;	×
467	case 1: goto label5;	×
468	case 2: goto label4;	×
469	case 3: goto label3;	×
470	case 4: goto label2;	×
471	case 5: goto label1;	×
472	}
473	}
474	}
475	else if ( *t == TANFUNCTION ) {	12✔
476	pimul = GetPiArgument(BHEAD t+FUNHEAD);	×
477	if ( pimul >= 0 && pimul < 24 ) {	×
478	if ( pimul == 6 \|\| pimul == 18 ) goto nextfun;	×
479	if ( pimul == 0 \|\| pimul == 12 ) goto getout;	×
480	if ( pimul > 12 ) { pimul = 24-pimul; nsgn = -1; }	×
481	else nsgn = 1;
482	if ( pimul > 6 ) { pimul = 12-pimul; nsgn = -nsgn; }	×
483	switch ( pimul ) {	×
484	case 1:	×
485	mpfr_sqrt_ui(auxr3,3L,RND);	×
486	mpfr_sub_ui(auxr3,auxr3,2L,RND);	×
487	if ( nsgn > 0 ) mpfr_neg(auxr3,auxr3,RND);	×
488	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
489	break;	×
490	case 2:	×
491	mpfr_sqrt_ui(auxr3,3L,RND);	×
492	mpfr_div_ui(auxr3,auxr3,3L,RND);	×
493	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
494	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
495	break;	×
496	case 3:
497	break;
498	case 4:	×
499	mpfr_sqrt_ui(auxr3,3L,RND);	×
500	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
501	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
502	break;	×
503	case 5:	×
504	mpfr_sqrt_ui(auxr3,3L,RND);	×
505	mpfr_add_ui(auxr3,auxr3,2L,RND);	×
506	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
507	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
508	break;	×
509	}
510	*t = 0;	×
511	goto nextfun;	×
512	}
513	}
514
515	if ( *t == AGMFUNCTION ) {	12✔
516	if ( GetFloatArgument(BHEAD auxr1,t,1) < 0 ) goto nextfun;	×
517	if ( GetFloatArgument(BHEAD auxr3,t,-2) < 0 ) goto nextfun;	×
518	}
519	else if ( GetFloatArgument(BHEAD auxr1,t,-1) < 0 ) goto nextfun;	12✔
520	nsgn = mpfr_sgn(auxr1);	12✔
521
522	switch ( *t ) {	12✔
523	case SQRTFUNCTION:	12✔
524	if ( nsgn < 0 ) goto nextfun;	12✔
525	if ( nsgn == 0 ) goto getout;	12✔
526	else mpfr_sqrt(auxr3,auxr1,RND);	12✔
527	mpfr_mul(auxr2,auxr2,auxr3,RND);	12✔
528	break;	12✔
529	case LNFUNCTION:	×
530	if ( nsgn <= 0 ) goto nextfun;	×
531	else mpfr_log(auxr3,auxr1,RND);	×
532	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
533	break;	×
534	case LI2FUNCTION: /* should be between -1 and +1 */	×
535	mpfr_abs(auxr3,auxr1,RND);	×
536	if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;	×
537	mpfr_li2(auxr3,auxr1,RND);	×
538	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
539	break;	×
540	case GAMMAFUN:	×
541	/*
542	We cannot do this when the argument is a non-positive integer
543	*/
544	if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] <= 0 ) goto nextfun;	×
545	if ( t[FUNHEAD] == t[1]-FUNHEAD && ABS(t[t[1]-1]) ==	×
546	t[FUNHEAD]-ARGHEAD-1 && t[t[1]-1] < 0 ) {	×
547	nsize = (-t[t[1]-1]-1)/2;	×
548	if ( t[t[1]-1-nsize] == 1 ) {	×
549	for ( i = 1; i < nsize; i++ ) {	×
550	if ( t[t[1]-1-nsize+i] != 0 ) break;	×
551	}
552	if ( i >= nsize ) goto nextfun;	×
553	}
554	}
555	mpfr_gamma(auxr3,auxr1,RND);	×
556	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
557	break;	×
558	case AGMFUNCTION:	×
559	mpfr_agm(auxr3,auxr1,auxr3,RND);	×
560	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
561	break;	×
562	case SINHFUNCTION:	×
563	mpfr_sinh(auxr3,auxr1,RND);	×
564	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
565	break;	×
566	case COSHFUNCTION:	×
567	mpfr_cosh(auxr3,auxr1,RND);	×
568	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
569	break;	×
570	case TANHFUNCTION:	×
571	mpfr_tanh(auxr3,auxr1,RND);	×
572	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
573	break;	×
574	case ASINHFUNCTION:	×
575	mpfr_asinh(auxr3,auxr1,RND);	×
576	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
577	break;	×
578	case ACOSHFUNCTION:	×
579	mpfr_acosh(auxr3,auxr1,RND);	×
580	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
581	break;	×
582	case ATANHFUNCTION:	×
583	if ( nsgn == 0 ) goto getout;	×
584	mpfr_abs(auxr3,auxr1,RND);	×
585	if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;	×
586	mpfr_atanh(auxr3,auxr1,RND);	×
587	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
588	break;	×
589	case ASINFUNCTION:	×
590	if ( nsgn == 0 ) goto getout;	×
591	mpfr_abs(auxr3,auxr1,RND);	×
592	if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;	×
593	mpfr_asin(auxr3,auxr1,RND);	×
594	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
595	break;	×
596	case ACOSFUNCTION:	×
597	if ( nsgn == 0 ) goto getout;	×
598	mpfr_abs(auxr3,auxr1,RND);	×
599	if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;	×
600	mpfr_acos(auxr3,auxr1,RND);	×
601	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
602	break;	×
603	case ATANFUNCTION:	×
604	mpfr_atan(auxr3,auxr1,RND);	×
605	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
606	break;	×
607	case SINFUNCTION:	×
608	mpfr_sin(auxr3,auxr1,RND);	×
609	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
610	break;	×
611	case COSFUNCTION:	×
612	mpfr_cos(auxr3,auxr1,RND);	×
613	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
614	break;	×
615	case TANFUNCTION:	×
616	mpfr_tan(auxr3,auxr1,RND);	×
617	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
618	break;	×
619	default:	×
620	goto nextfun;	×
621	break;	12✔
622	}
623	first = 0;	12✔
624	*t = 0;	12✔
625	goto nextfun;	12✔
626	}
627	else if ( pars[2] == ALLFUNCTIONS ) {	180✔
628	/*
629	Now we have to test whether this is one of our functions
630	*/
631	if ( t[1] == FUNHEAD ) goto nextfun;	126✔
632	switch ( *t ) {	126✔
633	case SQRTFUNCTION:	66✔
634	case LNFUNCTION:
635	case LI2FUNCTION:
636	case LINFUNCTION:
637	case ASINFUNCTION:
638	case ACOSFUNCTION:
639	case ATANFUNCTION:
640	case SINHFUNCTION:
641	case COSHFUNCTION:
642	case TANHFUNCTION:
643	case ASINHFUNCTION:
644	case ACOSHFUNCTION:
645	case ATANHFUNCTION:
646	case SINFUNCTION:
647	case COSFUNCTION:
648	case TANFUNCTION:
649	case AGMFUNCTION:
650	case SYMBOL:
651	goto TestArgument;	66✔
652	case MZV:	×
653	case EULER:
654	case MZVHALF:
655	goto nextfun;	×
656	default:	60✔
657	MLOCK(ErrorMessageLock);	60✔
658	MesPrint("Function in evaluate statement not yet implemented.");	60✔
659	MUNLOCK(ErrorMessageLock);	60✔
660	break;	60✔
661	}
662	}
663	else goto nextfun;	54✔
664	nextfun:	234✔
665	t += t[1];	234✔
666	}
667	if ( first == 1 ) return(Generator(BHEAD term,level));	114✔
668	mpfr_get_f(aux4,auxr2,RND);	114✔
669	/*
670	Step 3:
671	Now the regular coefficient, if it is not 1/1.
672	We have two cases: size +- 3, or bigger.
673	*/
674
675	nsize = term[*term-1];	114✔
676	if ( nsize < 0 ) { nsize = -nsize; nsgn = -1; }	114✔
677	else nsgn = 1;
678	if ( aux4->_mp_size < 0 ) {	114✔
679	aux4->_mp_size = -aux4->_mp_size;	×
680	nsgn = -nsgn;	×
681	}
682
683	if ( nsize == 3 ) {	114✔
684	if ( tstop[0] != 1 ) {	114✔
685	mpf_mul_ui(aux4,aux4,(ULONG)((UWORD)tstop[0]));	×
686	}
687	if ( tstop[1] != 1 ) {	114✔
688	mpf_div_ui(aux4,aux4,(ULONG)((UWORD)tstop[1]));	×
689	}
690	}
691	else {
692	RatToFloat(aux5,(UWORD *)tstop,nsize);	×
693	mpf_mul(aux4,aux4,aux5);	×
694	}
695	/*
696	Now we have to locate possible other float_ functions.
697	Note possible incompatibilities between the mpf and mpfr formats.
698	*/
699	t = term+1;
700	while ( t < tstop ) {	348✔
701	if ( *t == FLOATFUN ) {	234✔
702	UnpackFloat(aux5,t);	114✔
703	mpf_mul(aux4,aux4,aux5);	114✔
704	}
705	t += t[1];	234✔
706	}
707	/*
708	Now we should compose the new term in the WorkSpace.
709	*/
710	t = term+1;	114✔
711	newterm = AT.WorkPointer;	114✔
712	tt = newterm+1;	114✔
713	while ( t < tstop ) {	114✔
714	if ( t == 0 \|\| t == FLOATFUN ) t += t[1];	234✔
715	else {
716	i = t[1]; NCOPY(tt,t,i);	936✔
717	}
718	}
719	PackFloat(tt,aux4);	114✔
720	tt += tt[1];	114✔
721	tt++ = 1; tt++ = 1; tt++ = 3nsgn;	114✔
722	*newterm = tt-newterm;	114✔
723	AT.WorkPointer = tt;	114✔
724	retval = Generator(BHEAD newterm,level);	114✔
725	getout:	114✔
726	AT.WorkPointer = oldworkpointer;	114✔
727	return(retval);	114✔
728	}
729
730	/*
731	#] EvaluateFun :
732	#] mpfr_ :
733	*/

form-dev / form / 15701338753

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous