Commits (1)
......@@ -7,9 +7,11 @@ ctypedef fused ulong_or_object:
object
cpdef generic_power(a, n)
cdef generic_power_long(a, long n)
cpdef generic_power(a, n, mul=?, inv=?)
cdef generic_power_long(a, long n, mul, inv)
cdef generic_power_pos(a, ulong_or_object n) # n > 0
cdef generic_power_pos2(a, ulong_or_object n, mul)
cdef inline invert(a):
......
......@@ -7,6 +7,7 @@ square-and-multiply algorithm.
#*****************************************************************************
# Copyright (C) 2017 Jeroen Demeyer <J.Demeyer@UGent.be>
# 2018 Vincent Delecroix <20100.delecroix@gmail.com>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
......@@ -20,7 +21,7 @@ from cysignals.signals cimport sig_check
from .long cimport integer_check_long
cpdef generic_power(a, n):
cpdef generic_power(a, n, mul=None, inv=None):
"""
Return `a^n`.
......@@ -33,6 +34,13 @@ cpdef generic_power(a, n):
- ``n`` -- any integer (in the duck typing sense)
- ``mul`` -- an optional function for multiplication. If not
provided, simply use the multiplication implementation of
the underlying object.
- ``invert`` -- an optional function for inversion (only used
if ``mul`` is not ``None`` and exponent is negative)
EXAMPLES::
sage: from sage.arith.power import generic_power
......@@ -71,6 +79,13 @@ cpdef generic_power(a, n):
sage: x = SymbolicMul("x")
sage: print(generic_power(x, 7))
(((x*x)*(x*x))*((x*x)*x))
::
sage: generic_power(4, 6, mul=operator.add)
24
sage: generic_power(4, -6, mul=operator.add, inv=operator.neg)
-24
"""
if not n:
return one(a)
......@@ -80,15 +95,22 @@ cpdef generic_power(a, n):
if not integer_check_long(n, &value, &err):
raise NotImplementedError("non-integral exponents not supported")
if not err:
return generic_power_long(a, value)
return generic_power_long(a, value, mul, inv)
if n < 0:
n = -n
a = invert(a)
return generic_power_pos(a, n)
if mul is not None:
a = inv(a)
else:
a = invert(a)
if mul is None:
return generic_power_pos(a, n)
else:
return generic_power_pos2(a, n, mul)
cdef generic_power_long(a, long n):
cdef generic_power_long(a, long n, mul, inv):
"""
As ``generic_power`` but where ``n`` is a C long.
"""
......@@ -96,10 +118,16 @@ cdef generic_power_long(a, long n):
return one(a)
cdef unsigned long u = <unsigned long>n
if n < 0:
u = -u
a = invert(a)
return generic_power_pos(a, u)
if mul is None:
if n < 0:
u = -u
a = invert(a)
return generic_power_pos(a, u)
else:
if n < 0:
u = -u
a = inv(a)
return generic_power_pos2(a, u, mul)
cdef generic_power_pos(a, ulong_or_object n):
......@@ -124,3 +152,26 @@ cdef generic_power_pos(a, ulong_or_object n):
n >>= 1
return res
cdef generic_power_pos2(a, ulong_or_object n, mul):
"""
Return `a^n` where `n > 0`.
"""
# Find least significant set bit as starting point
apow = a
while not (n & 1):
sig_check()
apow = mul(apow, apow)
n >>= 1
# Now multiply together the correct factors a^(2^i)
res = apow
n >>= 1
while n:
sig_check()
apow = mul(apow, apow)
if n & 1:
res = mul(apow, res)
n >>= 1
return res