diff options
| -rw-r--r-- | arith.c | 88 | ||||
| -rw-r--r-- | arith.h | 2 | ||||
| -rw-r--r-- | txr.1 | 66 |
3 files changed, 156 insertions, 0 deletions
@@ -2970,6 +2970,92 @@ val digits(val n, val base) return digcommon(0, lit("digits"), n, base); } +val poly(val x, val seq) +{ + val self = lit("rpoly"); + val acc = zero; + seq_info_t si = seq_info(seq); + + if (!numberp(x)) + uw_throwf(error_s, lit("~a: bad argument ~s; number required"), + self, x, nao); + + switch (si.kind) { + case SEQ_NIL: + return acc; + case SEQ_LISTLIKE: + while (consp(seq)) { + val coeff = pop(&seq); + acc = plus(mul(acc, x), coeff); + } + + (void) endp(seq); + return acc; + case SEQ_VECLIKE: + { + cnum len = c_num(length(seq)), i = 0; + + while (i < len) { + val coeff = ref(seq, num(i++)); + acc = plus(mul(acc, x), coeff); + } + + return acc; + } + default: + uw_throwf(error_s, lit("~a: bad argument ~s; poly wants a sequence!"), + self, seq, nao); + + } +} + +val rpoly(val x, val seq) +{ + val self = lit("poly"); + val acc = zero; + val pow = x; + seq_info_t si = seq_info(seq); + + if (!numberp(x)) + uw_throwf(error_s, lit("~a: bad argument ~s; poly wants a number!"), + self, x, nao); + + switch (si.kind) { + case SEQ_NIL: + return acc; + case SEQ_LISTLIKE: + if (consp(seq)) + acc = pop(&seq); + + while (consp(seq)) { + val coeff = pop(&seq); + acc = plus(acc, mul(pow, coeff)); + if (seq) + pow = mul(pow, x); + } + + (void) endp(seq); + return acc; + case SEQ_VECLIKE: + { + cnum len = c_num(length(seq)), i = len; + + while (i > 0) { + val coeff = ref(seq, num(--i)); + acc = plus(mul(acc, x), coeff); + } + + return acc; + } + default: + uw_throwf(error_s, lit("~a: bad argument ~s; poly wants a sequence!"), + self, seq, nao); + + } + + return acc; +} + void arith_init(void) { mp_init(&NUM_MAX_MP); @@ -3007,6 +3093,8 @@ void arith_init(void) reg_fun(intern(lit("bits"), system_package), func_n1(bits)); reg_fun(intern(lit("digpow"), user_package), func_n2o(digpow, 1)); reg_fun(intern(lit("digits"), user_package), func_n2o(digits, 1)); + reg_fun(intern(lit("poly"), user_package), func_n2(poly)); + reg_fun(intern(lit("rpoly"), user_package), func_n2(rpoly)); } void arith_free_all(void) @@ -49,6 +49,8 @@ val width(val num); val bits(val obj); val digpow(val n, val base); val digits(val n, val base); +val poly(val x, val seq); +val rpoly(val x, val seq); noreturn void do_mp_error(val self, mp_err code); void arith_init(void); @@ -34851,6 +34851,72 @@ is zero. In that case, a one-element list containing zero is returned. (digpow 0) -> (0) .cble +.coNP Functions @ poly and @ rpoly +.synb +.mets (poly < arg << coeffs ) +.mets (rpoly < arg << coeffs ) +.syne +.desc +The +.code poly +and +.code rpoly +functions evaluate a polynomial, for the given numeric argument value +.meta arg +and the coefficients given by +.metn coeffs , +a sequence of numbers. + +If +.meta coeffs +is an empty sequence, it denotes the zero polynomial, whose value +is zero everywhere; the functions return zero in this case. + +Otherwise, the +.code poly +function considers +.meta coeffs +to hold the coefficients in the conventional order, namely in order +of decreasing degree of polynomial term. The first element of +.meta coeffs +is the leading coefficient, and the constant term appears as the last element. + +The +.code rpoly +function takes the coefficients in opposite order: the first element of +.meta coeffs +gives the constant term coefficient, and the last element gives the +leading coefficient. + +Note: except in the case of +.code rpoly +operating on a list or list-like sequence of coefficients, +Horner's method of evaluation is +used: a single result accumulator is initialized with zero, and then for each +successive coefficient, in order of decreasing term degree, the accumulator is +multiplied by the argument, and the coefficient is added. When +.code rpoly +operates on a list or list-like sequence, it makes a single +pass through the coefficients in order, thus taking them in increasing +term degree. It maintains two accumulators: one for successive powers of +.meta arg +and one for the resulting value. For each coefficient, the power +accumulator is updated by a multiplication by +.meta arg +and then this value is multiplied by the coefficient, and +that value is then added to the result accumulator. + +.TP* Examples: + +.cblk + ;; 2 + ;; evaluate x + 2x + 3 for x = 10. + (poly 10 '(1 2 3)) -> 123 + + ;; 2 + ;; evaluate 3x + 2x + 1 for x = 10. + (rpoly 10 '(1 2 3)) -> 321 +.cble .SS* Bit Operations In \*(TL, similarly to Common Lisp, bit operations on integers are based |