> restart;
> with(plots):
> 
> 
> 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:

Warning, the name changecoords has been redefined
> 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):
>   printf("\nThe Riemann sum is   %f\n",evalf(rrr));
>   G:=PLOT(seq(POLYGONS( evalf( op(i,intset) ), COLOR(RGB,0.0,1.0,0.0) ),i=1..nops(intset))):  
> #  with(plots):
>   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; 
> #  with(plots):
>     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:
> 
> 
> f1:=proc(x)
> 
> local t;
> if floor(sqrt(x))=ceil(sqrt(x)) then 
>   t:=evalf(floor(sqrt(x)));
>   return(abs(cos(Pi*x/50))*t);
> end if;
> return(-1);
> end:
> 
> dd1:=proc(x)
> local t,z;
> if f1(x)=-1 then
>   t:=ceil(sqrt(x));
>   t:=t*t;
>   z:=floor(sqrt(x));
>   z:=z*z;
>   return(min(0.9*(t-x),0.9*(x-z)));
> end if;
> return( evalf(min(1/2, 1/(1+ f1(x)) ) ));
> end:
> 
> 
> khbox(f1,0,10,dd1,colour=blue,thickness=2); 

The Riemann sum is   -3.501973

> ddel:=proc(x)
> local ddel;
> if (x=0) then
>   ddel:=1/200;
> else
>   ddel:=0.2*abs(x);
> end if;
> end:
> 
> f:=proc(x)
> local f;
> if (x=0) then
>   f:=0;
> else
>   f:=1/sqrt(abs(x));
> end if;
> end:
> 
> khbox(f,-0.2,0.2,ddel,y=0..30,colour=blue,thickness=2);


The Riemann sum is   1.484081

> ff:=x->piecewise(x<=0,1/(1+10*x^2),0<x,0.5);
>  ff1:=x->abs(D(ff)(x));
>  tdel:=x->piecewise(x<0,min(1.5,.15/(1+ff(x)+ff1(x))),x=0,.01,x>0,abs(x));
>  khbox(ff,-2,1,tdel,colour=blue,thickness=2);
                          /            1                \
      ff := x -> piecewise|x <= 0, ---------, 0 < x, 0.5|
                          |                2            |
                          \        1 + 10 x             /
                     ff1 := x -> |D(ff)(x)|
                      /          /            0.15       \               
tdel := x -> piecewise|x < 0, min|1.5, ------------------|, x = 0, 0.01, 
                      \          \     1 + ff(x) + ff1(x)/               

            \
  0 < x, |x||
            /

The Riemann sum is   0.926181

# 
# 
> del2:=proc(x)
> local del2,n;
> if (x=0) then
>   del2:=1/30;
> else 
>   n:=evalf(floor(1/x));
>   if (x=evalf(1/n)) then
>     del2:=4^(-n);
>   else
>     del2:=evalf(min(1/n-x,x-1/(n+1)));
>   end if;
> end if;
> # printf("x=%lf and delta=%lf\n",x,del2);
> end:
> 
> f2:=proc(x)
> local f2,n;
> if (x=0) then
>   f2:=0;
> else
>   n:=evalf(floor(1/x));
>   f2:=n*(-1)^n;
> end if;
> end:
> 
> 
> f4:=x->piecewise(x<=0,0,floor(1/x)*(-1)^(floor(1/x)));
> khbox(f4,0.,.5,del2,y=-30..30,colour=blue,thickness=2);
                         /                             /1\\
                         |                        floor|-||
                         |                /1\          \x/|
     f4 := x -> piecewise|x <= 0, 0, floor|-| (-1)        |
                         \                \x/             /

The Riemann sum is   0.089011

> 
> 
>