% FILE:     pquadn.erl
% PURPOSE:  recursive integral for N cores
% EXAMPLE:                integrand        A  B   Tol    N
%  1> pquadn:run(fun(X)->math:sqrt(X) end, 1, 2, 0.0001, 3).
%  0.6666666666660148
%  (2/3)x^(3/2) = (2/3)[1-0] = .666...
% METHOD:   Assuming that about N processors are available for this
% computation, spawn N processes and apportion part of the interval of
% integration to each.  Sum the results of the subintegrals.
% COMPILE:  c(pquadn).  
% DEBUG:    c(pquadn, {d,debug}).
% AUTHOR:   mckeeman@cs.dartmouth.edu

-module(pquadn).
-export([run/5]).

% set up integration subintervals, start workers, receive results.
run(F, A, B, Tol, N) ->
  Dx = (B-A)/N,                           % initial intervals size
  Xs = [A+I*Dx || I <- lists:seq(0, N-1)],% N sample points
  Master = self(),
  Vals = [receive {Worker,Val} -> Val end 
    || Worker <- [spawn(fun() -> 
      Master ! {self(), qVal(F, inf, X, F(X), X+Dx, F(X+Dx), Tol/N, 1)} end)
       || X <- Xs]],
  lists:sum(Vals).

% sequential computation of quadrature
qVal(F, Old, A, FA, B, FB, Tol, Lvl) ->
    M  = (A+B)/2,                         % midpoint
    FM = F(M),
    R  = (FA + 4*FM + FB)*(B-A)/6,        % Simpson's approx
    Diff = case Old of
      inf -> 2*Tol;                       % force once
      _   -> abs(R-Old/2)                 % compare new vs. old
    end,
    case Lvl<30 andalso Diff>Tol of  
      true ->                             % more work to do
          qVal(F, R, A, FA, M, FM, Tol/2, Lvl+1) 
        + qVal(F, R, M, FM, B, FB, Tol/2, Lvl+1);
      false -> R                          % good enough answer
    end.