Re: Another way to do x^n

"Julian V. Noble" <jvn@xxxxxxxxxxxx> wrote in message
news:e4jclh$8qn$1@xxxxxxxxxxxxxxxxxxxxxxxxxxx
John Doty wrote:
Julian V. Noble wrote:
I haven't seen this one (or if I have, I've forgotten :-( !!)
but I like it because it exposes the underlying algorithm
transparently:

\ raise float to positive integer power
\ algorithm: x^n = (x^2)^[n/2] * x^[n mod 2]

: f^2 FDUP F* ;

: (power) ( f: x -- x^n) ( n --) \ n assumed >= 0
DUP 0= IF DROP FDROP 1e0 EXIT THEN
DUP 1 = IF DROP EXIT THEN
DUP FDUP 2 MOD ( f: x x ) ( n n mod 2) RECURSE
FSWAP f^2 2/ ( f: x^[n mod 2] x^2) ( n/2) RECURSE
( f: x^[n mod 2] [x^2]^[n/2] ) F*
;

: power ( f: x -- x^n) ( n--)
FDUP F0= IF DROP FDROP 0e0 EXIT THEN
(power)
;

Anyway, it's for the people who always claim we never post
any code.


But it doesn't address the more serious problem: the code posted is
often excessively difficult for any but an ANS Forth expert to read.

1. 12 (!) stack adjustments. Yikes!

2. (power) is an incomprehensible 25 words long!

3. Recursion unnecessarily obfuscates the algorithm.

Here's how I would do it in LSE64:

# raise float to positive integer power
# algorithm: x^n = (x^2)^[n/2] * x^[n mod 2]

# Use variables to clarify and avoid stack gymnastics

variables: x n

# y multx yields z
# conditionally multiply stack top by x, if current x power needed

multx : n @ 2 % 0> && x @ *.

# square x

squarex : x @ dup *. x !

# divide n by 2

halven : n @ 2 / n !

# multiply stack top by x until we've exhausted n

(power) : n @ 0> && multx squarex halven repeat

# x n power yields x^n

power : n ! x ! 1.0 (power)

I think this is much simpler and clearer. Flow control by && (exit if
false) and repeat (tail recursion) is easier to follow. All definitions
are 7+-2 words: nothing excessively short or excessively long.
Definitions are not interupted by comments, making flow easier to read.
Few formal comments: strip down code until it speaks for itself.
Comments are for intentions, not mechanics. Clearer style in a better
Forth dialect.

De gustibus...

If you factor out the tests as

: test_exponent
S" DUP 0= IF DROP FDROP 1e0 EXIT THEN" EVALUATE
S" DUP 1 = IF DROP EXIT THEN" EVALUATE
; IMMEDIATE

the word power is only 9 words long. Gimme a break! Also, look
at the algorithm: basically it is

power(x,n) = power(x*x,n/2) * power(x, n mod 2 )

It's just MADE for recursion!

I think this depends on whether you are a mathematician or software
engineer. Recursion is often clear and obvious to mathematicians to not so
to the engineer. Mind you, if I were an engineer and not a Mathematician I'd
not use either method due to rounding errors (which are the same since the
same arithmatic is done in the same order but with different contol
structures).


Peter


.

Relevant Pages

  • Re: Another way to do x^n
    ... \ raise float to positive integer power ... DUP 0= IF DROP FDROP 1e0 EXIT THEN ... FDUP F0= IF DROP FDROP 0e0 EXIT THEN ... Recursion unnecessarily obfuscates the algorithm. ...
    (comp.lang.forth)
  • Re: Another way to do x^n
    ... \ raise float to positive integer power ... FDUP F0= IF DROP FDROP 0e0 EXIT THEN ... Recursion unnecessarily obfuscates the algorithm. ...
    (comp.lang.forth)
  • Re: Another way to do x^n
    ... \ raise float to positive integer power ... FDUP F0= IF DROP FDROP 0e0 EXIT THEN ... Recursion unnecessarily obfuscates the algorithm. ... # conditionally multiply stack top by x, ...
    (comp.lang.forth)
  • Re: to Matt,
    ... algorithm that can run BWT in linear time (with respect to the input ... Quantum computers are thing of past, ... I can use Jules' random compression ... drill a power cord from WB's super nuclear power plant to Sachin's ...
    (comp.compression)
  • Re: Powers of 5
    ... This should give you the modified algorithm almost immediately. ... Courses with greater emphasis on programming have entered ... require the argument to be a power of 2. ...
    (sci.math)