(* ::Package:: *)

   function [Q,fEvals] = PeriodicSimpson(f,a,b,P,M4,tol)
 % f is a handle to a function f that is defined on [a,b] and whose 
 % fourth derivative is bounded by M4. f has period P, i.e. f(x+P) = f(x) for all x.
 % Q is an estimate of I, the exact integral of f from a to b.
 % Q satisfies
 %                   | I - Q | <= tol where tol>0.
 % f-evals is the number of required f-evaluations
 
 % Let I(c,d) = integral of f from c to d
 
 % Determine k (integer) and fraction tau (0<=tau<1) so b = a + (k+tau)P
 k = floor((b-a)/P);
 tau = (b-(a+k*P))/P;
 % Thus b = a + (k+tau)P
 
 if k==0
     % Need |S(a,b,n) - I| <= tol. This is guaranteed if
     %   (b-a)^5*M4/(2880*n^4) <= tol
     % The best choice...
     n = ceil(((b-a)^5*M4/(2880*tol))^.25);
     Q = CompSimpson(f,a,b,n);
     fEvals = 2*n+1;
 else
    % I(a,b) = I(a,a+kP)  + I(a+kP,a+(k+tau)P)
    %        = k*I(a,a+P) + I(a,a+tau*P)
    %        = k*(I(a,a+tau*P) + I(a+tau*P,a+P)) + I(a,a+tau*P)
    %        = (k+1)*I(a,a+tau*P) + k*I(a+tau*P,a+P)
    
    % Our plan is to apply  composite simpson:
    
    % Q =  (k+1)S(a,a+tau*P,n1) + k*S(a+tau*P,a+P,n2)
    
    % A bound for the error in Q is therefore
    %
    %    (k+1)((M4/2880)*(tau*P)^5)/n1^4 + k*(M4/2880)*((1-tau)*P)^5/n2^4 
    % 
    % We must choose n1 and n2 so that this is less than tol, i.e.,
    
    %  (k+1)*tau^5/n1^4  + k (1-tau)^5 <= tol*2880/(M4*p^5)
    
    % Define the upper bound B by
    B = tol*2880/(M4*P^5);
    % Let lambda satisfy 0<lambda<1. If
    %           (k+1)*tau^5/n1^4 <= lambda*B
    % and 
    %            k (1-tau)^5/n2^4  <= (1-lambda)*B
    % then
    %          (k+1)*tau^5/n1^4  + k (1-tau)^5/n2^4 <= B
    % What is the best lambda? 
    % Since the total number of f-evals is (2n1+1)+(2n2+1)
    % We can estimate a good choice by trying a bunch of different
    % lambdas...
   
    lambdaVals = linspace(0,1,101);
    for j=1:101   
        lambda = lambdaVals(j);
        % For this lambda, figure out the optimal n1 and n2....
        % (k+1)*tau^5/n1^4 <= lambda*B   
        n1 = ceil(((k+1)*tau^5/(lambda*B))^.25 );
        %  k (1-tau)^5/n2^4  <= (1-lambda)*B
        n2 = ceil( ((k)*(1-tau)^5/((1-lambda)*B))^.25);
        % A call to CompSimpson with n subinterbvals involves 2n+1 f-evals
        work(j) = 2*(n1+n2+1);
    end
    [minWork,jmin] = min(work);
    lambda = lambdaVals(jmin);
    n1 = max(ceil(((k+1)*tau^5/(lambda*B))^.25 ),1);
    n2 = max(ceil( ((k)*(1-tau)^5/((1-lambda)*B))^.25),1);
   
    Q = (k+1)*CompSimpson(f,a,a+tau*P,n1) + k*CompSimpson(f,a+tau*P,a+P,n2);
    fEvals = 2*(n1+n2+1);
 end

% -2 if the required accuracy is not achieved for k = 0.
% -1 if the required accuracy is not achieved for k = 10 and k = 1000.
% -1 if more than 2000 evaluations for any k = 10 example.
% -1 if more than 7000 evaluations for any k = 1000 example.
