public
Last active — forked from pichi/bignum_root.erl

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

  • Download Gist
bignum_root.erl
Erlang
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
%% Power function (X ^ Y) and root function (X ^ (1/Y)) for
%% integers in Erlang
%% by Kenji Rikitake <kenji.rikitake@acm.org> 26-JAN-2010
%% modified by Hynek Vychodil <vychodil.hynek@gmail.com> 2-FEB-2010
%% modified by Kenji Rikitake <kenji.rikitake@acm.org> 3-FEB-2010
%% 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, Y) when is_integer(X), is_integer(Y), Y >= 0 ->
power(X, Y, 1).
 
power(_, 0, Acc) ->
Acc;
power(X, Y, Acc) when Y rem 2 =:= 1 ->
power(X, Y - 1, Acc * X);
power(X, Y, Acc) ->
power(X * X, Y div 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.

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.