[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Guile-commits] 03/03: Re-rewrite integer-expt in C
|
From: |
Andy Wingo |
|
Subject: |
[Guile-commits] 03/03: Re-rewrite integer-expt in C |
|
Date: |
Fri, 7 Jan 2022 16:22:33 -0500 (EST) |
wingo pushed a commit to branch wip-inline-digits
in repository guile.
commit 6c13cd098fcce03eb962d2d08a9d6e8eb7fec336
Author: Andy Wingo <wingo@pobox.com>
AuthorDate: Fri Jan 7 22:15:55 2022 +0100
Re-rewrite integer-expt in C
Calling out to Scheme was a performance regression.
* libguile/integers.h:
* libguile/integers.c (scm_integer_expt_ii, scm_integer_expt_zi): New
internal functions.
* libguile/numbers.c (scm_integer_expt): Go back to C. But, include
fast cases for inums and doubles.
* module/ice-9/boot-9.scm: Revert addition of integer-expt.
---
libguile/integers.c | 24 ++++++++++
libguile/integers.h | 3 ++
libguile/numbers.c | 120 ++++++++++++++++++++++++++++++++++++++++++------
module/ice-9/boot-9.scm | 40 +---------------
4 files changed, 135 insertions(+), 52 deletions(-)
diff --git a/libguile/integers.c b/libguile/integers.c
index cc8c2f4fc..143683738 100644
--- a/libguile/integers.c
+++ b/libguile/integers.c
@@ -2240,6 +2240,30 @@ scm_integer_lognot_z (struct scm_bignum *n)
return take_mpz (result);
}
+SCM
+scm_integer_expt_ii (scm_t_inum n, scm_t_inum k)
+{
+ ASSERT (k >= 0);
+ mpz_t res;
+ mpz_init (res);
+ mpz_ui_pow_ui (res, inum_magnitude (n), k);
+ if (n < 0 && (k & 1))
+ mpz_neg (res, res);
+ return take_mpz (res);
+}
+
+SCM
+scm_integer_expt_zi (struct scm_bignum *n, scm_t_inum k)
+{
+ ASSERT (k >= 0);
+ mpz_t res, zn;
+ mpz_init (res);
+ alias_bignum_to_mpz (n, zn);
+ mpz_pow_ui (res, zn, k);
+ scm_remember_upto_here_1 (n);
+ return take_mpz (res);
+}
+
static void
integer_init_mpz (mpz_ptr z, SCM n)
{
diff --git a/libguile/integers.h b/libguile/integers.h
index a232eb8cc..3fc53f019 100644
--- a/libguile/integers.h
+++ b/libguile/integers.h
@@ -124,6 +124,9 @@ SCM_INTERNAL int scm_integer_logbit_uz (unsigned long bit,
SCM_INTERNAL SCM scm_integer_lognot_i (scm_t_inum n);
SCM_INTERNAL SCM scm_integer_lognot_z (struct scm_bignum *n);
+SCM_INTERNAL SCM scm_integer_expt_ii (scm_t_inum n, scm_t_inum k);
+SCM_INTERNAL SCM scm_integer_expt_zi (struct scm_bignum *n, scm_t_inum k);
+
SCM_INTERNAL SCM scm_integer_modulo_expt_nnn (SCM n, SCM k, SCM m);
SCM_INTERNAL SCM scm_integer_lsh_iu (scm_t_inum n, unsigned long count);
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 20d525408..9ff870aa7 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -60,8 +60,8 @@
#include "bdw-gc.h"
#include "boolean.h"
#include "deprecation.h"
+#include "dynwind.h"
#include "eq.h"
-#include "eval.h"
#include "feature.h"
#include "finalizers.h"
#include "goops.h"
@@ -73,8 +73,6 @@
#include "simpos.h"
#include "smob.h"
#include "strings.h"
-#include "threads.h"
-#include "variable.h"
#include "values.h"
#include "numbers.h"
@@ -2923,23 +2921,119 @@ SCM_DEFINE (scm_modulo_expt, "modulo-expt", 3, 0, 0,
}
#undef FUNC_NAME
-static SCM integer_expt_var;
-
static void
-init_integer_expt_var (void)
+mpz_clear_on_unwind (void *mpz)
{
- integer_expt_var = scm_c_module_lookup (scm_the_root_module (),
- "integer-expt");
+ mpz_clear (mpz);
}
-SCM
-scm_integer_expt (SCM n, SCM k)
+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_i_pthread_once_t once = SCM_I_PTHREAD_ONCE_INIT;
- scm_i_pthread_once (&once, init_integer_expt_var);
+ // Fast cases first.
+ if (SCM_I_INUMP (k))
+ {
+ if (SCM_I_INUM (k) < 0)
+ {
+ if (SCM_NUMBERP (n) && scm_is_true (scm_zero_p (n)))
+ return scm_nan ();
+ k = scm_integer_negate_i (SCM_I_INUM (k));
+ n = scm_divide (n, SCM_UNDEFINED);
+ }
+ if (SCM_I_INUMP (n))
+ return scm_integer_expt_ii (SCM_I_INUM (n), SCM_I_INUM (k));
+ else if (SCM_BIGP (n))
+ return scm_integer_expt_zi (scm_bignum (n), SCM_I_INUM (k));
+ }
+ else if (SCM_BIGP (k))
+ {
+ if (scm_is_integer_negative_z (scm_bignum (k)))
+ {
+ if (SCM_NUMBERP (n) && scm_is_true (scm_zero_p (n)))
+ return scm_nan ();
+ k = scm_integer_negate_z (scm_bignum (k));
+ n = scm_divide (n, SCM_UNDEFINED);
+ }
+ if (scm_is_eq (n, SCM_INUM0) || scm_is_eq (n, SCM_INUM1))
+ return n;
+ else if (scm_is_eq (n, SCM_I_MAKINUM (-1)))
+ return scm_is_integer_odd_z (scm_bignum (k)) ? n : SCM_INUM1;
+ else if (scm_is_exact_integer (n))
+ scm_num_overflow ("integer-expt");
+ }
+ else
+ SCM_WRONG_TYPE_ARG (2, k);
+
+ // The general case.
+ if (scm_is_eq (k, SCM_INUM0))
+ return SCM_INUM1; /* n^(exact0) is exact 1, regardless of n */
+
+ 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));
+ }
+ }
+
+ mpz_t zk;
+ mpz_init (zk);
+ scm_to_mpz (k, zk);
- return scm_call_2 (scm_variable_ref (integer_expt_var), n, k);
+ scm_dynwind_begin (0);
+ scm_dynwind_unwind_handler (mpz_clear_on_unwind, zk, SCM_F_WIND_EXPLICITLY);
+ if (mpz_sgn (zk) == -1)
+ {
+ mpz_neg (zk, zk);
+ n = scm_divide (n, SCM_UNDEFINED);
+ }
+ SCM acc = SCM_INUM1;
+ while (1)
+ {
+ if (mpz_sgn (zk) == 0)
+ break;
+ if (mpz_cmp_ui(zk, 1) == 0)
+ {
+ acc = scm_product (acc, n);
+ break;
+ }
+ if (mpz_tstbit(zk, 0))
+ acc = scm_product (acc, n);
+ n = scm_product (n, n);
+ mpz_fdiv_q_2exp (zk, zk, 1);
+ }
+ scm_dynwind_end ();
+ return acc;
}
+#undef FUNC_NAME
static SCM
lsh (SCM n, SCM count, const char *fn)
diff --git a/module/ice-9/boot-9.scm b/module/ice-9/boot-9.scm
index e52352962..2323b1ec5 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-2022 Free Software Foundation, Inc.
+;;;; Copyright (C) 1995-2014, 2016-2021 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,44 +4618,6 @@ 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}
;;;