float.c

Go to the documentation of this file.

 
/* #[ License : */
/*
 *   Copyright (C) 1984-2026 J.A.M. Vermaseren
 *   When using this file you are requested to refer to the publication
 *   J.A.M.Vermaseren "New features of FORM" math-ph/0010025
 *   This is considered a matter of courtesy as the development was paid
 *   for by FOM the Dutch physics granting agency and we would like to
 *   be able to track its scientific use to convince FOM of its value
 *   for the community.
 *
 *   This file is part of FORM.
 *
 *   FORM is free software: you can redistribute it and/or modify it under the
 *   terms of the GNU General Public License as published by the Free Software
 *   Foundation, either version 3 of the License, or (at your option) any later
 *   version.
 *
 *   FORM is distributed in the hope that it will be useful, but WITHOUT ANY
 *   WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 *   FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 *   details.
 *
 *   You should have received a copy of the GNU General Public License along
 *   with FORM.  If not, see <http://www.gnu.org/licenses/>.
 */
/* #] License : */ 
/*
    #[ Includes : float.c
*/
 
#include "form3.h"
#include <math.h>
#include <gmp.h>
 
#define WITHCUTOFF
 
#define GMPSPREAD (GMP_LIMB_BITS/BITSINWORD)
 
 
void Form_mpf_init(mpf_t t);
void Form_mpf_clear(mpf_t t);
void Form_mpf_set_prec_raw(mpf_t t,ULONG newprec);
void FormtoZ(mpz_t z,UWORD *a,WORD na);
void ZtoForm(UWORD *a,WORD *na,mpz_t z);
long FloatToInteger(UWORD *out, mpf_t floatin, long *bitsused);
void IntegerToFloat(mpf_t result, UWORD *formlong, int longsize);
int FloatToRat(UWORD *ratout, WORD *nratout, mpf_t floatin);
int AddFloats(PHEAD WORD *fun3, WORD *fun1, WORD *fun2);
int MulFloats(PHEAD WORD *fun3, WORD *fun1, WORD *fun2);
int DivFloats(PHEAD WORD *fun3, WORD *fun1, WORD *fun2);
int AddRatToFloat(PHEAD WORD *outfun, WORD *infun, UWORD *formrat, WORD nrat);
int MulRatToFloat(PHEAD WORD *outfun, WORD *infun, UWORD *formrat, WORD nrat);
void SimpleDelta(mpf_t sum, int m);
void SimpleDeltaC(mpf_t sum, int m);
void SingleTable(mpf_t *tabl, int N, int m, int pow);
void DoubleTable(mpf_t *tabout, mpf_t *tabin, int N, int m, int pow);
void EndTable(mpf_t sum, mpf_t *tabin, int N, int m, int pow);
void deltaMZV(mpf_t, WORD *, int);
void deltaEuler(mpf_t, WORD *, int);
void deltaEulerC(mpf_t, WORD *, int);
void CalculateMZVhalf(mpf_t, WORD *, int);
void CalculateMZV(mpf_t, WORD *, int);
void CalculateEuler(mpf_t, WORD *, int);
int ExpandMZV(WORD *term, WORD level);
int ExpandEuler(WORD *term, WORD level);
int PackFloat(WORD *,mpf_t);
int UnpackFloat(mpf_t, WORD *);
void RatToFloat(mpf_t result, UWORD *formrat, int ratsize);
 
/*
    #] Includes : 
    #[ Some GMP float info :
 
    From gmp.h:
 
    typedef struct
    {
      int _mp_prec;     Max precision, in number of `mp_limb_t's.
                        Set by mpf_init and modified by
                        mpf_set_prec.  The area pointed to by the
                        _mp_d field contains `prec' + 1 limbs.
      int _mp_size;     abs(_mp_size) is the number of limbs the
                        last field points to.  If _mp_size is
                        negative this is a negative number.
      mp_exp_t _mp_exp; Exponent, in the base of `mp_limb_t'.
      mp_limb_t *_mp_d; Pointer to the limbs.
    } __mpf_struct;
 
    typedef __mpf_struct mpf_t[1];
 
    Currently:
        sizeof(int) = 4
        sizeof(mpf_t) = 24
        sizeof(mp_limb_t) = 8,
        sizeof(mp_exp_t) = 8,
        sizeof a pointer is 8.
    If any of this changes in the future we would have to change the
    PackFloat and UnpackFloat routines.
 
    Our float_ function is packed with the 4 arguments:
        -SNUMBER _mp_prec
        -SNUMBER _mp_size
        exponent which can be -SNUMBER or a regular term
        the limbs as n/1 in regular term format, or just an -SNUMBER.
    During normalization the sign is taken away from the second argument
    and put as the sign of the complete term. This is easier for the
    WriteInnerTerm routine in sch.c
 
    Wildcarding of functions excludes matches with float_
    Float_ always ends up inside the bracket, just like poly(rat)fun.
 
    From mpfr.h
 
    The formats here are different. This has, among others, to do with
    the rounding. The result is that we need different aux variables
    when working with mpfr and we need a different conversion routine
    to float. It also means that we need to treat the mzv_ etc functions
    completely separately. For now we try to develop that in the file
    Evaluate.c. At a later stage we may merge the two files.
    Originally we had the sqrt in float.c but it seems better to put
    it in Evaluate.c
 
    #] Some GMP float info : 
    #[ Low Level :
 
        In the low level routines we interact directly with the content
        of the GMP structs. This can be done safely, because their
        layout is documented. We pay particular attention to the init
        and clear routines, because they involve malloc/free calls.
 
        #[ Explanations :
 
    We need a number of facilities inside Form to deal with floating point
    numbers, taking into account that they are only meant for sporadic use.
    1: We need a special function float_ who's arguments contain a limb
       representation of the mpf_t of the gmp.
    2: Some functions may return a floating point number. This is then in
       the form of an occurrence of the function float_.
    3: We need a print routine to print this float_.
    4: We need routines for conversions:
        a: (long) int to float.
        b: rat to float.
        c: float to rat (if possible)
        d: float to integer (with rounding toward zero).
    5: We need calculational operations for add, sub, mul, div, exp.
    6: We need service routines to pack and unpack the function float_.
    7: The coefficient should be pulled inside float_.
    8: Float_ cannot occur inside PolyRatFun.
    9: We need to be able to treat float_ as the coefficient during sorting.
    10: For now, we need the functions mzvf_ and eulerf_ for the evaluation
        of mzv's and euler sums.
    11: We need a setting for the default precision.
    12: We need a special flag for the dirty field of arguments to avoid
        the normalization of the argument with the limbs.
        Alternatively, the float_ should be not a regular function, but
        get a status < FUNCTION. Just like SNUMBER, LNUMBER etc.
 
    Note that we can keep this thread-safe. The only routine that is not
    thread-safe is the print routine, but that is in accordance with all
    other print routines in Form. When called from MesPrint it needs to
    be protected by a lock, as is done standard inside Form.
 
        #] Explanations : 
        #[ Form_mpf_init :
*/
 
void Form_mpf_init(mpf_t t)
{
    mp_limb_t *d;
    int i, prec;
    if ( t->_mp_d ) { M_free(t->_mp_d,"Form_mpf_init"); t->_mp_d = 0; }
    prec = (AC.DefaultPrecision + 8*sizeof(mp_limb_t)-1)/(8*sizeof(mp_limb_t)) + 1;
    t->_mp_prec = prec;
    t->_mp_size = 0;
    t->_mp_exp = 0;
    d = (mp_limb_t *)Malloc1(prec*sizeof(mp_limb_t),"Form_mpf_init");
    if ( d == 0 ) {
        MesPrint("Fatal error in Malloc1 call in Form_mpf_init. Trying to allocate %ld bytes.",
            prec*sizeof(mp_limb_t));
        Terminate(-1);
    }
    t->_mp_d = d;
    for ( i = 0; i < prec; i++ ) d[i] = 0;
}
 
/*
        #] Form_mpf_init : 
        #[ Form_mpf_clear :
*/
 
void Form_mpf_clear(mpf_t t)
{
    if ( t->_mp_d ) { M_free(t->_mp_d,"Form_mpf_init"); t->_mp_d = 0; }
    t->_mp_prec = 0;
    t->_mp_size = 0;
    t->_mp_exp = 0;
}
 
/*
        #] Form_mpf_clear : 
        #[ Form_mpf_empty :
*/
 
void Form_mpf_empty(mpf_t t)
{
    int i;
    for ( i = 0; i < t->_mp_prec; i++ ) t->_mp_d[i] = 0;
    t->_mp_size = 0;
    t->_mp_exp = 0;
}
 
/*
        #] Form_mpf_empty : 
        #[ Form_mpf_set_prec_raw :
*/
 
void Form_mpf_set_prec_raw(mpf_t t,ULONG newprec)
{
    ULONG newpr = (newprec + 8*sizeof(mp_limb_t)-1)/(8*sizeof(mp_limb_t)) + 1;
    if ( newpr <= (ULONG)(t->_mp_prec) ) {
        int i, used = ABS(t->_mp_size) ;
        if ( newpr > (ULONG)used ) {
            for ( i = used; i < (int)newpr; i++ ) t->_mp_d[i] = 0;
        }
        if ( t->_mp_size < 0 ) newpr = -newpr;
        t->_mp_size = newpr;
    }
    else {
/*
        We can choose here between "forget about it" or making a new malloc.
        If the program is well designed, we should never get here though.
        For now we choose the "forget it" option.
*/
        MesPrint("Trying too big a precision %ld in Form_mpf_set_prec_raw.",newprec);
        MesPrint("Old maximal precision is %ld.",
            (size_t)(t->_mp_prec-1)*sizeof(mp_limb_t)*8);
        Terminate(-1);
    }
}
 
/*
        #] Form_mpf_set_prec_raw : 
        #[ PackFloat :
 
    Puts an object of type mpf_t inside a function float_
    float_(prec, nlimbs, exp, limbs);
    Complication: the limbs and the exponent are long's. That means
    two times the size of a Form (U)WORD.
    To save space we could split the limbs in two because Form stores
    coefficients as rationals with equal space for numerator and denominator.
    Hence for each limb we could put the most significant UWORD in the 
    numerator and the least significant limb in the denominator.
    This gives a problem with the compilation and subsequent potential
    normalization of the 'fraction'. Hence for now we take the wasteful
    approach. At a later stage we can try to veto normalization of FLOATFUN.
    Caution: gmp/fmp is little endian and Form is big endian.
 
    Zero is represented by zero limbs (and zero exponent).
*/
 
int PackFloat(WORD *fun,mpf_t infloat)
{
    WORD *t, nlimbs, num, nnum;
    mp_limb_t *d = infloat->_mp_d; /* Pointer to the limbs.  */
    int i;
    long e = infloat->_mp_exp;
    t = fun;
    *t++ = FLOATFUN;
    t++;
    FILLFUN(t);
    *t++ = -SNUMBER;
    *t++ = infloat->_mp_prec;
    *t++ = -SNUMBER;
    *t++ = infloat->_mp_size;
/*
    Now the exponent which is a signed long
*/
    if ( e < 0 ) {
        e = -e;
        if ( ( e >> (BITSINWORD-1) ) == 0 ) {
            *t++ = -SNUMBER; *t++ = -e;
        }
        else {
            *t++ = ARGHEAD+6; *t++ = 0; FILLARG(t);
            *t++ = 6;
            *t++ = (UWORD)e;
            *t++ = (UWORD)(e >> BITSINWORD);
            *t++ = 1; *t++ = 0; *t++ = -5;
        }
    }
    else {
        if ( ( e >> (BITSINWORD-1) ) == 0 ) {
            *t++ = -SNUMBER; *t++ = e;
        }
        else {
            *t++ = ARGHEAD+6; *t++ = 0; FILLARG(t);
            *t++ = 6;
            *t++ = (UWORD)e;
            *t++ = (UWORD)(e >> BITSINWORD);
            *t++ = 1; *t++ = 0; *t++ = 5;
        }
    }
/*
    and now the limbs. Note that the number of active limbs could be less
    than prec+1 in which case we copy the smaller number.
*/
    nlimbs = infloat->_mp_size < 0 ? -infloat->_mp_size: infloat->_mp_size;
    if ( nlimbs == 0 ) {
        *t++ = -SNUMBER;
        *t++ = 0;
    }
    else if ( nlimbs == 1 && (ULONG)(*d) < ((ULONG)1)<<(BITSINWORD-1) ) {
        *t++ = -SNUMBER;
        *t++ = (UWORD)(*d);
    }
    else {
      if ( d[nlimbs-1] < ((ULONG)1)<<BITSINWORD ) {
        num = 2*nlimbs-1; /* number of UWORDS needed */
      }
      else {
        num = 2*nlimbs;   /* number of UWORDS needed */
      }
      nnum = num;
      *t++ = 2*num+2+ARGHEAD;
      *t++ = 0;
      FILLARG(t);
      *t++ = 2*num+2;
      while ( num > 1 ) {
        *t++ = (UWORD)(*d); *t++ = (UWORD)(((ULONG)(*d))>>BITSINWORD); d++;
        num -= 2;
      }
      if ( num == 1 ) { *t++ = (UWORD)(*d); }
      *t++ = 1;
      for ( i = 1; i < nnum; i++ ) *t++ = 0;
      *t++ = 2*nnum+1; /* the sign is hidden in infloat->_mp_size */
    }
    fun[1] = t-fun;
    return(fun[1]);
}
 
/*
        #] PackFloat : 
        #[ UnpackFloat :
 
    Takes the arguments of a function float_ and puts it inside an mpf_t.
    For notation see commentary for PutInFloat.
 
    We should make a new allocation for the limbs if the precision exceeds
    the present precision of outfloat, or when outfloat has not been
    initialized yet.
 
    The return value is -1 if the contents of the FLOATFUN cannot be converted.
    Otherwise the return value is the number of limbs.
*/
 
int UnpackFloat(mpf_t outfloat,WORD *fun)
{
    UWORD *t;
    WORD *f;
    int num, i;
    mp_limb_t *d;
/*
    Very first step: check whether we can use float_:
*/
    GETIDENTITY
    if ( AT.aux_ == 0 ) {
        MLOCK(ErrorMessageLock);
        MesPrint("Illegal attempt at using a float_ function without proper startup.");
        MesPrint("Please use %#StartFloat <options> first.");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
/*
    Check first the number and types of the arguments
    There should be 4. -SNUMBER,-SNUMBER,-SNUMBER or a long number.
    + the limbs, either -SNUMBER or Long number in the form of a Form rational.
*/
    if ( TestFloat(fun) == 0 ) goto Incorrect;
    f = fun + FUNHEAD + 2;
    if ( ABS(f[1]) > f[-1]+1 ) goto Incorrect;
    if ( f[1] > outfloat->_mp_prec+1
      || outfloat->_mp_d == 0 ) {
        /*
            We need to make a new allocation.
            Unfortunately we cannot use Malloc1 because that is not
            recognised by gmp and hence we cannot clear with mpf_clear.
            Maybe we can do something about it by making our own
            mpf_init and mpf_clear?
        */                     
        if ( outfloat->_mp_d != 0 ) free(outfloat->_mp_d);
        outfloat->_mp_d = (mp_limb_t *)malloc((f[1]+1)*sizeof(mp_limb_t));
    }
    num = f[1];
    outfloat->_mp_size = num;
    if ( num < 0 ) { num = -num; }
    f += 2;
    if ( *f == -SNUMBER ) {
        outfloat->_mp_exp = (mp_exp_t)(f[1]);
        f += 2;
    }
    else {
        f += ARGHEAD+6;
        if ( f[-1] == -5 ) {
            outfloat->_mp_exp =
                -(mp_exp_t)((((ULONG)(f[-4]))<<BITSINWORD)+(UWORD)f[-5]);
        }
        else if ( f[-1] == 5 ) {
            outfloat->_mp_exp =
                 (mp_exp_t)((((ULONG)(f[-4]))<<BITSINWORD)+(UWORD)f[-5]);
        }
    }
/*
    And now the limbs if needed
*/
    d = outfloat->_mp_d;
    if ( outfloat->_mp_size == 0 ) {
        for ( i = 0; i < outfloat->_mp_prec; i++ ) *d++ = 0;
        return(0);
    }
    else if ( *f == -SNUMBER ) {
        *d++ = (mp_limb_t)f[1];
        for ( i = 0; i < outfloat->_mp_prec; i++ ) *d++ = 0;
        return(0);
    }
    num = (*f-ARGHEAD-2)/2; /* 2*number of limbs in the argument */
    t = (UWORD *)(f+ARGHEAD+1);
    while ( num > 1 ) {
        *d++ = (mp_limb_t)((((ULONG)(t[1]))<<BITSINWORD)+t[0]);
        t += 2; num -= 2;
    }
    if ( num == 1 ) *d++ = (mp_limb_t)((UWORD)(*t));
    for ( i = d-outfloat->_mp_d; i <= outfloat->_mp_prec; i++ ) *d++ = 0;
    return(0);
Incorrect:
    return(-1);
}
 
/*
        #] UnpackFloat : 
        #[ TestFloat :
 
        Tells whether the function (with its arguments) is a legal float_.
        We assume, as in the whole float library that one limb is 2 WORDs.
*/
 
int TestFloat(WORD *fun)
{
    WORD *fstop, *f, num, nnum, i;
    fstop = fun+fun[1];
    f = fun + FUNHEAD;
    num = 0;
/*
    Count arguments
*/
    while ( f < fstop ) { num++; NEXTARG(f); }
    if ( num != 4 ) return(0);
    f = fun + FUNHEAD;
    if ( *f != -SNUMBER ) return(0);
    if ( f[1] < 0 ) return(0);
    f += 2;
    if ( *f != -SNUMBER ) return(0);
    num = ABS(f[1]); /* number of limbs */
    f += 2;
    if ( *f == -SNUMBER ) { f += 2; }
    else {
        if ( *f != ARGHEAD+6 ) return(0);
        if ( ABS(f[ARGHEAD+5]) != 5 ) return(0);
        if ( f[ARGHEAD+3] != 1 ) return(0);
        if ( f[ARGHEAD+4] != 0 ) return(0);
        f += *f;
    }
    if ( num == 0 ) return(1);
    if ( *f == -SNUMBER ) { if ( num != 1 ) return(0); }
    else {
        nnum = (ABS(f[*f-1])-1)/2;
        if ( ( *f != 4*num+2+ARGHEAD ) && ( *f != 4*num+ARGHEAD ) ) return(0);
        if ( ( nnum != 2*num ) && ( nnum != 2*num-1 ) ) return(0);
        f += ARGHEAD+nnum+1;
        if ( f[0] != 1 ) return(0);
        for ( i = 1; i < nnum; i++ ) {
            if ( f[i] != 0 ) return(0);
        }
    }
    return(1);
}
 
/*
        #] TestFloat : 
        #[ FormtoZ :
 
        Converts a Form long integer to a GMP long integer.
 
        typedef struct
        {
          int _mp_alloc;   Number of *limbs* allocated and pointed
                           to by the _mp_d field.
          int _mp_size;    abs(_mp_size) is the number of limbs the
                           last field points to.  If _mp_size is
                           negative this is a negative number.
          mp_limb_t *_mp_d;Pointer to the limbs.
        } __mpz_struct;
        typedef __mpz_struct mpz_t[1];
*/
 
void FormtoZ(mpz_t z,UWORD *a,WORD na)
{
    int nlimbs, sgn = 1;
    mp_limb_t *d;
    UWORD *b = a;
    WORD nb = na;
 
    if ( nb == 0 ) {
        z->_mp_size = 0;
        z->_mp_d[0] = 0;
        return;
    }
    if ( nb < 0 ) { sgn = -1; nb = -nb; }
    nlimbs = (nb+1)/2;
    if ( mpz_cmp_si(z,0L) <= 1 ) {
        mpz_set_ui(z,10000000);
    }
    if ( z->_mp_alloc <= nlimbs ) {
        mpz_t zz;
        mpz_init(zz);
        while ( z->_mp_alloc <= nlimbs ) {
            mpz_mul(zz,z,z);
            mpz_set(z,zz);
        }
        mpz_clear(zz);
    }
    z->_mp_size = sgn*nlimbs;
    d = z->_mp_d;
    while ( nb > 1 ) {
        *d++ = (((mp_limb_t)(b[1]))<<BITSINWORD)+(mp_limb_t)(b[0]);
        b += 2; nb -= 2;
    }
    if ( nb == 1 ) { *d = (mp_limb_t)(*b); }
}
 
/*
        #] FormtoZ : 
        #[ ZtoForm :
 
        Converts a GMP long integer to a Form long integer.
        Mainly an exercise of going from little endian ULONGs.
        to big endian UWORDs
*/
 
void ZtoForm(UWORD *a,WORD *na,mpz_t z)
{
    mp_limb_t *d = z->_mp_d;
    int nlimbs = ABS(z->_mp_size), i;
    UWORD j;
    if ( nlimbs == 0 ) { *na = 0; return; }
    *na = 0;
    for ( i = 0; i < nlimbs-1; i++ ) {
        *a++ = (UWORD)(*d);
        *a++ = (UWORD)((*d++)>>BITSINWORD);
        *na += 2;
    }
    *na += 1;
    *a++ = (UWORD)(*d);
    j = (UWORD)((*d)>>BITSINWORD);
    if ( j != 0 ) { *a = j; *na += 1; }
    if ( z->_mp_size < 0 ) *na = -*na;
}
 
/*
        #] ZtoForm : 
        #[ FloatToInteger :
 
        Converts a GMP float to a long Form integer.
*/
 
long FloatToInteger(UWORD *out, mpf_t floatin, long *bitsused)
{
    mpz_t z;
    WORD nout, nx;
    UWORD x;
    mpz_init(z);
    mpz_set_f(z,floatin);
    ZtoForm(out,&nout,z);
    mpz_clear(z);
    x = out[nout-1]; nx = 0;
    while ( x ) { nx++; x >>= 1; }
    *bitsused = (nout-1)*BITSINWORD + nx;
    return(nout);
}
 
/*
        #] FloatToInteger : 
        #[ IntegerToFloat :
 
        Converts a Form long integer to a gmp float of default size.
        We assume that sizeof(unsigned long int) = 2*sizeof(UWORD).
*/
 
void IntegerToFloat(mpf_t result, UWORD *formlong, int longsize)
{
    mpz_t z;
    mpz_init(z);
    FormtoZ(z,formlong,longsize);
    mpf_set_z(result,z);
    mpz_clear(z);
}
 
/*
        #] IntegerToFloat : 
        #[ RatToFloat :
 
        Converts a Form rational to a gmp float of default size.
*/
 
void RatToFloat(mpf_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 ( num[nnum-1] == 0 ) { nnum--; }
    while ( den[nden-1] == 0 ) { nden--; }
    IntegerToFloat(aux2,num,nnum);
    IntegerToFloat(aux3,den,nden);
    mpf_div(result,aux2,aux3);
    if ( sgn > 0 ) mpf_neg(result,result);
}
 
/*
        #] RatToFloat : 
        #[ FloatFunToRat :
*/
 
WORD FloatFunToRat(PHEAD UWORD *ratout,WORD *in)
{
    WORD oldin = in[0], nratout;
    in[0] = FLOATFUN;
    UnpackFloat(aux4,in);
    FloatToRat(ratout,&nratout,aux4);
    in[0] = oldin;
    return(nratout);
}
 
/*
        #] FloatFunToRat : 
        #[ FloatToRat :
 
    Converts a gmp float to a Form rational.
    The algorithm we use is 'repeated fractions' of the style
        n0+1/(n1+1/(n2+1/(n3+.....)))
    The moment we get an n that is very large we should have the proper fraction
    Main problem: what if the sequence keeps going?
                 (there are always rounding errors)
    We need a cutoff with an error condition.
    Basically the product n0*n1*n2*... is an underlimit on the number of
    digits in the answer. We can put a limit on this.
    A return value of zero means that everything went fine.
*/
 
int FloatToRat(UWORD *ratout, WORD *nratout, mpf_t floatin)
{
    GETIDENTITY
    WORD *out = AT.WorkPointer;
    WORD s; /* the sign. */
    UWORD *a, *b, *c, *d, *mc;
    WORD na, nb, nc, nd, i;
    int nout;
    LONG oldpWorkPointer = AT.pWorkPointer;
    long bitsused = 0, totalbitsused = 0, totalbits = AC.DefaultPrecision-AC.MaxWeight-1;
    int retval = 0, startnul;
    WantAddPointers(AC.DefaultPrecision);  /* Horrible overkill */
    AT.pWorkSpace[AT.pWorkPointer++] = out;
    a = NumberMalloc("FloatToRat");
    b = NumberMalloc("FloatToRat");
 
    s = mpf_sgn(floatin);
    if ( s < 0 ) mpf_neg(floatin,floatin);
 
    Form_mpf_empty(aux1);
    Form_mpf_empty(aux2);
    Form_mpf_empty(aux3);
 
    mpf_trunc(aux1,floatin);     /* This should now be an integer */
    mpf_sub(aux2,floatin,aux1);  /* and this >= 0 */
    if ( mpf_sgn(aux1) == 0 ) { *out++ = 0; startnul = 1; }
    else {
        nout = FloatToInteger((UWORD *)out,aux1,&totalbitsused);
        out += nout;
        startnul = 0;
    }
    AT.pWorkSpace[AT.pWorkPointer++] = out;
    if ( mpf_sgn(aux2) ) {
      for(;;) {
        mpf_ui_div(aux3,1,aux2);
        mpf_trunc(aux1,aux3);
        mpf_sub(aux2,aux3,aux1);
        if ( mpf_sgn(aux1) == 0 ) { *out++ = 0; }
        else {
            nout = FloatToInteger((UWORD *)out,aux1,&bitsused);
            out += nout;
        }
        if ( bitsused > (totalbits-totalbitsused)/2 ) { break; }
        if ( mpf_sgn(aux2) == 0 ) {
            /*if ( startnul == 1 )*/ AT.pWorkSpace[AT.pWorkPointer++] = out;
            break;
        }
        totalbitsused += bitsused;
        AT.pWorkSpace[AT.pWorkPointer++] = out;
      }
      /* 
        For |floatin| << 1, we have startnul = 1 and hit the precision guard 
        already on the first continued-fraction step. The resulting float is 
        therefore close to the rational 0/1 and we can immediately return.
       */
      if ( startnul == 1 && AT.pWorkPointer == oldpWorkPointer + 2 ) {
        *ratout++ = 0; *ratout++ = 1; *ratout = *nratout = 3;
        goto ret;
      }
/*
      At this point we have the function with the repeated fraction.
      Now we should work out the fraction. Form code would be:
        repeat id   dum_(?a,x1?,x2?) = dum_(?a,x1+1/x2);
        id  dum_(x?) = x;
      We have to work backwards. This is why we made a list of pointers
      in AT.pWorkSpace
 
      Now we need the long integers a and b.
      Note that by construction there should never be a GCD!
*/
      out = (WORD *)(AT.pWorkSpace[--AT.pWorkPointer]);
      na = 1; a[0] = 1;
      c = (UWORD *)(AT.pWorkSpace[--AT.pWorkPointer]);
      nc = nb = ((UWORD *)out)-c;
      for ( i = 0; i < nb; i++ ) b[i] = c[i];
      mc = c = NumberMalloc("FloatToRat");
      while ( AT.pWorkPointer > oldpWorkPointer ) {
        d = (UWORD *)(AT.pWorkSpace[--AT.pWorkPointer]);
        nd = (UWORD *)(AT.pWorkSpace[AT.pWorkPointer+1])-d; /* This is the x1 in the formula */
        if ( nd == 1 && *d == 0 ) break;
        c = a; a = b; b = c; nc = na; na = nb; nb = nc; /* 1/x2 */
        MulLong(d,nd,a,na,mc,&nc);
        AddLong(mc,nc,b,nb,b,&nb);
      }
      NumberFree(mc,"FloatToRat");
      if ( startnul == 0 ) {
        c = a; a = b; b = c; nc = na; na = nb; nb = nc;
      }
    }
    else {
        c = (UWORD *)(AT.pWorkSpace[oldpWorkPointer]);
        na = (UWORD *)out - c;
        for ( i = 0; i < na; i++ ) a[i] = c[i];
        b[0] = 1; nb = 1;
    }
/*
    Now we have the fraction in a/b. Create the rational and add the sign.
*/
    d = ratout;
    if ( na == nb ) {
        *nratout = (2*na+1)*s;
        for ( i = 0; i < na; i++ ) *d++ = a[i];
        for ( i = 0; i < nb; i++ ) *d++ = b[i];
    }
    else if ( na < nb ) {
        *nratout = (2*nb+1)*s;
        for ( i = 0; i < na; i++ ) *d++ = a[i];
        for ( ; i < nb; i++ ) *d++ = 0;
        for ( i = 0; i < nb; i++ ) *d++ = b[i];
    }
    else {
        *nratout = (2*na+1)*s;
        for ( i = 0; i < na; i++ ) *d++ = a[i];
        for ( i = 0; i < nb; i++ ) *d++ = b[i];
        for ( ; i < na; i++ ) *d++ = 0;
    }
    *d = (UWORD)(*nratout);
    NumberFree(b,"FloatToRat");
    NumberFree(a,"FloatToRat");
    AT.pWorkPointer = oldpWorkPointer;
/*
    Just for check we go back to float
*/
    if ( s < 0 ) mpf_neg(floatin,floatin);
/*
    {
        WORD n = *nratout;
        RatToFloat(aux1,ratout,n);
        mpf_sub(aux2,floatin,aux1);
        gmp_printf("Diff is %.*Fe\n",40,aux2);
    }
*/
ret:
    return(retval);
}
 
/*
        #] FloatToRat : 
        #[ ZeroTable :
*/
        
void ZeroTable(mpf_t *tab, int N)
{
    int i, j;
    for ( i = 0; i < N; i++ ) {
        for ( j = 0; j < tab[i]->_mp_prec; j++ ) tab[i]->_mp_d[j] = 0;
    }
}
            
/*
        #] ZeroTable : 
        #[ ReadFloat :
 
    Used by the compiler. It reads a floating point number and
    outputs it at AT.WorkPointer as if it were a float_ function
    with its arguments.
    The call enters with s[-1] == TFLOAT.
*/
 
SBYTE *ReadFloat(SBYTE *s)
{
    GETIDENTITY
    SBYTE *ss, c;
    ss = s;
    while ( *ss > 0 ) ss++;
    c = *ss; *ss = 0;
    gmp_sscanf((char *)s,"%Ff",aux1);
    if ( AT.WorkPointer > AT.WorkTop ) {
        MLOCK(ErrorMessageLock);
        MesWork();
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
    PackFloat(AT.WorkPointer,aux1);
    *ss = c;
    return(ss);
}
 
/*
        #] ReadFloat : 
        #[ CheckFloat :
 
    For the tokenizer. Tests whether the input string can be a float.
*/
 
UBYTE *CheckFloat(UBYTE *ss, int *spec)
{
    GETIDENTITY
    UBYTE *s = ss;
    int gotdot = 0;
    /* A single dot is not a valid float */
    if ( *s == '.' && FG.cTable[s[-1]] != 1 && FG.cTable[s[1]] != 1 ) return(ss);
    while ( FG.cTable[s[-1]] == 1 ) s--;
    /* This cannot be a valid float */
    if ( *s != '.' && FG.cTable[*s] != 1 ) return(ss);
    *spec = 0;
    while ( FG.cTable[*s] == 1 ) s++;
    if ( *s == '.' ) {
        gotdot = 1;
        s++;
        while ( FG.cTable[*s] == 1 ) s++;
    }
/*
    Now we have had the mantissa part, which may be zero.
    Check for an exponent.
*/
    if ( *s == 'e' || *s == 'E' ) {
        s++;
        if ( *s == '-' || *s == '+' ) s++;
        if ( FG.cTable[*s] == 1 ) {
            s++;
            while ( FG.cTable[*s] == 1 ) s++;
        }
        else { return(ss); }
    }
    else if ( gotdot == 0 ) return(ss); /* No radix point and no exponent */
    if ( AT.aux_ == 0 ) { /* no floating point system */
        *spec = -1;
        return(s);
    }
    return(s);
}
 
/*
        #] CheckFloat : 
    #] Low Level : 
    #[ Float Routines :
        #[ SetFloatPrecision :
 
        Sets the default precision (in bits) of the floats and allocate
        buffer space for an output string.
        The buffer is used by PrintFloat (decimal output) and Strictrounding 
        (binary or decimal output), so it must accommodate the larger space
        requirement:
        exponent: up to 12 chars.
        mantissa: prec + a little bit extra
*/
 
int SetFloatPrecision(WORD prec)
{
    if ( prec <= 0 ) {
        MesPrint("&Illegal value for number of bits for floating point constants: %d",prec);
        return(-1);
    }
    else {
        AC.DefaultPrecision = prec;
        if ( AO.floatspace != 0 ) M_free(AO.floatspace,"floatspace");
        AO.floatsize = (prec+20)*sizeof(char);
        AO.floatspace = (UBYTE *)Malloc1(AO.floatsize,"floatspace");
        mpf_set_default_prec(prec);
        return(0);
    }
}
 
/*
        #] SetFloatPrecision : 
        #[ PrintFloat :
 
    Print the float_ function with its arguments as a floating point
    number with numdigits digits.
    If however we use rawfloat as a format option it prints it as a
    regular Form function and it should never come here because the
    print routines in sch.c should intercept it.
    The buffer AO.floatspace is allocated when the precision of the
    floats is set in the routine SetFloatPrecision.
*/
 
int PrintFloat(WORD *fun,int numdigits)
{
    GETIDENTITY
    UBYTE *s1, *s2;
    int n = 0;
    int prec = (AC.DefaultPrecision-AC.MaxWeight-1)*log10(2.0);
    if ( numdigits > prec || numdigits == 0 ) {
        numdigits = prec;
    }
/* 
    GMP's gmp_snprintf always prints a non-zero number before the decimal point, so we 
    ask for one digit less. 
*/
    if ( UnpackFloat(aux4,fun) == 0 )
        n = gmp_snprintf((char *)(AO.floatspace),AO.floatsize,"%.*Fe",numdigits-1,aux4);
    if ( n > 0 ) {
        int n1, n2 = n;
        s1 = AO.floatspace+n;
        while ( s1 > AO.floatspace && s1[-1] != 'e'
         && s1[-1] != 'E' && s1[-1] != '.' ) { s1--; n2--; }
        if ( s1 > AO.floatspace && s1[-1] != '.' ) {
            s1--; n2--;
            s2 = s1; n1 = n2;
            while ( s1[-1] == '0' ) { s1--; n1--; }
            if ( s1[-1] == '.' ) { s1++; n1++; }
            n -= (n2-n1);
            while ( n1 < n ) { *s1++ = *s2++; n1++; }
            *s1 = 0;
        }
        if ( AC.OutputMode == FORTRANMODE ) {
            s1 = AO.floatspace+n;
            while ( s1 > AO.floatspace && *s1 != 'e' && *s1 != 'E' ) { 
                s1--; 
            }
            if ( ( AO.DoubleFlag & 2 ) == 2 ) { *s1 = 'Q'; } /* Quadruple precision fortran */
            else if ( ( AO.DoubleFlag & 1 ) == 1 ) { *s1 = 'D'; } /* Double precision fortran */
            else { *s1 = 'E'; } /* Single precision fortran */
        }
        if ( AC.OutputMode == MATHEMATICAMODE ) {
            s1 = AO.floatspace+n;
            s2 = s1+2; *s2-- = 0;
            while ( s1 > AO.floatspace && *s1 != 'e' && *s1 != 'E' ) { 
                *s2-- = *s1--; 
            }
            *s2-- = '^'; *s2 = '*'; /* Replace 'e' by '^*' */
            n++;
        }
    }
    return(n);
}
 
/*
        #] PrintFloat : 
        #[ AddFloats :
*/
 
int AddFloats(PHEAD WORD *fun3, WORD *fun1, WORD *fun2)
{
    int retval = 0;
    if ( UnpackFloat(aux1,fun1) == 0 && UnpackFloat(aux2,fun2) == 0 ) {
        mpf_add(aux1,aux1,aux2);
        PackFloat(fun3,aux1);
    }
    else { retval = -1; }
    return(retval);
}
 
/*
        #] AddFloats : 
        #[ MulFloats :
*/
 
int MulFloats(PHEAD WORD *fun3, WORD *fun1, WORD *fun2)
{
    int retval = 0;
    if ( UnpackFloat(aux1,fun1) == 0 && UnpackFloat(aux2,fun2) == 0 ) {
        mpf_mul(aux1,aux1,aux2);
        PackFloat(fun3,aux1);
    }
    else { retval = -1; }
    return(retval);
}
 
/*
        #] MulFloats : 
        #[ DivFloats :
*/
 
int DivFloats(PHEAD WORD *fun3, WORD *fun1, WORD *fun2)
{
    int retval = 0;
    if ( UnpackFloat(aux1,fun1) == 0 && UnpackFloat(aux2,fun2) == 0 ) {
        mpf_div(aux1,aux1,aux2);
        PackFloat(fun3,aux1);
    }
    else { retval = -1; }
    return(retval);
}
 
/*
        #] DivFloats : 
        #[ AddRatToFloat :
 
        Note: this can be optimized, because often the rat is rather simple.
*/
 
int AddRatToFloat(PHEAD WORD *outfun, WORD *infun, UWORD *formrat, WORD nrat)
{
    int retval = 0;
    if ( UnpackFloat(aux2,infun) == 0 ) {
        RatToFloat(aux1,formrat,nrat);
        mpf_add(aux2,aux2,aux1);
        PackFloat(outfun,aux2);
    }
    else retval = -1;
    return(retval);
}
 
/*
        #] AddRatToFloat : 
        #[ MulRatToFloat :
 
        Note: this can be optimized, because often the rat is rather simple.
*/
 
int MulRatToFloat(PHEAD WORD *outfun, WORD *infun, UWORD *formrat, WORD nrat)
{
    int retval = 0;
    if ( UnpackFloat(aux2,infun) == 0 ) {
        RatToFloat(aux1,formrat,nrat);
        mpf_mul(aux2,aux2,aux1);
        PackFloat(outfun,aux2);
    }
    else retval = -1;
    return(retval);
}
 
/*
        #] MulRatToFloat : 
        #[ SetupMZVTables :
*/
 
void SetupMZVTables(void)
{
/*
    Sets up a table of N+1 mpf_t floats with variable precision.
    It assumes that each next element needs one bit less.
    Initializes all of them.
    We take some extra space, to have one limb overflow everywhere.
    Because the depth of the MZV's can get close to the weight
    and each deeper sum goes one higher, we make the tablesize a bit bigger.
    This may not be needed if we fiddle with the sum boundaries.
*/
#ifdef WITHPTHREADS
    int i, Nw, id, totnum;
    size_t N;
    mpf_t *a;
    Nw = AC.DefaultPrecision;
    N = (size_t)Nw;
    totnum = AM.totalnumberofthreads;
    for ( id = 0; id < totnum; id++ ) {
        AB[id]->T.mpf_tab1 = (void *)Malloc1((N+2)*sizeof(mpf_t),"mpftab1");
        a = (mpf_t *)AB[id]->T.mpf_tab1;
        for ( i = 0; i <=Nw; i++ ) {
/*
    As explained in the comment above, we could make this variable precision
    using mpf_init2.
*/
            mpf_init(a[i]);
        }
        AB[id]->T.mpf_tab2 = (void *)Malloc1((N+2)*sizeof(mpf_t),"mpftab2");
        a = (mpf_t *)AB[id]->T.mpf_tab2;
        for ( i = 0; i <=Nw; i++ ) {
            mpf_init(a[i]);
        }
    }
#else
    int i, Nw;
    size_t N;
    Nw = AC.DefaultPrecision;
    N = (size_t)Nw;
    AT.mpf_tab1 = (void *)Malloc1((N+2)*sizeof(mpf_t),"mpftab1");
    for ( i = 0; i <= Nw; i++ ) {
/*
    As explained in the comment above, we could make this variable precision
    using mpf_init2.
*/
        mpf_init(mpftab1[i]);
    }
    AT.mpf_tab2 = (void *)Malloc1((N+2)*sizeof(mpf_t),"mpftab2");
    for ( i = 0; i <= Nw; i++ ) {
        mpf_init(mpftab2[i]);
    }
#endif
    AS.delta_1 = (void *)Malloc1(sizeof(mpf_t),"delta1");
    mpf_init(mpfdelta1);
    SimpleDelta(mpfdelta1,1); /* this can speed up things. delta1 = ln(2) */
}
 
/*
        #] SetupMZVTables : 
        #[ SetupMPFTables :
    
        Allocates the aux variables
*/
 
void SetupMPFTables(void)
{
#ifdef WITHPTHREADS
    int id, totnum;
    mpf_t *a;
#ifdef WITHSORTBOTS
    totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads);
#endif
    for ( id = 0; id < totnum; id++ ) {
/*
        We work here with a[0] etc because the aux1 etc contain B which
        in the current routine would be AB[0] only
*/
        AB[id]->T.aux_ = (void *)Malloc1(sizeof(mpf_t)*8,"AB[id]->T.aux_");
        a = (mpf_t *)AB[id]->T.aux_;
        mpf_inits(a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7],(mpf_ptr)0);
        if ( AB[id]->T.indi1 ) M_free(AB[id]->T.indi1,"indi1");
        AB[id]->T.indi1 = (WORD *)Malloc1(sizeof(WORD)*AC.MaxWeight*2,"indi1");
        AB[id]->T.indi2 = AB[id]->T.indi1 + AC.MaxWeight;
    }
#else
    AT.aux_ = (void *)Malloc1(sizeof(mpf_t)*8,"AT.aux_");
    mpf_inits(aux1,aux2,aux3,aux4,aux5,auxjm,auxjjm,auxsum,(mpf_ptr)0);
    if ( AT.indi1 ) M_free(AT.indi1,"indi1");
    AT.indi1 = (WORD *)Malloc1(sizeof(WORD)*AC.MaxWeight*2,"indi1");
    AT.indi2 = AT.indi1 + AC.MaxWeight;
#endif
}
 
/*
        #] SetupMPFTables : 
        #[ ClearMZVTables :
*/
 
void ClearMZVTables(void)
{
#ifdef WITHPTHREADS
    int i, id, totnum;
    mpf_t *a;
    totnum = AM.totalnumberofthreads;
    for ( id = 0; id < totnum; id++ ) {
        if ( AB[id]->T.mpf_tab1 ) { 
/*
            We work here with a[0] etc because the aux1, mpftab1 etc contain B 
            which in the current routine would be AB[0] only
*/
            a = (mpf_t *)AB[id]->T.mpf_tab1;
            for ( i = 0; i <=AC.DefaultPrecision; i++ ) {
                mpf_clear(a[i]);
            }
            M_free(AB[id]->T.mpf_tab1,"mpftab1"); 
            AB[id]->T.mpf_tab1 = 0; 
        }
        if ( AB[id]->T.mpf_tab2 ) { 
            a = (mpf_t *)AB[id]->T.mpf_tab2;
            for ( i = 0; i <=AC.DefaultPrecision; i++ ) {
                mpf_clear(a[i]);
            }
            M_free(AB[id]->T.mpf_tab2,"mpftab2"); 
            AB[id]->T.mpf_tab2 = 0; 
        }
    }
#ifdef WITHSORTBOTS
    totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads);
#endif
    for ( id = 0; id < totnum; id++ ) {
        if ( AB[id]->T.aux_ ) { 
            a = (mpf_t *)AB[id]->T.aux_;
            mpf_clears(a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7],(mpf_ptr)0);
            M_free(AB[id]->T.aux_,"AB[id]->T.aux_");
            AB[id]->T.aux_ = 0; 
        }
        if ( AB[id]->T.indi1 ) { M_free(AB[id]->T.indi1,"indi1"); AB[id]->T.indi1 = 0; }
    }
#else
    int i;
    if ( AT.mpf_tab1 ) { 
        for ( i = 0; i <= AC.DefaultPrecision; i++ ) {
            mpf_clear(mpftab1[i]);
        }
        M_free(AT.mpf_tab1,"mpftab1"); 
        AT.mpf_tab1 = 0; 
    }
    if ( AT.mpf_tab2 ) { 
        for ( i = 0; i <= AC.DefaultPrecision; i++ ) {
            mpf_clear(mpftab2[i]);
        }
        M_free(AT.mpf_tab2,"mpftab2"); 
        AT.mpf_tab2 = 0; 
    }
    if ( AT.aux_ ) { 
        mpf_clears(aux1,aux2,aux3,aux4,aux5,auxjm,auxjjm,auxsum,(mpf_ptr)0);
        M_free(AT.aux_,"AT.aux_"); 
        AT.aux_ = 0; 
    }
    if ( AT.indi1 ) { M_free(AT.indi1,"indi1"); AT.indi1 = 0; }
#endif
    if ( AO.floatspace ) { M_free(AO.floatspace,"floatspace"); AO.floatspace = 0;
        AO.floatsize = 0; }
    if ( AS.delta_1 ) { 
        mpf_clear(mpfdelta1);
        M_free(AS.delta_1,"delta1"); 
        AS.delta_1 = 0; 
    }
}
 
/*
        #] ClearMZVTables : 
        #[ CoToFloat :
*/
 
int CoToFloat(UBYTE *s)
{
    GETIDENTITY
    if ( AT.aux_ == 0 ) {
        MesPrint("&Illegal attempt to convert to float_ without activating floating point numbers.");
        MesPrint("&Forgotten %#startfloat instruction?");
        return(1);
    }
    while ( *s == ' ' || *s == ',' || *s == '\t' ) s++;
    if ( *s ) {
        MesPrint("&Illegal argument(s) in ToFloat statement: '%s'",s);
        return(1);
    }
    Add2Com(TYPETOFLOAT);
    return(0);
}
 
/*
        #] CoToFloat : 
        #[ CoToRat :
*/
 
int CoToRat(UBYTE *s)
{
    GETIDENTITY
    if ( AT.aux_ == 0 ) {
        MesPrint("&Illegal attempt to convert from float_ without activating floating point numbers.");
        MesPrint("&Forgotten %#startfloat instruction?");
        return(1);
    }
    while ( *s == ' ' || *s == ',' || *s == '\t' ) s++;
    if ( *s ) {
        MesPrint("&Illegal argument(s) in ToRat statement: '%s'",s);
        return(1);
    }
    Add2Com(TYPETORAT);
    return(0);
}
 
/*
        #] CoToRat : 
        #[ CoStrictRounding : 
 
        Syntax: StrictRounding [precision][base]
        - precision: number of digits to round to (optional)
        - base: 'd' for decimal (base 10) or 'b' for binary (base 2)
        
        If no arguments are provided, uses default precision with binary base.
*/
int CoStrictRounding(UBYTE *s)
{
    GETIDENTITY
    WORD x;
    int base;
    if ( AT.aux_ == 0 ) {
        MesPrint("&Illegal attempt for strict rounding without activating floating point numbers.");
        MesPrint("&Forgotten %#startfloat instruction?");
        return(1);
    }
    while ( *s == ' ' || *s == ',' || *s == '\t' ) s++;
    if ( *s == 0 ) {
        /* No subkey, which means round to default precision */
        x = AC.DefaultPrecision - AC.MaxWeight - 1;
        base = 2;
    }
    else if ( FG.cTable[*s] == 1 ) { /* number */
        ParseNumber(x,s)
        if ( tolower(*s) == 'd' ) { base = 10; s++; }      /* decimal base */
        else if ( tolower(*s) == 'b' ){ base = 2; s++; }  /* binary base */
        else goto IllPar;  /* invalid base specification */
    }
    else {
        goto IllPar;
    }
    while ( *s == ' ' || *s == ',' || *s == '\t' ) s++;
    
    /* Check for invalid arguments */
    if ( *s ) {
IllPar:
        MesPrint("&Illegal argument(s) in StrictRounding statement: '%s'",s);
        return(1);
    }
    Add4Com(TYPESTRICTROUNDING,x,base);
    return(0);
} 
/*
        #] CoStrictRounding : 
        #[ CoChop :
 
        LHS notation of the chop statement: 
            TYPECHOP, length, FLOATFUN, ... 
        where FLOATFUN, ... represents the threshold of the chop statement in 
        the notation of a float_ function with its arguments. 
 
*/
 
int CoChop(UBYTE *s)
{
    GETIDENTITY
    UBYTE *ss, c;
    WORD *w, *OldWork;
    int spec, pow = 1;
    unsigned long x;
    if ( AT.aux_ == 0 ) {
        MesPrint("&Illegal attempt to chop a float_ without activating floating point numbers.");
        MesPrint("&Forgotten %#startfloat instruction?");
        return(1);
    }
    if ( *s == 0 ) {
        MesPrint("&Chop needs a number (float, rational or power) as an argument.");
        return(1);
    }
    /* Create TYPECHOP header */
    w = OldWork = AT.WorkPointer;
    *w++ = TYPECHOP; 
    w++;
 
    while ( *s == ' ' || *s == ',' || *s == '\t' ) s++;
    
    /* 
        The argument of chop can be 
        1: a floating-point number
        2: an integer, rational number or power
     */
    if ( FG.cTable[*s] == 1 || *s == '.' ) { 
        /* 1: Attempt to parse as floating-point number */
        ss = CheckFloat(s, &spec);
        if ( ss > s ) {
            /* CheckFloat found a valid float */
            AT.WorkPointer = w;
            /* 
                Reads the floating point number and outputs it at AT.WorkPointer as if it were a float_ 
                function with its arguments.
             */
            ReadFloat((SBYTE *)s); 
            s = ss;
            w += w[1];
        }
        else {
            /* 2: CheckFloat didn't find a float, we now try for rationals and powers */
            /* Parse the integer part (numerator for rationals) */
            if ( FG.cTable[*s] == 1 ) {
                ParseNumber(x,s)
                mpf_set_ui(aux1,x);
            }
            while ( *s == ' ' || *s == '\t' ) s++;
            /* Check for rational number or power*/
            if ( *s == '/' || *s == '^' ) {
                c = *s; s++; 
                while ( *s == ' ' || *s == '\t' ) s++;
                if ( *s == '-' ) { s++; pow = -1; } /* negative power */
                /* Parse the denominator or power */
                if ( FG.cTable[*s] == 1 ) {
                    ParseNumber(x,s)
                    if ( c == '/' ) { /* rational */
                        if ( x == 0 ) {
                            MesPrint("&Division by zero in chop statement.");
                            return(1);
                        }
                        /* Perform the division */
                        mpf_div_ui(aux1, aux1,x);
                    }
                    else { /* Power */
                        mpf_pow_ui(aux1,aux1,x);
                        if ( pow == -1 ) {
                            mpf_ui_div(aux1,(unsigned long) 1, aux1);
                        }
                    }
                }
            }
            /* Put aux1 in the notation of a float_ function */
            PackFloat(w, aux1);
            w += w[1];
        }
    }
    if ( *s ) {
        MesPrint("&Illegal argument(s) in Chop statement: '%s'.",s);
        return(1);
    }
    AT.WorkPointer = OldWork;
    AT.WorkPointer[1] = w - AT.WorkPointer;  /* Set total length */
    AddNtoL(AT.WorkPointer[1],AT.WorkPointer); /* Add the LHS to the compiler buffer */
    return(0);
}
 
/*
        #] CoChop : 
        #[ ToFloat :
 
        Converts the coefficient to floating point if it is still a rat.
*/
 
int ToFloat(PHEAD WORD *term, WORD level)
{
    WORD *t, *tstop, nsize, nsign = 3;
    t = term+*term;
    nsize = ABS(t[-1]);
    tstop = t - nsize;
    if ( t[-1] < 0 ) { nsign = -nsign; }
    if ( nsize == 3 && t[-2] == 1 && t[-3] == 1 ) { /* Could be float */
        t = term + 1;
        while ( t < tstop ) {
            if ( ( *t == FLOATFUN ) && ( t+t[1] == tstop ) ) {
                return(Generator(BHEAD term,level));
            }
            t += t[1];
        }
    }
    RatToFloat(aux4,(UWORD *)tstop,nsize);
    PackFloat(tstop,aux4);
    tstop += tstop[1];
    *tstop++ = 1; *tstop++ = 1; *tstop++ = nsign;
    *term = tstop - term;
    AT.WorkPointer = tstop;
    return(Generator(BHEAD term,level));
}
 
/*
        #] ToFloat : 
        #[ ToRat :
 
        Tries to convert the coefficient to rational if it is still a float.
*/
 
int ToRat(PHEAD WORD *term, WORD level)
{
    WORD *tstop, *t, nsize, nsign, ncoef;
/*
    1: find the float which should be at the end.
*/
    tstop = term + *term; nsize = ABS(tstop[-1]);
    nsign = tstop[-1] < 0 ? -1: 1; tstop -= nsize;
    t = term+1;
    while ( t < tstop ) {
        if ( *t == FLOATFUN && t + t[1] == tstop && TestFloat(t) &&
        nsize == 3 && tstop[0] == 1 && tstop[1] == 1 ) break;
        t += t[1];
    }
    if ( t < tstop ) {
/*
        Now t points at the float_ function and everything is correct.
        The result can go straight over the float_ function.
*/
        UnpackFloat(aux4,t);
        // If aux4 is zero, the term vanishes
        if ( mpf_sgn(aux4) == 0 ) return(0);
        if ( FloatToRat((UWORD *)t,&ncoef,aux4) == 0 ) {
            // Check if the resulting rational is zero
            if ( t[0] == 0 && t[1] == 1 && ncoef == 3 ) return(0);
            t += ABS(ncoef);
            t[-1] = ncoef*nsign;
            *term = t - term;
        }
    }
    return(Generator(BHEAD term,level));
}
 
/*
        #] ToRat : 
        #[ StrictRounding : 
 
        Rounds floating point numbers to a specified precision
        in a given base (decimal or binary).
*/
int StrictRounding(PHEAD WORD *term, WORD level, WORD prec, WORD base) {
    WORD *t,*tstop;
    int sign,size,maxprec = AC.DefaultPrecision-AC.MaxWeight-1;
    /* maxprec is in bits */
    if ( base == 2 && prec > maxprec ) {
        prec = maxprec;
    }
    if ( base == 10 && prec > (int)(maxprec*log10(2.0)) ) {
        prec = maxprec*log10(2.0);
    }
    /* Find the float which should be at the end. */
    tstop = term + *term; size = ABS(tstop[-1]);
    sign = tstop[-1] < 0 ? -1: 1; tstop -= size;
    t = term+1;
    while ( t < tstop ) {
        if ( *t == FLOATFUN && t + t[1] == tstop && TestFloat(t) &&
        size == 3 && tstop[0] == 1 && tstop[1] == 1) {
            break;
        }
        t += t[1];
    }
    if ( t < tstop ) {
/*
        Now t points at the float_ function and everything is correct.
        The result can go straight over the float_ function.
*/
        char *s;
        mp_exp_t exp;
        /* Extract the floating point value */
        UnpackFloat(aux4,t);
        /* Convert to string: 
           - Format as MeN with M the mantissa and N the exponent 
           - the generated string by mpf_get_str is the fraction/mantissa with 
           an implicit radix point immediately to the left of the first digit. 
           The applicable exponent is written in exp. */
        s = (char *)AO.floatspace;
        *s++ = '.';
        mpf_get_str(s,&exp, base, prec, aux4);
        while ( *s != 0 ) s++;
        *s++ = 'e';
        snprintf(s,AO.floatsize-(s-(char *)AO.floatspace),"%ld",exp);
        /* Negative base values are used to specify that the exponent is in decimal */
        mpf_set_str(aux4,(char *)AO.floatspace,-base);
        /* Pack the rounded floating point value back into the term */
        PackFloat(t,aux4);
        t+=t[1];
        *t++ = 1; *t++ = 1; *t++ = 3*sign;
        *term = t - term;
    }
    return(Generator(BHEAD term,level));
}
/*
        #] StrictRounding : 
        #[ Chop :
 
    Removes terms with a floating point number smaller than a given threshold. 
 
    Search for a FLOATFUN and compares its absolute value against the threshold 
    specified in the chop statement. This threshold can be obtained from the 
    LHS of the chop statement in the compiler buffer.
 
 */
int Chop(PHEAD WORD *term, WORD level)
{
    WORD *tstop, *t, nsize, *threshold; 
    CBUF *C = cbuf+AM.rbufnum;
    /* Find the float which should be at the end. */
    tstop = term + *term; 
    nsize = ABS(tstop[-1]); tstop -= nsize;
    t = term+1;
    while ( t < tstop ) {
        if ( *t == FLOATFUN && t + t[1] == tstop && TestFloat(t) &&
        nsize == 3 && tstop[0] == 1 && tstop[1] == 1 ) break;
        t += t[1];
    }
    if ( t < tstop ) {
        /* Get threshold value from compiler buffer */
        threshold = C->lhs[level];
        threshold += 2;  /* Skip TYPECHOP header */
        UnpackFloat(aux5, threshold); 
        
        /* Extract float and compute its absolute value */
        UnpackFloat(aux4, t);
        mpf_abs(aux4, aux4);
        
        /* Remove if < threshold */
        if ( mpf_cmp(aux4, aux5) < 0 ) return(0);
    }
    return(Generator(BHEAD term,level));
}
 
/*
        #] Chop : 
    #] Float Routines : 
    #[ Sorting :
 
        We need a method to see whether two terms need addition that could
        involve floating point numbers.
        1: if PolyWise is active, we do not need anything, because
           Poly(Rat)Fun and floating point are mutually exclusive.
        2: This means that there should be only interference in AddCoef.
           and the AddRat parts in MergePatches, MasterMerge, SortBotMerge
           and PF_GetLoser.
        The compare routine(s) should recognise the float_ and give off
        a code (SortFloatMode) inside the AT struct:
        0: no float_
        1: float_ in term1 only
        2: float_ in term2 only
        3: float_ in both terms
        Be careful: we have several compare routines, because various codes
        use their own routines and we just set a variable with its address.
        Currently we have Compare1, CompareSymbols and CompareHSymbols.
        Only Compare1 should be active for this. The others should make sure
        that the proper variable is always zero.
        Remember: the sign of the coefficient is in the last word of the term.
 
        We need two routines:
        1: AddWithFloat for SplitMerge which sorts by pointer.
        2: MergeWithFloat for MergePatches etc which keeps terms as much
           as possible in their location.
        The routines start out the same, because AT.SortFloatMode, set in
        Compare1, tells more or less what should be done.
        The difference is in where we leave the result.
 
        In the future we may want to add a feature that makes the result
        zero if the sum comes within a certain distance of the numerical
        accuracy, like for instance 5% of the number of bits.
 
        #[ AddWithFloat :
*/
 
int AddWithFloat(PHEAD WORD **ps1, WORD **ps2)
{
    GETBIDENTITY
    SORTING *S = AT.SS;
    WORD *coef1, *coef2, size1, size2, *fun1, *fun2, *fun3;
    WORD *s1,*s2,sign3,j,jj, *t1, *t2, i;
    s1 = *ps1;
    s2 = *ps2;
    coef1 = s1+*s1; size1 = coef1[-1]; coef1 -= ABS(coef1[-1]);
    coef2 = s2+*s2; size2 = coef2[-1]; coef2 -= ABS(coef2[-1]);
    if ( AT.SortFloatMode == 3 ) {
        fun1 = s1+1; while ( fun1 < coef1 && fun1[0] != FLOATFUN ) fun1 += fun1[1];
        fun2 = s2+1; while ( fun2 < coef2 && fun2[0] != FLOATFUN ) fun2 += fun2[1];
        UnpackFloat(aux1,fun1);
        if ( size1 < 0 ) mpf_neg(aux1,aux1);
        UnpackFloat(aux2,fun2);
        if ( size2 < 0 ) mpf_neg(aux2,aux2);
    }
    else if ( AT.SortFloatMode == 1 ) {
        fun1 = s1+1; while ( fun1 < coef1 && fun1[0] != FLOATFUN ) fun1 += fun1[1];
        UnpackFloat(aux1,fun1);
        if ( size1 < 0 ) mpf_neg(aux1,aux1);
        fun2 = coef2;
        RatToFloat(aux2,(UWORD *)coef2,size2);
    }
    else if ( AT.SortFloatMode == 2 ) {
        fun2 = s2+1; while ( fun2 < coef2 && fun2[0] != FLOATFUN ) fun2 += fun2[1];
        fun1 = coef1;
        RatToFloat(aux1,(UWORD *)coef1,size1);
        UnpackFloat(aux2,fun2);
        if ( size2 < 0 ) mpf_neg(aux2,aux2);
    }
    else {
        MLOCK(ErrorMessageLock);
        MesPrint("Illegal value %d for AT.SortFloatMode in AddWithFloat.",AT.SortFloatMode);
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        return(0);
    }
    mpf_add(aux3,aux1,aux2);
    sign3 = mpf_sgn(aux3);
    if ( sign3 < 0 ) mpf_neg(aux3,aux3);
    fun3 = TermMalloc("AddWithFloat");
    PackFloat(fun3,aux3);
/*
    Now we have to calculate where the sumterm fits.
    If we are sloppy, we can be faster, but run the risk to need the
    extension space, even when it is not needed. At slightly lower speed
    (ie first creating the result in scribble space) we are better off.
    This is why we use TermMalloc.
 
    The new term will have a rational coefficient of 1,1,+-3.
    The size will be (fun1 or fun2 - term) + fun3 +3;
*/
    if ( AT.SortFloatMode == 3 ) {
        if ( fun1[1] == fun3[1]  ) {
Over1:
            i = fun3[1]; t1 = fun1; t2 = fun3; NCOPY(t1,t2,i);
            *t1++ = 1; *t1++ = 1; *t1++ = sign3 < 0 ? -3: 3;
            *s1 = t1-s1; goto Finished;
        }
        else if ( fun2[1] == fun3[1]  ) {
Over2:
            i = fun3[1]; t1 = fun2; t2 = fun3; NCOPY(t1,t2,i);
            *t1++ = 1; *t1++ = 1; *t1++ = sign3 < 0 ? -3: 3;
            *s2 = t1-s2; *ps1 = s2; goto Finished;
        }
        else if ( fun1[1] >= fun3[1]  ) goto Over1;
        else if ( fun2[1] >= fun3[1]  ) goto Over2;
    }
    else if ( AT.SortFloatMode == 1 ) {
        if ( fun1[1] >= fun3[1]  ) goto Over1;
        else if ( fun3[1]+3 <= ABS(size2) ) goto Over2;
    }
    else if ( AT.SortFloatMode == 2 ) {
        if ( fun2[1] >= fun3[1]  ) goto Over2;
        else if ( fun3[1]+3 <= ABS(size1) ) goto Over1;
    }
/*
    Does not fit. Go to extension space.
*/
    jj = fun1-s1;
    j = jj+fun3[1]+3; /* space needed */
    if ( (S->sFill + j) >= S->sTop2 ) {
        GarbHand();
        s1 = *ps1; /* new position of the term after the garbage collection */
        fun1 = s1+jj;
    }
    t1 = S->sFill;
    for ( i = 0; i < jj; i++ ) *t1++ = s1[i];
    i = fun3[1]; s1 = fun3; NCOPY(t1,s1,i);
    *t1++ = 1; *t1++ = 1; *t1++ = sign3 < 0 ? -3: 3;
    *ps1 = S->sFill;
    **ps1 = t1-*ps1;
    S->sFill = t1;
Finished:
    *ps2 = 0;
    TermFree(fun3,"AddWithFloat");
    AT.SortFloatMode = 0;
    if ( **ps1 > AM.MaxTer/((LONG)(sizeof(WORD))) ) {
        MLOCK(ErrorMessageLock);
        MesPrint("Term too complex after addition in sort. MaxTermSize = %10l",
        AM.MaxTer/sizeof(WORD));
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
    return(1);
}
 
/*
        #] AddWithFloat : 
        #[ MergeWithFloat :
 
        Note that we always park the result on top of term1.
        This makes life easier, because the original AddRat in MergePatches
        does this as well.
*/
 
int MergeWithFloat(PHEAD WORD **interm1, WORD **interm2)
{
    GETBIDENTITY
    WORD *coef1, *coef2, size1, size2, *fun1, *fun2, *fun3, *tt;
    WORD sign3,jj, *t1, *t2, i, *term1 = *interm1, *term2 = *interm2;
    int retval = 0;
    coef1 = term1+*term1; size1 = coef1[-1]; coef1 -= ABS(size1);
    coef2 = term2+*term2; size2 = coef2[-1]; coef2 -= ABS(size2);
    if ( AT.SortFloatMode == 3 ) {
        fun1 = term1+1; while ( fun1 < coef1 && fun1[0] != FLOATFUN ) fun1 += fun1[1];
        fun2 = term2+1; while ( fun2 < coef2 && fun2[0] != FLOATFUN ) fun2 += fun2[1];
        UnpackFloat(aux1,fun1);
        if ( size1 < 0 ) mpf_neg(aux1,aux1);
        UnpackFloat(aux2,fun2);
        if ( size2 < 0 ) mpf_neg(aux2,aux2);
    }
    else if ( AT.SortFloatMode == 1 ) {
        fun1 = term1+1; while ( fun1 < coef1 && fun1[0] != FLOATFUN ) fun1 += fun1[1];
        UnpackFloat(aux1,fun1);
        if ( size1 < 0 ) mpf_neg(aux1,aux1);
        fun2 = coef2;
        RatToFloat(aux2,(UWORD *)coef2,size2);
    }
    else if ( AT.SortFloatMode == 2 ) {
        fun2 = term2+1; while ( fun2 < coef2 && fun2[0] != FLOATFUN ) fun2 += fun2[1];
        fun1 = coef1;
        RatToFloat(aux1,(UWORD *)coef1,size1);
        UnpackFloat(aux2,fun2);
        if ( size2 < 0 ) mpf_neg(aux2,aux2);
    }
    else {
        MLOCK(ErrorMessageLock);
        MesPrint("Illegal value %d for AT.SortFloatMode in MergeWithFloat.",AT.SortFloatMode);
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        return(0);
    }
    mpf_add(aux3,aux1,aux2);
    sign3 = mpf_sgn(aux3);
/*
    Now check whether we can park the result on top of one of the input terms.
*/
    if ( sign3 < 0 ) mpf_neg(aux3,aux3);
    fun3 = TermMalloc("MergeWithFloat");
    PackFloat(fun3,aux3);
    if ( AT.SortFloatMode == 3 ) {
        if ( fun1[1] + ABS(size1) == fun3[1] + 3 ) {
OnTopOf1:   t1 = fun3; t2 = fun1;
            for ( i = 0; i < fun3[1]; i++ ) *t2++ = *t1++;
            *t2++ = 1; *t2++ = 1; *t2++ = sign3 < 0 ? -3: 3;
            retval = 1;
        }
        else if ( fun1[1] + ABS(size1) > fun3[1] + 3 ) {
Shift1:     t2 = term1 + *term1; tt = t2;
            *--t2 = sign3 < 0 ? -3: 3; *--t2 = 1; *--t2 = 1;
            t1 = fun3 + fun3[1]; for ( i = 0; i < fun3[1]; i++ ) *--t2 = *--t1;
            t1 = fun1;
            while ( t1 > term1 ) *--t2 = *--t1;
            *t2 = tt-t2; term1 = t2;
            retval = 1;
        }
        else { /* Here we have to move term1 to the left to make room. */
            jj = fun3[1]-fun1[1]+3-ABS(size1); /* This is positive */
Over1:      t2 = term1-jj; t1 = term1;
            while ( t1 < fun1 ) *t2++ = *t1++;
            term1 -= jj;
            *term1 += jj;
            for ( i = 0; i < fun3[1]; i++ ) *t2++ = fun3[i];
            *t2++ = 1; *t2++ = 1;  *t2++ = sign3 < 0 ? -3: 3;
            retval = 1;
        }
    }
    else if ( AT.SortFloatMode == 1 ) {
        if ( fun1[1] + ABS(size1) == fun3[1] + 3 ) goto OnTopOf1;
        else if ( fun1[1] + ABS(size1) > fun3[1] + 3 ) goto Shift1;
        else {
            jj = fun3[1]-fun1[1]+3-ABS(size1); /* This is positive */
            goto Over1;
        }
    }
    else { /* Can only be 2, based on previous tests */
        if ( fun3[1] + 3 == ABS(size1) ) goto OnTopOf1;
        else if ( fun3[1] + 3 < ABS(size1) ) goto Shift1;
        else {
            jj = fun3[1]+3-ABS(size1); /* This is positive */
            goto Over1;
        }
    }
    *interm1 = term1;
    TermFree(fun3,"MergeWithFloat");
    AT.SortFloatMode = 0;
    return(retval);
}
 
/*
        #] MergeWithFloat : 
    #] Sorting : 
    #[ MZV :
 
        The functions here allow the arbitrary precision calculation
        of MZV's and Euler sums.
        They are called when the functions mzv_ and/or euler_ are used
        and the evaluate statement is applied.
        The output is in a float_ function.
        The expand statement tries to express the functions in terms of a basis.
        The bases are the 'standard basis for the euler sums and the
        pushdown bases from the datamine, unless otherwise specified.
 
        #[ SimpleDelta :
*/
 
void SimpleDelta(mpf_t sum, int m)
{
    long s = 1;
    long prec = AC.DefaultPrecision;
    unsigned long xprec = (unsigned long)prec, x;
    int j, jmax, n;
    mpf_t jm,jjm;
    mpf_init(jm); mpf_init(jjm);
    if ( m < 0 ) { s = -1; m = -m; }
 
/*
    We will loop until 1/2^j/j^m is smaller than the default precision.
    Just running to prec is however overkill, specially when m is large.
    We try to estimate a better value.
    jmax = prec - ((log2(prec)-1)*m).
    Hence we need the leading bit of prec.
    We are still overshooting a bit.
*/
    n = 0; x = xprec;
    while ( x ) { x >>= 1; n++; }
/* 
    We have now n = floor(log2(x))+1. 
*/
    n--;
    jmax = (int)((int)xprec - (n-1)*m);
/*
    For small prec and large m, the estimate can be wrong and even be negative, 
    so we increase jmax until jmax + m*log2(jmax) > prec
*/
    if ( jmax < 0 ) jmax = 1;
    do {
        n = 0;
        x = (unsigned long)jmax;
        while (x) { x >>= 1; n++; }
        n--; // floor(log2(jmax))
        jmax++; 
    } while ( jmax + m * n <= prec );
    mpf_set_ui(sum,0);
    for ( j = 1; j <= jmax; j++ ) {
#ifdef WITHCUTOFF
        xprec--;
        mpf_set_prec_raw(jm,xprec);
        mpf_set_prec_raw(jjm,xprec);
#endif
        mpf_set_ui(jm,1L);
        mpf_div_ui(jjm,jm,(unsigned long)j);
        mpf_pow_ui(jm,jjm,(unsigned long)m);
        mpf_div_2exp(jjm,jm,(unsigned long)j);
        if ( s == -1 && j%2 == 1 ) mpf_sub(sum,sum,jjm);
        else                       mpf_add(sum,sum,jjm);
    }
    mpf_clear(jjm); mpf_clear(jm);
}
 
/*
        #] SimpleDelta : 
        #[ SimpleDeltaC :
*/
 
void SimpleDeltaC(mpf_t sum, int m)
{
    long s = 1;
    long prec = AC.DefaultPrecision;
    unsigned long xprec = (unsigned long)prec, x;
    int j, jmax, n;
    mpf_t jm,jjm;
    mpf_init(jm); mpf_init(jjm);
    if ( m < 0 ) { s = -1; m = -m; }
/*
    We will loop until 1/2^j/j^m is smaller than the default precision.
    Just running to prec is however overkill, specially when m is large.
    We try to estimate a better value.
    jmax = prec - ((log2(prec)-1)*m).
    Hence we need the leading bit of prec.
    We are still overshooting a bit.
*/
    n = 0; x = xprec;
    while ( x ) { x >>= 1; n++; }
/* 
    We have now n = floor(log2(x))+1. 
*/
    n--;
    jmax = (int)((int)xprec - (n-1)*m);
/*
    For small prec and large m, the estimate can be wrong and even be negative, 
    so we increase jmax until jmax + m*log2(jmax) > prec
*/
    if ( jmax < 0 ) jmax = 1;
    do {
        n = 0;
        x = (unsigned long)jmax;
        while (x) { x >>= 1; n++; }
        n--; // floor(log2(jmax))
        jmax++; 
    } while ( jmax + m * n <= prec );
    if ( s < 0 ) jmax /= 2;
    mpf_set_si(sum,0L);
    for ( j = 1; j <= jmax; j++ ) {
#ifdef WITHCUTOFF
        xprec--;
        mpf_set_prec_raw(jm,xprec);
        mpf_set_prec_raw(jjm,xprec);
#endif
        mpf_set_ui(jm,1L);
        mpf_div_ui(jjm,jm,(unsigned long)j);
        mpf_pow_ui(jm,jjm,(unsigned long)m);
        if ( s == -1 ) mpf_div_2exp(jjm,jm,2*(unsigned long)j);
        else           mpf_div_2exp(jjm,jm,(unsigned long)j);
        mpf_add(sum,sum,jjm);
    }
    mpf_clear(jjm); mpf_clear(jm);
}
 
/*
        #] SimpleDeltaC : 
        #[ SingleTable :
*/
 
void SingleTable(mpf_t *tabl, int N, int m, int pow)
{
/*
    Creates a table T(1,...,N) with
        T(i) = sum_(j,i,N,[sign_(j)]/2^j/j^m)
    To make this table we have two options:
    1: run the sum backwards with the potential problem that the 
       precision is difficult to manage.
    2: run the sum forwards. Take sum_(j,1,i-1,...) and later subtract.
    When doing Euler sums we may need also 1/4^j
    pow: 1 -> 1/2^j
         2 -> 1/4^j
*/
    GETIDENTITY
    long s = 1;
    int j;
#ifdef WITHCUTOFF
    long prec = mpf_get_default_prec();
#endif
    mpf_t jm,jjm;
    mpf_init(jm); mpf_init(jjm);
    if ( pow < 1 || pow > 2 ) {
        MLOCK(ErrorMessageLock);
        MesPrint("Wrong parameter pow in SingleTable: %d\n",pow);
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
    if ( m < 0 ) { m = -m; s = -1; }
    mpf_set_si(auxsum,0L);
    for ( j = N; j > 0; j-- ) {
#ifdef WITHCUTOFF
        mpf_set_prec_raw(jm,prec-j);
        mpf_set_prec_raw(jjm,prec-j);
#endif
        mpf_set_ui(jm,1L);
        mpf_div_ui(jjm,jm,(unsigned long)j);
        mpf_pow_ui(jm,jjm,(unsigned long)m);
        mpf_div_2exp(jjm,jm,pow*(unsigned long)j);
        if ( pow == 2 ) mpf_add(auxsum,auxsum,jjm);
        else if ( s == -1 && j%2 == 1 ) mpf_sub(auxsum,auxsum,jjm);
        else                       mpf_add(auxsum,auxsum,jjm);
/*
        And now copy auxsum to tablelement j
*/
        mpf_set(tabl[j],auxsum);
    }
    mpf_clear(jjm); mpf_clear(jm);
}
 
/*
        #] SingleTable : 
        #[ DoubleTable :
*/
 
void DoubleTable(mpf_t *tabout, mpf_t *tabin, int N, int m, int pow)
{
    GETIDENTITY
    long s = 1;
#ifdef WITHCUTOFF
    long prec = mpf_get_default_prec();
#endif
    int j;
    mpf_t jm,jjm;
    mpf_init(jm); mpf_init(jjm);
    if ( pow < -1 || pow > 2 ) {
        MLOCK(ErrorMessageLock);
        MesPrint("Wrong parameter pow in DoubleTable: %d\n",pow);
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
    if ( m < 0 ) { m = -m; s = -1; }
    mpf_set_ui(auxsum,0L);
    for ( j = N; j > 0; j-- ) {
#ifdef WITHCUTOFF
        mpf_set_prec_raw(jm,prec-j);
        mpf_set_prec_raw(jjm,prec-j);
#endif
        mpf_set_ui(jm,1L);
        mpf_div_ui(jjm,jm,(unsigned long)j);
        mpf_pow_ui(jm,jjm,(unsigned long)m);
        if ( pow == -1 ) {
            mpf_mul_2exp(jjm,jm,(unsigned long)j);
            mpf_mul(jm,jjm,tabin[j+1]);
        }
        else if ( pow > 0 ) {
            mpf_div_2exp(jjm,jm,pow*(unsigned long)j);
            mpf_mul(jm,jjm,tabin[j+1]);
        }
        else {
            mpf_mul(jm,jm,tabin[j+1]);
        }
        if ( pow == 2 ) mpf_add(auxsum,auxsum,jm);
        else if ( s == -1 && j%2 == 1 ) mpf_sub(auxsum,auxsum,jm);
        else                       mpf_add(auxsum,auxsum,jm);
/*
        And now copy auxsum to tablelement j
*/
        mpf_set(tabout[j],auxsum);
    }
    mpf_clear(jjm); mpf_clear(jm);
}
 
/*
        #] DoubleTable : 
        #[ EndTable :
*/
 
void EndTable(mpf_t sum, mpf_t *tabin, int N, int m, int pow)
{
    long s = 1;
#ifdef WITHCUTOFF
    long prec = mpf_get_default_prec();
#endif
    int j;
    mpf_t jm,jjm;
    mpf_init(jm); mpf_init(jjm);
    if ( pow < -1 || pow > 2 ) {
        MLOCK(ErrorMessageLock);
        MesPrint("Wrong parameter pow in EndTable: %d\n",pow);
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
    if ( m < 0 ) { m = -m; s = -1; }
    mpf_set_si(sum,0L);
    for ( j = N; j > 0; j-- ) {
#ifdef WITHCUTOFF
        mpf_set_prec_raw(jm,prec-j);
        mpf_set_prec_raw(jjm,prec-j);
#endif
        mpf_set_ui(jm,1L);
        mpf_div_ui(jjm,jm,(unsigned long)j);
        mpf_pow_ui(jm,jjm,(unsigned long)m);
        if ( pow == -1 ) {
            mpf_mul_2exp(jjm,jm,(unsigned long)j);
            mpf_mul(jm,jjm,tabin[j+1]);
        }
        else if ( pow > 0 ) {
            mpf_div_2exp(jjm,jm,pow*(unsigned long)j);
            mpf_mul(jm,jjm,tabin[j+1]);
        }
        else {
            mpf_mul(jm,jm,tabin[j+1]);
        }
        if ( s == -1 && j%2 == 1 ) mpf_sub(sum,sum,jm);
        else                       mpf_add(sum,sum,jm);
    }
    mpf_clear(jjm); mpf_clear(jm);
}
 
/*
        #] EndTable : 
        #[ deltaMZV :
*/
 
void deltaMZV(mpf_t result, WORD *indexes, int depth)
{
    GETIDENTITY
/*
    Because all sums go roughly like 1/2^j we need about depth steps.
*/
/*  MesPrint("deltaMZV(%a)",depth,indexes); */
    if ( depth == 1 ) {
        if ( indexes[0] == 1 ) {
            mpf_set(result,mpfdelta1);
            return;
        }
        SimpleDelta(result,indexes[0]);
    }
    else if ( depth == 2 ) {
        if ( indexes[0] == 1 && indexes[1] == 1 ) {
            mpf_pow_ui(result,mpfdelta1,2);
            mpf_div_2exp(result,result,1);
        }
        else {
            SingleTable(mpftab1,AC.DefaultPrecision-AC.MaxWeight+1,indexes[0],1);
            EndTable(result,mpftab1,AC.DefaultPrecision-AC.MaxWeight,indexes[1],0);
        };
    }
    else if ( depth > 2 ) {
        mpf_t *mpftab3;
        int d;
/*
        Check first whether this is a power of delta1 = ln(2)
*/
        for ( d = 0; d < depth; d++ ) {
            if ( indexes[d] != 1 ) break;
        }
        if ( d == depth ) { /* divide by fac_(depth) */
            mpf_pow_ui(result,mpfdelta1,depth);
            for ( d = 2; d <= depth; d++ ) {
                mpf_div_ui(result,result,(unsigned long)d);
            }
        }
        else {
            d = depth-1;
            SingleTable(mpftab1,AC.DefaultPrecision-AC.MaxWeight+d,*indexes,1);
            d--; indexes++;
            for(;;) {
                DoubleTable(mpftab2,mpftab1,AC.DefaultPrecision-AC.MaxWeight+d,*indexes,0);
                d--; indexes++;
                if ( d == 0 ) {
                    EndTable(result,mpftab2,AC.DefaultPrecision-AC.MaxWeight,*indexes,0);
                    break;
                }
                mpftab3 = (mpf_t *)AT.mpf_tab1; AT.mpf_tab1 = AT.mpf_tab2;
                AT.mpf_tab2 = (void *)mpftab3;
            }
        }
    }
    else if ( depth == 0 ) {
        mpf_set_ui(result,1L);
    }
}
 
/*
        #] deltaMZV : 
        #[ deltaEuler :
 
        Regular Euler delta with - signs, but everywhere 1/2^j
*/
 
void deltaEuler(mpf_t result, WORD *indexes, int depth)
{
    GETIDENTITY
    int m;
#ifdef DEBUG
    int i;
    printf(" deltaEuler(");
    for ( i = 0; i < depth; i++ ) {
        if ( i != 0 ) printf(",");
        printf("%d",indexes[i]);
    }
    printf(") = ");
#endif
    if ( depth == 1 ) {
        SimpleDelta(result,indexes[0]);
    }
    else if ( depth == 2 ) {
        SingleTable(mpftab1,AC.DefaultPrecision-AC.MaxWeight+1,indexes[0],1);
        m = indexes[1]; if ( indexes[0] < 0 ) m = -m;
        EndTable(result,mpftab1,AC.DefaultPrecision-AC.MaxWeight,m,0);
    }
    else if ( depth > 2 ) {
        int d;
        mpf_t *mpftab3;
        d = depth-1;
        SingleTable(mpftab1,AC.DefaultPrecision-AC.MaxWeight+d,*indexes,1);
        d--; indexes++;
        m = *indexes; if ( indexes[-1] < 0 ) m = -m;
        for(;;) {
            DoubleTable(mpftab2,mpftab1,AC.DefaultPrecision-AC.MaxWeight+d,m,0);
            d--; indexes++;
            m = *indexes; if ( indexes[-1] < 0 ) m = -m;
            if ( d == 0 ) {
                EndTable(result,mpftab2,AC.DefaultPrecision-AC.MaxWeight,m,0);
                break;
            }
            mpftab3 = (mpf_t *)AT.mpf_tab1; AT.mpf_tab1 = AT.mpf_tab2;
            AT.mpf_tab2 = (void *)mpftab3;
        }
    }
    else if ( depth == 0 ) {
        mpf_set_ui(result,1L);
    }
#ifdef DEBUG
    gmp_printf("%.*Fe\n",40,result);
#endif
}
 
/*
        #] deltaEuler : 
        #[ deltaEulerC :
 
        Conjugate Euler delta without - signs, but possibly 1/4^j
        When there is a - in the string we have 1/4.
*/
 
void deltaEulerC(mpf_t result, WORD *indexes, int depth)
{
    GETIDENTITY
    int m;
#ifdef DEBUG
    int i;
    printf(" deltaEulerC(");
    for ( i = 0; i < depth; i++ ) {
        if ( i != 0 ) printf(",");
        printf("%d",indexes[i]);
    }
    printf(") = ");
#endif
    mpf_set_ui(result,0);
    if ( depth == 1 ) {
        if ( indexes[0] == 0 ) {
            MLOCK(ErrorMessageLock);
            MesPrint("Illegal index in depth=1 deltaEulerC: %d\n",indexes[0]);
            MUNLOCK(ErrorMessageLock);
            Terminate(-1);
        }
        if ( indexes[0] < 0 ) SimpleDeltaC(result,indexes[0]);
        else                  SimpleDelta(result,indexes[0]);
    }
    else if ( depth == 2 ) {
        int par;
        m = indexes[0];
        if ( m < 0 ) SingleTable(mpftab1,AC.DefaultPrecision-AC.MaxWeight+depth,-m,2);
        else         SingleTable(mpftab1,AC.DefaultPrecision-AC.MaxWeight+depth, m,1);
        m = indexes[1];
        if ( m < 0 ) { m = -m; par = indexes[0] < 0 ? 0: 1; }
        else { par = indexes[0] < 0 ? -1: 0; }
        EndTable(result,mpftab1,AC.DefaultPrecision-AC.MaxWeight,m,par);
    }
    else if ( depth > 2 ) {
        int d, par;
        mpf_t *mpftab3;
        d = depth-1;
        m = indexes[0];
        ZeroTable(mpftab1,AC.DefaultPrecision+1);
        if ( m < 0 ) SingleTable(mpftab1,AC.DefaultPrecision-AC.MaxWeight+d+1,-m,2);
        else         SingleTable(mpftab1,AC.DefaultPrecision-AC.MaxWeight+d+1, m,1);
        d--; indexes++; m = indexes[0];
        if ( m < 0 ) { m = -m; par = indexes[-1] < 0 ? 0: 1; }
        else { par = indexes[-1] < 0 ? -1: 0; }
        for(;;) {
            ZeroTable(mpftab2,AC.DefaultPrecision-AC.MaxWeight+d);
            DoubleTable(mpftab2,mpftab1,AC.DefaultPrecision-AC.MaxWeight+d,m,par);
            d--; indexes++; m = indexes[0];
            if ( m < 0 ) { m = -m; par = indexes[-1] < 0 ? 0: 1; }
            else { par = indexes[-1] < 0 ? -1: 0; }
            if ( d == 0 ) {
                mpf_set_ui(result,0);
                EndTable(result,mpftab2,AC.DefaultPrecision-AC.MaxWeight,m,par);
                break;
            }
            mpftab3 = (mpf_t *)AT.mpf_tab1; AT.mpf_tab1 = AT.mpf_tab2;
            AT.mpf_tab2 = (void *)mpftab3;
        }
    }
    else if ( depth == 0 ) {
        mpf_set_ui(result,1L);
    }
#ifdef DEBUG
    gmp_printf("%.*Fe\n",40,result);
#endif
}
 
/*
        #] deltaEulerC : 
        #[ CalculateMZVhalf :
 
        This routine is mainly for support when large numbers of
        MZV's have to be calculated at the same time.
*/
 
void CalculateMZVhalf(mpf_t result, WORD *indexes, int depth)
{
    int i;
    if ( depth < 0 ) {
        MLOCK(ErrorMessageLock);
        MesPrint("Illegal depth in CalculateMZVhalf: %d",depth);
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
    for ( i = 0; i < depth; i++ ) {
        if ( indexes[i] <= 0 ) {
            MLOCK(ErrorMessageLock);
            MesPrint("Illegal index[%d] in CalculateMZVhalf: %d",i,indexes[i]);
            MUNLOCK(ErrorMessageLock);
            Terminate(-1);
        }
    }
    deltaMZV(result,indexes,depth);
}
 
/*
        #] CalculateMZVhalf : 
        #[ CalculateMZV :
*/
 
void CalculateMZV(mpf_t result, WORD *indexes, int depth)
{
    GETIDENTITY
    int num1, num2 = 0, i, j = 0;
    if ( depth < 0 ) {
        MLOCK(ErrorMessageLock);
        MesPrint("Illegal depth in CalculateMZV: %d",depth);
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
    if ( indexes[0] == 1 ) {
        MLOCK(ErrorMessageLock);
        MesPrint("Divergent MZV in CalculateMZV");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
/*  MesPrint("calculateMZV(%a)",depth,indexes); */
    for ( i = 0; i < depth; i++ ) {
        if ( indexes[i] <= 0 ) {
            MLOCK(ErrorMessageLock);
            MesPrint("Illegal index[%d] in CalculateMZV: %d",i,indexes[i]);
            MUNLOCK(ErrorMessageLock);
            Terminate(-1);
        }
        AT.indi1[i] = indexes[i];
    }
    num1 = depth; i = 0;
    deltaMZV(result,indexes,depth);
    for(;;) {
        if ( AT.indi1[0] > 1 ) {
            AT.indi1[0] -= 1;
            for ( j = num2; j > 0; j-- ) AT.indi2[j] = AT.indi2[j-1];
            AT.indi2[0] = 1; num2++;
        }
        else {
            AT.indi2[0] += 1;
            for ( j = 1; j < num1; j++ ) AT.indi1[j-1] = AT.indi1[j];
            num1--;
            if ( num1 == 0 ) break;
        }
        deltaMZV(aux1,AT.indi1,num1);
        deltaMZV(aux2,AT.indi2,num2);
        mpf_mul(aux3,aux1,aux2);
        mpf_add(result,result,aux3);
    }
    deltaMZV(aux3,AT.indi2,num2);
    mpf_add(result,result,aux3);
}
 
/*
        #] CalculateMZV : 
        #[ CalculateEuler :
 
        There is a problem of notation here.
        Basically there are two notations.
        1: Z(m1,m2,m3) = sum_(j3,1,inf,sign_(m3)/j3^abs_(m3)*
                         sum_(j2,j3+1,inf,sign_(m2)/j2^abs_(m2)*
                         sum_(j1,j2+1,inf,sign_(m1)/j1^abs_(m1))))
           See Broadhurst,1996
        2: L(m1,m2,m3) = sum_(j1,1,inf,sign_(m1)*
                         sum_(j2,1,inf,sign_(m2)*
                         sum_(j3,1,inf,sign_(m3)
                            /j3^abs_(m3)
                            /(j2+j3)^abs_(m2)
                            /(j1+j2+j3)^abs_(m1) )))
           See Borwein, Bradley, Broadhurst and Lisonek, 1999
        The L corresponds with the H of the datamine up to sign_(m1*m2*m3)
        The algorithm follows 2, but the more common notation is 1.
        Hence we start with a conversion.
*/
 
void CalculateEuler(mpf_t result, WORD *Zindexes, int depth)
{
    GETIDENTITY
    int s1, num1, num2, i, j;
 
    WORD *indexes = AT.WorkPointer;
    for ( i = 0; i < depth; i++ ) indexes[i] = Zindexes[i];
    for ( i = 0; i < depth-1; i++ ) {
        if ( Zindexes[i] < 0 ) {
            for ( j = i+1; j < depth; j++ ) indexes[j] = -indexes[j];
        }
    }
 
    if ( depth < 0 ) {
        MLOCK(ErrorMessageLock);
        MesPrint("Illegal depth in CalculateEuler: %d\n",depth);
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
    if ( indexes[0] == 1 ) {
        MLOCK(ErrorMessageLock);
        MesPrint("Divergent Euler sum in CalculateEuler\n");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
    for ( i = 0, j = 0; i < depth; i++ ) {
        if ( indexes[i] == 0 ) {
            MLOCK(ErrorMessageLock);
            MesPrint("Illegal index[%d] in CalculateEuler: %d\n",i,indexes[i]);
            MUNLOCK(ErrorMessageLock);
            Terminate(-1);
        }
        if ( indexes[i] < 0 ) j = 1;
        AT.indi1[i] = indexes[i];
    }
    if ( j == 0 ) {
        CalculateMZV(result,indexes,depth);
        return;
    }
    num1 = depth; AT.indi2[0] = 0; num2 = 0;
    deltaEuler(result,AT.indi1,depth);
    s1 = 0;
    for(;;) {
        s1++;
        if ( AT.indi1[0] > 1 ) {
            AT.indi1[0] -= 1;
            for ( j = num2; j > 0; j-- ) AT.indi2[j] = AT.indi2[j-1];
            AT.indi2[0] = 1; num2++;
        }
        else if ( AT.indi1[0] < -1 ) {
            AT.indi1[0] += 1;
            for ( j = num2; j > 0; j-- ) AT.indi2[j] = AT.indi2[j-1];
            AT.indi2[0] = 1; num2++;
        }
        else if ( AT.indi1[0] == -1 ) {
            for ( j = num2; j > 0; j-- ) AT.indi2[j] = AT.indi2[j-1];
            AT.indi2[0] = -1; num2++;
            for ( j = 1; j < num1; j++ ) AT.indi1[j-1] = AT.indi1[j];
            num1--;
        }
        else {
            AT.indi2[0] = AT.indi2[0] < 0 ? AT.indi2[0]-1: AT.indi2[0]+1;
            for ( j = 1; j < num1; j++ ) AT.indi1[j-1] = AT.indi1[j];
            num1--;
        }
        if ( num1 == 0 ) break;
        deltaEuler(aux1,AT.indi1,num1);
        deltaEulerC(aux2,AT.indi2,num2);
        mpf_mul(aux3,aux1,aux2);
        if ( (depth+num1+num2+s1)%2 == 0 ) mpf_add(result,result,aux3);
        else                               mpf_sub(result,result,aux3);
    }
    deltaEulerC(aux3,AT.indi2,num2);
    if ( (depth+num2+s1)%2 == 0 ) mpf_add(result,result,aux3);
    else                          mpf_sub(result,result,aux3);
}
 
/*
        #] CalculateEuler : 
        #[ ExpandMZV :
*/
 
int ExpandMZV(WORD *term, WORD level)
{
    DUMMYUSE(term);
    DUMMYUSE(level);
    return(0);
}
 
/*
        #] ExpandMZV : 
        #[ ExpandEuler :
*/
 
int ExpandEuler(WORD *term, WORD level)
{
    DUMMYUSE(term);
    DUMMYUSE(level);
    return(0);
}
 
/*
        #] ExpandEuler : 
        #[ EvaluateEuler :
 
    There are various steps here:
    1: locate an Euler function.
    2: evaluate it into a float.
    3: convert the coefficient to a float and multiply.
    4: Put the new term together.
    5: Restart to see whether there are more Euler functions.
    The compiler should check that there is no conflict between
    the evaluate command and a potential polyfun or polyratfun.
    Yet, when reading in expressions there can be a conflict.
    Hence the Normalize routine should check as well.
 
    par == MZV: MZV
    par == EULER: Euler
    par == MZVHALF: MZV sums in 1/2 rather than 1. -> deltaMZV.
    par == ALLMZVFUNCTIONS: all of the above.
*/
 
int EvaluateEuler(PHEAD WORD *term, WORD level, WORD par)
{
    WORD *t, *tstop, *tt, *tnext, sumweight, depth, first = 1, *newterm, i;
    WORD *oldworkpointer = AT.WorkPointer, nsize, *indexes;
    int retval;
    tstop = term + *term; tstop -= ABS(tstop[-1]);
    if ( AT.WorkPointer < term+*term ) AT.WorkPointer = term + *term;
/*
    Step 1. We also need to verify that the argument we find is legal
*/
    t = term+1;
    while ( t < tstop ) {
        if ( ( *t == par ) || ( ( par == ALLMZVFUNCTIONS ) && (
            *t == MZV || *t == EULER || *t == MZVHALF ) ) ) {
    /* all arguments should be small integers */
            indexes = AT.WorkPointer;
            sumweight = 0; depth = 0;
            tnext = t + t[1]; tt = t + FUNHEAD;
            while ( tt < tnext ) {
                if ( *tt != -SNUMBER ) goto nextfun;
                if ( tt[1] == 0 ) goto nextfun;
                indexes[depth] = tt[1];
                sumweight += ABS(tt[1]); depth++;
                tt += 2;
            }
            /* euler sum without arguments, i.e. mzv_(), euler_() or mzvhalf_() */
            if ( depth == 0) goto nextfun;
            if ( sumweight > AC.MaxWeight ) {
                MLOCK(ErrorMessageLock);
                MesPrint("Error: Weight of Euler/MZV sum greater than %d.",AC.MaxWeight);
                MesPrint("Please increase the maximum weight in %#startfloat.");
                MUNLOCK(ErrorMessageLock);
                Terminate(-1);
            }
/*
            Step 2: evaluate:
*/
            AT.WorkPointer += depth;
            if ( first ) {
                if ( *t == MZV ) CalculateMZV(aux4,indexes,depth);
                else if ( *t == EULER ) CalculateEuler(aux4,indexes,depth);
                else if ( *t == MZVHALF ) CalculateMZVhalf(aux4,indexes,depth);
                first = 0;
            }
            else {
                if ( *t == MZV ) CalculateMZV(aux5,indexes,depth);
                else if ( *t == EULER ) CalculateEuler(aux5,indexes,depth);
                else if ( *t == MZVHALF ) CalculateMZVhalf(aux5,indexes,depth);
                mpf_mul(aux4,aux4,aux5);
            }
            *t = 0;
        }
nextfun:
        t += t[1];
    }
    if ( first == 1 ) return(Generator(BHEAD term,level));
/*
    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 ) {
        mpf_neg(aux4,aux4);
        nsize = -nsize;
    }
    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.
*/
    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;
    *newterm = tt-newterm;
    AT.WorkPointer = tt;
    retval = Generator(BHEAD newterm,level);
    AT.WorkPointer = oldworkpointer;
    return(retval);
}
 
/*
        #] EvaluateEuler : 
    #] MZV : 
    #[ Functions :
        #[ CoEvaluate :
 
        Current subkeys: mzv_, euler_ and sqrt_.
*/
 
int CoEvaluate(UBYTE *s)
{
    GETIDENTITY
    UBYTE *subkey, c;
    WORD numfun, type;
    int error = 0;
    if ( AT.aux_ == 0 ) {
        MesPrint("&Illegal attempt to evaluate a function without activating floating point numbers.");
        MesPrint("&Forgotten %#startfloat instruction?");
        return(1);
    }
    while ( *s == ' ' || *s == ',' || *s == '\t' ) s++;
    if ( *s == 0 ) {
/*
        No subkey, which means all functions.
*/
        Add3Com(TYPEEVALUATE,ALLFUNCTIONS);
/*
        The MZV, EULER and MZVHALF are done separately
*/
        Add3Com(TYPEEVALUATE,ALLMZVFUNCTIONS);
        return(0);  
    }
    while ( *s ) {
    subkey = s;
    while ( FG.cTable[*s] == 0 ) s++;
      if ( *s == '2' ) s++; /* cases li2_ and atan2_ */
      if ( *s == '_' ) s++;
      c = *s; *s = 0;
/*
        We still need provisions for pi_ and possibly other constants.
*/
      if ( ( ( type = GetName(AC.varnames,subkey,&numfun,NOAUTO) ) != CFUNCTION )
            || ( functions[numfun].spec != 0 ) ) {
 
        if ( type == CSYMBOL ) {
            Add4Com(TYPEEVALUATE,SYMBOL,numfun);
            break;
        }
/*
            This cannot work.
*/
        MesPrint("&%s should be a built in function that can be evaluated numerically.",s);
        return(1);
      }
      else {
        switch ( numfun+FUNCTION ) {
            case MZV:
            case EULER:
            case MZVHALF:
/*
            The following functions are treated in evaluate.c
*/
            case SQRTFUNCTION:
            case LNFUNCTION:
            case EXPFUNCTION:
            case SINFUNCTION:
            case COSFUNCTION:
            case TANFUNCTION:
            case ASINFUNCTION:
            case ACOSFUNCTION:
            case ATANFUNCTION:
            case ATAN2FUNCTION:
            case SINHFUNCTION:
            case COSHFUNCTION:
            case TANHFUNCTION:
            case ASINHFUNCTION:
            case ACOSHFUNCTION:
            case ATANHFUNCTION:
            case LI2FUNCTION:
            case LINFUNCTION:
            case AGMFUNCTION:
            case GAMMAFUN:
/*
            At a later stage we can add more functions from mpfr here
                mpfr_(number,arg(s))
*/
                Add3Com(TYPEEVALUATE,numfun+FUNCTION);
                break;
            default:
                MesPrint("&%s should be a built in function that can be evaluated numerically.",s);
                error = 1;
                break;
        }
      }
      *s = c;
      while ( *s == ' ' || *s == ',' || *s == '\t' ) s++;
    }
    return(error);
}
 
/*
        #] CoEvaluate : 
    #] Functions : 
*/