From 3ad3ac740fc004ec24b4ddefa0425e32d8590d98 Mon Sep 17 00:00:00 2001 From: Andy Wingo Date: Mon, 3 Jan 2022 16:19:44 +0100 Subject: [PATCH] Reimplement integer-expt in Scheme * libguile/numbers.c (integer_expt_var): New static variable. (init_integer_expt_var): New helper. (scm_integer_expt): Delegate to Scheme. * module/ice-9/boot-9.scm (integer-expt): Reimplement in Scheme. Misses some optimizations for fractions but that is probably OK! --- libguile/numbers.c | 138 ++++++---------------------------------- module/ice-9/boot-9.scm | 40 +++++++++++- 2 files changed, 57 insertions(+), 121 deletions(-) diff --git a/libguile/numbers.c b/libguile/numbers.c index 0afa83b79..ab7a659c8 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -61,6 +61,7 @@ #include "boolean.h" #include "deprecation.h" #include "eq.h" +#include "eval.h" #include "feature.h" #include "finalizers.h" #include "goops.h" @@ -72,6 +73,8 @@ #include "simpos.h" #include "smob.h" #include "strings.h" +#include "threads.h" +#include "variable.h" #include "values.h" #include "numbers.h" @@ -3208,128 +3211,23 @@ SCM_DEFINE (scm_modulo_expt, "modulo-expt", 3, 0, 0, } #undef FUNC_NAME -SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0, - (SCM n, SCM k), - "Return @var{n} raised to the power @var{k}. @var{k} must be an\n" - "exact integer, @var{n} can be any number.\n" - "\n" - "Negative @var{k} is supported, and results in\n" - "@math{1/@var{n}^abs(@var{k})} in the usual way.\n" - "@math{@var{n}^0} is 1, as usual, and that\n" - "includes @math{0^0} is 1.\n" - "\n" - "@lisp\n" - "(integer-expt 2 5) @result{} 32\n" - "(integer-expt -3 3) @result{} -27\n" - "(integer-expt 5 -3) @result{} 1/125\n" - "(integer-expt 0 0) @result{} 1\n" - "@end lisp") -#define FUNC_NAME s_scm_integer_expt +static SCM integer_expt_var; + +static void +init_integer_expt_var (void) { - scm_t_inum i2 = 0; - SCM z_i2 = SCM_BOOL_F; - int i2_is_big = 0; - SCM acc = SCM_I_MAKINUM (1L); - - /* Specifically refrain from checking the type of the first argument. - This allows us to exponentiate any object that can be multiplied. - If we must raise to a negative power, we must also be able to - take its reciprocal. */ - if (!SCM_LIKELY (SCM_I_INUMP (k)) && !SCM_LIKELY (SCM_BIGP (k))) - SCM_WRONG_TYPE_ARG (2, k); - - if (SCM_UNLIKELY (scm_is_eq (k, SCM_INUM0))) - return SCM_INUM1; /* n^(exact0) is exact 1, regardless of n */ - else if (SCM_UNLIKELY (scm_is_eq (n, SCM_I_MAKINUM (-1L)))) - return scm_is_false (scm_even_p (k)) ? n : SCM_INUM1; - /* The next check is necessary only because R6RS specifies different - behavior for 0^(-k) than for (/ 0). If n is not a scheme number, - we simply skip this case and move on. */ - else if (SCM_NUMBERP (n) && scm_is_true (scm_zero_p (n))) - { - /* k cannot be 0 at this point, because we - have already checked for that case above */ - if (scm_is_true (scm_positive_p (k))) - return n; - else /* return NaN for (0 ^ k) for negative k per R6RS */ - return scm_nan (); - } - else if (SCM_FRACTIONP (n)) - { - /* Optimize the fraction case by (a/b)^k ==> (a^k)/(b^k), to avoid - needless reduction of intermediate products to lowest terms. - If a and b have no common factors, then a^k and b^k have no - common factors. Use 'scm_i_make_ratio_already_reduced' to - construct the final result, so that no gcd computations are - needed to exponentiate a fraction. */ - if (scm_is_true (scm_positive_p (k))) - return scm_i_make_ratio_already_reduced - (scm_integer_expt (SCM_FRACTION_NUMERATOR (n), k), - scm_integer_expt (SCM_FRACTION_DENOMINATOR (n), k)); - else - { - k = scm_difference (k, SCM_UNDEFINED); - return scm_i_make_ratio_already_reduced - (scm_integer_expt (SCM_FRACTION_DENOMINATOR (n), k), - scm_integer_expt (SCM_FRACTION_NUMERATOR (n), k)); - } - } - - if (SCM_I_INUMP (k)) - i2 = SCM_I_INUM (k); - else if (SCM_BIGP (k)) - { - z_i2 = scm_i_clonebig (k, 1); - scm_remember_upto_here_1 (k); - i2_is_big = 1; - } - else - SCM_WRONG_TYPE_ARG (2, k); - - if (i2_is_big) - { - if (mpz_sgn(SCM_I_BIG_MPZ (z_i2)) == -1) - { - mpz_neg (SCM_I_BIG_MPZ (z_i2), SCM_I_BIG_MPZ (z_i2)); - n = scm_divide (n, SCM_UNDEFINED); - } - while (1) - { - if (mpz_sgn(SCM_I_BIG_MPZ (z_i2)) == 0) - { - return acc; - } - if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2), 1) == 0) - { - return scm_product (acc, n); - } - if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2), 0)) - acc = scm_product (acc, n); - n = scm_product (n, n); - mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2), SCM_I_BIG_MPZ (z_i2), 1); - } - } - else - { - if (i2 < 0) - { - i2 = -i2; - n = scm_divide (n, SCM_UNDEFINED); - } - while (1) - { - if (0 == i2) - return acc; - if (1 == i2) - return scm_product (acc, n); - if (i2 & 1) - acc = scm_product (acc, n); - n = scm_product (n, n); - i2 >>= 1; - } - } + integer_expt_var = scm_c_module_lookup (scm_the_root_module (), + "integer-expt"); +} + +SCM +scm_integer_expt (SCM n, SCM k) +{ + static scm_i_pthread_once_t once = SCM_I_PTHREAD_ONCE_INIT; + scm_i_pthread_once (&once, init_integer_expt_var); + + return scm_call_2 (scm_variable_ref (integer_expt_var), n, k); } -#undef FUNC_NAME /* Efficiently compute (N * 2^COUNT), where N is an exact integer, and COUNT > 0. */ diff --git a/module/ice-9/boot-9.scm b/module/ice-9/boot-9.scm index 2323b1ec5..e52352962 100644 --- a/module/ice-9/boot-9.scm +++ b/module/ice-9/boot-9.scm @@ -1,6 +1,6 @@ ;;; -*- mode: scheme; coding: utf-8; -*- -;;;; Copyright (C) 1995-2014, 2016-2021 Free Software Foundation, Inc. +;;;; Copyright (C) 1995-2014, 2016-2022 Free Software Foundation, Inc. ;;;; ;;;; This library is free software; you can redistribute it and/or ;;;; modify it under the terms of the GNU Lesser General Public @@ -4618,6 +4618,44 @@ when none is available, reading FILE-NAME with READER." +;;; {Math helpers} +;;; + +(define (integer-expt n k) + "Return @var{n} raised to the power @var{k}. @var{k} must be an exact +integer, @var{n} can be any number. + +Negative @var{k} is supported, and results in +@math{1/@var{n}^abs(@var{k})} in the usual way. @math{@var{n}^0} is 1, +as usual, and that includes @math{0^0} is 1. + +@lisp +(integer-expt 2 5) @result{} 32 +(integer-expt -3 3) @result{} -27 +(integer-expt 5 -3) @result{} 1/125 +(integer-expt 0 0) @result{} 1 +@end lisp" + (cond + ((not (exact-integer? k)) + (scm-error 'wrong-type-arg "integer-expt" + "Wrong type (expected an exact integer): ~S" + (list k) #f)) + ((negative? k) + (if (and (number? n) (zero? n)) + +nan.0 + (integer-expt (/ n) (- k)))) + (else + (let lp ((acc 1) (k k) (n n)) + (cond + ((eqv? k 0) acc) + ((eqv? k 1) (* acc n)) + (else + (lp (if (odd? k) (* acc n) acc) + (ash k -1) + (* n n)))))))) + + + ;;; {R6RS and R7RS} ;;;