Power function (X ^ Y) and root function (X ^ (1/Y)) for integers in Erlang

bignum_root.erl

%% Power function (X ^ Y) and root function (X ^ (1/Y)) for
%% integers in Erlang
%% by Kenji Rikitake <[email protected]> 26-JAN-2010
%% modified by Hynek Vychodil <[email protected]> 2-FEB-2010
%% modified by Kenji Rikitake <[email protected]> 3-FEB-2010
%% modified by Hynek Vychodil <[email protected]> 24-JUL-2016
%% Distributed under MIT license at the end of the source code.

-module(bignum_root).

-export([power/2, root/2, sqrt/1]).

%% if an integer exceeds DIGITS, then do the significand/exponent split
%% computation. The value 300 is decided by the overflowing test of
%% math:log/1 and math:sqrt/1.
-define(DIGITS, 300).

%% computing X ^ Y

power(X, 0) when is_integer(X) -> 1;
power(X, Y) when is_integer(X), is_integer(Y), Y > 0 ->
    Bits = bits(Y, []),
    power(X, Bits, X).

power(_, [], Acc) -> Acc;
power(X, [0|Bits], Acc) ->
    power(X, Bits, Acc*Acc);
power(X, [1|Bits], Acc) ->
    power(X, Bits, Acc*Acc*X).

bits(1, Acc) -> Acc;
bits(Y, Acc) ->
    bits(Y div 2, [Y rem 2 | Acc]).

%% sqrt/1 as a special case of root/2
sqrt(X) -> root(X, 2).

%% computing bignum root estimation
estimate_root(X, Y) when is_integer(X), Y >= 2 ->
    % estimation on the larger side
    L = integer_to_list(X),
    Len = length(L),
    case Len > ?DIGITS of
	% prevent math module function overflow
	% by splitting the significand and exponent
        true ->
            Exp = Len - ?DIGITS,
            M = list_to_integer(lists:sublist(L, ?DIGITS)),
            RM = math:exp(math:log(M) / Y) * math:pow(10, (Exp rem Y) / Y),
	    R = list_to_integer(integer_to_list(trunc(RM) + 1) ++ 
				 lists:duplicate(Exp div Y, $0));
        false ->
	    % estimation on the larger side
	    case Y =:= 2 of
		true ->
		    R = trunc(math:sqrt(X)) + 1;
		false ->
		    R = trunc(math:exp(math:log(X) / Y)) + 1
	    end
    end,
    R.

%% computing X ^ (1/Y)
%% with Newton-Raphson method

root(X, 1) when is_integer(X) -> X;
root(X, 2) when is_integer(X), X > 0 -> sqrt(X, estimate_root(X, 2));
root(X, Y) when is_integer(X), X > 0, is_integer(Y), Y >= 3 ->
    % estimation of the root by the float calc
    root(X, Y, estimate_root(X, Y)).

root(X, Y, E) ->
    % Newton-Raphson method of solving (Y)th-root
    EP = power(E, Y - 1),
    Err = E * EP - X,
    E2 = E - Err div (Y * EP),
    if
	E2 =/= E -> root(X, Y, E2);
	Err > 0 -> E - 1;	% Solve rounding error
	true -> E
    end.

sqrt(X, E) ->
    % Newton-Raphson method of solving square root
    Err = E * E - X,
    E2 = E - Err div (2 * E),
    if
	E2 =/= E -> sqrt(X, E2);
	Err > 0 -> E - 1;	% Solve rounding error
	true -> E
    end.

%%  MIT License:
%%  Copyright (c) 2010 by Kenji Rikitake.
%%  Copyright (c) 2010 by Hynek Vychodil.
%%
%%  Permission is hereby granted, free of charge, to any person
%%  obtaining a copy of this software and associated documentation
%%  files (the "Software"), to deal in the Software without
%%  restriction, including without limitation the rights to use,
%%  copy, modify, merge, publish, distribute, sublicense, and/or sell
%%  copies of the Software, and to permit persons to whom the
%%  Software is furnished to do so, subject to the following
%%  conditions:
%%
%%  The above copyright notice and this permission notice shall be
%%  included in all copies or substantial portions of the Software.
%%
%%  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
%%  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
%%  OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
%%  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
%%  HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
%%  WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
%%  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
%%  OTHER DEALINGS IN THE SOFTWARE.