[m-rev.] for review: float__pow: accept negative exponent; Russian peasants alg

[Peter Moulder]

On Thu, Feb 28, 2002 at 11:27:50AM +1100, Ralph Becket wrote:

> In the interests of clarity, I'd rather see you use Exp // 2 than
> unchecked_right_shift(Exp, 1) and let the compiler do the obvious
> optimization.

I've changed it to div.  `//' requires the compiler to work out that Exp
is non-negative if it's to implement with a shift.  (I haven't checked
whether the compiler does "do the obvious optimization" for div
<constant_power_of_two>, but it's reasonable to expect that the
optimization will be present in future.  Same goes for `mod
<constant_power_of_two>'.)

> Similarly, I'd like to see even/1 go in int.m rather than use the idiom
> Exp /\ 1 = 0.

Done.  (Have also changed int__pow to russian peasants alg.  Presumably
smaller range of exponent for int__pow, so less performance difference,
but can't hurt.  Have tried to keep the float and int versions as
similar as possible.)

> My last (very minor) niggle is that I feel the compiler should perform
> accumulator introduction automatically...  Since it *can* do such a
> thing, I'd like to see it happen in the library at least since the
> head recursive code is easier to read.

It's a bit tougher in this case, as it requires knowing that (EXPR) = 1.0
* (EXPR).

I've left in the explicit accumulator version, though have changed the
name & documentation to maximize understandability.

I'd like to commit this version unless there are any major issues with
the below patch.  If there are non-functional changes that someone wants
to make, such as not using explicit accumulator or renaming or changing
comments or whatever, then I'd prefer for someone else to do that in a
separate commit.


Testing done:

  $ ./int_pow_tst
  2^6 = 64:       PASS
  5^0 = 1:        PASS
  -2^0 = 1:       PASS
  -2^3 = -8:      PASS
  3^4 = 81:       PASS
  4^15 = 1073741824:      PASS
  $ ./float_pow_tst
  2.00000000000000^6 = 64.0000000000000:  PASS
  0.500000000000000^-2 = 4.00000000000000:        PASS
  5.12340000000000^0 = 1.00000000000000:  PASS
  -2.34320000000000^0 = 1.00000000000000: PASS
  -2.00000000000000^3 = -8.00000000000000:        PASS
  -2.00000000000000^-3 = -0.125000000000000:      PASS
  3.00000000000000^4 = 81.0000000000000:  PASS
  2.00000000000000^49 = 562949953421312.: PASS
  4.00000000000000^15 = 1073741824.00000: PASS
  8.00000000000000^-2 = 0.0156250000000000:       PASS

Bugs in Mercury's constant outputting (e.g. writing only 15 digits,
which isn't enough to distinguish every float (double) value, and
mishandling infinity) restrict the range of float tests that can be
done.

pjm.


Estimated hours taken: 4

Let float__pow handle negative exponents.
Use "Russian peasants" algorithm for float__pow and int__pow.
New predicates int__even/1, int__odd/1.

library/float.m:
        float__pow: now works with negative Exponent too.
        Also change to use Russian peasants algorithm.
library/int.m:
	int__pow: Change to use Russian peasants algorithm.
NEWS:
        Document these changes.
        (Also s/float/float.m/ in an unrelated item.)


Index: NEWS
===================================================================
RCS file: /home/mercury1/repository/mercury/NEWS,v
retrieving revision 1.252
diff -d -u -r1.252 NEWS
--- NEWS	13 Mar 2002 12:19:26 -0000	1.252
+++ NEWS	15 Mar 2002 04:19:52 -0000
@@ -139,7 +139,7 @@
   `ops__init_op_table' has been renamed `ops__init_mercury_op_table'.
   `ops__max_priority' is now a function taking an operator table argument.
 
-* The predicates and functions in int.m, float, math.m and array.m now
+* The predicates and functions in int.m, float.m, math.m and array.m now
   generate exceptions rather than program aborts on domain errors and
   out-of-bounds array accesses. There are new functions
   `float__unchecked_quotient/2', `int__unchecked_quotient/2' and
@@ -150,6 +150,9 @@
   float.m and extras/complex_numbers.  (Because of rounding errors,
   the functions aren't actually reversible.)
 
+* float__pow now works for negative exponents, and runs much faster
+  for large exponents.
+
 * We've removed the destructive update modes of string__set_char,
   string__set_char_det and string__unsafe_set_char. The compiler
   currently always stores constant strings in static data, even
@@ -206,6 +209,8 @@
   representations of term components as strings.
 
 * We've made the outputs of the string concatenation primitives unique.
+
+* New convenience/readability predicates `int__even/1' and `int__odd/1'.
 
 * We've removed the long obsolete `int__builtin_*' and
   `float__builtin_float_*' predicates, which were synonyms
Index: int.m
===================================================================
RCS file: /home/mercury1/repository/mercury/library/int.m,v
retrieving revision 1.83
diff -d -u -r1.83 int.m
--- int.m	18 Feb 2002 07:01:04 -0000	1.83
+++ int.m	14 Mar 2002 06:46:29 -0000
@@ -173,6 +173,14 @@
 :- func unchecked_right_shift(int, int) = int.
 :- mode unchecked_right_shift(in, in) = uo is det.
 
+	% even(X) is equivalent to (X mod 2 = 0).
+:- pred even(int).
+:- mode even(in) is semidet.
+
+	% odd(X) is equivalent to (not even(X)), i.e. (X mod 2 = 1).
+:- pred odd(int).
+:- mode odd(in) is semidet.
+
 	% bitwise and
 :- func int /\ int = int.
 :- mode in  /\ in  = uo  is det.
@@ -381,6 +389,14 @@
 		)
 	).
 
+:- pragma inline(even/1).
+even(X):-
+	(X /\ 1) = 0.
+
+:- pragma inline(odd/1).
+odd(X):-
+	(X /\ 1) \= 0.
+
 int__abs(Num) = Abs :-
 	int__abs(Num, Abs).
 
@@ -417,26 +433,29 @@
 		Min = Y
 	).
 
-int__pow(Val, Exp) = Result :-
-	int__pow(Val, Exp, Result).
+int__pow(Base, Exp) = Result :-
+	int__pow(Base, Exp, Result).
 
-int__pow(Val, Exp, Result) :-
+int__pow(Base, Exp, Result) :-
 	( domain_checks, Exp < 0 ->
 		throw(math__domain_error("int__pow"))
 	;
-		int__pow_2(Val, Exp, 1, Result)
+		Result = int__multiply_by_pow(1, Base, Exp)
 	).
 
-:- pred int__pow_2(int, int, int, int).
-:- mode int__pow_2(in, in, in, out) is det.
-
-int__pow_2(Val, Exp, Result0, Result) :-
+:- func int__multiply_by_pow(int, int, int) = int.
+	% Returns Scale0 * (Base ** Exp).
+	% Requires that Exp >= 0.
+int__multiply_by_pow(Scale0, Base, Exp) = Result :-
 	( Exp = 0 ->
-		Result = Result0
+		Result = Scale0
 	;
-		Exp1 is Exp - 1,
-		Result1 is Result0 * Val,
-		int__pow_2(Val, Exp1, Result1, Result)
+		( odd(Exp) ->
+			Scale1 = Scale0 * Base
+		;
+			Scale1 = Scale0
+		),
+		Result = int__multiply_by_pow(Scale1, Base * Base, Exp div 2)
 	).
 
 int__log2(X) = N :-
Index: float.m
===================================================================
RCS file: /home/mercury1/repository/mercury/library/float.m,v
retrieving revision 1.43
diff -d -u -r1.43 float.m
--- float.m	18 Feb 2002 07:01:03 -0000	1.43
+++ float.m	14 Mar 2002 06:46:29 -0000
@@ -59,7 +59,7 @@
 :- mode in    * in    = uo  is det.
 
 	% division
-	% Throws an `math__domain_error' exception if the right
+	% Throws a `math__domain_error' exception if the right
 	% operand is zero. See the comments at the top of math.m
 	% to find out how to disable this check.
 :- func float / float = float.
@@ -135,9 +135,9 @@
 :- func min(float, float) = float.
 
 	% pow(Base, Exponent) returns Base raised to the power Exponent.
-	% The exponent must be an integer greater or equal to 0.
-	% Currently this function runs at O(n), where n is the value
-	% of the exponent.
+	% Fewer domain restrictions than math__pow: works for negative Base,
+	% and float__pow(B, 0) = 1.0 for all B, even B=0.0.
+	% Only pow(0, <negative>) throws a `math__domain_error' exception.
 :- func pow(float, int) = float.
 
 	% Compute a non-negative integer hash value for a float.
@@ -337,20 +337,57 @@
 		Min = Y
 	).
 
-% float_pow(Base, Exponent) = Answer.
-%	XXXX This function could be more efficient, with an int_mod pred, to
-%	reduce O(N) to O(logN) of the exponent.
-float__pow(X, Exp) = Ans :-
-	( Exp < 0 ->
-		throw(math__domain_error("float__pow"))
-	; Exp = 1 ->
-		Ans =  X
-	; Exp = 0 ->
-		Ans = 1.0
+
+float__pow(Base, Exp) = Ans :-
+	( Exp >= 0 ->
+		Ans = float__multiply_by_pow(1.0, Base, Exp)
 	;
-		New_e is Exp - 1,
-		Ans is X * float__pow(X, New_e)
+		( domain_checks, Base = 0.0 ->
+			throw(math__domain_error("float:pow"))
+		;
+			Ans = unchecked_quotient(1.0,
+				float__multiply_by_pow(1.0, Base, -Exp))
+			% See below re use of unchecked_quotient.
+		)
 	).
+
+:- func float__multiply_by_pow(float, float, int) = float.
+	% Returns Scale0 * (Base ** Exp) (where X ** 0 == 1.0 for all X).
+	% Requires that Exp >= 0.
+	% Uses a simple "Russian peasants" algorithm.  O(lg(Exp+1)).
+float__multiply_by_pow(Scale0, Base, Exp) = Result :-
+	( Exp = 0 ->
+		Result = Scale0
+	;
+		( odd(Exp) ->
+			Scale1 = Scale0 * Base
+		;
+			Scale1 = Scale0
+		),
+		Result = float__multiply_by_pow(Scale1, Base * Base, Exp div 2)
+	).
+
+	% The reason for using unchecked_quotient in float__pow is so
+	% that float__pow(+/-0.5, -1111) gives +/-infinity rather than
+	% a domain error.  (N.B. This relies on unchecked_quotient(1.0,
+	% +/-0.0) giving +/-infinity, whereas the documentation in
+	% float.m says that the results are undefined.)
+	% Using Result = float__multiply_by_pow(1.0, 1.0 / Base, -Exp)
+	% would give the right behaviour for underflow, but isn't
+	% generally as accurate.
+
+	% (Efficiency note: An optimization used by `power' in SGI's STL
+	%  implementation is to test for Exp=0 and (for non-zero Exp) handle
+	%  low zero bits in Exp before calling this loop: the loop for the low
+	%  zero bits needs only square Base, it needn't update Acc until the
+	%  end of that loop at which point Acc can be simply assigned from the
+	%  then-current value of Base.  This optimization would be especially
+	%  valuable for expensive `*' operations; maybe provide a
+	%  std_util__monoid_pow(func(T,T)=T MonoidOperator, T Identity, int
+	%  Exp, T Base) = T Result function to complement the existing
+	%  std_util__pow function.)
+
+%---------------------------------------------------------------------------%
 
 :- pragma foreign_proc("C", float__hash(F::in) = (H::out),
 	[will_not_call_mercury, promise_pure, thread_safe],
--------------------------------------------------------------------------
mercury-reviews mailing list
post:  mercury-reviews at cs.mu.oz.au
administrative address: owner-mercury-reviews at cs.mu.oz.au
unsubscribe: Address: mercury-reviews-request at cs.mu.oz.au Message: unsubscribe
subscribe:   Address: mercury-reviews-request at cs.mu.oz.au Message: subscribe
--------------------------------------------------------------------------