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

[Peter Moulder]

This diff has two parts: one is just to implement the russian peasants
algorithm (O(lg Exp) instead of current O(Exp); comments in the code
suggest that int__mod etc. weren't available when float__pow was written);
the other change involves an interface change, in that I'm proposing
accepting negative exponents instead of either throwing an exception
(current CVS version) or calling error/1 (0.10.1 version).

Discussion of this proposed interface change is part of the patch.  In
the last numbered release (0.10.1), float__pow called error for
negative exponents, so backwards compatibility is only even mildly an
issue for people using snapshots since 2 Sep when it was changed to
throw an exception instead of calling error.

Can people look at that discussion and then comment on the proposed
change if you can see something I've missed?

In any case, the change of algorithm is orthogonal to any change in
interface, and I don't expect any issues with that part.  I've done
some private testing of the change (source code available on request),
though I haven't added anything to tests/general/float_test.m, which
looks as if it's in a state of flux (most of it being commented out).


pjm.


Estimated hours taken: 2

library/float.m:
	float__pow: now works with negative Exponent too.
	Also 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.245
diff -d -u -r1.245 NEWS
--- NEWS	24 Feb 2002 11:53:16 -0000	1.245
+++ NEWS	27 Feb 2002 07:15:29 -0000
@@ -108,7 +108,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
@@ -118,6 +118,9 @@
 * We've removed the buggy reverse modes of the arithmetic functions in
   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
Index: library/float.m
===================================================================
RCS file: /home/mercury1/repository/mercury/library/float.m,v
retrieving revision 1.43
diff -d -u -r1.43 float.m
--- library/float.m	18 Feb 2002 07:01:03 -0000	1.43
+++ library/float.m	27 Feb 2002 07:15:29 -0000
@@ -135,9 +135,11 @@
 :- 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 (including B=0.0).
+	% The only domain restriction is introduced by the `/' operation
+	% for negative Exponents, for which Base must not be zero (nor
+	% so small that Base**-Exponent underflows to zero).
 :- func pow(float, int) = float.
 
 	% Compute a non-negative integer hash value for a float.
@@ -337,20 +339,80 @@
 		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
-	;
-		New_e is Exp - 1,
-		Ans is X * float__pow(X, New_e)
+
+% float__pow(Base, Exponent) = Answer.
+
+float__pow(Base, Exp) = Ans :-
+	( Exp < 0 ->
+		% Old behaviour (since 2 September 2001 in CVS):
+		%throw(math__domain_error("float__pow"))
+		% Behaviour before 2Sep, including in version 0.10.1:
+		%error("float__pow taken with exponent < 0\n")
+
+		% Arguments in favour of retaining old behaviour:
+		%   - Almost none that I can think of.  Even backwards compatibility
+		%     isn't really an issue given that the program simply aborted
+		%     in the last numbered-release version of the library.
+		%
+		% Arguments in favour of allowing negative exponent:
+		%   - There's no mathematical reason not to allow negative
+		%     exponent (except perhaps for the Base=+/-0.0 case,
+		%     see above).
+		%   - If we're doing a check for negative Exponent anyway
+		%     then there's no efficiency reason not to handle
+		%     negative exponents.
+		%     (Note that at least until now the code hasn't
+		%     even called domain_checks before checking for
+		%     negative exponent.)
+		%   - It's useful.  (Especially for Base = 10.0 or 2.0.)
+		%   - It seems strange that `<<' in Mercury accepts a negative
+		%     right operand when float__pow doesn't.
+		%
+		% Alternatives are to use the old behaviour except also
+		% test domain_checks; and/or to delay making a decision so
+		% that more people can see and add to this discussion before
+		% we change the behaviour.
+		%
+		% (If we decide to reject or delay allowing float__pow for
+		% negative exponents, then consider adding a ldexp function
+		% to math.m or equivalently a `<<' operator to float.m; these
+		% at least address the case where Base is a constant 2.0.)
+
+		Ans = 1.0 / float__acc_pow(1.0, Base, -Exp)
+	;
+		Ans = float__acc_pow(1.0, Base, Exp)
 	).
+
+
+:- func float__acc_pow(float, float, int) = float.
+
+% Returns Acc0 * (Base ** Exp) (where X ** 0 == 1.0 for all X).
+% Here Exp must be >= 0.
+% Uses a simple "Russian peasants" algorithm.  O(lg(Exp+1)).
+%
+% (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.)
+
+float__acc_pow(Acc0, Base, Exp) = Ans :-
+	( Exp = 0 ->
+		Ans = Acc0
+	;
+		( Exp /\ 1 \= 0 ->
+			Acc1 = Base * Acc0
+		;
+			Acc1 = Acc0
+		),
+		Ans = float__acc_pow(Acc1, Base * Base,
+				     unchecked_right_shift(Exp, 1))
+	).
+
 
 :- 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
--------------------------------------------------------------------------