% --------------------------------------------------------------
% Erlang module 'series'
%
% Approximates the infinite sum:
%   sum( (n^a)/(b^n) ) where n <- 1 .. infinity
% --------------------------------------------------------------

-module( series).
-export( [
   sum_na_over_bn/2,
   sum_na_over_bn_numerator/2,
   close_sum_na_over_bn_numerator/2,
   list_close_sum_na_over_bn_numerator_for_a/1,
   list_close_sum_na_over_bn_numerator_for_b/1,
   list_exact_numerator_when_a_is_1/0,
   list_exact_numerator_when_a_is_2/0,
   exact_numerator_when_a_is_1/1,
   exact_numerator_when_a_is_2/1,
   exact_when_a_is_1/1,
   exact_when_a_is_2/1
]).

% --------------------------------------------------------------
% Math helper functions
% --------------------------------------------------------------

% --------------------------------------------------------------
% get_close_int( A )
%
%   If A is a float that is very close to an integer, returns
%   the integer. This doesn't just round, A has to be very close
%   to the integer value. Otherwise returns 'not_close'.

get_close_int( A ) when is_integer( A) -> A;
get_close_int( A ) when is_float( A) ->
  Delta = 0.0001,
  Close = trunc( if A < 0 -> A - Delta; true -> A + Delta end),
  Diff  = abs( A - Close),
  case Diff =< (2 * Delta) of
    true  -> Close;
    false -> not_close
  end.

% --------------------------------------------------------------
% a_to_the_nth_power( X, N )
%
%   Simple function to calculate (X^N). Uses logs like a
%   sliderule, so the answer is not exact even when A and N
%   are both integers.
%
%   Works for X > 0.
%   Works for any value of N.
%
%   Improvement: You could define this for (X == 0), as long as
%   (N /= 0) (don't forget 0.0 =/= 0).
%
%   Improvement: You could define this for (X < 0) as long as
%   N is an integer.
%
%   Improvment: If both A and N are integers maybe this should
%   always return an integer? It's probably what people expect.
%   And the following can overflow, while overflow is unlikely
%   if we stay in integers whenever possible.

a_to_the_nth_power( A, N ) ->
  math:exp( math:log( A) * N).

% --------------------------------------------------------------
% sum_a_thru_b( Fn, A, B )
%
%   Sums Fn( N) where N goes from A to B inclusive.
%   A and B should be integers, A <= B.
%   Fn( N) takes one integer arg and returns something
%   (probably a number) that works with lists:sum(..).
%
%   sum_a_thru_b( fun( N ) -> N end, 1, 10) returns 55
%   (which is 1+2+3+4+5+6+7+8+9+10).

sum_a_thru_b( Fn, A, B ) ->
  lists:sum( lists:map( Fn, lists:seq( A, B))).

% --------------------------------------------------------------
% Functions to calculate results without summing
%
%   Calculates the formulas:
%     B/((B-1)^(A+1))
%     (B*(B+1))/((B-1)^(A+1))
% --------------------------------------------------------------

% --------------------------------------------------------------

exact_denominator( A, B ) ->
  a_to_the_nth_power( B - 1, A + 1).

% --------------------------------------------------------------

exact_numerator_when_a_is_1( B ) -> B.
exact_numerator_when_a_is_2( B ) -> B * (B + 1).

% Returns the list:
%   [2,3,4,5,6,7,8,9,10,11]
list_exact_numerator_when_a_is_1( ) ->
  lists:map(
    fun exact_numerator_when_a_is_1/1,
    lists:seq( 2, 11)).

% Returns the list:
%   [6,12,20,30,42,56,72,90,110,132]
list_exact_numerator_when_a_is_2( ) ->
  lists:map(
    fun exact_numerator_when_a_is_2/1,
    lists:seq( 2, 11)).

% --------------------------------------------------------------

exact_when_a_is_1( B ) ->
  exact_numerator_when_a_is_1( B) / exact_denominator( 1, B).

exact_when_a_is_2( B ) ->
  exact_numerator_when_a_is_2( B) / exact_denominator( 2, B).

% --------------------------------------------------------------
% Functions to calculate terms in the series
%
%   The series looks like this:
%     sum( (N^A)/(B^N) ) where N <- 1 .. infinity
%
%   A can be any number.
%   B must be > 1.
%     -1 < B =< 1 clearly does not converge.
%     B =< 0 does not work because a_to_the_nth_power( X, N)
%       does not work when X <= 0 (log(0) is negative infinity).
%       Although we could change a_to_the_nth_power( X, N) to
%       work when (X < 0) and N is an integer.
% --------------------------------------------------------------

% --------------------------------------------------------------
% na_over_bn( A, B, N )
%
%   Calculates: (N^A)/(B^N)

na_over_bn( A, B, N ) ->
  a_to_the_nth_power( N, A) /
  a_to_the_nth_power( B, N).

% --------------------------------------------------------------
% sum_na_over_bn( A, B )
%
%   Sums: (N^A)/(B^N)
%     for N <- 1 .. 100 (inclusive)

sum_na_over_bn( A, B ) ->
  sum_a_thru_b(
    fun( N ) -> na_over_bn( A, B, N) end,
    1,
    100).

% --------------------------------------------------------------
% sum_na_over_bn_numerator( A, B )
%
%   Returns sum( (N^A)/(B^N) ) (N <- 1..100) multiplied by
%   (B-1)^(A+1). We do this because:
%
%     sum( (N^0)/(B^N) ) is 1/((B-1)^1)
%     sum( (N^1)/(B^N) ) is B/((B-1)^2)
%     sum( (N^2)/(B^N) ) is (B*(B+1))/((B-1)^3)
%
%   We are interested in the numerator series:
%     1   B   B(B+1)  ..?..
%   The next term is NOT B(B+1)(B+2) like you might expect.
%   When B is 2 the series looks like this:
%     1   2   6   26   150   1082   9366   94586   1091670

sum_na_over_bn_numerator( A, B ) ->
  sum_na_over_bn( A, B ) *
  a_to_the_nth_power( B - 1, A + 1).

% --------------------------------------------------------------
% close_sum_na_over_bn_numerator( A, B )
%
%   Pretties the result from sum_na_over_bn_numerator( A, B )
%   by rounding the float numbers to nearby integers.
%   The floats have to be very close to integer values or they
%   will not come thru.

close_sum_na_over_bn_numerator( A, B ) ->
  get_close_int(
    sum_na_over_bn_numerator( A, B )).

% --------------------------------------------------------------
% list_close_sum_na_over_bn_numerator_for_a( A )
%
%   Generates the series:
%     sum( (N^A)/(B^N) ) for N <- 1..100 (infinity)
%   Multiplied by:
%     (B-1)^(A+1)
%   For values of B <- 2 .. 11 (inclusive)
%
%   For A <- 2, result is:
%     [6,12,20,30,42,56,72,90,110,132]
%   See http://www.research.att.com/~njas/sequences/A002378
%   See http://mathworld.wolfram.com/PronicNumber.html
%
%   For A <- 3, result is:
%     [26,66,132,230,366,546,776,1062,1410,1826]
%
%   For A <- 4, result is:
%     [150,480,1140,2280,4074,6720,10440,15480,22110,30624]
%
%   For A <- 1, result is:
%     [2,3,4,5,6,7,8,9,10,11]
%
%   For A <- 0, result is:
%     [1,1,1,1,1,1,1,1,1,1]
%
%   Result is not integer when A is not an integer,
%   or when A < 0.

list_close_sum_na_over_bn_numerator_for_a( A ) ->
  lists:map(
    fun( B ) -> close_sum_na_over_bn_numerator( A, B) end,
    lists:seq( 2, 11)).

% --------------------------------------------------------------
% list_close_sum_na_over_bn_numerator_for_b( B )
%
%   Generates the series:
%     sum( (N^A)/(B^N) ) for N <- 1..100 (infinity)
%   Multiplied by:
%     (B-1)^(A+1)
%   For values of A <- 1 .. 10 (inclusive)
%
%   For B <- 2, result is:
%     [2,6,26,150,1082,9366,94586,1091670,14174522,204495126]
%   See http://mathworld.wolfram.com/
%         StirlingNumberoftheSecondKind.html
%   See http://www.research.att.com/~njas/sequences/A000629
%
%   For B <- 3, result is:
%     [3,12,66,480,4368,47712,608016,8855040,145083648,
%      2641216512]
%   See http://www.research.att.com/~njas/sequences/A009362
%
%   For B <- 4, result is:
%     [4,20,132,1140,12324,160020,2424132,41967540,817374564,
%      17688328020]
%
%   For B <- 5, result is:
%     [5,30,230,2280,28280,421680,7336880,145879680,3263031680,
%      81097294080]

list_close_sum_na_over_bn_numerator_for_b( B ) ->
  lists:map(
    fun( A ) -> close_sum_na_over_bn_numerator( A, B) end,
    lists:seq( 1, 10)).