9364948935

Committed 04 Jun 2024 09:49AM UTC coverage: 49.979% (-0.02%) from 49.999%

Build # 9364948935

Build Type

Pull #526

github

Committed by

web-flow

Commit Message

Merge 7062bd769 into 83e3d4185

Pull Request Pull Request #526: RFC: better debugging

Run Details

52 of 415 new or added lines in 46 files covered. (12.53%)

32 existing lines in 2 files now uncovered.

41391 of 82816 relevant lines covered (49.98%)

878690.77 hits per line

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

0.0

/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] == AM.numpi ) ) {
                        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] == AM.numpi && 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] == AM.numpi ) 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] != AM.numpi || 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 ( 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:
                                        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)	×
80	{
81	/*
82	mpfr_prec_t fprec = prec * mp_bits_per_limb - EXTRAPRECISION;
83	*/
84	mpfr_prec_t fprec = prec + 1;	×
85	mpfr_set_default_prec(fprec);	×
86	#ifdef WITHPTHREADS
87	int totnum = AM.totalnumberofthreads, id;
88	mpfr_t *a;
89	#ifdef WITHSORTBOTS
90	totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads);
91	#endif
92	if ( AB[0]->T.auxr_ ) {
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++ ) {
101	AB[id]->T.auxr_ = (void )Malloc1(sizeof(mpfr_t)5,"AB[id]->T.auxr_");
102	a = (mpfr_t *)AB[id]->T.auxr_;
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);
108	}
109	#else
110	if ( AT.auxr_ ) {
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_");
116	mpfr_inits2(fprec,auxr1,auxr2,auxr3,auxr4,auxr5,(mpfr_ptr)0);
117	#endif
118	}	×
119
120	/*
121	#] SetfFloatPrecision :
122	#[ ClearfFloat :
123	*/
124
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)	×
175	{
176	WORD term, tn, t, tstop, first, ncoef, arg, argp, abspar;	×
177	arg = fun+FUNHEAD;	×
178	abspar = ABS(par);	×
179	while ( arg < fun+fun[1] && abspar > 1 ) { abspar--; NEXTARG(arg) }	×
180	if ( arg >= fun+fun[1] \|\| abspar!= 1 ) return(-1);	×
181	if ( par < 0 ) {	×
182	argp = arg; NEXTARG(argp); if ( argp != fun+fun[1] ) return(-1);	×
183	}
184	if ( *arg < 0 ) {	×
185	if ( *arg == -SNUMBER ) {	×
186	mpfr_set_si(f_out,(LONG)(arg[1]),RND);	×
187	}
188	else if ( *arg == -SYMBOL && ( arg[1] == AM.numpi ) ) {	×
189	mpfr_const_pi(f_out,RND);	×
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);	×
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);	×
NEW 238	TERMINATE(-1);	×
239	}
240	else if ( t[1] == SYMBOL && t[2] == 4 && t[3] == AM.numpi && 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] == AM.numpi ) 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] != AM.numpi \|\| 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)	×
364	{
365	WORD t, tstop, tt, tnext, *newterm, i;	×
366	WORD *oldworkpointer = AT.WorkPointer, nsize, nsgn;	×
367	int retval = 0, first = 1, pimul;	×
368
369	tstop = term + *term; tstop -= ABS(tstop[-1]);	×
370	if ( AT.WorkPointer < term+term ) AT.WorkPointer = term + term;	×
371	t = term+1;	×
372	mpfr_set_ui(auxr2,1L,RND);	×
373	while ( t < tstop ) {	×
374	if ( pars[2] == t ) { / have to do this one if possible */	×
375	TestArgument:	×
376	/*
377	There must be a single argument, except for the AGM function
378	*/
379	tnext = t+t[1]; tt = t+FUNHEAD; NEXTARG(tt);	×
380	if ( tt != tnext && *t != AGMFUNCTION ) goto nextfun;	×
381	if ( *t == SINFUNCTION ) {	×
382	pimul = GetPiArgument(BHEAD t+FUNHEAD);	×
383	if ( pimul >= 0 && pimul < 24 ) {	×
384	if ( pimul == 0 \|\| pimul == 12 ) goto getout;	×
385	if ( pimul > 12 ) { pimul = 24-pimul; nsgn = -1; }	×
386	else nsgn = 1;
387	if ( pimul > 6 ) pimul = 12-pimul;	×
388	switch ( pimul ) {	×
389	case 1:
390	label1:	×
391	mpfr_sqrt_ui(auxr3,3L,RND);	×
392	mpfr_sub_ui(auxr3,auxr3,1L,RND);	×
393	mpfr_sqrt_ui(auxr1,2L,RND);	×
394	mpfr_div_ui(auxr1,auxr1,4L,RND);	×
395	mpfr_mul(auxr3,auxr3,auxr1,RND);	×
396	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
397	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
398	break;	×
399	case 2:
400	label2:	×
401	mpfr_div_ui(auxr2,auxr2,2L,RND);	×
402	if ( nsgn < 0 ) mpfr_neg(auxr2,auxr2,RND);	×
403	break;
404	case 3:
405	label3:	×
406	mpfr_sqrt_ui(auxr3,2L,RND);	×
407	mpfr_div_ui(auxr3,auxr3,2L,RND);	×
408	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
409	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
410	break;	×
411	case 4:
412	label4:	×
413	mpfr_sqrt_ui(auxr3,3L,RND);	×
414	mpfr_div_ui(auxr3,auxr3,2L,RND);	×
415	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
416	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
417	break;	×
418	case 5:
419	label5:	×
420	mpfr_sqrt_ui(auxr3,3L,RND);	×
421	mpfr_add_ui(auxr3,auxr3,1L,RND);	×
422	mpfr_sqrt_ui(auxr1,2L,RND);	×
423	mpfr_div_ui(auxr1,auxr1,4L,RND);	×
424	mpfr_mul(auxr3,auxr3,auxr1,RND);	×
425	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
426	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
427	break;	×
428	case 6:
429	label6:	×
430	if ( nsgn < 0 ) mpfr_neg(auxr2,auxr2,RND);	×
431	break;
432	}
433	*t = 0; first = 0;	×
434	goto nextfun;	×
435	}
436	}
437	else if ( *t == COSFUNCTION ) {	×
438	pimul = GetPiArgument(BHEAD t+FUNHEAD);	×
439	if ( pimul >= 0 && pimul < 24 ) {	×
440	if ( pimul > 12 ) pimul = 24-12;	×
441	if ( pimul > 6 ) { pimul = 12-pimul; nsgn = -1; }	×
442	else nsgn = 1;
443	if ( pimul == 6 ) goto getout;	×
444	switch ( pimul ) {	×
445	case 0: goto label6;	×
446	case 1: goto label5;	×
447	case 2: goto label4;	×
448	case 3: goto label3;	×
449	case 4: goto label2;	×
450	case 5: goto label1;	×
451	}
452	}
453	}
454	else if ( *t == TANFUNCTION ) {	×
455	pimul = GetPiArgument(BHEAD t+FUNHEAD);	×
456	if ( pimul >= 0 && pimul < 24 ) {	×
457	if ( pimul == 6 \|\| pimul == 18 ) goto nextfun;	×
458	if ( pimul == 0 \|\| pimul == 12 ) goto getout;	×
459	if ( pimul > 12 ) { pimul = 24-pimul; nsgn = -1; }	×
460	else nsgn = 1;
461	if ( pimul > 6 ) { pimul = 12-pimul; nsgn = -nsgn; }	×
462	switch ( pimul ) {	×
463	case 1:	×
464	mpfr_sqrt_ui(auxr3,3L,RND);	×
465	mpfr_sub_ui(auxr3,auxr3,2L,RND);	×
466	if ( nsgn > 0 ) mpfr_neg(auxr3,auxr3,RND);	×
467	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
468	break;	×
469	case 2:	×
470	mpfr_sqrt_ui(auxr3,3L,RND);	×
471	mpfr_div_ui(auxr3,auxr3,3L,RND);	×
472	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
473	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
474	break;	×
475	case 3:
476	break;
477	case 4:	×
478	mpfr_sqrt_ui(auxr3,3L,RND);	×
479	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
480	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
481	break;	×
482	case 5:	×
483	mpfr_sqrt_ui(auxr3,3L,RND);	×
484	mpfr_add_ui(auxr3,auxr3,2L,RND);	×
485	if ( nsgn < 0 ) mpfr_neg(auxr3,auxr3,RND);	×
486	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
487	break;	×
488	}
489	*t = 0;	×
490	goto nextfun;	×
491	}
492	}
493
494	if ( *t == AGMFUNCTION ) {	×
495	if ( GetFloatArgument(BHEAD auxr1,t,1) < 0 ) goto nextfun;	×
496	if ( GetFloatArgument(BHEAD auxr3,t,-2) < 0 ) goto nextfun;	×
497	}
498	else if ( GetFloatArgument(BHEAD auxr1,t,-1) < 0 ) goto nextfun;	×
499	nsgn = mpfr_sgn(auxr1);	×
500
501	switch ( *t ) {	×
502	case SQRTFUNCTION:	×
503	if ( nsgn < 0 ) goto nextfun;	×
504	if ( nsgn == 0 ) goto getout;	×
505	else mpfr_sqrt(auxr3,auxr1,RND);	×
506	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
507	break;	×
508	case LNFUNCTION:	×
509	if ( nsgn <= 0 ) goto nextfun;	×
510	else mpfr_log(auxr3,auxr1,RND);	×
511	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
512	break;	×
513	case LI2FUNCTION: /* should be between -1 and +1 */	×
514	mpfr_abs(auxr3,auxr1,RND);	×
515	if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;	×
516	mpfr_li2(auxr3,auxr1,RND);	×
517	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
518	break;	×
519	case GAMMAFUN:	×
520	/*
521	We cannot do this when the argument is a non-positive integer
522	*/
523	if ( t[FUNHEAD] == -SNUMBER && t[FUNHEAD+1] <= 0 ) goto nextfun;	×
524	if ( t[FUNHEAD] == t[1]-FUNHEAD && ABS(t[t[1]-1]) ==	×
525	t[FUNHEAD]-ARGHEAD-1 && t[t[1]-1] < 0 ) {	×
526	nsize = (-t[t[1]-1]-1)/2;	×
527	if ( t[t[1]-1-nsize] == 1 ) {	×
528	for ( i = 1; i < nsize; i++ ) {	×
529	if ( t[t[1]-1-nsize+i] != 0 ) break;	×
530	}
531	if ( i >= nsize ) goto nextfun;	×
532	}
533	}
534	mpfr_gamma(auxr3,auxr1,RND);	×
535	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
536	break;	×
537	case AGMFUNCTION:	×
538	mpfr_agm(auxr3,auxr1,auxr3,RND);	×
539	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
540	break;	×
541	case SINHFUNCTION:	×
542	mpfr_sinh(auxr3,auxr1,RND);	×
543	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
544	break;	×
545	case COSHFUNCTION:	×
546	mpfr_cosh(auxr3,auxr1,RND);	×
547	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
548	break;	×
549	case TANHFUNCTION:	×
550	mpfr_tanh(auxr3,auxr1,RND);	×
551	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
552	break;	×
553	case ASINHFUNCTION:	×
554	mpfr_asinh(auxr3,auxr1,RND);	×
555	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
556	break;	×
557	case ACOSHFUNCTION:	×
558	mpfr_acosh(auxr3,auxr1,RND);	×
559	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
560	break;	×
561	case ATANHFUNCTION:	×
562	if ( nsgn == 0 ) goto getout;	×
563	mpfr_abs(auxr3,auxr1,RND);	×
564	if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;	×
565	mpfr_atanh(auxr3,auxr1,RND);	×
566	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
567	break;	×
568	case ASINFUNCTION:	×
569	if ( nsgn == 0 ) goto getout;	×
570	mpfr_abs(auxr3,auxr1,RND);	×
571	if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;	×
572	mpfr_asin(auxr3,auxr1,RND);	×
573	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
574	break;	×
575	case ACOSFUNCTION:	×
576	if ( nsgn == 0 ) goto getout;	×
577	mpfr_abs(auxr3,auxr1,RND);	×
578	if ( mpfr_cmp_ui(auxr3,1L) > 0 ) goto nextfun;	×
579	mpfr_acos(auxr3,auxr1,RND);	×
580	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
581	break;	×
582	case ATANFUNCTION:	×
583	mpfr_atan(auxr3,auxr1,RND);	×
584	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
585	break;	×
586	case SINFUNCTION:	×
587	mpfr_sin(auxr3,auxr1,RND);	×
588	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
589	break;	×
590	case COSFUNCTION:	×
591	mpfr_cos(auxr3,auxr1,RND);	×
592	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
593	break;	×
594	case TANFUNCTION:	×
595	mpfr_tan(auxr3,auxr1,RND);	×
596	mpfr_mul(auxr2,auxr2,auxr3,RND);	×
597	break;	×
598	default:	×
599	goto nextfun;	×
600	break;	×
601	}
602	first = 0;	×
603	*t = 0;	×
604	goto nextfun;	×
605	}
606	else if ( pars[2] == ALLFUNCTIONS ) {	×
607	/*
608	Now we have to test whether this is one of our functions
609	*/
610	if ( t[1] == FUNHEAD ) goto nextfun;	×
611	switch ( *t ) {	×
612	case SQRTFUNCTION:	×
613	case LNFUNCTION:
614	case LI2FUNCTION:
615	case LINFUNCTION:
616	case ASINFUNCTION:
617	case ACOSFUNCTION:
618	case ATANFUNCTION:
619	case SINHFUNCTION:
620	case COSHFUNCTION:
621	case TANHFUNCTION:
622	case ASINHFUNCTION:
623	case ACOSHFUNCTION:
624	case ATANHFUNCTION:
625	case SINFUNCTION:
626	case COSFUNCTION:
627	case TANFUNCTION:
628	case AGMFUNCTION:
629	goto TestArgument;	×
630	case MZV:	×
631	case EULER:
632	case MZVHALF:
633	goto nextfun;	×
634	default:	×
635	MLOCK(ErrorMessageLock);	×
636	MesPrint("Function in evaluate statement not yet implemented.");	×
637	MUNLOCK(ErrorMessageLock);	×
638	break;
639	}
640	}
641	else goto nextfun;	×
642	nextfun:	×
643	t += t[1];	×
644	}
645	if ( first == 1 ) return(Generator(BHEAD term,level));	×
646	mpfr_get_f(aux4,auxr2,RND);	×
647	/*
648	Step 3:
649	Now the regular coefficient, if it is not 1/1.
650	We have two cases: size +- 3, or bigger.
651	*/
652
653	nsize = term[*term-1];	×
654	if ( nsize < 0 ) { nsize = -nsize; nsgn = -1; }	×
655	else nsgn = 1;
656	if ( aux4->_mp_size < 0 ) {	×
657	aux4->_mp_size = -aux4->_mp_size;	×
658	nsgn = -nsgn;	×
659	}
660
661	if ( nsize == 3 ) {	×
662	if ( tstop[0] != 1 ) {	×
663	mpf_mul_ui(aux4,aux4,(ULONG)((UWORD)tstop[0]));	×
664	}
665	if ( tstop[1] != 1 ) {	×
666	mpf_div_ui(aux4,aux4,(ULONG)((UWORD)tstop[1]));	×
667	}
668	}
669	else {
670	RatToFloat(aux5,(UWORD *)tstop,nsize);	×
671	mpf_mul(aux4,aux4,aux5);	×
672	}
673	/*
674	Now we have to locate possible other float_ functions.
675	Note possible incompatibilities between the mpf and mpfr formats.
676	*/
677	t = term+1;
678	while ( t < tstop ) {	×
679	if ( *t == FLOATFUN ) {	×
680	UnpackFloat(aux5,t);	×
681	mpf_mul(aux4,aux4,aux5);	×
682	}
683	t += t[1];	×
684	}
685	/*
686	Now we should compose the new term in the WorkSpace.
687	*/
688	t = term+1;	×
689	newterm = AT.WorkPointer;	×
690	tt = newterm+1;	×
691	while ( t < tstop ) {	×
692	if ( t == 0 \|\| t == FLOATFUN ) t += t[1];	×
693	else {
694	i = t[1]; NCOPY(tt,t,i);	×
695	}
696	}
697	PackFloat(tt,aux4);	×
698	tt += tt[1];	×
699	tt++ = 1; tt++ = 1; tt++ = 3nsgn;	×
700	*newterm = tt-newterm;	×
701	AT.WorkPointer = tt;	×
702	retval = Generator(BHEAD newterm,level);	×
703	getout:	×
704	AT.WorkPointer = oldworkpointer;	×
705	return(retval);	×
706	}
707
708	/*
709	#] EvaluateFun :
710	#] mpfr_ :
711	*/

vermaseren / form / 9364948935

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