Exponentiation operator

From Rosetta Code

Most all programming languages have a built-in implementation of exponentiation. Re-implement integer exponentiation for both int^int and real^int as both a procedure, and an operator (if your language supports operator definition).

If the language supports operator (or procedure) overloading, then an overloaded form should be provided for both int^int and real^int variants.

Contents 1 Ada 2 ALGOL 68 3 C++ 4 Forth 5 Haskell 6 J 7 Java 8 Logo 9 Lucid

[edit] Ada

First we declare the specifications of the two procedures and the two corresponding operators (written as functions with quoted operators as their names):

package Integer_Exponentiation is
   --  int^int
   procedure Exponentiate (Argument : in     Integer;
                           Exponent : in     Natural;
                           Result   :    out Integer);
   function "**" (Left  : Integer;
                  Right : Natural) return Integer;

   --  real^int
   procedure Exponentiate (Argument : in     Float;
                           Exponent : in     Integer;
                           Result   :    out Float);
   function "**" (Left  : Float;
                  Right : Integer) return Float;
end Integer_Exponentiation;

Now we can create a test program:

with Ada.Float_Text_IO, Ada.Integer_Text_IO, Ada.Text_IO;
with Integer_Exponentiation;

procedure Test_Integer_Exponentiation is
   use Ada.Float_Text_IO, Ada.Integer_Text_IO, Ada.Text_IO;
   use Integer_Exponentiation;
   R : Float;
   I : Integer;
begin
   Exponentiate (Argument => 2.5, Exponent => 3, Result => R);
   Put ("2.5 ^ 3 = ");
   Put (R, Fore => 2, Aft => 4, Exp => 0);
   New_Line;

   Exponentiate (Argument => -12, Exponent => 3, Result => I);
   Put ("-12 ^ 3 = ");
   Put (I, Width => 7);
   New_Line;
end Test_Integer_Exponentiation;

Finally we can implement the procedures and operations:

package body Integer_Exponentiation is
   --  int^int
   procedure Exponentiate (Argument : in     Integer;
                           Exponent : in     Natural;
                           Result   :    out Integer) is
   begin
      Result := 1;
      for Counter in 1 .. Exponent loop
         Result := Result * Argument;
      end loop;
   end Exponentiate;

   function "**" (Left  : Integer;
                  Right : Natural) return Integer is
      Result : Integer;
   begin
      Exponentiate (Argument => Left,
                    Exponent => Right,
                    Result   => Result);
      return Result;
   end "**";

   --  real^int
   procedure Exponentiate (Argument : in     Float;
                           Exponent : in     Integer;
                           Result   :    out Float) is
   begin
      Result := 1.0;
      if Exponent < 0 then
         for Counter in Exponent .. -1 loop
            Result := Result / Argument;
         end loop;
      else
         for Counter in 1 .. Exponent loop
            Result := Result * Argument;
         end loop;
      end if;
   end Exponentiate;

   function "**" (Left  : Float;
                  Right : Integer) return Float is
      Result : Float;
   begin
      Exponentiate (Argument => Left,
                    Exponent => Right,
                    Result   => Result);
      return Result;
   end "**";
end Integer_Exponentiation;

[edit] ALGOL 68

main:(
  INT two=2, thirty=30; # test constants #
  PROC VOID undefined;

# First implement exponentiation using a rather slow but sure FOR loop #
  PROC int pow = (INT base, exponent)INT: ( # PROC cannot be over loaded #
    IF exponent<0 THEN undefined FI;
    INT out:=( exponent=0 | 1 | base );
    FROM 2 TO exponent DO out*:=base OD;
    out
  );

  printf(($" One Gibi-unit is: int pow("g(0)","g(0)")="g(0)" - (cost: "g(0)
           " INT multiplications)"l$,two, thirty, int pow(two,thirty),thirty-1));

# implement exponentiation using a faster binary technique and WHILE LOOP #
  OP ** = (INT base, exponent)INT: (
    BITS binary exponent:=BIN exponent ; # do exponent arithmetic in binary #
    INT out := IF bits width ELEM binary exponent THEN base ELSE 1 FI;
    INT sq := IF exponent < 0 THEN undefined; ~ ELSE base FI;

    WHILE
      binary exponent := binary exponent SHR 1;
      binary exponent /= BIN 0
    DO
      sq *:= sq;
      IF bits width ELEM binary exponent THEN out *:= sq FI
    OD;
    out
  );

  printf(($" One Gibi-unit is: "g(0)"**"g(0)"="g(0)" - (cost: "g(0)
           " INT multiplications)"l$,two, thirty, two ** thirty,8));

  OP ** = (REAL in base, INT in exponent)REAL: ( # ** INT Operator can be overloaded #
    REAL base := ( in exponent<0 | 1/in base | in base);
    INT exponent := ABS in exponent;
    BITS binary exponent:=BIN exponent ; # do exponent arithmetic in binary #
    REAL out := IF bits width ELEM binary exponent THEN base ELSE 1 FI;
    REAL sq := base;

    WHILE
      binary exponent := binary exponent SHR 1;
      binary exponent /= BIN 0
    DO
      sq *:= sq;
      IF bits width ELEM binary exponent THEN out *:= sq FI
    OD;
    out
  );

  printf(($" One Gibi-unit is: "g(0,1)"**"g(0)"="g(0,1)" - (cost: "g(0)
           " REAL multiplications)"l$, 2.0, thirty, 2.0 ** thirty,8));

  OP ** = (REAL base, REAL exponent)REAL: ( # ** REAL Operator can be overloaded #
    exp(ln(base)*exponent)
  );

  printf(($" One Gibi-unit is: "g(0,1)"**"g(0,1)"="g(0,1)" - (cost: "
           "depends on precision)"l$, 2.0, 30.0, 2.0 ** 30.0))
)

Output

One Gibi-unit is: int pow(2,30)=1073741824 - (cost: 29 INT multiplications)
One Gibi-unit is: 2**30=1073741824 - (cost: 8 INT multiplications)
One Gibi-unit is: 2.0**30=1073741824.0 - (cost: 8 REAL multiplications)
One Gibi-unit is: 2.0**30.0=1073741824.0 - (cost: depends on precision)

[edit] C++

While C++ does allow operator overloading, it does not have an exponentiation operator, therefore only a function definition is given. For non-negative exponents the integer and floating point versions are exactly the same, for obvious reasons. For negative exponents, the integer exponentiation would not give integer results; therefore there are several possibilities:

1. Use floating point results even for integer exponents.
2. Use integer results for integer exponents and give an error for negative exponents.
3. Use integer results for integer exponents and return just the integer part (i.e. return 0 if the base is larger than one and the exponent is negative).

The third option somewhat resembles the integer division rules, and has the nice property that it can use the exact same algorithm as the floating point version. Therefore this option is chosen here. Actually the template can be used with any type which supports multiplication, division and explicit initialization from int. Note that there are several aspects about int which are not portably defined; most notably it is not guaranteed

• that the negative of a valid int is again a valid int; indeed for most implementations, the minimal value doesn't have a positive counterpart,
• whether the result of a%b is positive or negative if a is negative, and in which direction the corresponding division is rounded (however, it is guaranteed that (a/b)*b + a%b == a)

The code below tries to avoid those platform dependencies. Note that bitwise operations wouldn't help here either, because the representation of negative numbers can vary as well.

 
template<typename Number>
 Number power(Number base, int exponent)
{
  int zerodir;
  Number factor;
  if (exponent < 0)
  {
    zerodir = 1;
    factor = Number(1)/base;
  }
  else
  {
    zerodir = -1;
    factor = base;
  }
 
  Number result(1);
  while (exponent != 0)
  {
    if (exponent % 2 != 0)
    {
      result *= factor;
      exponent += zerodir;
    }
    else
    {
      factor *= factor;
      exponent /= 2;
    }
  }
  return result;
}
 

[edit] Forth

: ** ( n m -- n^m )
  1 swap  0 ?do over * loop  nip ;

: f**n ( f n -- f^n )
  dup 0= if
    drop fdrop 1e
  else dup 1 and if
    1- fdup recurse f*
  else
    2/ fdup f* recurse
  then then ;

[edit] Haskell

Here's the exponentiation operator from the Prelude:

(^) :: (Num a, Integral b) => a -> b -> a
_ ^ 0           =  1
x ^ n | n > 0   =  f x (n-1) x where
  f _ 0 y = y
  f a d y = g a d  where
    g b i | even i  = g (b*b) (i `quot` 2)
          | otherwise = f b (i-1) (b*y)
_ ^ _           = error "Prelude.^: negative exponent"

There's no difference in Haskell between a procedure (or function) and an operator, other than the infix notation. This routine is overloaded for any integral exponent (which includes the arbitrarily large Integer type) and any numeric type for the bases (including, for example, Complex). It uses the fast "binary" exponentiation algorithm. For a negative exponent, the type of the base must support division (and hence reciprocals):

(^^) :: (Fractional a, Integral b) => a -> b -> a
x ^^ n = if n >= 0 then x^n else recip (x^(negate n))

This rules out e.g. the integer types as base values in this case. Haskell also has a third exponentiation operator,

(**) :: Floating a => a -> a -> a
x ** y = exp (log x * y)

which is used for floating point arithmetic.

[edit] J

J is concretely specified, which makes it easy to define primitives in terms of other primitives (this is especially true of mathematical primitives, given the language's mathematical emphasis).

So we have any number of options. Here's the simplest, equivalent to the for each number, product = product * number of other languages. The base may be any number, and the exponent may be any non-negative integer (including zero):

   exp  =.  */@:#~ 
   
   10 exp 3
1000
   
   10 exp 0
1

We can make this more general by allowing the exponent to be any integer (including negatives), at the cost of a slight increase in complexity:

   exp  =.  *@:] %: */@:(#~|)
   
   10 exp _3
0.001

J also provides for a calculus of functions, permitting us to define exponentiation in its full generality, as the inverse of log (i.e. exp = log^-1):

   exp  =.  ^.^:_1
   
   81 exp 0.5   
  9

Note that the definition above does not use the primitive exponentiation function ^ . The carets in it represent different (but related) things . The function may be parsed as follows: ^. ^: _1 . The first part, ^., is the primitive logarithm operator (e.g. 10 = 1000^.3) .

The second part, ^: , is interesting: it is a "meta operator". It takes two arguments: a function f on its left, and a number N on its right. It produces a new function, which, when given an argument, applies f to that argument N times. For example, if we had a function increment, then increment^:3 X would increment X three times, so the result would be X+3.

In the case of ^. ^: _1 , f is ^. (i.e. logarithm) and N is -1. Therefore we apply log negative one times or the inverse of log once (precisely as in log^-1).

Similarly, we can define log as the reverse of the inverse of root. That is, x pow y = y root^-1 x:

   exp  =.  %:^:_1~
   
   81 exp 0.5
  9

Compare this with the previous definition: it is the same, except that %: , root, has been substituted for ^. , logarithm, and the arguments have been reversed (or reflected) with ~.

That is, J is telling us that power is the same as the reflex of the inverse of root, exactly as we'd expect.

One last note: we said these definitions are the same as ^ in its full generality. What is meant by that? Well, in the context of this puzzle, it means both the base and exponent may be any real number. But J goes further than that: it also permits complex numbers.

Let's use Euler's famous formula, e^pi*i = -1 as an example:

   pi =. 3.14159265358979323846
   e  =. 2.71828182845904523536
   i  =. 2 %: _1                  NB.  Square root of -1
   
   e^(pi*i)
_1

And, as stated, our redefinition is equivalent:

   exp =. %:^:_1~
   
   e exp (pi*i)
_1

[edit] Java

Java does not support operator definition. This example is unoptimized, but will handle negative exponents as well. It is unnecessary to show int^int since an int in Java will be cast as a double.

public class Exp{
   public static void main(String[] args){
      System.out.println(pow(2,30));
      System.out.println(pow(2.0,30)); //tests
      System.out.println(pow(2.0,-2));
   }

   public static double pow(double base, int exp){
      if(exp < 0) return 1 / pow(base, -exp);
      double ans = base;
      for(exp -= 1;exp > 0;--exp) ans *= base;
      return ans;
   }
}

Output:

1.073741824E9
1.073741824E9
0.25

[edit] Logo

to int_power :n :m
  if equal? 0 :m [output 1]
  if equal? 0 modulo :m 2 [output int_power :n*:n :m/2]
  output :n * int_power :n :m-1
end

[edit] Lucid

Some misconceptions about Lucid

pow(n,x)
   k = n fby k div 2;
   p = x fby p*p;
   y =1 fby if even(k) then y else y*p;
   result y asa k eq 0;
end