diff options
Diffstat (limited to 'ic-reals-6.3/math-lib')
54 files changed, 5070 insertions, 0 deletions
diff --git a/ic-reals-6.3/math-lib/Makefile b/ic-reals-6.3/math-lib/Makefile new file mode 100644 index 0000000..b370191 --- /dev/null +++ b/ic-reals-6.3/math-lib/Makefile @@ -0,0 +1,51 @@ +OBJS = init.o \ + tan_R.o \ + tan_Q.o \ + atan_R.o \ + atan_Q.o \ + tanh_R.o \ + tanh_Q.o \ + atanh_R.o \ + atanh_Q.o \ + sin_R.o \ + sin_Q.o \ + asin_R.o \ + asin_Q.o \ + sinh_R.o \ + sinh_Q.o \ + asinh_R.o \ + asinh_Q.o \ + cos_R.o \ + cos_Q.o \ + acos_R.o \ + acos_Q.o \ + cosh_R.o \ + cosh_Q.o \ + acosh_R.o \ + acosh_Q.o \ + secant.o \ + cosecant.o \ + cotangent.o \ + log_R.o \ + log_Q.o \ + exp_R.o \ + exp_Q.o \ + sqrt_R.o \ + sqrt_Q.o \ + pow_R_R.o \ + pi.o \ + abs_R.o \ + neg_R.o \ + stdTensorCont.o + +force: $(OBJS) +$(OBJS): ../real.h ../real-impl.h + +stdTensor.o : math-lib.h +tan_R.o : math-lib.h +atan_R.o : math-lib.h +exp_R.o : math-lib.h +log_R.o : math-lib.h + +clean: + rm -f $(OBJS) diff --git a/ic-reals-6.3/math-lib/abs_R.c b/ic-reals-6.3/math-lib/abs_R.c new file mode 100644 index 0000000..cc2507d --- /dev/null +++ b/ic-reals-6.3/math-lib/abs_R.c @@ -0,0 +1,467 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" + +/* + * This computes the absolute value of a real. For now the argument + * must be either a vector or a stream. + */ +Real +abs_R(Real x) +{ + Real r; + Cls *cls; + SignX *signX; + int sign_a, sign_b; + void force_To_SignX_From_Abs_Entry(); + void force_To_SignX_From_Abs_Cont(); + void force_To_DigsX_From_Abs_Entry(); + void force_To_DigsX_From_Abs_Cont(); + void force_To_DigsX_From_Abs_Copy_Entry(); + void force_To_DigsX_From_Abs_Copy_Cont(); + void force_To_DigsX_From_Abs_Copy_and_Negate(); + void force_To_DigsX_From_Abs_Copy_and_Negate_Entry(); + void force_To_DigsX_From_Abs_Copy_and_Negate_Cont(); + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(force_To_SignX_From_Abs_Entry, + "force_To_SignX_From_Abs_Entry", 2); + registerForceFunc(force_To_SignX_From_Abs_Cont, + "force_To_SignX_From_Abs_Cont", 2); + registerForceFunc(force_To_DigsX_From_Abs_Entry, + "force_To_DigsX_From_Abs_Entry", 3); + registerForceFunc(force_To_DigsX_From_Abs_Cont, + "force_To_DigsX_From_Abs_Cont", 3); + registerForceFunc(force_To_DigsX_From_Abs_Copy_Entry, + "force_To_DigsX_From_Abs_Copy_Entry", 3); + registerForceFunc(force_To_DigsX_From_Abs_Copy_Cont, + "force_To_DigsX_From_Abs_Copy_Cont", 3); + registerForceFunc(force_To_DigsX_From_Abs_Copy_and_Negate_Entry, + "force_To_DigsX_From_Abs_Copy_and_Negate_Entry", 3); + registerForceFunc(force_To_DigsX_From_Abs_Copy_and_Negate_Cont, + "force_To_DigsX_From_Abs_Copy_and_Negate_Cont", 3); + doneInit++; + } + + /* + * A slight optimization. In the case of a vector we just inspect the + * entries and then make them both the same sign (in a copy of the vector). + */ + if (x->gen.tag.type == VECTOR) { + sign_a = mpz_sgn(x->vec.vec[0]); + sign_b = mpz_sgn(x->vec.vec[1]); + if (sign_a * sign_b < 0) { + r = vector_Z(x->vec.vec[0], x->vec.vec[1]); + mpz_neg(r->vec.vec[0], r->vec.vec[0]); + return r; + } + else + return x; + } + + /* + * For all other argument types, we simply work with the stream + * representation of the argument. + */ + if (x->gen.tag.isSigned) { + r = makeStream(x); /* this will make a signed stream */ + + /* + * All the work is done in the closure. + */ + cls = allocCls(force_To_SignX_From_Abs_Entry, (void *) r); +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(cls, r); + endGraphUpdate(); +#endif + signX = allocSignX((Real) cls, SIGN_UNKN); + signX->force = force_To_SignX_From_Abs_Entry; + } + else + return x; + return NULL; +} + +/* + * In spite of the fact that the output is positive, we don't + * really know the sign of abs_R(x) until we get the sign of x. + */ +void +force_To_SignX_From_Abs_Entry() +{ + SignX *signX, *source; + Cls *cls; + void force_To_SignX_From_Abs_Cont(); + + signX = (SignX *) POP; + cls = (Cls *) signX->x; + source = (SignX *) cls->userData; + + PUSH_2(force_To_SignX_From_Abs_Cont, signX); + + if (source->tag.value == SIGN_UNKN) + PUSH_2(source->force, source); +} + +/* + * When this gets activated, then we know the sign of the argument + * has been determined. The sign of the output is the same as the argument + * unless the argument is SNEG in which case we switch to SPOS. + */ +void +force_To_SignX_From_Abs_Cont() +{ + SignX *signX; + Cls *cls; + SignX *source; + void force_To_DigsX_From_Abs_Entry(); + + signX = (SignX *) POP; + cls = (Cls *) signX->x; + source = (SignX *) cls->userData; + + /* + * This allocates an empty DigsX structure between the sign and the + * closure. That's where we will copy digits to eventually. + */ + introDigsX(signX); + signX->x->digsX.force = force_To_DigsX_From_Abs_Entry; + + switch (source->tag.value) { + case SPOS : + case SZERO : + case SINF : + signX->tag.value = source->tag.value; + break; + case SNEG : + signX->tag.value = SPOS; /* the only case where sign is changed */ + break; + case SIGN_UNKN : + Error(FATAL, E_INT, "force_To_SignX_From_Abs_Cont", "sign not set"); + return; + default : + Error(FATAL, E_INT, "force_To_SignX_From_Abs_Cont", "bad sign"); + return; + } + + /* + * Now we skip over the sign of the argument. It has been ``absorbed''. + * But we keep a copy of the value for future use. + */ + cls->tag.value = source->tag.value; + cls->userData = (void *) source->x; + +#ifdef DAVINCI + beginGraphUpdate(); + deleteOnlyEdge(cls, source); + newEdgeToOnlyChild(cls, source->x); + endGraphUpdate(); +#endif +} + +void +force_To_DigsX_From_Abs_Entry() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + void force_To_DigsX_From_Abs_Cont(); + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + + PUSH_3(force_To_DigsX_From_Abs_Cont, target, digitsNeeded); + + /* + * Now see if the source has the number of digits we need. If not, + * then force the remaining. + */ + if (source->count < (unsigned int)digitsNeeded) + PUSH_3(source->force, source, digitsNeeded - source->count); +} + +/* + * The is the function which translates one stream of digits into + * another. The rewrite rules are as follows depending on the sign of + * the argument x. + * SZERO : if the word in the characteristic pair is negative, then + * the argument is negative and we must negate the word and every + * word in subsequent characteristic pairs. + * SINF : if the word in the characteristic pair is positive, then + * the argument is negative and we must negate the word and every + * word in subsequent characteristic pairs. + * SNEG : the number is negative, so we must negate the word of every + * characteristic pair. + * SPOS : just copy characteristic pairs from the source to the target. + * There *might* be an opportunity to optimize things by just + * short circuiting around the closure. Needs thought. + * The code is more or less copied from the force methods for DigsX structures + * with the negation stuff thrown in along the way. + */ +void +force_To_DigsX_From_Abs_Cont() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + void force_To_DigsX_From_Abs_Copy_Entry(); + void force_To_DigsX_From_Abs_Copy_and_Negate_Entry(); + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + + if (source->count > 0) { + switch (cls->tag.value) { + case SZERO : +#ifdef OPTIM + if (source->count <= DIGITS_PER_WORD) { + if (source->word.small < 0) + cls->force = force_To_DigsX_From_Abs_Copy_and_Negate_Entry; + else + cls->force = force_To_DigsX_From_Abs_Copy_Entry; + } + else { +#endif + if (mpz_sgn(source->word.big) < 0) + cls->force = force_To_DigsX_From_Abs_Copy_and_Negate_Entry; + else + cls->force = force_To_DigsX_From_Abs_Copy_Entry; +#ifdef OPTIM + } +#endif + break; + + case SINF: +#ifdef OPTIM + if (source->count <= DIGITS_PER_WORD) { + if (source->word.small > 0) + cls->force = force_To_DigsX_From_Abs_Copy_and_Negate_Entry; + else + cls->force = force_To_DigsX_From_Abs_Copy_Entry; + } + else { +#endif + if (mpz_sgn(source->word.big) > 0) + cls->force = force_To_DigsX_From_Abs_Copy_and_Negate_Entry; + else + cls->force = force_To_DigsX_From_Abs_Copy_Entry; +#ifdef OPTIM + } +#endif + break; + + case SNEG: + cls->force = force_To_DigsX_From_Abs_Copy_and_Negate_Entry; + break; + + case SPOS: + cls->force = force_To_DigsX_From_Abs_Copy_Entry; + break; + + default : + Error(FATAL, E_INT, "force_To_DigsX_From_Abs_Cont", + "bad sign"); + } + target->force = cls->force; + PUSH_3(target->force, target, digitsNeeded); + } +} + +void +force_To_DigsX_From_Abs_Copy_and_Negate_Entry() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + void force_To_DigsX_From_Abs_Copy_and_Negate_Cont(); + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + + PUSH_3(force_To_DigsX_From_Abs_Copy_and_Negate_Cont, target, digitsNeeded); + + /* + * See if the source has the number of digits we need. If not, + * then force the remaining. + */ + if (source->count < (unsigned int)digitsNeeded) + PUSH_3(source->force, source, digitsNeeded - source->count); +} + +void +force_To_DigsX_From_Abs_Copy_and_Negate_Cont() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + +#ifdef OPTIM + if (target->count + source->count <= DIGITS_PER_WORD) + target->word.small = (target->word.small << source->count) + - source->word.small; + else { + if (target->count <= DIGITS_PER_WORD) + mpz_init_set_si(target->word.big, target->word.small); +#endif + mpz_mul_2exp(target->word.big, target->word.big, source->count); +#ifdef OPTIM + if (source->count <= DIGITS_PER_WORD) { + if (source->word.small >= 0) + mpz_sub_ui(target->word.big, target->word.big, + source->word.small); + else + mpz_add_ui(target->word.big, target->word.big, + -(source->word.small)); + } + else +#endif + mpz_sub(target->word.big, target->word.big, source->word.big); +#ifdef OPTIM + } +#endif + + target->count += source->count; + +#ifdef TRACE + if (TRACE) { + debugp("force_To_DigsX_From_Abs_Copy_and_Negate", + "%x %x emitted=%d\n", + (unsigned) target, + (unsigned) source, + source->count); + } +#endif + /* + * We've consumed the source so advance to the next possible source + * of information + */ + cls->userData = (void *) source->x; + +#ifdef DAVINCI + beginGraphUpdate(); + deleteOnlyEdge(cls, source); + newEdgeToOnlyChild(cls, source->x); + endGraphUpdate(); +#endif + newDigsX(target); +} + +void +force_To_DigsX_From_Abs_Copy_Entry() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + void force_To_DigsX_From_Abs_Copy_Cont(); + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + + PUSH_3(force_To_DigsX_From_Abs_Copy_Cont, target, digitsNeeded); + + /* + * See if the source has the number of digits we need. If not, + * then force the remaining. + */ + if (source->count < (unsigned int)digitsNeeded) + PUSH_3(source->force, source, digitsNeeded - source->count); +} + +void +force_To_DigsX_From_Abs_Copy_Cont() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + +#ifdef OPTIM + if (target->count + source->count <= DIGITS_PER_WORD) + target->word.small = (target->word.small << source->count) + + source->word.small; + else { + if (target->count <= DIGITS_PER_WORD) + mpz_init_set_si(target->word.big, target->word.small); +#endif + mpz_mul_2exp(target->word.big, target->word.big, source->count); +#ifdef OPTIM + if (source->count <= DIGITS_PER_WORD) { + if (source->word.small >= 0) { + mpz_add_ui(target->word.big, target->word.big, + source->word.small); + } + else { + mpz_sub_ui(target->word.big, target->word.big, + -(source->word.small)); + } + } + else +#endif + mpz_add(target->word.big, target->word.big, source->word.big); +#ifdef OPTIM + } +#endif + + target->count += source->count; + +#ifdef TRACE + if (TRACE) { + debugp("force_To_DigsX_From_Abs_Copy", + "%x %x emitted=%d\n", + (unsigned) target, + (unsigned) source, + source->count); + } +#endif + /* + * We've consumed the source so advance to the next possible source + * of information + */ + cls->userData = (void *) source->x; + +#ifdef DAVINCI + beginGraphUpdate(); + deleteOnlyEdge(cls, source); + newEdgeToOnlyChild(cls, source->x); + endGraphUpdate(); +#endif + newDigsX(target); +} + +#ifdef JUNK +In some cases it would be better to use the matices to +negate things. But to do this for absolute value +we must know the sign of the argument. + +As a rule we want to use a matrix for reciprocal and negation +except at the top level when we might prefer to use algorithms +in Peters thesis. This is also the case when we real is +the argument to a predicate. +#endif diff --git a/ic-reals-6.3/math-lib/abs_R.c~ b/ic-reals-6.3/math-lib/abs_R.c~ new file mode 100644 index 0000000..e9a7fad --- /dev/null +++ b/ic-reals-6.3/math-lib/abs_R.c~ @@ -0,0 +1,462 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" + +/* + * This computes the absolute value of a real. For now the argument + * must be either a vector or a stream. + */ +Real +abs_R(Real x) +{ + Real r; + Cls *cls; + SignX *signX; + int sign_a, sign_b; + void force_To_SignX_From_Abs_Entry(); + void force_To_SignX_From_Abs_Cont(); + void force_To_DigsX_From_Abs_Entry(); + void force_To_DigsX_From_Abs_Cont(); + void force_To_DigsX_From_Abs_Copy_Entry(); + void force_To_DigsX_From_Abs_Copy_Cont(); + void force_To_DigsX_From_Abs_Copy_and_Negate(); + void force_To_DigsX_From_Abs_Copy_and_Negate_Entry(); + void force_To_DigsX_From_Abs_Copy_and_Negate_Cont(); + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(force_To_SignX_From_Abs_Entry, + "force_To_SignX_From_Abs_Entry", 2); + registerForceFunc(force_To_SignX_From_Abs_Cont, + "force_To_SignX_From_Abs_Cont", 2); + registerForceFunc(force_To_DigsX_From_Abs_Entry, + "force_To_DigsX_From_Abs_Entry", 3); + registerForceFunc(force_To_DigsX_From_Abs_Cont, + "force_To_DigsX_From_Abs_Cont", 3); + registerForceFunc(force_To_DigsX_From_Abs_Copy_Entry, + "force_To_DigsX_From_Abs_Copy_Entry", 3); + registerForceFunc(force_To_DigsX_From_Abs_Copy_Cont, + "force_To_DigsX_From_Abs_Copy_Cont", 3); + registerForceFunc(force_To_DigsX_From_Abs_Copy_and_Negate_Entry, + "force_To_DigsX_From_Abs_Copy_and_Negate_Entry", 3); + registerForceFunc(force_To_DigsX_From_Abs_Copy_and_Negate_Cont, + "force_To_DigsX_From_Abs_Copy_and_Negate_Cont", 3); + doneInit++; + } + + /* + * A slight optimization. In the case of a vector we just inspect the + * entries and then make them both the same sign (in a copy of the vector). + */ + if (x->gen.tag.type == VECTOR) { + sign_a = mpz_sgn(x->vec.vec[0]); + sign_b = mpz_sgn(x->vec.vec[1]); + if (sign_a * sign_b < 0) { + r = vector_Z(x->vec.vec[0], x->vec.vec[1]); + mpz_neg(r->vec.vec[0], r->vec.vec[0]); + return r; + } + else + return x; + } + + /* + * For all other argument types, we simply work with the stream + * representation of the argument. + */ + if (x->gen.tag.isSigned) { + r = makeStream(x); /* this will make a signed stream */ + + /* + * All the work is done in the closure. + */ + cls = allocCls(force_To_SignX_From_Abs_Entry, (void *) r); + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(cls, r); + endGraphUpdate(); + } + signX = allocSignX((Real) cls, SIGN_UNKN); + signX->force = force_To_SignX_From_Abs_Entry; + } + else + return x; +} + +/* + * In spite of the fact that the output is positive, we don't + * really know the sign of abs_R(x) until we get the sign of x. + */ +void +force_To_SignX_From_Abs_Entry() +{ + SignX *signX, *source; + Cls *cls; + void force_To_SignX_From_Abs_Cont(); + + signX = (SignX *) POP; + cls = (Cls *) signX->x; + source = (SignX *) cls->userData; + + PUSH_2(force_To_SignX_From_Abs_Cont, signX); + + if (source->tag.value == SIGN_UNKN) + PUSH_2(source->force, source); +} + +/* + * When this gets activated, then we know the sign of the argument + * has been determined. The sign of the output is the same as the argument + * unless the argument is SNEG in which case we switch to SPOS. + */ +void +force_To_SignX_From_Abs_Cont() +{ + SignX *signX; + Cls *cls; + SignX *source; + void force_To_DigsX_From_Abs_Entry(); + + signX = (SignX *) POP; + cls = (Cls *) signX->x; + source = (SignX *) cls->userData; + + /* + * This allocates an empty DigsX structure between the sign and the + * closure. That's where we will copy digits to eventually. + */ + introDigsX(signX); + signX->x->digsX.force = force_To_DigsX_From_Abs_Entry; + + switch (source->tag.value) { + case SPOS : + case SZERO : + case SINF : + signX->tag.value = source->tag.value; + break; + case SNEG : + signX->tag.value = SPOS; /* the only case where sign is changed */ + break; + case SIGN_UNKN : + Error(FATAL, E_INT, "force_To_SignX_From_Abs_Cont", "sign not set"); + default : + Error(FATAL, E_INT, "force_To_SignX_From_Abs_Cont", "bad sign"); + } + + /* + * Now we skip over the sign of the argument. It has been ``absorbed''. + * But we keep a copy of the value for future use. + */ + cls->tag.value = source->tag.value; + cls->userData = (void *) source->x; + + if (DAVINCI) { + beginGraphUpdate(); + deleteOnlyEdge(cls, source); + newEdgeToOnlyChild(cls, source->x); + endGraphUpdate(); + } +} + +void +force_To_DigsX_From_Abs_Entry() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + void force_To_DigsX_From_Abs_Cont(); + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + + PUSH_3(force_To_DigsX_From_Abs_Cont, target, digitsNeeded); + + /* + * Now see if the source has the number of digits we need. If not, + * then force the remaining. + */ + if (source->count < digitsNeeded) + PUSH_3(source->force, source, digitsNeeded - source->count); +} + +/* + * The is the function which translates one stream of digits into + * another. The rewrite rules are as follows depending on the sign of + * the argument x. + * SZERO : if the word in the characteristic pair is negative, then + * the argument is negative and we must negate the word and every + * word in subsequent characteristic pairs. + * SINF : if the word in the characteristic pair is positive, then + * the argument is negative and we must negate the word and every + * word in subsequent characteristic pairs. + * SNEG : the number is negative, so we must negate the word of every + * characteristic pair. + * SPOS : just copy characteristic pairs from the source to the target. + * There *might* be an opportunity to optimize things by just + * short circuiting around the closure. Needs thought. + * The code is more or less copied from the force methods for DigsX structures + * with the negation stuff thrown in along the way. + */ +void +force_To_DigsX_From_Abs_Cont() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + void force_To_DigsX_From_Abs_Copy_Entry(); + void force_To_DigsX_From_Abs_Copy_and_Negate_Entry(); + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + + if (source->count > 0) { + switch (cls->tag.value) { + case SZERO : +#ifdef OPTIM + if (source->count <= DIGITS_PER_WORD) { + if (source->word.small < 0) + cls->force = force_To_DigsX_From_Abs_Copy_and_Negate_Entry; + else + cls->force = force_To_DigsX_From_Abs_Copy_Entry; + } + else { +#endif + if (mpz_sgn(source->word.big) < 0) + cls->force = force_To_DigsX_From_Abs_Copy_and_Negate_Entry; + else + cls->force = force_To_DigsX_From_Abs_Copy_Entry; +#ifdef OPTIM + } +#endif + break; + + case SINF: +#ifdef OPTIM + if (source->count <= DIGITS_PER_WORD) { + if (source->word.small > 0) + cls->force = force_To_DigsX_From_Abs_Copy_and_Negate_Entry; + else + cls->force = force_To_DigsX_From_Abs_Copy_Entry; + } + else { +#endif + if (mpz_sgn(source->word.big) > 0) + cls->force = force_To_DigsX_From_Abs_Copy_and_Negate_Entry; + else + cls->force = force_To_DigsX_From_Abs_Copy_Entry; +#ifdef OPTIM + } +#endif + break; + + case SNEG: + cls->force = force_To_DigsX_From_Abs_Copy_and_Negate_Entry; + break; + + case SPOS: + cls->force = force_To_DigsX_From_Abs_Copy_Entry; + break; + + default : + Error(FATAL, E_INT, "force_To_DigsX_From_Abs_Cont", + "bad sign"); + } + target->force = cls->force; + PUSH_3(target->force, target, digitsNeeded); + } +} + +void +force_To_DigsX_From_Abs_Copy_and_Negate_Entry() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + void force_To_DigsX_From_Abs_Copy_and_Negate_Cont(); + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + + PUSH_3(force_To_DigsX_From_Abs_Copy_and_Negate_Cont, target, digitsNeeded); + + /* + * See if the source has the number of digits we need. If not, + * then force the remaining. + */ + if (source->count < digitsNeeded) + PUSH_3(source->force, source, digitsNeeded - source->count); +} + +void +force_To_DigsX_From_Abs_Copy_and_Negate_Cont() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + +#ifdef OPTIM + if (target->count + source->count <= DIGITS_PER_WORD) + target->word.small = (target->word.small << source->count) + - source->word.small; + else { + if (target->count <= DIGITS_PER_WORD) + mpz_init_set_si(target->word.big, target->word.small); +#endif + mpz_mul_2exp(target->word.big, target->word.big, source->count); +#ifdef OPTIM + if (source->count <= DIGITS_PER_WORD) { + if (source->word.small >= 0) + mpz_sub_ui(target->word.big, target->word.big, + source->word.small); + else + mpz_add_ui(target->word.big, target->word.big, + -(source->word.small)); + } + else +#endif + mpz_sub(target->word.big, target->word.big, source->word.big); +#ifdef OPTIM + } +#endif + + target->count += source->count; + + if (TRACE) { + debugp("force_To_DigsX_From_Abs_Copy_and_Negate", + "%x %x emitted=%d\n", + (unsigned) target, + (unsigned) source, + source->count); + } + + /* + * We've consumed the source so advance to the next possible source + * of information + */ + cls->userData = (void *) source->x; + + if (DAVINCI) { + beginGraphUpdate(); + deleteOnlyEdge(cls, source); + newEdgeToOnlyChild(cls, source->x); + endGraphUpdate(); + } + newDigsX(target); +} + +void +force_To_DigsX_From_Abs_Copy_Entry() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + void force_To_DigsX_From_Abs_Copy_Cont(); + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + + PUSH_3(force_To_DigsX_From_Abs_Copy_Cont, target, digitsNeeded); + + /* + * See if the source has the number of digits we need. If not, + * then force the remaining. + */ + if (source->count < digitsNeeded) + PUSH_3(source->force, source, digitsNeeded - source->count); +} + +void +force_To_DigsX_From_Abs_Copy_Cont() +{ + DigsX *target, *source; + Cls *cls; + int digitsNeeded; + + target = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) target->x; + source = (DigsX *) cls->userData; + +#ifdef OPTIM + if (target->count + source->count <= DIGITS_PER_WORD) + target->word.small = (target->word.small << source->count) + + source->word.small; + else { + if (target->count <= DIGITS_PER_WORD) + mpz_init_set_si(target->word.big, target->word.small); +#endif + mpz_mul_2exp(target->word.big, target->word.big, source->count); +#ifdef OPTIM + if (source->count <= DIGITS_PER_WORD) { + if (source->word.small >= 0) { + mpz_add_ui(target->word.big, target->word.big, + source->word.small); + } + else { + mpz_sub_ui(target->word.big, target->word.big, + -(source->word.small)); + } + } + else +#endif + mpz_add(target->word.big, target->word.big, source->word.big); +#ifdef OPTIM + } +#endif + + target->count += source->count; + + if (TRACE) { + debugp("force_To_DigsX_From_Abs_Copy", + "%x %x emitted=%d\n", + (unsigned) target, + (unsigned) source, + source->count); + } + + /* + * We've consumed the source so advance to the next possible source + * of information + */ + cls->userData = (void *) source->x; + + if (DAVINCI) { + beginGraphUpdate(); + deleteOnlyEdge(cls, source); + newEdgeToOnlyChild(cls, source->x); + endGraphUpdate(); + } + newDigsX(target); +} + +#ifdef JUNK +In some cases it would be better to use the matices to +negate things. But to do this for absolute value +we must know the sign of the argument. + +As a rule we want to use a matrix for reciprocal and negation +except at the top level when we might prefer to use algorithms +in Peters thesis. This is also the case when we real is +the argument to a predicate. +#endif diff --git a/ic-reals-6.3/math-lib/acos_Q.c b/ic-reals-6.3/math-lib/acos_Q.c new file mode 100644 index 0000000..f584c93 --- /dev/null +++ b/ic-reals-6.3/math-lib/acos_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +acos_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return acos_R(r); +} + +Real +acos_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return acos_R(r); +} diff --git a/ic-reals-6.3/math-lib/acos_R.c b/ic-reals-6.3/math-lib/acos_R.c new file mode 100644 index 0000000..0f62f41 --- /dev/null +++ b/ic-reals-6.3/math-lib/acos_R.c @@ -0,0 +1,26 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * acos(x) = atan((sqrt(1-x^2))/x) + */ +Real +acos_R(Real x) +{ + Real r; + + r = mul_R_R(x, x); + r = sub_Int_R(1, r); + r = sqrt_R(r); + r = div_R_R(r, x); + r = atan_R(r); + return r; +} diff --git a/ic-reals-6.3/math-lib/acos_R.c~ b/ic-reals-6.3/math-lib/acos_R.c~ new file mode 100644 index 0000000..4de1e7a --- /dev/null +++ b/ic-reals-6.3/math-lib/acos_R.c~ @@ -0,0 +1,26 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * acos(x) = atan((sqrt(1-x^2))/x) + */ +Real +acos_R(Real x) +{ + Real r, s; + + r = mul_R_R(x, x); + r = sub_Int_R(1, r); + r = sqrt_R(r); + r = div_R_R(r, x); + r = atan_R(r); + return r; +} diff --git a/ic-reals-6.3/math-lib/acosh_Q.c b/ic-reals-6.3/math-lib/acosh_Q.c new file mode 100644 index 0000000..9f0846f --- /dev/null +++ b/ic-reals-6.3/math-lib/acosh_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +acosh_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return acosh_R(r); +} + +Real +acosh_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return acosh_R(r); +} diff --git a/ic-reals-6.3/math-lib/acosh_R.c b/ic-reals-6.3/math-lib/acosh_R.c new file mode 100644 index 0000000..467539b --- /dev/null +++ b/ic-reals-6.3/math-lib/acosh_R.c @@ -0,0 +1,22 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +acosh_R(Real x) +{ + Real r; + + r = tensor_Int(x, x, 1, 0, 0, 0, 0, 0, -1, 1); + r = sqrt_R(r); + r = add_R_R(x, r); + r = log_R(r); + return r; +} diff --git a/ic-reals-6.3/math-lib/asin_Q.c b/ic-reals-6.3/math-lib/asin_Q.c new file mode 100644 index 0000000..b51eeaf --- /dev/null +++ b/ic-reals-6.3/math-lib/asin_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +asin_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return asin_R(r); +} + +Real +asin_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return asin_R(r); +} diff --git a/ic-reals-6.3/math-lib/asin_R.c b/ic-reals-6.3/math-lib/asin_R.c new file mode 100644 index 0000000..c12b7c4 --- /dev/null +++ b/ic-reals-6.3/math-lib/asin_R.c @@ -0,0 +1,26 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * asin(x) = atan(x/sqrt(1-x^2)) + */ +Real +asin_R(Real x) +{ + Real r; + + r = mul_R_R(x, x); + r = sub_Int_R(1, r); + r = sqrt_R(r); + r = div_R_R(x, r); + r = atan_R(r); + return r; +} diff --git a/ic-reals-6.3/math-lib/asin_R.c~ b/ic-reals-6.3/math-lib/asin_R.c~ new file mode 100644 index 0000000..009bcf2 --- /dev/null +++ b/ic-reals-6.3/math-lib/asin_R.c~ @@ -0,0 +1,26 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * asin(x) = atan(x/sqrt(1-x^2)) + */ +Real +asin_R(Real x) +{ + Real r, s; + + r = mul_R_R(x, x); + r = sub_Int_R(1, r); + r = sqrt_R(r); + r = div_R_R(x, r); + r = atan_R(r); + return r; +} diff --git a/ic-reals-6.3/math-lib/asinh_Q.c b/ic-reals-6.3/math-lib/asinh_Q.c new file mode 100644 index 0000000..4200b8d --- /dev/null +++ b/ic-reals-6.3/math-lib/asinh_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +asinh_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return asinh_R(r); +} + +Real +asinh_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return asinh_R(r); +} diff --git a/ic-reals-6.3/math-lib/asinh_R.c b/ic-reals-6.3/math-lib/asinh_R.c new file mode 100644 index 0000000..7c20839 --- /dev/null +++ b/ic-reals-6.3/math-lib/asinh_R.c @@ -0,0 +1,22 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +asinh_R(Real x) +{ + Real r; + + r = tensor_Int(x, x, 1, 0, 0, 0, 0, 0, 1, 1); + r = sqrt_R(r); + r = add_R_R(x, r); + r = log_R(r); + return r; +} diff --git a/ic-reals-6.3/math-lib/atan_Q.c b/ic-reals-6.3/math-lib/atan_Q.c new file mode 100644 index 0000000..f23a945 --- /dev/null +++ b/ic-reals-6.3/math-lib/atan_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +atan_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return atan_R(r); +} + +Real +atan_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return atan_R(r); +} diff --git a/ic-reals-6.3/math-lib/atan_R.c b/ic-reals-6.3/math-lib/atan_R.c new file mode 100644 index 0000000..c79efab --- /dev/null +++ b/ic-reals-6.3/math-lib/atan_R.c @@ -0,0 +1,135 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +static TenXY *nextTensor(Real, Real, int); +static void atanInside(); +static void atanPositive(); +static void atanNegative(); + +Real +atan_R(Real x) +{ + Real u, v, w; + Bool xLtEq1, xGtEqNeg1, xLtEq1_and_GtEqNeg1; + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(atanInside, "atanInside", 2); + registerForceFunc(atanPositive, "atanPositive", 2); + registerForceFunc(atanNegative, "atanNegative", 2); + doneInit++; + } + + xLtEq1 = ltEq_R_QInt(x, 1, 1); + xGtEqNeg1 = gtEq_R_QInt(x, -1, 1); + xLtEq1_and_GtEqNeg1 = and_B_B(xLtEq1, xGtEqNeg1); + + u = (Real) allocCls(atanInside, (void *) x); + u->cls.tag.isSigned = TRUE; + v = (Real) allocCls(atanPositive, (void *) x); + v->cls.tag.isSigned = TRUE; + w = (Real) allocCls(atanNegative, (void *) x); + w->cls.tag.isSigned = TRUE; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(u, x); + newEdgeToOnlyChild(v, x); + newEdgeToOnlyChild(w, x); + endGraphUpdate(); +#endif + + /* + * The order of the tests in the alt is not semantically + * significant. The tests are applied in order, so there is a + * very modest performance improvement by putting the -1,1 tests + * before the 999/1000,-999/1000 tests. Also there is a very + * small win by doing the negative test before the positive test. + */ + return realIf(3, + xLtEq1_and_GtEqNeg1, u, + gtEq_R_QInt(x, 999, 1000), v, + ltEq_R_QInt(x, -999, 1000), w); +} + +static TenXY * +nextTensor(Real x, Real y, int n) +{ + return (TenXY *) tensor_Int(x, y, 2*n+1, 1, -1, -(2*n+1), + 2*n+1, 2*n+1, 2*n+1, 2*n+1); +} + +static void +atanInside() +{ + Cls *cls, *newCls; + Real x; + ClsData *data; + void stdTensorCont(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + if ((data = (ClsData *) malloc(sizeof(ClsData))) == NULL) + Error(FATAL, E_INT, "atan_R", "malloc failed"); + + data->n = 1; + data->x = x; + data->nextTensor = nextTensor; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + cls->redirect = tensor_Int(x, (Real) newCls, 1, 1, 1, -1, 0, 1, 0, 1); + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif +} + +/* + * This implements the identity: + * x > 0 implies atan(x) = atan(szero * x) + \pi/4 + */ +static void +atanPositive() +{ + Cls *cls; + Real x, y; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + y = atan_R(matrix_Int(x, 1, 1, -1, 1)); + cls->redirect = add_R_R(y, div_R_Int(Pi, 4)); +} + +/* + * This implements the identity: + * x < 0 implies atan(x) = atan(szero * (sneg^(-1) * x) - \pi/4 + */ +static void +atanNegative() +{ + Cls *cls; + Real x, y; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + y = atan_R(matrix_Int(x, 1, -1, 1, 1)); + cls->redirect = sub_R_R(y, div_R_Int(Pi, 4)); +} diff --git a/ic-reals-6.3/math-lib/atan_R.c~ b/ic-reals-6.3/math-lib/atan_R.c~ new file mode 100644 index 0000000..a6c1154 --- /dev/null +++ b/ic-reals-6.3/math-lib/atan_R.c~ @@ -0,0 +1,135 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +static TenXY *nextTensor(Real, Real, int); +static void atanInside(); +static void atanPositive(); +static void atanNegative(); + +Real +atan_R(Real x) +{ + Real u, v, w; + Bool xLtEq1, xGtEqNeg1, xLtEq1_and_GtEqNeg1; + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(atanInside, "atanInside", 2); + registerForceFunc(atanPositive, "atanPositive", 2); + registerForceFunc(atanNegative, "atanNegative", 2); + doneInit++; + } + + xLtEq1 = ltEq_R_QInt(x, 1, 1); + xGtEqNeg1 = gtEq_R_QInt(x, -1, 1); + xLtEq1_and_GtEqNeg1 = and_B_B(xLtEq1, xGtEqNeg1); + + u = (Real) allocCls(atanInside, (void *) x); + u->cls.tag.isSigned = TRUE; + v = (Real) allocCls(atanPositive, (void *) x); + v->cls.tag.isSigned = TRUE; + w = (Real) allocCls(atanNegative, (void *) x); + w->cls.tag.isSigned = TRUE; + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(u, x); + newEdgeToOnlyChild(v, x); + newEdgeToOnlyChild(w, x); + endGraphUpdate(); + } + + /* + * The order of the tests in the alt is not semantically + * significant. The tests are applied in order, so there is a + * very modest performance improvement by putting the -1,1 tests + * before the 999/1000,-999/1000 tests. Also there is a very + * small win by doing the negative test before the positive test. + */ + return realIf(3, + xLtEq1_and_GtEqNeg1, u, + gtEq_R_QInt(x, 999, 1000), v, + ltEq_R_QInt(x, -999, 1000), w); +} + +static TenXY * +nextTensor(Real x, Real y, int n) +{ + return (TenXY *) tensor_Int(x, y, 2*n+1, 1, -1, -(2*n+1), + 2*n+1, 2*n+1, 2*n+1, 2*n+1); +} + +static void +atanInside() +{ + Cls *cls, *newCls; + Real x, atan_x; + ClsData *data; + void stdTensorCont(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + if ((data = (ClsData *) malloc(sizeof(ClsData))) == NULL) + Error(FATAL, E_INT, "atan_R", "malloc failed"); + + data->n = 1; + data->x = x; + data->nextTensor = nextTensor; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + cls->redirect = tensor_Int(x, (Real) newCls, 1, 1, 1, -1, 0, 1, 0, 1); + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} + +/* + * This implements the identity: + * x > 0 implies atan(x) = atan(szero * x) + \pi/4 + */ +static void +atanPositive() +{ + Cls *cls, *newCls; + Real x, y; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + y = atan_R(matrix_Int(x, 1, 1, -1, 1)); + cls->redirect = add_R_R(y, div_R_Int(Pi, 4)); +} + +/* + * This implements the identity: + * x < 0 implies atan(x) = atan(szero * (sneg^(-1) * x) - \pi/4 + */ +static void +atanNegative() +{ + Cls *cls, *newCls; + Real x, y; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + y = atan_R(matrix_Int(x, 1, -1, 1, 1)); + cls->redirect = sub_R_R(y, div_R_Int(Pi, 4)); +} diff --git a/ic-reals-6.3/math-lib/atanh_Q.c b/ic-reals-6.3/math-lib/atanh_Q.c new file mode 100644 index 0000000..d66cf14 --- /dev/null +++ b/ic-reals-6.3/math-lib/atanh_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +atanh_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return atanh_R(r); +} + +Real +atanh_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return atanh_R(r); +} diff --git a/ic-reals-6.3/math-lib/atanh_R.c b/ic-reals-6.3/math-lib/atanh_R.c new file mode 100644 index 0000000..c88e2d4 --- /dev/null +++ b/ic-reals-6.3/math-lib/atanh_R.c @@ -0,0 +1,24 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * One could do this more efficiently, but it will have to wait. + */ +Real +atanh_R(Real x) +{ + Real r; + + r = matrix_Int(x, 1, -1, 1, 1); + r = log_R(r); + r = div_R_Int(r, 2); + return r; +} diff --git a/ic-reals-6.3/math-lib/cos_Q.c b/ic-reals-6.3/math-lib/cos_Q.c new file mode 100644 index 0000000..b2c472a --- /dev/null +++ b/ic-reals-6.3/math-lib/cos_Q.c @@ -0,0 +1,46 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +cos_QZ(mpz_t a, mpz_t b) +{ + mpz_t x; + Real r; + + mpz_init(x); + mpz_mul_ui(x, b, 2); + r = tan_QZ(a, x); + mpz_clear(x); + r = tensor_Int(r, r, -1, 1, 0, 0, 0, 0, 1, 1); + return r; +} + +Real +cos_QInt(int a, int b) +{ + Real r; + mpz_t ap, bp; + + /* check for overflow */ + if ((unsigned int)b > 0x3FFFFFFF || (unsigned int)b < 0xC0000000) { + mpz_init_set_si(ap, a); + mpz_init_set_si(bp, b); + mpz_mul_ui(bp, bp, 2); + r = tan_QZ(ap, bp); + mpz_clear(ap); + mpz_clear(bp); + } + else + r = tan_QInt(a, b * 2); + + r = tensor_Int(r, r, -1, 1, 0, 0, 0, 0, 1, 1); + return r; +} diff --git a/ic-reals-6.3/math-lib/cos_Q.c~ b/ic-reals-6.3/math-lib/cos_Q.c~ new file mode 100644 index 0000000..b1a4a69 --- /dev/null +++ b/ic-reals-6.3/math-lib/cos_Q.c~ @@ -0,0 +1,46 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +cos_QZ(mpz_t a, mpz_t b) +{ + mpz_t x; + Real r; + + mpz_init(x); + mpz_mul_ui(x, b, 2); + r = tan_QZ(a, x); + mpz_clear(x); + r = tensor_Int(r, r, -1, 1, 0, 0, 0, 0, 1, 1); + return r; +} + +Real +cos_QInt(int a, int b) +{ + Real r; + mpz_t ap, bp; + + /* check for overflow */ + if (b > 0x3FFFFFFF || b < 0xC0000000) { + mpz_init_set_si(ap, a); + mpz_init_set_si(bp, b); + mpz_mul_ui(bp, bp, 2); + r = tan_QZ(ap, bp); + mpz_clear(ap); + mpz_clear(bp); + } + else + r = tan_QInt(a, b * 2); + + r = tensor_Int(r, r, -1, 1, 0, 0, 0, 0, 1, 1); + return r; +} diff --git a/ic-reals-6.3/math-lib/cos_R.c b/ic-reals-6.3/math-lib/cos_R.c new file mode 100644 index 0000000..9eafe58 --- /dev/null +++ b/ic-reals-6.3/math-lib/cos_R.c @@ -0,0 +1,21 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +cos_R(Real x) +{ + Real r; + + r = div_R_Int(x, 2); + r = tan_R(r); + r = tensor_Int(r, r, -1, 1, 0, 0, 0, 0, 1, 1); + return r; +} diff --git a/ic-reals-6.3/math-lib/cosecant.c b/ic-reals-6.3/math-lib/cosecant.c new file mode 100644 index 0000000..207974b --- /dev/null +++ b/ic-reals-6.3/math-lib/cosecant.c @@ -0,0 +1,98 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These functions should probably be broken into separate files. + */ + +/* + * cosec x = 1 / (sin x) + */ +Real +cosec_R(Real x) +{ + return div_Int_R(1, sin_R(x)); +} + +Real +cosec_QInt(int a, int b) +{ + return div_Int_R(1, sin_QInt(a, b)); +} + +Real +cosec_QZ(mpz_t a, mpz_t b) +{ + return div_Int_R(1, sin_QZ(a, b)); +} + +/* + * acosec x = asin (1/x) + */ +Real +acosec_R(Real x) +{ + return asin_R(div_Int_R(1, x)); +} + +Real +acosec_QInt(int a, int b) +{ + return asin_QInt(b, a); +} + +Real +acosec_QZ(mpz_t a, mpz_t b) +{ + return asin_QZ(b, a); +} + +/* + * cosech x = 1 / (sinh x) + */ +Real +cosech_R(Real x) +{ + return div_Int_R(1, sinh_R(x)); +} + +Real +cosech_QInt(int a, int b) +{ + return div_Int_R(1, sinh_QInt(a, b)); +} + +Real +cosech_QZ(mpz_t a, mpz_t b) +{ + return div_Int_R(1, sinh_QZ(a, b)); +} + +/* + * acosech x = asinh (1/x) + */ +Real +acosech_R(Real x) +{ + return asinh_R(div_Int_R(1, x)); +} + +Real +acosech_QInt(int a, int b) +{ + return asinh_QInt(b, a); +} + +Real +acosech_QZ(mpz_t a, mpz_t b) +{ + return asinh_QZ(b, a); +} diff --git a/ic-reals-6.3/math-lib/cosh_Q.c b/ic-reals-6.3/math-lib/cosh_Q.c new file mode 100644 index 0000000..53ea74d --- /dev/null +++ b/ic-reals-6.3/math-lib/cosh_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +cosh_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return cosh_R(r); +} + +Real +cosh_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return cosh_R(r); +} diff --git a/ic-reals-6.3/math-lib/cosh_R.c b/ic-reals-6.3/math-lib/cosh_R.c new file mode 100644 index 0000000..6dd28b1 --- /dev/null +++ b/ic-reals-6.3/math-lib/cosh_R.c @@ -0,0 +1,19 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +cosh_R(Real x) +{ + Real r; + + r = exp_R(x); + return tensor_Int(r, r, 1, 0, 0, 1, 0, 1, 1, 0); +} diff --git a/ic-reals-6.3/math-lib/cotangent.c b/ic-reals-6.3/math-lib/cotangent.c new file mode 100644 index 0000000..de6ec44 --- /dev/null +++ b/ic-reals-6.3/math-lib/cotangent.c @@ -0,0 +1,98 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These functions should probably be broken into separate files. + */ + +/* + * cotan x = 1 / (tan x) + */ +Real +cotan_R(Real x) +{ + return div_Int_R(1, tan_R(x)); +} + +Real +cotan_QInt(int a, int b) +{ + return div_Int_R(1, tan_QInt(a, b)); +} + +Real +cotan_QZ(mpz_t a, mpz_t b) +{ + return div_Int_R(1, tan_QZ(a, b)); +} + +/* + * acotan x = atan (1/x) + */ +Real +acotan_R(Real x) +{ + return atan_R(div_Int_R(1, x)); +} + +Real +acotan_QInt(int a, int b) +{ + return atan_QInt(b, a); +} + +Real +acotan_QZ(mpz_t a, mpz_t b) +{ + return atan_QZ(b, a); +} + +/* + * cotanh x = 1 / (tanh x) + */ +Real +cotanh_R(Real x) +{ + return div_Int_R(1, tanh_R(x)); +} + +Real +cotanh_QInt(int a, int b) +{ + return div_Int_R(1, tanh_QInt(a, b)); +} + +Real +cotanh_QZ(mpz_t a, mpz_t b) +{ + return div_Int_R(1, tanh_QZ(a, b)); +} + +/* + * acotanh x = atanh (1/x) + */ +Real +acotanh_R(Real x) +{ + return atanh_R(div_Int_R(1, x)); +} + +Real +acotanh_QInt(int a, int b) +{ + return atanh_QInt(b, a); +} + +Real +acotanh_QZ(mpz_t a, mpz_t b) +{ + return atanh_QZ(b, a); +} diff --git a/ic-reals-6.3/math-lib/exp_Q.c b/ic-reals-6.3/math-lib/exp_Q.c new file mode 100644 index 0000000..2e7e183 --- /dev/null +++ b/ic-reals-6.3/math-lib/exp_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +exp_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return exp_R(r); +} + +Real +exp_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return exp_R(r); +} diff --git a/ic-reals-6.3/math-lib/exp_R.c b/ic-reals-6.3/math-lib/exp_R.c new file mode 100644 index 0000000..2a0d9ac --- /dev/null +++ b/ic-reals-6.3/math-lib/exp_R.c @@ -0,0 +1,122 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +static TenXY *nextTensor(Real, Real, int); +static void expInside(); +static void expOutside(); + +Real +exp_R(Real x) +{ + Real u, v; + Bool xLtEq1, xGtEqNeg1, xLtEq1_and_GtEqNeg1; + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(expInside, "expInside", 2); + registerForceFunc(expOutside, "expOutside", 2); + doneInit++; + } + + /* + if (x->gen.tag.type == VECTOR) + return exp_QZ(x->vec.vec[0], x->vec.vec[1]); + */ + + xLtEq1 = ltEq_R_QInt(x, 1, 1); + xGtEqNeg1 = gtEq_R_QInt(x, -1, 1); + xLtEq1_and_GtEqNeg1 = and_B_B(xLtEq1, xGtEqNeg1); + + u = (Real) allocCls(expInside, (void *) x); + u->cls.tag.isSigned = TRUE; + v = (Real) allocCls(expOutside, (void *) x); + v->cls.tag.isSigned = TRUE; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(u, x); + newEdgeToOnlyChild(v, x); + endGraphUpdate(); +#endif + + /* + * The order of the tests in the alt is not semantically + * significant. The tests are applied in order, so there is a + * very modest performance improvement by putting the -1,1 tests + * before the 999/1000,-999/1000 tests. Also there is a very + * small win by doing the negative test before the positive test. + */ + return realIf(4, + not_B(xLtEq1_and_GtEqNeg1), v, + xLtEq1_and_GtEqNeg1, u, + gtEq_R_QInt(x, 999, 1000), v, + ltEq_R_QInt(x, -999, 1000), v); +} + +static TenXY * +nextTensor(Real x, Real y, int n) +{ + return (TenXY *) tensor_Int(x, y, -1, 0, 0, 1, 2*n+1, 2*n+1, 2*n+1, 2*n+1); +} + +static void +expInside() +{ + Cls *cls, *newCls; + Real x; + ClsData *data; + void stdTensorCont(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + if ((data = (ClsData *) malloc(sizeof(ClsData))) == NULL) + Error(FATAL, E_INT, "exp_R", "malloc failed"); + + data->n = 1; + data->x = x; + data->nextTensor = nextTensor; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + cls->redirect = tensor_Int(x, (Real) newCls, 0, -1, 1, 0, 1, 1, 1, 1); + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif +} + +static void +expOutside() +{ + Cls *cls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + x = div_R_Int(x, 2); + + w = exp_R(x); + + cls->redirect = mul_R_R(w, w); + +#ifdef DAVINCI + beginGraphUpdate(); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif +} diff --git a/ic-reals-6.3/math-lib/exp_R.c~ b/ic-reals-6.3/math-lib/exp_R.c~ new file mode 100644 index 0000000..9723dec --- /dev/null +++ b/ic-reals-6.3/math-lib/exp_R.c~ @@ -0,0 +1,123 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +static TenXY *nextTensor(Real, Real, int); +static void expInside(); +static void expOutside(); + +Real +exp_R(Real x) +{ + Real u, v; + Bool xLtEq1, xGtEqNeg1, xLtEq1_and_GtEqNeg1; + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(expInside, "expInside", 2); + registerForceFunc(expOutside, "expOutside", 2); + doneInit++; + } + + /* + if (x->gen.tag.type == VECTOR) + return exp_QZ(x->vec.vec[0], x->vec.vec[1]); + */ + + xLtEq1 = ltEq_R_QInt(x, 1, 1); + xGtEqNeg1 = gtEq_R_QInt(x, -1, 1); + xLtEq1_and_GtEqNeg1 = and_B_B(xLtEq1, xGtEqNeg1); + + u = (Real) allocCls(expInside, (void *) x); + u->cls.tag.isSigned = TRUE; + v = (Real) allocCls(expOutside, (void *) x); + v->cls.tag.isSigned = TRUE; + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(u, x); + newEdgeToOnlyChild(v, x); + endGraphUpdate(); + } + + /* + * The order of the tests in the alt is not semantically + * significant. The tests are applied in order, so there is a + * very modest performance improvement by putting the -1,1 tests + * before the 999/1000,-999/1000 tests. Also there is a very + * small win by doing the negative test before the positive test. + */ + return realIf(4, + not_B(xLtEq1_and_GtEqNeg1), v, + xLtEq1_and_GtEqNeg1, u, + gtEq_R_QInt(x, 999, 1000), v, + ltEq_R_QInt(x, -999, 1000), v); +} + +static TenXY * +nextTensor(Real x, Real y, int n) +{ + return (TenXY *) tensor_Int(x, y, -1, 0, 0, 1, 2*n+1, 2*n+1, 2*n+1, 2*n+1); +} + +static void +expInside() +{ + Cls *cls, *newCls; + Real x, exp_x; + ClsData *data; + void stdTensorCont(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + if ((data = (ClsData *) malloc(sizeof(ClsData))) == NULL) + Error(FATAL, E_INT, "exp_R", "malloc failed"); + + data->n = 1; + data->x = x; + data->nextTensor = nextTensor; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + cls->redirect = tensor_Int(x, (Real) newCls, 0, -1, 1, 0, 1, 1, 1, 1); + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} + +static void +expOutside() +{ + Cls *cls; + Real u, v, w, x; + Bool xLtEq1, xGtEqNeg1, xLtEq1_and_GtEqNeg1; + + cls = (Cls *) POP; + x = (Real) cls->userData; + x = div_R_Int(x, 2); + + w = exp_R(x); + + cls->redirect = mul_R_R(w, w); + + if (DAVINCI) { + beginGraphUpdate(); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} diff --git a/ic-reals-6.3/math-lib/init.c b/ic-reals-6.3/math-lib/init.c new file mode 100644 index 0000000..be3e6ec --- /dev/null +++ b/ic-reals-6.3/math-lib/init.c @@ -0,0 +1,67 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +/* + * Temporary storage for use in the math library + */ +mpz_t tmpw_z, tmpx_z, tmpy_z, tmpz_z; + +Real E; + +void +initReals() +{ + void initRealBase(); + void stdTensorCont(); + void failCls(); + + initRealBase(); + + mpz_init(tmpw_z); + mpz_init(tmpx_z); + mpz_init(tmpy_z); + mpz_init(tmpz_z); + + registerForceFunc(stdTensorCont, "stdTensorCont", 2); + registerForceFunc(failCls, "failCls", 2); + + initPi(); + E = exp_QInt(1, 1); /* there is better way but this will come later */ +} + +Real +realError(char *p) +{ + void failCls(); + Cls *cls; + + cls = allocCls(failCls, (void *) p); + cls->tag.isSigned = FALSE; + + return (Real) cls; +} + +void +failCls() +{ + Cls *cls; + char *p; + + cls = (Cls *) POP; + p = (char *) cls->userData; + + if (p != NULL) + Error(FATAL, E_INT, "failCls", p); + else + Error(FATAL, E_INT, "failCls", "failure closure activated"); +} diff --git a/ic-reals-6.3/math-lib/log_Q.c b/ic-reals-6.3/math-lib/log_Q.c new file mode 100644 index 0000000..0e8d1b7 --- /dev/null +++ b/ic-reals-6.3/math-lib/log_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +log_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return log_R(r); +} + +Real +log_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return log_R(r); +} diff --git a/ic-reals-6.3/math-lib/log_R.c b/ic-reals-6.3/math-lib/log_R.c new file mode 100644 index 0000000..d7266a9 --- /dev/null +++ b/ic-reals-6.3/math-lib/log_R.c @@ -0,0 +1,227 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +void log2Cont(); +void logInside(); +void logLow(); +void logHigh(); + +/* + * This diverges when the argument is zero. + */ +Real +log_R(Real x) +{ + Bool xGtEq0, xLtEq2, xGtEqOneHalf, xLtEq2_and_GtEqOneHalf; + Bool xLtEqOverOneHalf_and_GtEq0; + Real in, low, high, ltZero; + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(logInside, "logInside", 2); + registerForceFunc(logLow, "logLow", 2); + registerForceFunc(logHigh, "logHigh", 2); + registerForceFunc(log2Cont, "log2Cont", 3); + doneInit++; + } + + xGtEq0 = gtEq_R_0(x); + xLtEq2 = ltEq_R_QInt(x, 2, 1); + xGtEqOneHalf = gtEq_R_QInt(x, 1, 2); + xLtEq2_and_GtEqOneHalf = and_B_B(xLtEq2, xGtEqOneHalf); + xLtEqOverOneHalf_and_GtEq0 = and_B_B(ltEq_R_QInt(x, 1001, 2000), xGtEq0); + + in = (Real) allocCls(logInside, (void *) x); + in->cls.tag.isSigned = TRUE; + low = (Real) allocCls(logLow, (void *) x); + low->cls.tag.isSigned = TRUE; + high = (Real) allocCls(logHigh, (void *) x); + high->cls.tag.isSigned = TRUE; + ltZero = realError("(log_R x) and x < 0"); + + return realIf(4, + xLtEq2_and_GtEqOneHalf, in, + gtEq_R_QInt(x, 1999, 1000), high, + xLtEqOverOneHalf_and_GtEq0, low, + not_B(xGtEq0), ltZero); +} + +static TenXY *nextTensor(Real, Real, int); + +void +logInside() +{ + Cls *cls, *newCls; + ClsData *data; + Real x; + void stdTensorCont(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + if ((data = (ClsData *) malloc(sizeof(ClsData))) == NULL) + Error(FATAL, E_INT, "logInside", "malloc failed"); + + data->n = 1; + data->x = x; + data->nextTensor = nextTensor; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + cls->redirect = tensor_Int(x, (Real) newCls, 1, 0, 1, 1, -1, 1, -1, 0); + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif +} + +void +logHigh() +{ + Cls *cls, *newCls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + w = div_R_Int(x, 2); + w = log_R(w); + newCls = allocCls(log2Cont, (void *) 0); + newCls->tag.isSigned = FALSE; + w = add_R_R(w, (Real) newCls); + w->tenXY.forceY = log2Cont; /* see the note below */ + + /* now guard tensor to prevent it being copied */ + cls->redirect = (Real) allocSignX(w, SIGN_UNKN); +} + +void +logLow() +{ + Cls *cls, *newCls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + w = mul_R_Int(x, 2); + w = log_R(w); + newCls = allocCls(log2Cont, (void *) 0); + newCls->tag.isSigned = FALSE; + w = sub_R_R(w, (Real) newCls); + w->tenXY.forceY = log2Cont; /* see the note below */ + + /* now guard tensor to prevent it being copied */ + cls->redirect = (Real) allocSignX(w, SIGN_UNKN); +} + +/* + * This closure works a little different from most. Most simply create + * an lft and set the redirect field of the closure to point to the new + * heap object. In this case, we can do a bit better. This continuation + * represents (log 2). One option would be to follow the usual strategy + * and to arrange for (log 2) to be shared whenever possible. However, + * we know that this is used internally to log_R only and only in + * those circumstances when the argument falls outside the working range + * of the log tensor chain. Moreover, we know that the consumer of this + * real is either an addition or substraction tensor. Whatever it is, + * it is certain to be a tensor. The other argument of the tensor will + * reduce at most to a matrix and hence the tensor itself will never + * reduce. + * + * This property allows us to adopt a different strategy. Rather than create + * the next matrix in the (log 2) chain, we put the matrix directly into + * the consuming tensor. This avoids creating garbage in the stack and + * avoids a separate reduction step. + * + * There is a slighlt better scheme than this. Rather than plunk in + * each successive + * matrix into the tensor (each providing 4 digits), it would be better + * to accumulate a sequence of matrices (up to the capacity of a small matrix) + * and then plunk it into the tensor. This would avoid using bignum stuff until + * it becomes necessary. Perhaps I'll do this later. + * + * For this it is useful to bear in mind that when a sequence of matrices + * starts from n = 0, then the largest entry is always d. When the sequence + * starts with n > 0, then the largest entry is c. To decide if there + * is going to be an overflow, it suffices to check if c (for example) + * will overflow. + */ +static void nextMatrix(Tensor, int); + +void +log2Cont() +{ + TenXY *tenXY; + Cls *cls; + int digitsNeeded; + int n; + + tenXY = (TenXY *) POP; + digitsNeeded = POP; + cls = (Cls *) tenXY->y; + + n = (int) cls->userData; + + while (digitsNeeded > 0) { + nextMatrix(tenXY->ten, n); + if (n == 0) + digitsNeeded -= 2; + else + digitsNeeded -= 4; + n += 1; + } + cls->userData = (void *) n; +} + +/* + * We set a rather arbitrary (but large) limit on the value of n + * for log 2. On a 32 bit machine it is 536,870,911. + * Larger than this and the matrix entries are no longer small. + */ +static void +nextMatrix(Tensor ten, int n) +{ + SmallMatrix smallMat; + + if (n == 0) { + smallMat[0][0] = 1; + smallMat[0][1] = 1; + smallMat[1][0] = 1; + smallMat[1][1] = 2; + } + else { + if (n <= (MAXINT - 2) / 4) { + smallMat[0][0] = 3 * n + 1; + smallMat[0][1] = 2 * n + 1; + smallMat[1][0] = 4 * n + 2; + smallMat[1][1] = 3 * n + 2; + } + else + Error(FATAL, E_INT, "nextMatrix (log2)", + "n > %d", (MAXINT - 2) / 4); + } + multVectorPairTimesSmallMatrix(ten[0], ten[1], smallMat); + multVectorPairTimesSmallMatrix(ten[2], ten[3], smallMat); + normalizeTensor(ten); +} + +static TenXY * +nextTensor(Real x, Real y, int n) +{ + return (TenXY *) tensor_Int(x, y, n, 0, 2*n+1, n+1, n+1, 2*n+1, 0, n); +} diff --git a/ic-reals-6.3/math-lib/log_R.c~ b/ic-reals-6.3/math-lib/log_R.c~ new file mode 100644 index 0000000..bc66045 --- /dev/null +++ b/ic-reals-6.3/math-lib/log_R.c~ @@ -0,0 +1,227 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +void log2Cont(); +void logInside(); +void logLow(); +void logHigh(); + +/* + * This diverges when the argument is zero. + */ +Real +log_R(Real x) +{ + Bool xGtEq0, xLtEq2, xGtEqOneHalf, xLtEq2_and_GtEqOneHalf; + Bool xLtEqOverOneHalf_and_GtEq0; + Real in, low, high, ltZero; + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(logInside, "logInside", 2); + registerForceFunc(logLow, "logLow", 2); + registerForceFunc(logHigh, "logHigh", 2); + registerForceFunc(log2Cont, "log2Cont", 3); + doneInit++; + } + + xGtEq0 = gtEq_R_0(x); + xLtEq2 = ltEq_R_QInt(x, 2, 1); + xGtEqOneHalf = gtEq_R_QInt(x, 1, 2); + xLtEq2_and_GtEqOneHalf = and_B_B(xLtEq2, xGtEqOneHalf); + xLtEqOverOneHalf_and_GtEq0 = and_B_B(ltEq_R_QInt(x, 1001, 2000), xGtEq0); + + in = (Real) allocCls(logInside, (void *) x); + in->cls.tag.isSigned = TRUE; + low = (Real) allocCls(logLow, (void *) x); + low->cls.tag.isSigned = TRUE; + high = (Real) allocCls(logHigh, (void *) x); + high->cls.tag.isSigned = TRUE; + ltZero = realError("(log_R x) and x < 0"); + + return realIf(4, + xLtEq2_and_GtEqOneHalf, in, + gtEq_R_QInt(x, 1999, 1000), high, + xLtEqOverOneHalf_and_GtEq0, low, + not_B(xGtEq0), ltZero); +} + +static TenXY *nextTensor(Real, Real, int); + +void +logInside() +{ + Cls *cls, *newCls; + ClsData *data; + Real x; + void stdTensorCont(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + if ((data = (ClsData *) malloc(sizeof(ClsData))) == NULL) + Error(FATAL, E_INT, "logInside", "malloc failed"); + + data->n = 1; + data->x = x; + data->nextTensor = nextTensor; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + cls->redirect = tensor_Int(x, (Real) newCls, 1, 0, 1, 1, -1, 1, -1, 0); + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} + +void +logHigh() +{ + Cls *cls, *newCls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + w = div_R_Int(x, 2); + w = log_R(w); + newCls = allocCls(log2Cont, (void *) 0); + newCls->tag.isSigned = FALSE; + w = add_R_R(w, (Real) newCls); + w->tenXY.forceY = log2Cont; /* see the note below */ + + /* now guard tensor to prevent it being copied */ + cls->redirect = (Real) allocSignX(w, SIGN_UNKN); +} + +void +logLow() +{ + Cls *cls, *newCls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + w = mul_R_Int(x, 2); + w = log_R(w); + newCls = allocCls(log2Cont, (void *) 0); + newCls->tag.isSigned = FALSE; + w = sub_R_R(w, (Real) newCls); + w->tenXY.forceY = log2Cont; /* see the note below */ + + /* now guard tensor to prevent it being copied */ + cls->redirect = (Real) allocSignX(w, SIGN_UNKN); +} + +/* + * This closure works a little different from most. Most simply create + * an lft and set the redirect field of the closure to point to the new + * heap object. In this case, we can do a bit better. This continuation + * represents (log 2). One option would be to follow the usual strategy + * and to arrange for (log 2) to be shared whenever possible. However, + * we know that this is used internally to log_R only and only in + * those circumstances when the argument falls outside the working range + * of the log tensor chain. Moreover, we know that the consumer of this + * real is either an addition or substraction tensor. Whatever it is, + * it is certain to be a tensor. The other argument of the tensor will + * reduce at most to a matrix and hence the tensor itself will never + * reduce. + * + * This property allows us to adopt a different strategy. Rather than create + * the next matrix in the (log 2) chain, we put the matrix directly into + * the consuming tensor. This avoids creating garbage in the stack and + * avoids a separate reduction step. + * + * There is a slighlt better scheme than this. Rather than plunk in + * each successive + * matrix into the tensor (each providing 4 digits), it would be better + * to accumulate a sequence of matrices (up to the capacity of a small matrix) + * and then plunk it into the tensor. This would avoid using bignum stuff until + * it becomes necessary. Perhaps I'll do this later. + * + * For this it is useful to bear in mind that when a sequence of matrices + * starts from n = 0, then the largest entry is always d. When the sequence + * starts with n > 0, then the largest entry is c. To decide if there + * is going to be an overflow, it suffices to check if c (for example) + * will overflow. + */ +static void nextMatrix(Tensor, int); + +void +log2Cont() +{ + TenXY *tenXY; + Cls *cls; + int digitsNeeded; + int n; + + tenXY = (TenXY *) POP; + digitsNeeded = POP; + cls = (Cls *) tenXY->y; + + n = (int) cls->userData; + + while (digitsNeeded > 0) { + nextMatrix(tenXY->ten, n); + if (n == 0) + digitsNeeded -= 2; + else + digitsNeeded -= 4; + n += 1; + } + cls->userData = (void *) n; +} + +/* + * We set a rather arbitrary (but large) limit on the value of n + * for log 2. On a 32 bit machine it is 536,870,911. + * Larger than this and the matrix entries are no longer small. + */ +static void +nextMatrix(Tensor ten, int n) +{ + SmallMatrix smallMat; + + if (n == 0) { + smallMat[0][0] = 1; + smallMat[0][1] = 1; + smallMat[1][0] = 1; + smallMat[1][1] = 2; + } + else { + if (n <= (MAXINT - 2) / 4) { + smallMat[0][0] = 3 * n + 1; + smallMat[0][1] = 2 * n + 1; + smallMat[1][0] = 4 * n + 2; + smallMat[1][1] = 3 * n + 2; + } + else + Error(FATAL, E_INT, "nextMatrix (log2)", + "n > %d", (MAXINT - 2) / 4); + } + multVectorPairTimesSmallMatrix(ten[0], ten[1], smallMat); + multVectorPairTimesSmallMatrix(ten[2], ten[3], smallMat); + normalizeTensor(ten); +} + +static TenXY * +nextTensor(Real x, Real y, int n) +{ + return (TenXY *) tensor_Int(x, y, n, 0, 2*n+1, n+1, n+1, 2*n+1, 0, n); +} diff --git a/ic-reals-6.3/math-lib/math-lib.h b/ic-reals-6.3/math-lib/math-lib.h new file mode 100644 index 0000000..51c25be --- /dev/null +++ b/ic-reals-6.3/math-lib/math-lib.h @@ -0,0 +1,18 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +/* + * This is the structure holding the closure data for those functions + * which use the standard tensor continuation. + */ + +typedef struct { + int n; + Real x; + TenXY *(*nextTensor)(Real, Real, int); +} ClsData; diff --git a/ic-reals-6.3/math-lib/neg_R.c b/ic-reals-6.3/math-lib/neg_R.c new file mode 100644 index 0000000..3ae644c --- /dev/null +++ b/ic-reals-6.3/math-lib/neg_R.c @@ -0,0 +1,16 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +neg_R(Real x) +{ + return mul_R_Int(x, -1); +} diff --git a/ic-reals-6.3/math-lib/pi.c b/ic-reals-6.3/math-lib/pi.c new file mode 100644 index 0000000..15c32d2 --- /dev/null +++ b/ic-reals-6.3/math-lib/pi.c @@ -0,0 +1,151 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" + +static mpz_t k; + +Real Pi; + +static void piCont(); +static void nextMatrix(Tensor, int); + +void +initPi() +{ + Real sq, cls; + DigsX *digsX; + TenXY *tenXY; + static int doneInit; + void force_To_DigsX_From_TenXY_Entry(); + + if (!doneInit) { + registerForceFunc(piCont, "piCont", 2); + doneInit++; + } + + mpz_init_set_str(k, "10939058860032000", 10); + + sq = sqrt_QInt(10005, 1); + cls = (Real) allocCls(piCont, (void *) 0); + cls->gen.tag.isSigned = FALSE; + + tenXY = (TenXY *) div_R_R(sq, cls); + tenXY->forceY = piCont; + + digsX = allocDigsX(); + digsX->x = (Real) tenXY; + digsX->force = force_To_DigsX_From_TenXY_Entry; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, tenXY); + endGraphUpdate(); +#endif + + Pi = (Real) digsX; +} + +/* + * It might be possible compute the next matrix in the sequence + * in terms of the previous matrix. This would mean having much more + * state but presumably much quicker. + */ +static void +piCont() +{ + TenXY *tenXY; + Cls *cls; + int digitsNeeded; + int n; + + tenXY = (TenXY *) POP; + digitsNeeded = (int) POP; + + cls = (Cls *) tenXY->y; + n = (int) cls->userData; + + /* + * The two constants below have been determined empirically. The second + * may be off by one (ie 47) in some circumstances. + */ + while (digitsNeeded > 0) { + nextMatrix(tenXY->ten, n); + if (n == 0) + digitsNeeded -= 27; + else + digitsNeeded -= 46; + n +=1; + } + cls->userData = (void *) n; +} + +static void +nextMatrix(Tensor ten, int n) +{ + SmallMatrix smallMat; + + if (n == 0) { + smallMat[0][0] = 6795705; + smallMat[0][1] = 213440; + smallMat[1][0] = 6795704; + smallMat[1][1] = 213440; + multVectorPairTimesSmallMatrix(ten[0], ten[1], smallMat); + multVectorPairTimesSmallMatrix(ten[2], ten[3], smallMat); + } + else { + /* b = 2n - 1 */ + mpz_set_ui(tmpb_z, n); + mpz_mul_2exp(tmpb_z, tmpb_z, 1); + mpz_sub_ui(tmpb_z, tmpb_z, 1); + + /* tmp = 6n - 5, tmpf_z = 6n */ + mpz_set_ui(tmpa_z, n); + mpz_mul_ui(tmpa_z, tmpa_z, 6); + mpz_set(tmpf_z, tmpa_z); + mpz_sub_ui(tmpa_z, tmpa_z, 5); + + /* tmpf_z = 6n - 1 */ + mpz_sub_ui(tmpf_z, tmpf_z, 1); + + /* b = (2n - 1) * (6n - 5) * (6n - 1) */ + mpz_mul(tmpb_z, tmpb_z, tmpa_z); + mpz_mul(tmpb_z, tmpb_z, tmpf_z); + + /* c = b*(545140134n + 13591409) */ + mpz_set_ui(tmpa_z, 545140134); + mpz_mul_ui(tmpa_z, tmpa_z, n); + mpz_add_ui(tmpa_z, tmpa_z, 13591409); + mpz_mul(tmpc_z, tmpb_z, tmpa_z); + + /* d = b * (n + 1) */ + mpz_set_ui(tmpd_z, n + 1); + mpz_mul(tmpd_z, tmpd_z, tmpb_z); + + /* e = 10939058860032000*n^4 */ + mpz_set_ui(tmpe_z, n); + mpz_mul(tmpe_z, tmpe_z, tmpe_z); + mpz_mul(tmpe_z, tmpe_z, tmpe_z); + mpz_mul(tmpe_z, tmpe_z, k); + + mpz_sub(bigTmpMat[0][0], tmpe_z, tmpd_z); + mpz_set(bigTmpMat[1][1], bigTmpMat[0][0]); + mpz_add(bigTmpMat[0][1], tmpe_z, tmpd_z); + mpz_set(bigTmpMat[1][0], bigTmpMat[0][1]); + + mpz_sub(bigTmpMat[0][0], bigTmpMat[0][0], tmpc_z); + mpz_add(bigTmpMat[0][1], bigTmpMat[0][1], tmpc_z); + mpz_sub(bigTmpMat[1][0], bigTmpMat[1][0], tmpc_z); + mpz_add(bigTmpMat[1][1], bigTmpMat[1][1], tmpc_z); + + multVectorPairTimesMatrix(ten[0], ten[1], bigTmpMat); + multVectorPairTimesMatrix(ten[2], ten[3], bigTmpMat); + } +} diff --git a/ic-reals-6.3/math-lib/pi.c~ b/ic-reals-6.3/math-lib/pi.c~ new file mode 100644 index 0000000..a95eef5 --- /dev/null +++ b/ic-reals-6.3/math-lib/pi.c~ @@ -0,0 +1,151 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" + +static mpz_t k; + +Real Pi; + +static void piCont(); +static void nextMatrix(Tensor, int); + +void +initPi() +{ + Real sq, cls; + DigsX *digsX; + TenXY *tenXY; + static int doneInit; + void force_To_DigsX_From_TenXY_Entry(); + + if (!doneInit) { + registerForceFunc(piCont, "piCont", 2); + doneInit++; + } + + mpz_init_set_str(k, "10939058860032000", 10); + + sq = sqrt_QInt(10005, 1); + cls = (Real) allocCls(piCont, (void *) 0); + cls->gen.tag.isSigned = FALSE; + + tenXY = (TenXY *) div_R_R(sq, cls); + tenXY->forceY = piCont; + + digsX = allocDigsX(); + digsX->x = (Real) tenXY; + digsX->force = force_To_DigsX_From_TenXY_Entry; + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, tenXY); + endGraphUpdate(); + } + + Pi = (Real) digsX; +} + +/* + * It might be possible compute the next matrix in the sequence + * in terms of the previous matrix. This would mean having much more + * state but presumably much quicker. + */ +static void +piCont() +{ + TenXY *tenXY; + Cls *cls; + int digitsNeeded; + int n; + + tenXY = (TenXY *) POP; + digitsNeeded = (int) POP; + + cls = (Cls *) tenXY->y; + n = (int) cls->userData; + + /* + * The two constants below have been determined empirically. The second + * may be off by one (ie 47) in some circumstances. + */ + while (digitsNeeded > 0) { + nextMatrix(tenXY->ten, n); + if (n == 0) + digitsNeeded -= 27; + else + digitsNeeded -= 46; + n +=1; + } + cls->userData = (void *) n; +} + +static void +nextMatrix(Tensor ten, int n) +{ + SmallMatrix smallMat; + + if (n == 0) { + smallMat[0][0] = 6795705; + smallMat[0][1] = 213440; + smallMat[1][0] = 6795704; + smallMat[1][1] = 213440; + multVectorPairTimesSmallMatrix(ten[0], ten[1], smallMat); + multVectorPairTimesSmallMatrix(ten[2], ten[3], smallMat); + } + else { + /* b = 2n - 1 */ + mpz_set_ui(tmpb_z, n); + mpz_mul_2exp(tmpb_z, tmpb_z, 1); + mpz_sub_ui(tmpb_z, tmpb_z, 1); + + /* tmp = 6n - 5, tmpf_z = 6n */ + mpz_set_ui(tmpa_z, n); + mpz_mul_ui(tmpa_z, tmpa_z, 6); + mpz_set(tmpf_z, tmpa_z); + mpz_sub_ui(tmpa_z, tmpa_z, 5); + + /* tmpf_z = 6n - 1 */ + mpz_sub_ui(tmpf_z, tmpf_z, 1); + + /* b = (2n - 1) * (6n - 5) * (6n - 1) */ + mpz_mul(tmpb_z, tmpb_z, tmpa_z); + mpz_mul(tmpb_z, tmpb_z, tmpf_z); + + /* c = b*(545140134n + 13591409) */ + mpz_set_ui(tmpa_z, 545140134); + mpz_mul_ui(tmpa_z, tmpa_z, n); + mpz_add_ui(tmpa_z, tmpa_z, 13591409); + mpz_mul(tmpc_z, tmpb_z, tmpa_z); + + /* d = b * (n + 1) */ + mpz_set_ui(tmpd_z, n + 1); + mpz_mul(tmpd_z, tmpd_z, tmpb_z); + + /* e = 10939058860032000*n^4 */ + mpz_set_ui(tmpe_z, n); + mpz_mul(tmpe_z, tmpe_z, tmpe_z); + mpz_mul(tmpe_z, tmpe_z, tmpe_z); + mpz_mul(tmpe_z, tmpe_z, k); + + mpz_sub(bigTmpMat[0][0], tmpe_z, tmpd_z); + mpz_set(bigTmpMat[1][1], bigTmpMat[0][0]); + mpz_add(bigTmpMat[0][1], tmpe_z, tmpd_z); + mpz_set(bigTmpMat[1][0], bigTmpMat[0][1]); + + mpz_sub(bigTmpMat[0][0], bigTmpMat[0][0], tmpc_z); + mpz_add(bigTmpMat[0][1], bigTmpMat[0][1], tmpc_z); + mpz_sub(bigTmpMat[1][0], bigTmpMat[1][0], tmpc_z); + mpz_add(bigTmpMat[1][1], bigTmpMat[1][1], tmpc_z); + + multVectorPairTimesMatrix(ten[0], ten[1], bigTmpMat); + multVectorPairTimesMatrix(ten[2], ten[3], bigTmpMat); + } +} diff --git a/ic-reals-6.3/math-lib/pow_R_R.c b/ic-reals-6.3/math-lib/pow_R_R.c new file mode 100644 index 0000000..9da7446 --- /dev/null +++ b/ic-reals-6.3/math-lib/pow_R_R.c @@ -0,0 +1,20 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +pow_R_R(Real x, Real y) +{ + Real r; + + r = log_R(x); + r = mul_R_R(y, r); + return exp_R(r); +} diff --git a/ic-reals-6.3/math-lib/secant.c b/ic-reals-6.3/math-lib/secant.c new file mode 100644 index 0000000..7261546 --- /dev/null +++ b/ic-reals-6.3/math-lib/secant.c @@ -0,0 +1,98 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These functions should probably be broken into separate files. + */ + +/* + * sec x = 1 / (cos x) + */ +Real +sec_R(Real x) +{ + return div_Int_R(1, cos_R(x)); +} + +Real +sec_QInt(int a, int b) +{ + return div_Int_R(1, cos_QInt(a, b)); +} + +Real +sec_QZ(mpz_t a, mpz_t b) +{ + return div_Int_R(1, cos_QZ(a, b)); +} + +/* + * asec x = acos (1/x) + */ +Real +asec_R(Real x) +{ + return acos_R(div_Int_R(1, x)); +} + +Real +asec_QInt(int a, int b) +{ + return acos_QInt(b, a); +} + +Real +asec_QZ(mpz_t a, mpz_t b) +{ + return acos_QZ(b, a); +} + +/* + * sech x = 1 / (cosh x) + */ +Real +sech_R(Real x) +{ + return div_Int_R(1, cosh_R(x)); +} + +Real +sech_QInt(int a, int b) +{ + return div_Int_R(1, cosh_QInt(a, b)); +} + +Real +sech_QZ(mpz_t a, mpz_t b) +{ + return div_Int_R(1, cosh_QZ(a, b)); +} + +/* + * asech x = acosh (1/x) + */ +Real +asech_R(Real x) +{ + return acosh_R(div_Int_R(1, x)); +} + +Real +asech_QInt(int a, int b) +{ + return acosh_QInt(b, a); +} + +Real +asech_QZ(mpz_t a, mpz_t b) +{ + return acosh_QZ(b, a); +} diff --git a/ic-reals-6.3/math-lib/sin_Q.c b/ic-reals-6.3/math-lib/sin_Q.c new file mode 100644 index 0000000..3e40ee7 --- /dev/null +++ b/ic-reals-6.3/math-lib/sin_Q.c @@ -0,0 +1,46 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +sin_QZ(mpz_t a, mpz_t b) +{ + mpz_t x; + Real r; + + mpz_init(x); + mpz_mul_ui(x, b, 2); + r = tan_QZ(a, x); + mpz_clear(x); + r = tensor_Int(r, r, 0, 1, 1, 0, 1, 0, 0, 1); + return r; +} + +Real +sin_QInt(int a, int b) +{ + Real r; + mpz_t ap, bp; + + /* check for overflow */ + if ((unsigned int)b > 0x3FFFFFFF || (unsigned int)b < 0xC0000000) { + mpz_init_set_si(ap, a); + mpz_init_set_si(bp, b); + mpz_mul_ui(bp, bp, 2); + r = tan_QZ(ap, bp); + mpz_clear(ap); + mpz_clear(bp); + } + else + r = tan_QInt(a, b * 2); + + r = tensor_Int(r, r, 0, 1, 1, 0, 1, 0, 0, 1); + return r; +} diff --git a/ic-reals-6.3/math-lib/sin_Q.c~ b/ic-reals-6.3/math-lib/sin_Q.c~ new file mode 100644 index 0000000..171c641 --- /dev/null +++ b/ic-reals-6.3/math-lib/sin_Q.c~ @@ -0,0 +1,46 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +sin_QZ(mpz_t a, mpz_t b) +{ + mpz_t x; + Real r; + + mpz_init(x); + mpz_mul_ui(x, b, 2); + r = tan_QZ(a, x); + mpz_clear(x); + r = tensor_Int(r, r, 0, 1, 1, 0, 1, 0, 0, 1); + return r; +} + +Real +sin_QInt(int a, int b) +{ + Real r; + mpz_t ap, bp; + + /* check for overflow */ + if (b > 0x3FFFFFFF || b < 0xC0000000) { + mpz_init_set_si(ap, a); + mpz_init_set_si(bp, b); + mpz_mul_ui(bp, bp, 2); + r = tan_QZ(ap, bp); + mpz_clear(ap); + mpz_clear(bp); + } + else + r = tan_QInt(a, b * 2); + + r = tensor_Int(r, r, 0, 1, 1, 0, 1, 0, 0, 1); + return r; +} diff --git a/ic-reals-6.3/math-lib/sin_R.c b/ic-reals-6.3/math-lib/sin_R.c new file mode 100644 index 0000000..6b68c6d --- /dev/null +++ b/ic-reals-6.3/math-lib/sin_R.c @@ -0,0 +1,21 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +sin_R(Real x) +{ + Real r; + + r = div_R_Int(x, 2); + r = tan_R(r); + r = tensor_Int(r, r, 0, 1, 1, 0, 1, 0, 0, 1); + return r; +} diff --git a/ic-reals-6.3/math-lib/sinh_Q.c b/ic-reals-6.3/math-lib/sinh_Q.c new file mode 100644 index 0000000..3870b80 --- /dev/null +++ b/ic-reals-6.3/math-lib/sinh_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +sinh_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return sinh_R(r); +} + +Real +sinh_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return sinh_R(r); +} diff --git a/ic-reals-6.3/math-lib/sinh_R.c b/ic-reals-6.3/math-lib/sinh_R.c new file mode 100644 index 0000000..c5198f1 --- /dev/null +++ b/ic-reals-6.3/math-lib/sinh_R.c @@ -0,0 +1,19 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +sinh_R(Real x) +{ + Real r; + + r = exp_R(x); + return tensor_Int(r, r, 1, 0, 0, 1, 0, 1, -1, 0); +} diff --git a/ic-reals-6.3/math-lib/sqrt_Q.c b/ic-reals-6.3/math-lib/sqrt_Q.c new file mode 100644 index 0000000..688e146 --- /dev/null +++ b/ic-reals-6.3/math-lib/sqrt_Q.c @@ -0,0 +1,140 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" + +/* + * It would be better to use a union of three int's and three mpz_t's and + * then go to the bignumns when we overflow a machine word. + */ +typedef struct { + mpz_t a, b, c; +} RatSqrtData; + +static int doneInit = 0; + +static void ratSqrtCont(); +static bool emitDigitFromRatSqrt(RatSqrtData *, Digit *); + +Real +sqrt_QInt(int a, int b) +{ + RatSqrtData *rsd; + Cls *cls; + DigsX *digsX; + + if (!doneInit) { + registerForceFunc(ratSqrtCont, "ratSqrtCont", 3); + doneInit++; + } + + if ((rsd = (RatSqrtData *) malloc(sizeof(RatSqrtData))) == NULL) + Error(FATAL, E_INT, "sqrt_QInt", "malloc failed"); + + mpz_init_set_ui(rsd->a, a); + mpz_init_set_ui(rsd->b, b); + mpz_init_set_ui(rsd->c, a - b); + + cls = allocCls(ratSqrtCont, (void *) rsd); + cls->tag.isSigned = FALSE; + + digsX = allocDigsX(); + digsX->x = (Real) cls; + digsX->force = ratSqrtCont; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, cls); + endGraphUpdate(); +#endif + + return (Real) digsX; +} + +Real +sqrt_QZ(mpz_t a, mpz_t b) +{ + RatSqrtData *rsd; + Cls *cls; + DigsX *digsX; + + if (!doneInit) { + registerForceFunc(ratSqrtCont, "ratSqrtCont", 3); + doneInit++; + } + + if ((rsd = (RatSqrtData *) malloc(sizeof(RatSqrtData))) == NULL) + Error(FATAL, E_INT, "sqrt_QZ", "malloc failed"); + + mpz_init_set(rsd->a, a); + mpz_init_set(rsd->b, b); + mpz_init(rsd->c); + mpz_sub(rsd->c, rsd->a, rsd->b); + + cls = allocCls(ratSqrtCont, (void *) rsd); + cls->tag.isSigned = FALSE; + + digsX = allocDigsX(); + digsX->x = (Real) cls; + digsX->force = ratSqrtCont; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, cls); + endGraphUpdate(); +#endif + + return (Real) digsX; +} + +static void +ratSqrtCont() +{ + DigsX *digsX; + Cls *cls; + int digitsNeeded; + RatSqrtData *rsd; + + digsX = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) digsX->x; + + rsd = (RatSqrtData *) cls->userData; + + emitDigits(digsX, (edf) emitDigitFromRatSqrt, (void *) rsd, digitsNeeded); + newDigsX(digsX); +} + +static bool +emitDigitFromRatSqrt(RatSqrtData *rsd, Digit *d) +{ + mpz_sub(tmpd_z, rsd->b, rsd->a); + mpz_mul_2exp(tmpd_z, tmpd_z, 1); + mpz_add(tmpd_z, tmpd_z, rsd->c); + + switch (mpz_sgn(tmpd_z)) { + case 0 : + case 1 : + *d = DNEG; + mpz_mul_2exp(rsd->a, rsd->a, 2); /* a = 4a */ + mpz_set(rsd->b, tmpd_z); + break; + case -1 : + *d = DPOS; + mpz_neg(rsd->a, tmpd_z); + mpz_mul_2exp(rsd->b, rsd->b, 2); /* b = 4b */ + break; + default : + Error(FATAL, E_INT, "emitDigitFromRatSqrt", + "bad value returned by mpz_sgn"); + break; + } + return TRUE; +} diff --git a/ic-reals-6.3/math-lib/sqrt_Q.c~ b/ic-reals-6.3/math-lib/sqrt_Q.c~ new file mode 100644 index 0000000..0a59773 --- /dev/null +++ b/ic-reals-6.3/math-lib/sqrt_Q.c~ @@ -0,0 +1,140 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" + +/* + * It would be better to use a union of three int's and three mpz_t's and + * then go to the bignumns when we overflow a machine word. + */ +typedef struct { + mpz_t a, b, c; +} RatSqrtData; + +static int doneInit = 0; + +static void ratSqrtCont(); +static bool emitDigitFromRatSqrt(RatSqrtData *, Digit *); + +Real +sqrt_QInt(int a, int b) +{ + RatSqrtData *rsd; + Cls *cls; + DigsX *digsX; + + if (!doneInit) { + registerForceFunc(ratSqrtCont, "ratSqrtCont", 3); + doneInit++; + } + + if ((rsd = (RatSqrtData *) malloc(sizeof(RatSqrtData))) == NULL) + Error(FATAL, E_INT, "sqrt_QInt", "malloc failed"); + + mpz_init_set_ui(rsd->a, a); + mpz_init_set_ui(rsd->b, b); + mpz_init_set_ui(rsd->c, a - b); + + cls = allocCls(ratSqrtCont, (void *) rsd); + cls->tag.isSigned = FALSE; + + digsX = allocDigsX(); + digsX->x = (Real) cls; + digsX->force = ratSqrtCont; + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, cls); + endGraphUpdate(); + } + + return (Real) digsX; +} + +Real +sqrt_QZ(mpz_t a, mpz_t b) +{ + RatSqrtData *rsd; + Cls *cls; + DigsX *digsX; + + if (!doneInit) { + registerForceFunc(ratSqrtCont, "ratSqrtCont", 3); + doneInit++; + } + + if ((rsd = (RatSqrtData *) malloc(sizeof(RatSqrtData))) == NULL) + Error(FATAL, E_INT, "sqrt_QZ", "malloc failed"); + + mpz_init_set(rsd->a, a); + mpz_init_set(rsd->b, b); + mpz_init(rsd->c); + mpz_sub(rsd->c, rsd->a, rsd->b); + + cls = allocCls(ratSqrtCont, (void *) rsd); + cls->tag.isSigned = FALSE; + + digsX = allocDigsX(); + digsX->x = (Real) cls; + digsX->force = ratSqrtCont; + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, cls); + endGraphUpdate(); + } + + return (Real) digsX; +} + +static void +ratSqrtCont() +{ + DigsX *digsX; + Cls *cls; + int digitsNeeded; + RatSqrtData *rsd; + + digsX = (DigsX *) POP; + digitsNeeded = (int) POP; + cls = (Cls *) digsX->x; + + rsd = (RatSqrtData *) cls->userData; + + emitDigits(digsX, (edf) emitDigitFromRatSqrt, (void *) rsd, digitsNeeded); + newDigsX(digsX); +} + +static bool +emitDigitFromRatSqrt(RatSqrtData *rsd, Digit *d) +{ + mpz_sub(tmpd_z, rsd->b, rsd->a); + mpz_mul_2exp(tmpd_z, tmpd_z, 1); + mpz_add(tmpd_z, tmpd_z, rsd->c); + + switch (mpz_sgn(tmpd_z)) { + case 0 : + case 1 : + *d = DNEG; + mpz_mul_2exp(rsd->a, rsd->a, 2); /* a = 4a */ + mpz_set(rsd->b, tmpd_z); + break; + case -1 : + *d = DPOS; + mpz_neg(rsd->a, tmpd_z); + mpz_mul_2exp(rsd->b, rsd->b, 2); /* b = 4b */ + break; + default : + Error(FATAL, E_INT, "emitDigitFromRatSqrt", + "bad value returned by mpz_sgn"); + break; + } + return TRUE; +} diff --git a/ic-reals-6.3/math-lib/sqrt_R.c b/ic-reals-6.3/math-lib/sqrt_R.c new file mode 100644 index 0000000..3333949 --- /dev/null +++ b/ic-reals-6.3/math-lib/sqrt_R.c @@ -0,0 +1,469 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" + +Real +sqrt_R(Real x) +{ + Bool xLtEq4, xGtEqOneFifth; + Cls *in, *gtEq3, *ltEqOneQuarter; + Real lt0; + void force_To_DigsX_From_Sqrt_TenXY_Entry(); + void sqrtInside(); + void sqrtGtEq3(); + void sqrtLtEqOneQuarter(); + void sqrtLtEqOneQuarterCont(); + void force_To_DigsX_From_Sqrt_TenXY_Entry(); + void force_To_DigsX_From_Sqrt_TenXY_Cont(); + void force_To_DigsX_From_Sqrt_MatX_Entry(); + void force_To_DigsX_From_Sqrt_MatX_Cont(); + void force_To_DigsX_From_Sqrt_Reduce(); + void force_To_MatX_Until_Refining(); + static int doneInit = 0; + Real sqrt_QZ(mpz_t, mpz_t); + + if (!doneInit) { + registerForceFunc(force_To_DigsX_From_Sqrt_TenXY_Entry, + "force_To_DigsX_From_Sqrt_TenXY_Entry", 3); + registerForceFunc(force_To_DigsX_From_Sqrt_TenXY_Cont, + "force_To_DigsX_From_Sqrt_TenXY_Cont", 3); + registerForceFunc(force_To_DigsX_From_Sqrt_MatX_Entry, + "force_To_DigsX_From_Sqrt_MatX_Entry", 3); + registerForceFunc(force_To_DigsX_From_Sqrt_MatX_Cont, + "force_To_DigsX_From_Sqrt_MatX_Cont", 3); + registerForceFunc(force_To_DigsX_From_Sqrt_Reduce, + "force_To_DigsX_From_Sqrt_Reduce", 2); + registerForceFunc(force_To_MatX_Until_Refining, + "force_To_MatX_Until_Refining", 2); + registerForceFunc(sqrtInside, "sqrtInside", 2); + registerForceFunc(sqrtGtEq3, "sqrtGtEq3", 2); + registerForceFunc(sqrtLtEqOneQuarter, "sqrtLtEqOneQuarter", 2); + registerForceFunc(sqrtLtEqOneQuarterCont, "sqrtLtEqOneQuarterCont", 2); + doneInit++; + } + + if (x->gen.tag.type == VECTOR) + return sqrt_QZ(x->vec.vec[0], x->vec.vec[1]); + + xLtEq4 = ltEq_R_QInt(x, 4, 1); + xGtEqOneFifth = gtEq_R_QInt(x, 1, 5); + + in = allocCls(sqrtInside, (void *) x); + in->tag.isSigned = FALSE; + + gtEq3 = allocCls(sqrtGtEq3, (void *) x); + gtEq3->tag.isSigned = FALSE; + + ltEqOneQuarter = allocCls(sqrtLtEqOneQuarter, (void *) x); + ltEqOneQuarter->tag.isSigned = FALSE; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(in, x); + newEdgeToOnlyChild(gtEq3, x); + newEdgeToOnlyChild(ltEqOneQuarter, x); + endGraphUpdate(); +#endif + + lt0 = realError("(sqrt_R x) and x < 0"); + + /* + * The order of the tests in the realAlt is not semantically + * significant. Note that the sqrt_R is always unsigned since + * each case is unsigned. + */ + return realIf(4, + and_B_B(xLtEq4, xGtEqOneFifth), (Real) in, + gtEq_R_QInt(x, 3, 1), (Real) gtEq3, + and_B_B(ltEq_R_QInt(x, 1, 4), gtEq_R_0(x)), (Real) ltEqOneQuarter, + lt_R_0(x), lt0); +} + +void +sqrtInside() +{ + Cls *cls; + TenXY *tenXY; + DigsX *digsX; + void force_To_TenXY_X_Until_Refining(); + void force_To_DigsX_From_Sqrt_TenXY_Entry(); + Real x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + digsX = allocDigsX(); + digsX->force = force_To_DigsX_From_Sqrt_TenXY_Entry; + + tenXY = (TenXY *) tensor_Int(x, (Real) digsX, 1, 0, 2, 1, 1, 2, 0, 1); + tenXY->tag.isSigned = FALSE; + + digsX->x = (Real) tenXY; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, tenXY); + endGraphUpdate(); +#endif + + /* + * We must still absorb the sign (if any) and enough information + * to ensure the tensor is refining. Note that, if we get to + * this point, then we know the argument is within the interval + * 1/5 >= x <= 4. I claim that the n'th tensor for n >= 2 + * can always be made refining by absorbing + * information from the left. + */ + PUSH_2(force_To_TenXY_X_Until_Refining, tenXY); + + if (tenXY->x->gen.tag.isSigned) + PUSH_2(tenXY->forceX, tenXY); + + cls->redirect = (Real) digsX; + cls->userData = NULL; + +#ifdef DAVINCI + beginGraphUpdate(); + deleteOnlyEdge(cls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif +} + +void +sqrtGtEq3() +{ + Cls *cls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + w = div_R_Int(x, 4); + w = sqrt_R(w); + cls->redirect = mul_R_Int(w, 2); + cls->userData = NULL; + +#ifdef DAVINCI + beginGraphUpdate(); + deleteOnlyEdge(cls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif +} + +void +sqrtLtEqOneQuarter() +{ + Cls *cls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + w = mul_R_Int(x, 4); + w = sqrt_R(w); + cls->redirect = div_R_Int(w, 2); + cls->userData = NULL; + +#ifdef DAVINCI + beginGraphUpdate(); + deleteOnlyEdge(cls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif +} + +#ifdef LATER +void +sqrtLtEqOneQuarter() +{ + Cls *cls, *newCls; + void sqrtLtEqOneQuarterCont(); + Real x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + newCls = allocCls(sqrtLtEqOneQuarterCont, (void *) x); + newCls->tag.isSigned = FALSE; + cls->redirect = matrix_Int((Real) newCls, 0, 1, 1, 2); + cls->userData = NULL; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + deleteOnlyEdge(cls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif +} +#endif + +void +sqrtLtEqOneQuarterCont() +{ + Cls *cls; + Real w, x; + void force_To_MatX_Until_Refining(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + /* w = matrix_Int(x, 0, 1, 1, 0); reciprocal of x */ + w = mul_R_Int(x, 4); + w = sqrt_R(w); + /* cls->redirect = matrix_Int(w, 1, 0, -2, 1); */ + cls->redirect = w; + cls->redirect->gen.tag.isSigned = FALSE; + cls->userData = NULL; + +#ifdef DAVINCI + beginGraphUpdate(); + deleteOnlyEdge(cls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif + + PUSH_2(force_To_MatX_Until_Refining, cls->redirect); +} + +void +force_To_DigsX_From_Sqrt_TenXY_Entry() +{ + DigsX *digsX; + int digitsNeeded; + void force_To_DigsX_From_Sqrt_TenXY_Cont(); + void force_To_DigsX_From_DigsX_Entry(); + void force_To_DigsX_From_Sqrt_Reduce(); + + digsX = (DigsX *) POP; + digitsNeeded = (int) POP; + + /* + * The current strategy is that when forced, we create a list + * of DigsX structures rather than a single struct with the + * requested number of digits. So, after we are done emitting, + * we arrange to reduce the list at the end. + * + * A better strategy would be to allocate two DigsX structures at the + * very start. The second one would be private and only visible to the sqrt + * closure. Then the cycle would be: + * emit from tensor into DigsX-2 + * absorb DigsX-2 digits into tensor + * absorb DigsX-2 digits into DigsX-1 + * clear DigsX-2 of all digits. + * + * This would avoid allocating a chain of DigsX structures and the + * need for reducing at the end. + PUSH_3(force_To_DigsX_From_DigsX_Entry, digsX, digitsNeeded); + */ + PUSH_2(force_To_DigsX_From_Sqrt_Reduce, digsX); + PUSH_3(force_To_DigsX_From_Sqrt_TenXY_Cont, digsX, digitsNeeded); +} + +void +force_To_DigsX_From_Sqrt_TenXY_Cont() +{ + DigsX *digsX; + TenXY *tenXY; + int digitsNeeded = 0; + int nX; + int digitsEmitted; + int epsDelTensorX(Tensor, int); + bool emitDigitFromTensor(Tensor, Digit *); + + digsX = (DigsX *) POP; + digitsNeeded = (int) POP; + tenXY = (TenXY *) digsX->x; + + digitsEmitted = emitDigits(digsX, + (edf) emitDigitFromTensor, + (void *) tenXY->ten, + digitsNeeded); + +#ifdef TRACE + int bitsShifted = 0; + if (digitsEmitted > 0) + bitsShifted = normalizeTensor(tenXY->ten); + + if (TRACE) { + debugp("force_To_DigsX_From_Sqrt_TenXY_Cont", + "%x %x emitted=%d shifted=%d\n", + (unsigned) digsX, + (unsigned) tenXY, + digitsEmitted, + bitsShifted); + } +#endif + + digitsNeeded -= digitsEmitted; + + if (digsX->count > 0) + newDigsX(digsX); + + if (digitsNeeded <= 0) + return; + + /* + * So now we emitted what we can but still need more. First arrange + * to come back and try to emit again after forcing the necessary + * number of digits from the the argument. + */ + if (digsX->count > 0) + PUSH_3(force_To_DigsX_From_Sqrt_TenXY_Cont, digsX->x, digitsNeeded); + else + PUSH_3(force_To_DigsX_From_Sqrt_TenXY_Cont, digsX, digitsNeeded); + + /* + * Now absorb everything emitted into the tensor. + */ + absorbDigsXIntoTenXY_Y(tenXY); + + nX = epsDelTensorX(tenXY->ten, digitsNeeded); +#ifdef TRACE + if (TRACE) { + debugp("force_To_DigsX_From_Sqrt_TenXY_Cont", + "%x %x nX=%d\n", + (unsigned) digsX, + (unsigned) tenXY, + nX); + } +#endif + + if (nX > 0) + PUSH_3(tenXY->forceX, tenXY, nX); + else + PUSH_3(tenXY->forceX, tenXY, 1); +} + +void +force_To_DigsX_From_Sqrt_MatX_Entry() +{ + DigsX *digsX; + int digitsNeeded; + void force_To_DigsX_From_Sqrt_MatX_Cont(); + void force_To_DigsX_From_DigsX_Entry(); + void force_To_DigsX_From_Sqrt_Reduce(); + + digsX = (DigsX *) POP; + digitsNeeded = (int) POP; + +/* + PUSH_3(force_To_DigsX_From_DigsX_Entry, digsX, digitsNeeded); +*/ + PUSH_2(force_To_DigsX_From_Sqrt_Reduce, digsX); + PUSH_3(force_To_DigsX_From_Sqrt_MatX_Cont, digsX, digitsNeeded); +} + +/* + * This is not used at this stage. This will be used when reduction is + * performed after the tensor is first allocated and when the tensor + * reduces to a matrix. + */ +void +force_To_DigsX_From_Sqrt_MatX_Cont() +{ + DigsX *digsX; + MatX *matX; + int digitsNeeded = 0; + int digitsEmitted; + bool emitDigitFromMatrix(Matrix, Digit *); + + digsX = (DigsX *) POP; + digitsNeeded = (int) POP; + matX = (MatX *) digsX->x; + + digitsEmitted = emitDigits(digsX, + (edf) emitDigitFromMatrix, + (void *) matX->mat, + digitsNeeded); + + +#ifdef TRACE + int bitsShifted = 0; + if (digitsEmitted > 0) + bitsShifted = normalizeMatrix(matX->mat); + + if (TRACE) { + debugp("force_To_DigsX_From_Sqrt_MatX_Cont", + "%x %x emitted=%d shifted=%d\n", + (unsigned) digsX, + (unsigned) matX, + digitsEmitted, + bitsShifted); + } +#endif + + digitsNeeded -= digitsEmitted; + + if (digsX->count > 0) + newDigsX(digsX); + + if (digitsNeeded <= 0) + return; + + /* + * Now absorb everything emitted into the matrix. + */ + absorbDigsXIntoMatX(matX); + + /* + * So now we emitted what we can but still need more. First arrange + * to come back and try to emit again after forcing the necessary + * number of digits from the the argument. + */ + if (digsX->count > 0) + PUSH_3(force_To_DigsX_From_Sqrt_MatX_Cont, digsX->x, digitsNeeded); + else + PUSH_3(force_To_DigsX_From_Sqrt_MatX_Cont, digsX, digitsNeeded); +} + +/* + * In some cases when we generate a matrix which is not refining, but + * where we no that the argument of the matrix is constrained in such a + * way that after a finite amount of absorption, it will become refining. + * What we do is force information from the argument until the matrix is + * refining. + */ +void +force_To_MatX_Until_Refining() +{ + MatX *matX; + int sgn; + + matX = (MatX *) POP; + + if (matX->tag.type != MATX) + return; + + sgn = matrixSign(matX->mat); + + if (sgn > 0) /* matrix is refining and entries positive */ + return; + + if (sgn < 0) { /* matrix is refining and entries negative */ + negateMatrix(matX->mat); + return; + } + + PUSH_2(force_To_MatX_Until_Refining, matX); + PUSH_3(matX->force, matX, 1); +} + +void +force_To_DigsX_From_Sqrt_Reduce() +{ + DigsX *digsX; + + digsX = (DigsX *) POP; + reduceDigsXList(digsX); +} diff --git a/ic-reals-6.3/math-lib/sqrt_R.c~ b/ic-reals-6.3/math-lib/sqrt_R.c~ new file mode 100644 index 0000000..b4e3b18 --- /dev/null +++ b/ic-reals-6.3/math-lib/sqrt_R.c~ @@ -0,0 +1,466 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" + +Real +sqrt_R(Real x) +{ + Bool xLtEq4, xGtEqOneFifth; + Cls *in, *gtEq3, *ltEqOneQuarter; + Real lt0; + void force_To_DigsX_From_Sqrt_TenXY_Entry(); + void sqrtInside(); + void sqrtGtEq3(); + void sqrtLtEqOneQuarter(); + void sqrtLtEqOneQuarterCont(); + void force_To_DigsX_From_Sqrt_TenXY_Entry(); + void force_To_DigsX_From_Sqrt_TenXY_Cont(); + void force_To_DigsX_From_Sqrt_MatX_Entry(); + void force_To_DigsX_From_Sqrt_MatX_Cont(); + void force_To_DigsX_From_Sqrt_Reduce(); + void force_To_MatX_Until_Refining(); + static int doneInit = 0; + Real sqrt_QZ(mpz_t, mpz_t); + + if (!doneInit) { + registerForceFunc(force_To_DigsX_From_Sqrt_TenXY_Entry, + "force_To_DigsX_From_Sqrt_TenXY_Entry", 3); + registerForceFunc(force_To_DigsX_From_Sqrt_TenXY_Cont, + "force_To_DigsX_From_Sqrt_TenXY_Cont", 3); + registerForceFunc(force_To_DigsX_From_Sqrt_MatX_Entry, + "force_To_DigsX_From_Sqrt_MatX_Entry", 3); + registerForceFunc(force_To_DigsX_From_Sqrt_MatX_Cont, + "force_To_DigsX_From_Sqrt_MatX_Cont", 3); + registerForceFunc(force_To_DigsX_From_Sqrt_Reduce, + "force_To_DigsX_From_Sqrt_Reduce", 2); + registerForceFunc(force_To_MatX_Until_Refining, + "force_To_MatX_Until_Refining", 2); + registerForceFunc(sqrtInside, "sqrtInside", 2); + registerForceFunc(sqrtGtEq3, "sqrtGtEq3", 2); + registerForceFunc(sqrtLtEqOneQuarter, "sqrtLtEqOneQuarter", 2); + registerForceFunc(sqrtLtEqOneQuarterCont, "sqrtLtEqOneQuarterCont", 2); + doneInit++; + } + + if (x->gen.tag.type == VECTOR) + return sqrt_QZ(x->vec.vec[0], x->vec.vec[1]); + + xLtEq4 = ltEq_R_QInt(x, 4, 1); + xGtEqOneFifth = gtEq_R_QInt(x, 1, 5); + + in = allocCls(sqrtInside, (void *) x); + in->tag.isSigned = FALSE; + + gtEq3 = allocCls(sqrtGtEq3, (void *) x); + gtEq3->tag.isSigned = FALSE; + + ltEqOneQuarter = allocCls(sqrtLtEqOneQuarter, (void *) x); + ltEqOneQuarter->tag.isSigned = FALSE; + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(in, x); + newEdgeToOnlyChild(gtEq3, x); + newEdgeToOnlyChild(ltEqOneQuarter, x); + endGraphUpdate(); + } + + lt0 = realError("(sqrt_R x) and x < 0"); + + /* + * The order of the tests in the realAlt is not semantically + * significant. Note that the sqrt_R is always unsigned since + * each case is unsigned. + */ + return realIf(4, + and_B_B(xLtEq4, xGtEqOneFifth), (Real) in, + gtEq_R_QInt(x, 3, 1), (Real) gtEq3, + and_B_B(ltEq_R_QInt(x, 1, 4), gtEq_R_0(x)), (Real) ltEqOneQuarter, + lt_R_0(x), lt0); +} + +void +sqrtInside() +{ + Cls *cls; + TenXY *tenXY; + DigsX *digsX; + void force_To_TenXY_X_Until_Refining(); + void force_To_DigsX_From_Sqrt_TenXY_Entry(); + Real x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + digsX = allocDigsX(); + digsX->force = force_To_DigsX_From_Sqrt_TenXY_Entry; + + tenXY = (TenXY *) tensor_Int(x, (Real) digsX, 1, 0, 2, 1, 1, 2, 0, 1); + tenXY->tag.isSigned = FALSE; + + digsX->x = (Real) tenXY; + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(digsX, tenXY); + endGraphUpdate(); + } + + /* + * We must still absorb the sign (if any) and enough information + * to ensure the tensor is refining. Note that, if we get to + * this point, then we know the argument is within the interval + * 1/5 >= x <= 4. I claim that the n'th tensor for n >= 2 + * can always be made refining by absorbing + * information from the left. + */ + PUSH_2(force_To_TenXY_X_Until_Refining, tenXY); + + if (tenXY->x->gen.tag.isSigned) + PUSH_2(tenXY->forceX, tenXY); + + cls->redirect = (Real) digsX; + cls->userData = NULL; + + if (DAVINCI) { + beginGraphUpdate(); + deleteOnlyEdge(cls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} + +void +sqrtGtEq3() +{ + Cls *cls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + w = div_R_Int(x, 4); + w = sqrt_R(w); + cls->redirect = mul_R_Int(w, 2); + cls->userData = NULL; + + if (DAVINCI) { + beginGraphUpdate(); + deleteOnlyEdge(cls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} + +void +sqrtLtEqOneQuarter() +{ + Cls *cls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + w = mul_R_Int(x, 4); + w = sqrt_R(w); + cls->redirect = div_R_Int(w, 2); + cls->userData = NULL; + + if (DAVINCI) { + beginGraphUpdate(); + deleteOnlyEdge(cls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} + +#ifdef LATER +void +sqrtLtEqOneQuarter() +{ + Cls *cls, *newCls; + void sqrtLtEqOneQuarterCont(); + Real x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + + newCls = allocCls(sqrtLtEqOneQuarterCont, (void *) x); + newCls->tag.isSigned = FALSE; + cls->redirect = matrix_Int((Real) newCls, 0, 1, 1, 2); + cls->userData = NULL; + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + deleteOnlyEdge(cls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} +#endif + +void +sqrtLtEqOneQuarterCont() +{ + Cls *cls; + Real w, x; + void force_To_MatX_Until_Refining(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + /* w = matrix_Int(x, 0, 1, 1, 0); reciprocal of x */ + w = mul_R_Int(x, 4); + w = sqrt_R(w); + /* cls->redirect = matrix_Int(w, 1, 0, -2, 1); */ + cls->redirect = w; + cls->redirect->gen.tag.isSigned = FALSE; + cls->userData = NULL; + + if (DAVINCI) { + beginGraphUpdate(); + deleteOnlyEdge(cls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } + + PUSH_2(force_To_MatX_Until_Refining, cls->redirect); +} + +void +force_To_DigsX_From_Sqrt_TenXY_Entry() +{ + DigsX *digsX; + TenXY *tenXY; + int digitsNeeded; + void force_To_DigsX_From_Sqrt_TenXY_Cont(); + void force_To_DigsX_From_DigsX_Entry(); + void force_To_DigsX_From_Sqrt_Reduce(); + + digsX = (DigsX *) POP; + digitsNeeded = (int) POP; + tenXY = (TenXY *) digsX->x; + + /* + * The current strategy is that when forced, we create a list + * of DigsX structures rather than a single struct with the + * requested number of digits. So, after we are done emitting, + * we arrange to reduce the list at the end. + * + * A better strategy would be to allocate two DigsX structures at the + * very start. The second one would be private and only visible to the sqrt + * closure. Then the cycle would be: + * emit from tensor into DigsX-2 + * absorb DigsX-2 digits into tensor + * absorb DigsX-2 digits into DigsX-1 + * clear DigsX-2 of all digits. + * + * This would avoid allocating a chain of DigsX structures and the + * need for reducing at the end. + PUSH_3(force_To_DigsX_From_DigsX_Entry, digsX, digitsNeeded); + */ + PUSH_2(force_To_DigsX_From_Sqrt_Reduce, digsX); + PUSH_3(force_To_DigsX_From_Sqrt_TenXY_Cont, digsX, digitsNeeded); +} + +void +force_To_DigsX_From_Sqrt_TenXY_Cont() +{ + DigsX *digsX; + TenXY *tenXY; + int digitsNeeded; + int nX; + int digitsEmitted, bitsShifted; + int epsDelTensorX(Tensor, int); + bool emitDigitFromTensor(Tensor, Digit *); + + digsX = (DigsX *) POP; + digitsNeeded = (int) POP; + tenXY = (TenXY *) digsX->x; + + digitsEmitted = emitDigits(digsX, + (edf) emitDigitFromTensor, + (void *) tenXY->ten, + digitsNeeded); + + if (digitsEmitted > 0) + bitsShifted = normalizeTensor(tenXY->ten); + + if (TRACE) { + debugp("force_To_DigsX_From_Sqrt_TenXY_Cont", + "%x %x emitted=%d shifted=%d\n", + (unsigned) digsX, + (unsigned) tenXY, + digitsEmitted, + bitsShifted); + } + + digitsNeeded -= digitsEmitted; + + if (digsX->count > 0) + newDigsX(digsX); + + if (digitsNeeded <= 0) + return; + + /* + * So now we emitted what we can but still need more. First arrange + * to come back and try to emit again after forcing the necessary + * number of digits from the the argument. + */ + if (digsX->count > 0) + PUSH_3(force_To_DigsX_From_Sqrt_TenXY_Cont, digsX->x, digitsNeeded); + else + PUSH_3(force_To_DigsX_From_Sqrt_TenXY_Cont, digsX, digitsNeeded); + + /* + * Now absorb everything emitted into the tensor. + */ + absorbDigsXIntoTenXY_Y(tenXY); + + nX = epsDelTensorX(tenXY->ten, digitsNeeded); + + if (TRACE) { + debugp("force_To_DigsX_From_Sqrt_TenXY_Cont", + "%x %x nX=%d\n", + (unsigned) digsX, + (unsigned) tenXY, + nX); + } + + if (nX > 0) + PUSH_3(tenXY->forceX, tenXY, nX); + else + PUSH_3(tenXY->forceX, tenXY, 1); +} + +void +force_To_DigsX_From_Sqrt_MatX_Entry() +{ + DigsX *digsX; + MatX *matX; + int digitsNeeded; + void force_To_DigsX_From_Sqrt_MatX_Cont(); + void force_To_DigsX_From_DigsX_Entry(); + void force_To_DigsX_From_Sqrt_Reduce(); + + digsX = (DigsX *) POP; + digitsNeeded = (int) POP; + matX = (MatX *) digsX->x; + +/* + PUSH_3(force_To_DigsX_From_DigsX_Entry, digsX, digitsNeeded); +*/ + PUSH_2(force_To_DigsX_From_Sqrt_Reduce, digsX); + PUSH_3(force_To_DigsX_From_Sqrt_MatX_Cont, digsX, digitsNeeded); +} + +/* + * This is not used at this stage. This will be used when reduction is + * performed after the tensor is first allocated and when the tensor + * reduces to a matrix. + */ +void +force_To_DigsX_From_Sqrt_MatX_Cont() +{ + DigsX *digsX; + MatX *matX; + int digitsNeeded; + int nX; + int digitsEmitted, bitsShifted; + bool emitDigitFromMatrix(Matrix, Digit *); + + digsX = (DigsX *) POP; + digitsNeeded = (int) POP; + matX = (MatX *) digsX->x; + + digitsEmitted = emitDigits(digsX, + (edf) emitDigitFromMatrix, + (void *) matX->mat, + digitsNeeded); + + if (digitsEmitted > 0) + bitsShifted = normalizeMatrix(matX->mat); + + if (TRACE) { + debugp("force_To_DigsX_From_Sqrt_MatX_Cont", + "%x %x emitted=%d shifted=%d\n", + (unsigned) digsX, + (unsigned) matX, + digitsEmitted, + bitsShifted); + } + + digitsNeeded -= digitsEmitted; + + if (digsX->count > 0) + newDigsX(digsX); + + if (digitsNeeded <= 0) + return; + + /* + * Now absorb everything emitted into the matrix. + */ + absorbDigsXIntoMatX(matX); + + /* + * So now we emitted what we can but still need more. First arrange + * to come back and try to emit again after forcing the necessary + * number of digits from the the argument. + */ + if (digsX->count > 0) + PUSH_3(force_To_DigsX_From_Sqrt_MatX_Cont, digsX->x, digitsNeeded); + else + PUSH_3(force_To_DigsX_From_Sqrt_MatX_Cont, digsX, digitsNeeded); +} + +/* + * In some cases when we generate a matrix which is not refining, but + * where we no that the argument of the matrix is constrained in such a + * way that after a finite amount of absorption, it will become refining. + * What we do is force information from the argument until the matrix is + * refining. + */ +void +force_To_MatX_Until_Refining() +{ + MatX *matX; + int sgn; + + matX = (MatX *) POP; + + if (matX->tag.type != MATX) + return; + + sgn = matrixSign(matX->mat); + + if (sgn > 0) /* matrix is refining and entries positive */ + return; + + if (sgn < 0) { /* matrix is refining and entries negative */ + negateMatrix(matX->mat); + return; + } + + PUSH_2(force_To_MatX_Until_Refining, matX); + PUSH_3(matX->force, matX, 1); +} + +void +force_To_DigsX_From_Sqrt_Reduce() +{ + DigsX *digsX; + + digsX = (DigsX *) POP; + reduceDigsXList(digsX); +} diff --git a/ic-reals-6.3/math-lib/stdTensorCont.c b/ic-reals-6.3/math-lib/stdTensorCont.c new file mode 100644 index 0000000..dcf9cd3 --- /dev/null +++ b/ic-reals-6.3/math-lib/stdTensorCont.c @@ -0,0 +1,59 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +/* + * A number of functions on reals behave in a uniform way in the sense + * that they create a chain of tensors lazily. The only difference are + * the tensors they create. This is the standard tensor continuation. + * It is the force method in a closure and the assumption is that + * the closure data contains a pointer to a function which when called + * produces the next tensor in the chain. + * + * If the argument is a stream terminated by a vector, then there is + * no reduction. As a rule, guarding a vector by a stream inhibits + * reduction of a matrix or tensor against the rational. + * + * When the tensor is created, there may be some information + * about the argument available. At present, we do not absorb this + * information. Perhaps we should. It is not clear what the best policy is. + * + * Note that in the code below, the closure (Cls) is never shared so + * we are free to clobber the data in it. + */ +void +stdTensorCont() +{ + Cls *cls, *newCls; + TenXY *newTenXY; + ClsData *data; + void force_To_TenXY_X_Until_Refining(); + + cls = (Cls *) POP; + + data = (ClsData *) cls->userData; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + newTenXY = (*(data->nextTensor))(data->x, (Real) newCls, data->n); + newTenXY->tag.isSigned = FALSE; + data->n++; + + PUSH_2(force_To_TenXY_X_Until_Refining, newTenXY); + + if (newTenXY->x->gen.tag.isSigned) + PUSH_2(newTenXY->forceX, newTenXY); + + cls->userData = NULL; + cls->redirect = (Real) newTenXY; +} diff --git a/ic-reals-6.3/math-lib/stdTensorCont.c~ b/ic-reals-6.3/math-lib/stdTensorCont.c~ new file mode 100644 index 0000000..3dc2443 --- /dev/null +++ b/ic-reals-6.3/math-lib/stdTensorCont.c~ @@ -0,0 +1,59 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +/* + * A number of functions on reals behave in a uniform way in the sense + * that they create a chain of tensors lazily. The only difference are + * the tensors they create. This is the standard tensor continuation. + * It is the force method in a closure and the assumption is that + * the closure data contains a pointer to a function which when called + * produces the next tensor in the chain. + * + * If the argument is a stream terminated by a vector, then there is + * no reduction. As a rule, guarding a vector by a stream inhibits + * reduction of a matrix or tensor against the rational. + * + * When the tensor is created, there may be some information + * about the argument available. At present, we do not absorb this + * information. Perhaps we should. It is not clear what the best policy is. + * + * Note that in the code below, the closure (Cls) is never shared so + * we are free to clobber the data in it. + */ +void +stdTensorCont() +{ + Cls *cls, *newCls; + TenXY *tenXY, *newTenXY; + ClsData *data; + void force_To_TenXY_X_Until_Refining(); + + cls = (Cls *) POP; + + data = (ClsData *) cls->userData; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + newTenXY = (*(data->nextTensor))(data->x, (Real) newCls, data->n); + newTenXY->tag.isSigned = FALSE; + data->n++; + + PUSH_2(force_To_TenXY_X_Until_Refining, newTenXY); + + if (newTenXY->x->gen.tag.isSigned) + PUSH_2(newTenXY->forceX, newTenXY); + + cls->userData = NULL; + cls->redirect = (Real) newTenXY; +} diff --git a/ic-reals-6.3/math-lib/tan_Q.c b/ic-reals-6.3/math-lib/tan_Q.c new file mode 100644 index 0000000..fe18828 --- /dev/null +++ b/ic-reals-6.3/math-lib/tan_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +tan_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return tan_R(r); +} + +Real +tan_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return tan_R(r); +} diff --git a/ic-reals-6.3/math-lib/tan_R.c b/ic-reals-6.3/math-lib/tan_R.c new file mode 100644 index 0000000..93e19ab --- /dev/null +++ b/ic-reals-6.3/math-lib/tan_R.c @@ -0,0 +1,122 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +static TenXY *nextTensor(Real, Real, int); +static void tanInside(); +static void tanOutside(); + +Real +tan_R(Real x) +{ + Real u, v; + Bool xLtEq1, xGtEqNeg1, xLtEq1_and_GtEqNeg1; + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(tanInside, "tanInside", 2); + registerForceFunc(tanOutside, "tanOutside", 2); + doneInit++; + } + + /* + if (x->gen.tag.type == VECTOR) + return tan_QZ(x->vec.vec[0], x->vec.vec[1]); + */ + + xLtEq1 = ltEq_R_QInt(x, 1, 1); + xGtEqNeg1 = gtEq_R_QInt(x, -1, 1); + xLtEq1_and_GtEqNeg1 = and_B_B(xLtEq1, xGtEqNeg1); + + u = (Real) allocCls(tanInside, (void *) x); + u->cls.tag.isSigned = TRUE; + v = (Real) allocCls(tanOutside, (void *) x); + v->cls.tag.isSigned = TRUE; + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(u, x); + newEdgeToOnlyChild(v, x); + endGraphUpdate(); +#endif + + /* + * The order of the tests in the alt is not semantically + * significant. The tests are applied in order, so there is a + * very modest performance improvement by putting the -1,1 tests + * before the 999/1000,-999/1000 tests. Also there is a very + * small win by doing the negative test before the positive test. + */ + return realIf(4, + not_B(xLtEq1_and_GtEqNeg1), v, + xLtEq1_and_GtEqNeg1, u, + gtEq_R_QInt(x, 999, 1000), v, + ltEq_R_QInt(x, -999, 1000), v); +} + +static TenXY * +nextTensor(Real x, Real y, int n) +{ + return (TenXY *) tensor_Int(x, y, 0, 2, -2, 0, 2*n+1, 2*n+1, 2*n+1, 2*n+1); +} + +static void +tanInside() +{ + Cls *cls, *newCls; + Real x; + ClsData *data; + void stdTensorCont(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + if ((data = (ClsData *) malloc(sizeof(ClsData))) == NULL) + Error(FATAL, E_INT, "tan_R", "malloc failed"); + + data->n = 1; + data->x = x; + data->nextTensor = nextTensor; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + cls->redirect = tensor_Int(x, (Real) newCls, 1, 1, 1, -1, 0, 1, 0, 1); + +#ifdef DAVINCI + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif +} + +static void +tanOutside() +{ + Cls *cls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + x = div_R_Int(x, 2); + + w = tan_R(x); + + cls->redirect = tensor_Int(w, w, 0, -1, 1, 0, 1, 0, 0, 1); + +#ifdef DAVINCI + beginGraphUpdate(); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); +#endif +} diff --git a/ic-reals-6.3/math-lib/tan_R.c~ b/ic-reals-6.3/math-lib/tan_R.c~ new file mode 100644 index 0000000..5969b00 --- /dev/null +++ b/ic-reals-6.3/math-lib/tan_R.c~ @@ -0,0 +1,122 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" +#include "real-impl.h" +#include "math-lib.h" + +static TenXY *nextTensor(Real, Real, int); +static void tanInside(); +static void tanOutside(); + +Real +tan_R(Real x) +{ + Real u, v; + Bool xLtEq1, xGtEqNeg1, xLtEq1_and_GtEqNeg1; + static int doneInit = 0; + + if (!doneInit) { + registerForceFunc(tanInside, "tanInside", 2); + registerForceFunc(tanOutside, "tanOutside", 2); + doneInit++; + } + + /* + if (x->gen.tag.type == VECTOR) + return tan_QZ(x->vec.vec[0], x->vec.vec[1]); + */ + + xLtEq1 = ltEq_R_QInt(x, 1, 1); + xGtEqNeg1 = gtEq_R_QInt(x, -1, 1); + xLtEq1_and_GtEqNeg1 = and_B_B(xLtEq1, xGtEqNeg1); + + u = (Real) allocCls(tanInside, (void *) x); + u->cls.tag.isSigned = TRUE; + v = (Real) allocCls(tanOutside, (void *) x); + v->cls.tag.isSigned = TRUE; + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(u, x); + newEdgeToOnlyChild(v, x); + endGraphUpdate(); + } + + /* + * The order of the tests in the alt is not semantically + * significant. The tests are applied in order, so there is a + * very modest performance improvement by putting the -1,1 tests + * before the 999/1000,-999/1000 tests. Also there is a very + * small win by doing the negative test before the positive test. + */ + return realIf(4, + not_B(xLtEq1_and_GtEqNeg1), v, + xLtEq1_and_GtEqNeg1, u, + gtEq_R_QInt(x, 999, 1000), v, + ltEq_R_QInt(x, -999, 1000), v); +} + +static TenXY * +nextTensor(Real x, Real y, int n) +{ + return (TenXY *) tensor_Int(x, y, 0, 2, -2, 0, 2*n+1, 2*n+1, 2*n+1, 2*n+1); +} + +static void +tanInside() +{ + Cls *cls, *newCls; + Real x, tan_x; + ClsData *data; + void stdTensorCont(); + + cls = (Cls *) POP; + x = (Real) cls->userData; + + if ((data = (ClsData *) malloc(sizeof(ClsData))) == NULL) + Error(FATAL, E_INT, "tan_R", "malloc failed"); + + data->n = 1; + data->x = x; + data->nextTensor = nextTensor; + + newCls = allocCls(stdTensorCont, (void *) data); + newCls->tag.isSigned = FALSE; + + cls->redirect = tensor_Int(x, (Real) newCls, 1, 1, 1, -1, 0, 1, 0, 1); + + if (DAVINCI) { + beginGraphUpdate(); + newEdgeToOnlyChild(newCls, x); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} + +static void +tanOutside() +{ + Cls *cls; + Real w, x; + + cls = (Cls *) POP; + x = (Real) cls->userData; + x = div_R_Int(x, 2); + + w = tan_R(x); + + cls->redirect = tensor_Int(w, w, 0, -1, 1, 0, 1, 0, 0, 1); + + if (DAVINCI) { + beginGraphUpdate(); + drawEqEdge(cls, cls->redirect); + endGraphUpdate(); + } +} diff --git a/ic-reals-6.3/math-lib/tanh_Q.c b/ic-reals-6.3/math-lib/tanh_Q.c new file mode 100644 index 0000000..159e3de --- /dev/null +++ b/ic-reals-6.3/math-lib/tanh_Q.c @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +/* + * These will be specialized in a later version. + */ +Real +tanh_QInt(int a, int b) +{ + Real r; + + r = vector_Int(a, b); + return tanh_R(r); +} + +Real +tanh_QZ(mpz_t a, mpz_t b) +{ + Real r; + + r = vector_Z(a, b); + return tanh_R(r); +} diff --git a/ic-reals-6.3/math-lib/tanh_R.c b/ic-reals-6.3/math-lib/tanh_R.c new file mode 100644 index 0000000..41be93b --- /dev/null +++ b/ic-reals-6.3/math-lib/tanh_R.c @@ -0,0 +1,19 @@ +/* + * Copyright (C) 2000, Imperial College + * + * This file is part of the Imperial College Exact Real Arithmetic Library. + * See the copyright notice included in the distribution for conditions + * of use. + */ + +#include <stdio.h> +#include "real.h" + +Real +tanh_R(Real x) +{ + Real r; + + r = exp_R(x); + return tensor_Int(r, r, 1, 1, 0, 0, 0, 0, -1, 1); +} |