18125652002

Committed 30 Sep 2025 09:38AM UTC coverage: 54.367% (+0.3%) from 54.112%

Build # 18125652002

Build Type

push

github

Committed by

tueda

Commit Message

refactor: remove workaround for FLINT >= 3.0.0 and < 3.2.0

We now require FLINT >= 3.2.0. The workaround is no longer necessary.

Run Details

44447 of 81754 relevant lines covered (54.37%)

2275204.07 hits per line

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

75.81

/sources/float.c

/** @file float.c
 * 
 *  This file contains numerical routines that deal with floating point numbers.
 *        We use the GMP for arbitrary precision floating point arithmetic.
 *        The functions of the type sin, cos, ln, sqrt etc are handled by the
 *        mpfr library. The reason that for MZV's and the general notation we use
 *        the GMP (mpf_) is because the contents of the structs have been frozen
 *        and can be used for storage in a Form function float_. With mpfr_ this
 *        is not safely possible. All mpfr_ related code is in the file evaluate.c.
 */
/* #[ License : */
/*
 *   Copyright (C) 1984-2023 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, int *, int);
void deltaEuler(mpf_t, int *, int);
void deltaEulerC(mpf_t, int *, int);
void CalculateMZVhalf(mpf_t, int *, int);
void CalculateMZV(mpf_t, int *, int);
void CalculateEuler(mpf_t, int *, 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 ) {
        }
        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 ) {
                MesPrint("Illegal attempt at using a float_ function without proper startup.");
                MesPrint("Please use %#StartFloat <options> first.");
                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)+f[-5]);
                }
                else if ( f[-1] == 5 ) {
                        outfloat->_mp_exp =
                                 (mp_exp_t)((((ULONG)(f[-4]))<<BITSINWORD)+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[0]; 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;
        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;
                        totalbitsused += bitsused;
                }
                if ( bitsused > (totalbits-totalbitsused)/2 ) { break; }
                if ( mpf_sgn(aux2) == 0 ) {
                        /*if ( startnul == 1 )*/ AT.pWorkSpace[AT.pWorkPointer++] = out;
                        break;
                }
                AT.pWorkSpace[AT.pWorkPointer++] = out;
          }
/*
          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;
          if ( nc > 10 ) {
                mc = c;
                c = (UWORD *)(AT.pWorkSpace[--AT.pWorkPointer]);
                nc = nb = ((UWORD *)mc)-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 = *ratout;
                RatToFloat(aux1,ratout,n);
                mpf_sub(aux2,floatin,aux1);
                gmp_printf("Diff is %.*Fe\n",40,aux2);
        }
*/
        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 zero = 1, gotdot = 0;
        while ( FG.cTable[s[-1]] == 1 ) s--;
        *spec = 0;
        if ( FG.cTable[*s] == 1 ) {
                while ( *s == '0' ) s++;
                if ( FG.cTable[*s] == 1 ) {
                        s++;
                        while ( FG.cTable[*s] == 1 ) s++;
                        zero = 0;
                }
                if ( *s == '.' ) { goto dot; }
        }
        else if ( *s == '.' ) {
dot:
                gotdot = 1;
                s++;
                if ( FG.cTable[*s] != 1 && zero == 1 ) return(ss);
                while ( *s == '0' ) s++;
                if ( FG.cTable[*s] == 1 ) {
                        s++;
                        while ( FG.cTable[*s] == 1 ) s++;
                        zero = 0;
                }
        }
        else return(ss);
/*
        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);
        if ( AT.aux_ == 0 ) { /* no floating point system */
                *spec = -1;
                return(s);
        }
        if ( zero ) *spec = 1;
        return(s);
}

/*
                 #] CheckFloat : 
          #] Low Level : 
          #[ Float Routines :
                 #[ SetFloatPrecision :

                We set the default precision of the floats and allocate
                space for an output string if we want to write the float.
                Space needed: exponent: up to 12 chars.
                mantissa 2+10*prec/33 + 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 = ((10*prec)/33+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++ ) {
                if ( AB[id]->T.mpf_tab1 ) {
                        a = (mpf_t *)AB[id]->T.mpf_tab1;
                        for ( i = 0; i <=Nw; i++ ) {
                                mpf_clear(a[i]);
                        }
                        M_free(AB[id]->T.mpf_tab1,"mpftab1");
                }
                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]);
                }
                if ( AB[id]->T.mpf_tab2 ) {
                        a = (mpf_t *)AB[id]->T.mpf_tab2;
                        for ( i = 0; i <=Nw; i++ ) {
                                mpf_clear(a[i]);
                        }
                        M_free(AB[id]->T.mpf_tab2,"mpftab2");
                }
                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;
        if ( AT.mpf_tab1 ) {
                for ( i = 0; i <= Nw; i++ ) {
                        mpf_clear(mpftab1[i]);
                }
                M_free(AT.mpf_tab1,"mpftab1");
        }
        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]);
        }
        if ( AT.mpf_tab2 ) {
                for ( i = 0; i <= Nw; i++ ) {
                        mpf_clear(mpftab2[i]);
                }
                M_free(AT.mpf_tab2,"mpftab2");
        }
        AT.mpf_tab2 = (void *)Malloc1((N+2)*sizeof(mpf_t),"mpftab2");
        for ( i = 0; i <= Nw; i++ ) {
                mpf_init(mpftab2[i]);
        }
#endif
        if ( AS.delta_1 ) {
                mpf_clear(mpfdelta1);
                M_free(AS.delta_1,"delta1");
        }
        AS.delta_1 = (void *)Malloc1(sizeof(mpf_t),"delta1");
        mpf_init(mpfdelta1);
        SimpleDelta(mpfdelta1,1); /* this can speed up things. delta1 = ln(2) */
/*
        Finally the character buffer for printing
        if ( AO.floatspace ) M_free(AO.floatspace,"floatspace");
        AO.floatspace = (UBYTE *)Malloc1(((10*AC.DefaultPrecision)/33+40)*sizeof(UBYTE),"floatspace");
*/
}

/*
                 #] SetupMZVTables : 
                 #[ SetupMPFTables :
*/

void SetupMPFTables(void)
{
#ifdef WITHPTHREADS
        int id, totnum;
        mpf_t *a;
/*
        Now the aux variables
*/
#ifdef WITHSORTBOTS
        totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads);
#endif
    for ( id = 0; id < totnum; id++ ) {
                if ( AB[id]->T.aux_ ) {
/*
                        We work here with a[0] etc because the aux1 etc contain B which
                        in the current routine would be AB[0] only
*/
                        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_ = (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 = (int *)Malloc1(sizeof(int)*AC.MaxWeight*2,"indi1");
                AB[id]->T.indi2 = AB[id]->T.indi1 + AC.MaxWeight;
        }
#else
/*
        Now the aux variables
*/
        if ( AT.aux_ ) {
                mpf_clears(aux1,aux2,aux3,aux4,aux5,auxjm,auxjjm,auxsum,(mpf_ptr)0);
                M_free(AT.aux_,"AT.aux");
        }
        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 = (int *)Malloc1(sizeof(int)*AC.MaxWeight*2,"indi1");
        AT.indi2 = AT.indi1 + AC.MaxWeight;
#endif
/*
        Finally the character buffer for printing
        if ( AO.floatspace ) M_free(AO.floatspace,"floatspace");
        AO.floatspace = (UBYTE *)Malloc1(((10*AC.DefaultPrecision)/33+40)*sizeof(UBYTE),"floatspace");
*/
}

/*
                 #] 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 : 
                 #[ 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 ( FloatToRat((UWORD *)t,&ncoef,aux4) == 0 ) {
                        t += ABS(ncoef);
                        t[-1] = ncoef*nsign;
                        *term = t - term;
                }
        }
        return(Generator(BHEAD term,level));
}

/*
                 #] ToRat : 
          #] 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 ) {        /* May be rare! */
                *ps1 = 0; *ps2 = 0; AT.SortFloatMode = 0; return(0);
        }
        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,j,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);
        if ( sign3 == 0 ) {        /* May be very rare! */
                AT.SortFloatMode = 0; return(0);
        }
/*
        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. */
Over1:                jj = fun3[1]-fun1[1]+3-ABS(size1); /* This is positive */
                        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 goto Over1;
        }
        else { /* Can only be 2, based on previous tests */
                if ( fun3[1] + 3 == ABS(size1) ) {
                        t2 = coef1; t1 = fun3;
                        for ( i = 0; i < fun3[1]; i++ ) *t2++ = *t1++;
                        *t2++ = 1; *t2++ = 1;  *t2++ = sign3 < 0 ? -3: 3;
                        retval = 1;
                }
                else if ( fun3[1] + 3 < ABS(size1) ) {
                        j = ABS(size1) - fun3[1] - 3;
                        t2 = term1 + *term1; tt = t2;
                        *--t2 = sign3 < 0 ? -3: 3; *--t2 = 1; *--t2 = 1;
                        t2 -= fun3[1]; t1 = t2-j;
                        while ( t2 > term1 ) *--t2 = *--t1;
                        *t2 = tt-t2; term1 = t2;
                        retval = 1;
                }
                else 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);
        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);
        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 ) {
                printf("Wrong parameter pow in SingleTable: %d\n",pow);
                exit(-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 ) {
                printf("Wrong parameter pow in SingleTable: %d\n",pow);
                exit(-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 ) {
                printf("Wrong parameter pow in SingleTable: %d\n",pow);
                exit(-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, int *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, int *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, int *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 ) {
                        printf("Illegal index in depth=1 deltaEulerC: %d\n",indexes[0]);
                }
                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, int *indexes, int depth)
{
        int i;
        if ( depth < 0 ) {
                printf("Illegal depth in CalculateMZVhalf: %d\n",depth);
                exit(-1);
        }
        for ( i = 0; i < depth; i++ ) {
                if ( indexes[i] <= 0 ) {
                        printf("Illegal index[%d] in CalculateMZVhalf: %d\n",i,indexes[i]);
                        exit(-1);
                }
        }
        deltaMZV(result,indexes,depth);
}

/*
                 #] CalculateMZVhalf : 
                 #[ CalculateMZV :
*/

void CalculateMZV(mpf_t result, int *indexes, int depth)
{
        GETIDENTITY
        int num1, num2 = 0, i, j = 0;
        if ( depth < 0 ) {
                printf("Illegal depth in CalculateMZV: %d\n",depth);
                exit(-1);
        }
        if ( indexes[0] == 1 ) {
                printf("Divergent MZV in CalculateMZV\n");
                exit(-1);
        }
/*        MesPrint("calculateMZV(%a)",depth,indexes); */
        for ( i = 0; i < depth; i++ ) {
                if ( indexes[i] <= 0 ) {
                        printf("Illegal index[%d] in CalculateMZV: %d\n",i,indexes[i]);
                        exit(-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, int *Zindexes, int depth)
{
        GETIDENTITY
        int s1, num1, num2, i, j;

        int *indexes = (int *)(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 ) {
                printf("Illegal depth in CalculateEuler: %d\n",depth);
                exit(-1);
        }
        if ( indexes[0] == 1 ) {
                printf("Divergent Euler sum in CalculateEuler\n");
                exit(-1);
        }
        for ( i = 0, j = 0; i < depth; i++ ) {
                if ( indexes[i] == 0 ) {
                        printf("Illegal index[%d] in CalculateEuler: %d\n",i,indexes[i]);
                        exit(-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 ) {
                                MesPrint("Error: Weight of Euler/MZV sum greater than %d",sumweight);
                                MesPrint("Please increase MaxWeight in form.set.");
                                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 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 : 
*/

tueda / form / 18125652002

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