aitkenDEL:=proc(x,del)
local dx,DEL;
dx:=del(x);
DEL:=-(dx^2)/(del(x+dx)-dx);
return(DEL);
end:


aitkenlmt:=proc(left,right,del,limprecis)
local x,tt;

Digits:=20;
x:=left:
tt:=aitkenDEL(x,del); 
while abs(tt)>limprecis and (x+tt)>left and x<right do  
   x:=x+tt;
   tt:=aitkenDEL(x,del);
end do:
if (x>=right) then
  x:=right;
else
  x:=evalf((10^(-4))*ceil((x)*10^4));
end;
return(x);
end:


lmt:=proc(left,right,tdel)
local x,td;
global intset,tagdiv,limprecis,intctr,lll;

Digits:=20;
lll:=0.0;
x:=left: 
td:=evalf(tdel(x));
while td>limprecis and x>=left and x<right do  
  x:=evalf(x+td);
  lll:=x;
  td:=evalf(tdel(x));
end do:

if (x>=right) then
  x:=right;
else
  x:=evalf((10^(-4))*ceil((x)*10^4));
end;

return(x);
end:

# 
weakkurz:=proc(a,b,f,maxdp,delta)
global intset,tagdiv,limprecis,intctr,lll;
local md,dm,r,width,ans,i,gotone,dl,dr,de,left,fv,
      right,khsum,depth,tmp,tmpintset,h,ndiv,thisb,L,mx,s,origintset;

left:=a;
khsum:=0;
origintset:=intset;
while (left<b) do
 dl:=evalf(delta(left));
 if (left+dl>b) then
   dl:=b-left;
 end if;
 fv:=f(left);
 khsum:=khsum+fv*dl;
 intset:={[[left,0],[left,fv],[left+dl,fv],[left+dl,0]]} union intset;
 tagdiv:=[op(tagdiv),[[left,left+dl],left]];
 intctr:=intctr+1;
 printf("WEAK FINE: [left,right] = [%lg,%lg]\n",left,left+dl); 
 left:=left+dl;
end do;
return(khsum); 
end:



semistrongkurz:=proc(a,b,f,maxdp,delta)
global intset,tagdiv,limprecis,intctr,lll;
local md,dm,r,width,ans,i,gotone,dl,dr,de,left,fv,
      right,khsum,depth,tmp,tmpintset,h,ndiv,thisb,L,mx,s,gotmax,origintset;

left:=a;
khsum:=0;
origintset:=intset;
while (left<b) do
 right:=b;
 while (left<right) do 
    gotone:=0;
    depth:=0;
    gotmax:=0;
    while (gotone=0) do
      width:=right-left;
      md:=left+right;
      md:=0.5*md;
      dl:=delta(left);
      dr:=delta(right);
      dm:=delta(md);
      L:=[left+dl,right+dr,md+dm];
      if ((dl+dr)<width) then 
        L[1]:=left-1.0;
        L[2]:=left-1.0;
      end if;  
      if((dm+dm) < width) then
        L[3]:=left-1.0;
      end if; 
      mx:=max(op(L));    
      if (mx>=left) then
        gotone:=1;
      elif (depth=maxdp) then
        printf("HIT MAX DEPTH!\n");
        gotmax:=1;
        gotone:=1;
       else
        right:=(left+right)/2.0;
        depth:=depth+1;
      end if;
    end do;

    r:=min(b,mx);

    if (gotmax=1) then
      r:=evalf((10^(-4))*ceil((right)*10^4));
      if (delta(r)>=(r-left)) then
        fv:=f(r);
        ans:=fv*(r-left);
        intset:={[[left,0],[left,fv],[r,fv],[r,0]]} union intset; 
        tagdiv:=[op(tagdiv),[[left,r],r]];
      else
        printf("WARNING: NOT FINE ON %lg to %lg\n",left,r);
      end if;
    elif (mx=L[3]) then
      fv:=f(md);
      ans:=fv*(r-left);
      intset:={[[left,0],[left,fv],[r,fv],[r,0]]} union intset; 
      tagdiv:=[op(tagdiv),[[left,r],md]];
    elif (mx=L[1]) then
      fv:=f(left);
      ans:=fv*(r-left);
      intset:={[[left,0],[left,fv],[r,fv],[r,0]]} union intset; 
      tagdiv:=[op(tagdiv),[[left,r],left]];
    else
      s:=min(left+dl,right);
      fv:=f(left);
      ans:=fv*(s-left);
      intset:={[[left,0],[left,fv],[s,fv],[s,0]]} union intset; 
      tagdiv:=[op(tagdiv),[[left,s],left]];
      fv:=f(right);
      ans:=ans+fv*(r-s);
      intset:={[[s,0],[s,fv],[r,fv],[r,0]]} union intset; 
      tagdiv:=[op(tagdiv),[[s,r],right]];
    end if;

    intctr:=intctr+1;
    right:=r;
    khsum:=khsum+ans;  
    left:=right;
  end do;
end do;
return(khsum); 
end:






strongkurz:=proc(a,b,f,maxdp,delta)
global intset,tagdiv,limprecis,intctr,lll;
local md,dm,r,width,ans,i,gotone,dl,dr,de,left,fv,
      right,khsum,depth,tmp,tmpintset,h,ndiv,thisb,L,mx,s,origintset;

left:=a;
khsum:=0;
origintset:=intset;
while (left<b) do
 right:=b;
 while (left<right) do 
    gotone:=0;
    depth:=0;
    while (gotone=0) do
      width:=right-left;
      md:=left+right;
      md:=0.5*md;
      dl:=delta(left);
      dr:=delta(right);
      dm:=delta(md);
      L:=[left+dl,right+dr,md+dm];
 #     printf("left,middle,right,width,delleft,delmid,delright=%f %f %f    %f   %f %f %f\n",left,md,right,width,dl,dm,dr);
      if ((dl+dr)<width) then 
        L[1]:=left-1.0;
        L[2]:=left-1.0;
      end if;  
      if((dm+dm) < width) then
        L[3]:=left-1.0;
      end if; 
      mx:=max(op(L));    
      if (mx>=left) then
        gotone:=1;
      elif (depth=maxdp) then
        printf("HIT MAX DEPTH!\n");
        intset:=origintset;
        gotone:=1;
        return(infinity);
      else
        right:=(left+right)/2.0;
        depth:=depth+1;
      end if;
    end do;

    r:=min(b,mx);
    if (mx=L[3]) then
      fv:=f(md);
      ans:=fv*(r-left);
      intset:={[[left,0],[left,fv],[r,fv],[r,0]]} union intset; 
      tagdiv:=[op(tagdiv),[[left,r],md]];
    elif (mx=L[1]) then
      fv:=f(left);
      ans:=fv*(r-left);
      intset:={[[left,0],[left,fv],[r,fv],[r,0]]} union intset;
      tagdiv:=[op(tagdiv),[[left,r],left]];
    else
      s:=min(left+dl,right);
      fv:=f(left);
      ans:=fv*(s-left);
      intset:={[[left,0],[left,fv],[s,fv],[s,0]]} union intset;
      tagdiv:=[op(tagdiv),[[left,s],left]];
      if (abs(r-s)>0.0000001) then
        fv:=f(right);
        ans:=ans+fv*(r-s);
        intset:={[[s,0],[s,fv],[r,fv],[r,0]]} union intset;
        tagdiv:=[op(tagdiv),[[s,r],right]];
      end if;
    end if;

    intctr:=intctr+1;
    right:=r;
    khsum:=khsum+ans;  
    left:=right;
  end do;
end do;
return(khsum); 
end:

split:=proc(a,b,f,delta)
global intset,tagdiv,limprecis,intctr,lll;
local rmdr,maxdp,ok,tmpintset,khsum,left,right,ans,s,t,fv;

tmpintset:={};
left:=a;
khsum:=0;
maxdp:=15;
while (left<b) do

  ok:=1;
  right:=left;
  try
    right:=timelimit(5,lmt(left,b,delta));
  catch:
    printf("TIME out! lll=%lf\n",lll);
    ok:=0;
  end;

  if (right<=left) then
    ok:=0;
  end if;
  if (ok=1) then  
    rmdr:=evalf(right-delta(right));
    if (rmdr>left) then
      ans:=strongkurz(left,rmdr,f,maxdp,delta);
      if (ans<>infinity) then
        khsum:=khsum+ans;
      else
        ans:=weakkurz(left,rmdr,f,maxdp,delta);
      end if;
    end if;
    s:=max(left,rmdr);
    t:=min(b,evalf(right+delta(right)));
    fv:=f(right);
    ans:=fv*(t-s);
    khsum:=khsum+ans;
    intset:={[[s,0],[s,fv],[t,fv],[t,0]]} union intset;
    tagdiv:=[op(tagdiv),[[s,t],right]];
    left:=t;
  else
    ok:=1;
    left:=evalf(left);
    try
  #right:=timelimit(5,aitkenlmt(left,b,delta,limprecis));
                     right:=timelimit(5,aitkenlmt(lll,b,delta,limprecis));
    catch:
      printf("Time out!\n");
      ok:=0;
    end;
    if (right<=left) then
      right:=b;
    end if;
    rmdr:=evalf(right-delta(right));    
    if (rmdr>left) then
      ans:=semistrongkurz(left,rmdr,f,maxdp,delta);      
    end if;
    s:=max(left,rmdr);
    t:=min(b,evalf(right+delta(right)));
    fv:=f(right);
    ans:=fv*(t-s);
    khsum:=khsum+ans;
    intset:={[[s,0],[s,fv],[t,fv],[t,0]]} union intset;
    tagdiv:=[op(tagdiv),[[s,t],right]];
    left:=evalf(t);
  end if;   
end do;
return(khsum); 
end:

# 
khbox:=proc()
  global intset,tagdiv,limprecis,intctr,lll;

  local rrr,G,H,f,a,b,delta,i,Z; 
  if nargs < 4 then 
    printf("This routine needs at least 4 arguments: the function f, the endpoints\n");
    printf("of the interval a,b and the delta function.\n");
    return;
  end if;
  intset:={};
  tagdiv:={};
  f:=args[1]; 
  a:=args[2]; 
  b:=args[3]; 
  delta:=args[4];  
  limprecis:=(2)*10^(-9);
  intctr:=0;
  rrr:=evalf(split(a,b,f,delta));
  intset:=evalf(intset):
  tagdiv:=evalf(tagdiv):
  G:=PLOT(seq(POLYGONS( evalf( op(i,intset) ), COLOR(RGB,0.0,1.0,0.0) ),i=1..nops(intset))):  
  setoptions(symbol=CIRCLE,symbolsize=25):
  if nargs>4 then
    H:=plot(f,a..b,args[5..nargs],discont=true);
  else  
    H:=plot(f,a..b,discont=true);
  end if; 
  #:
    display({G,H});
end:

khsum:=proc()
  global intset,tagdiv,limprecis,intctr,lll;

  local rrr,G,H,f,a,b,delta,i,Z; 
  if nargs < 4 then 
    printf("This routine needs at least 4 arguments: the function f, the endpoints\n");
    printf("of the interval a,b and the delta function.\n");
    return;
  end if;
  intset:={};
  tagdiv:={};
  f:=args[1]; 
  a:=args[2]; 
  b:=args[3]; 
  delta:=args[4];  
  limprecis:=(2)*10^(-9);
  intctr:=0;
  rrr:=evalf(split(a,b,f,delta));
  intset:=evalf(intset):
  tagdiv:=evalf(tagdiv):
  printf("\nThe Riemann sum is   %f\n\n",evalf(rrr));
end: