import graph;
import patterns;
usepackage("amsmath");
//usepackage("amsfonts");
//usepackage("MinionPro","fourierbb,mathlf,textosf");
usepackage("fourier");
add("hatch",hatch(4mm));

size(300,00);
//unitsize(2cm);

struct fctesc {
  real [] sub ;
  real [] vsub;
  real[] vpt;
}

fctesc operator init() {return new fctesc;}
fctesc dftesc (real[] x, real []y, real [] z)
{
  fctesc f;
  f.sub=x;
  f.vsub=y;
  f.vpt=z;
  return f;
  }

real[] submin (real [] s1, real[] s2)
{
  real [] fs;
  int n1=s1.length;
  int n2=s2.length;
  int i=0;
  int j=0;
  bool t2=true;
  bool t1=true;
     do {
      if(t1 && t2) {
        if(s1[i]<s2[j]) {fs.push(s1[i]); i=i+1;}
        else {
                  if(s1[i]>s2[j]) {fs.push(s2[j]); j=j+1;}
                  else {fs.push(s1[i]); i=i+1; j=j+1;}
                  }
        } else {
          if(t1) {fs.push(s1[i]); i=i+1;};
          if(t2)  {fs.push(s2[j]);  j=j+1;};
        }
      if(i==n1) {t1=false;}
      if(j==n2) {t2=false;}
      }
      while (t1 || t2);
     return fs;
};

fctesc fsubmin(fctesc f)
{
  real [] sub=f.sub;
  real [] vsub=f.vsub;
  real [] vpt=f.vpt;
  real [] x,y,z;
  int n=sub.length;
  x.push(sub[0]);
  y.push(vsub[0]);
  z.push(vpt[0]);
  int n1=sub.length;
  int k=0;
  for (int i=1;i<n-1;++i ) {
    if((y[k]!=vsub[i]) || (y[k]!=vpt[i]) || (vpt[i]!=vsub[i])) {
      x.push(sub[i]);
      y.push(vsub[i]);
      z.push(vpt[i]);
      ++k;
    }
  }
  x.push(sub[n-1]);
  z.push(vpt[n-1]);
  fctesc ff;
  ff.sub=x;ff.vsub=y;ff.vpt=z;
  return ff;
}
real vfctesc( fctesc f, real s)
{
  real r;
  int n=f.sub.length;
  if( (s<f.sub[0]) || (s>f.sub[n-1])) { abort("cavapas");};
  int i=search(f.sub,s);
  if(i==f.sub.length-1) {
    if(s==f.sub[i]) {r=f.vpt[i];} else {
      r=f.vsub[n-2];
    }
  };
  if(i<f.sub.length-1) {r=f.vsub[i];
    if(s==f.sub[i]) {r=f.vpt[i];};};
  return r;
}
fctesc somme(fctesc f1, fctesc f2)
{
  if(f1.sub[0]!=f2.sub[0]) {abort("cavapas");};
  if(f1.sub[f1.sub.length-1]!=f2.sub[f2.sub.length-1]) {abort("cavpas");};
  real []x,y,z;  
  x=submin(f1.sub,f2.sub);
  for (int i=0;i<x.length-1 ; ++i) {
    z.push(vfctesc(f1,x[i])+vfctesc(f2,x[i]));
    y.push(vfctesc(f1,(x[i]+x[i+1])/2)+vfctesc(f2,(x[i]+x[i+1])/2));
  }
  fctesc f;
  z.push(vfctesc(f1,x[x.length-1])+vfctesc(f2,x[x.length-1]));
  f.sub=x;
  f.vpt=z;
  f.vsub=y;
  return f;
}

fctesc soust(fctesc f1, fctesc f2)
{
  if(f1.sub[0]!=f2.sub[0]) {abort("cavapas");};
  if(f1.sub[f1.sub.length-1]!=f2.sub[f2.sub.length-1]) {abort("cavapas");};
  real []x,y,z;  
  x=submin(f1.sub,f2.sub);
  for (int i=0;i<x.length-1 ; ++i) {
    z.push(vfctesc(f1,x[i])-vfctesc(f2,x[i]));
    y.push(vfctesc(f1,(x[i]+x[i+1])/2)-vfctesc(f2,(x[i]+x[i+1])/2));
  }
  fctesc f;
  z.push(vfctesc(f1,x[x.length-1])-vfctesc(f2,x[x.length-1]));
  f.sub=x;
  f.vpt=z;
  f.vsub=y;
  return f;
}

fctesc produit(fctesc f1, fctesc f2)
{
  if(f1.sub[0]!=f2.sub[0]) {abort("cavapas");};
  if(f1.sub[f1.sub.length-1]!=f2.sub[f2.sub.length-1]) {abort("cavapas");};
  real []x,y,z;  
  x=submin(f1.sub,f2.sub);
  for (int i=0;i<x.length-1 ; ++i) {
    y.push(vfctesc(f1,x[i])*vfctesc(f2,x[i]));
    z.push(vfctesc(f1,(x[i]+x[i+1])/2)*vfctesc(f2,(x[i]+x[i+1])/2));
  }
  fctesc f;
  z.push(vfctesc(f1,x[x.length-1])*vfctesc(f2,x[x.length-1]));
  f.sub=x;
  f.vpt=z;
  f.vsub=y;
  return f;
}

void trace (picture pic=currentpicture, fctesc f, pen p=currentpen)
{
  fctesc fm;
  fm=fsubmin(f);
  int n=fm.sub.length;
  real [] sub=fm.sub;
  real [] vsub=fm.vsub;
  real [] vpt=fm.vpt;
  path mg=scale(1.4)*subpath(W..N..E..S..cycle,1,3);
  frame mmg; 
  draw(mmg,mg);
  path md=scale(1.4)*subpath(unitcircle,1,3);
  frame mmd;
  draw(mmd,md);
  pen ip=invisible;
  //pen mmg=makepen(scale(0.8mm)*mg);
  for (int i=0;i<n-1 ;++i ) {
    pair z1,z2;
    z1=(sub[i],vsub[i]);
    z2=(sub[i+1],vsub[i]);
    if((vpt[i]==vsub[i]) && (vpt[i+1]==vsub[i])) {
      draw(pic,z1--z2,p);
      dot(pic,z1,p);
      dot(pic,z2,p);
    };
    if((vpt[i]==vsub[i]) && (vpt[i+1]!=vsub[i])) {
      draw(pic,z1--z2,p,EndDotMargin);
      draw(pic,z2,ip,marker(mmd));
      dot(pic,(sub[i],vpt[i]),p);
    };
    if((vpt[i]!=vsub[i]) && (vpt[i+1]==vsub[i])) {
      draw(pic,z1--z2,p,BeginDotMargin);
      draw(pic,z1,ip,marker(mmg));
      dot(pic,(sub[i],vpt[i]),p);
    };
    if((vpt[i]!=vsub[i]) && (vpt[i+1]!=vsub[i])) {
      draw(pic,z1--z2,p,DotMargins);
      draw(pic,z1,ip,marker(mmg));
      draw(pic,z2,ip,marker(mmd));
      dot(pic,(sub[i],vpt[i]),p);
    };
  }
  dot(pic,(sub[n-1],vpt[n-1]),p);
}

int ttest(real x)
{
  int v;
  if(x<0) {v=-1;};
  if(x>0) {v=1;};
  if(x==0) {v=0;};
  return v;
}

void tracei (picture pic=currentpicture, fctesc f, pen p=currentpen)
{
  add("hatch",hatch(4mm,p));
  add("antihatch",hatch(4mm,NW,p));
  fctesc fm;
  fm=fsubmin(f);
  int v;
  int n=fm.sub.length;
  real [] sub=fm.sub;
  real [] vsub=fm.vsub;
  v=ttest(vsub[0]);
  path g;
  g=(sub[0],0)--(sub[0],vsub[0])--(sub[1],vsub[0]);
  int w;
  for (int i=1;i<n-1 ;++i ) {
    w=ttest(vsub[i]);
    if(w==v) {g=g--(sub[i],vsub[i])--(sub[i+1],vsub[i]);
      if(i==n-2) {g=g--(sub[i+1],0)--cycle;fill(pic,g,pattern("hatch"));
      }
    };
    if(w!=v) {g=g--(sub[i],0)--cycle;
      if(v!=0) {fill(pic,g,pattern("hatch"));};
      //if(v<0) {fill(pic,g,pattern("antihatch"));};
      g=(sub[i],0)--(sub[i],vsub[i])--(sub[i+1],vsub[i]);
      v=ttest(vsub[i]);
    }
  }
}

void encadesc(picture pic=currentpicture, fctesc f1, fctesc f2, pen p=currentpen, bool special=false)
{
  add("hatch",hatch(3mm,p));
  add("antihatch",hatch(3mm,NW,p));
  fctesc fm1;
  fctesc fm2;
  fm1=fsubmin(f1);
  fm2=fsubmin(f2);
  int v;
  int n1=fm1.sub.length;
  int n2=fm2.sub.length;
  real [] sub1=fm1.sub;
  real [] vsub1=fm1.vsub;
  real [] sub2=fm2.sub;
  real [] vsub2=fm2.vsub;
  path g;
  for (int i=0;i<n1-1 ;++i ) {
    g=g--(sub1[i],vsub1[i])--(sub1[i+1],vsub1[i]);
  }
  for (int i=n2-2;i>=0 ; i=i-1) {
    g=g--(sub2[i+1],vsub2[i])--(sub2[i],vsub2[i]);
  }
  g=g--cycle;
  if(special) {fill(pic,g,p);}
else {
  fill(pic,g,pattern("hatch"));
};
}

void tracesubd(picture pic=currentpicture, real [] a, fctesc f, string mt="x",pen p=currentpen)
{
  int n=a.length;
  string stt;
  draw((a[0],0)--(a[0],vfctesc(f,a[0]+.00001)),dashed);
  stt="$"+mt+"_"+string(0)+"$";
  label(stt,(a[0],0),N,Fill(white));
  for (int i=1; i<n-1; ++i) {
    draw((a[i],0)--(a[i],vfctesc(f,a[i]-.00001)),dashed);
    draw((a[i],0)--(a[i],vfctesc(f,a[i]+.00001)),dashed);
  }
  for (int i=1; i<n-1; ++i) {
    stt="$"+mt+"_"+string(i)+"$";
    label(stt,(a[i],0),N,Fill(white));
  }
  draw((a[n-1],0)--(a[n-1],vfctesc(f,a[n-1]-.00001)),dashed);
  stt="$"+mt+"_"+string(n-1)+"$";
  label(stt,(a[n-1],0),N,Fill(white));
}


fctesc fs,fs1,fs2,fs3;
//real [] a={1,2,3,4};
//real [] b={0,1,4};
//real []c={1,0,2,3};

//real [] a1={1,2.3,2.5,3.3,3.6,4}, b1={1,4,5,-1,0}, c1={0,4,-1,2,2,3};
//fs=dftesc(a, b,c);
//fs2=produit(fs,fs1);
//fs2=somme(fs,fs1);
//write(fs2.sub,fs2.vsub,fs2.vpt);
//fctesc nfs=fsubmin(fs2);
//write(nfs.sub,nfs.vsub,nfs.vpt);
//write(a,b);
//write(c);

//trace(fs2);
//xlimits(1,4);
//ylimits(-1,7);
//xaxis("$x$", YZero,Ticks);
//yaxis("$y$",Left,Ticks);
//trace(fs2,red);
//tracei(fs2,red);

//real [] b3={1.2,4.3,5.5,-0.3,0.2}, c3={0,4,-1,2,2,3};
//fs3=dftesc(a1,b3,c3);
//encadesc(fs1,fs3);
//trace(fs1,red);
//tracei(fs1,blue);



/// Exemple 1, fonction en escalier, subdivision adaptée

real [] a={1,2.5,3,3.4,4,4.7};
real [] b={-1,-0.5,-0.5,1,2};
real [] c={0,1,-0.5,1.5,-1,2};


fs=dftesc(a,b,c);
xlimits(0,5);
ylimits(-1.2,2.2);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",Left,Ticks);
trace(fs,red);
tracesubd(a,fs);
shipout("escalier1");
erase();


// Exemple 2, fonction en escalier, subdivision (minimale) adaptée


real [] a={1,2.5,3.4,4,4.7};
real [] b={-1,-0.5,1,2};
real [] c={0,1,1.5,-1,2};


fs=dftesc(a,b,c);
xlimits(0,5);
ylimits(-1.2,2.2);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",Left,Ticks);
trace(fs,red);
tracesubd(a,fs);
shipout("escalier2");
erase();



// Exemple 3, fonction en escalier, subdivision non adaptée


real [] a={1,2.5,3.4,4,4.7};
real [] b={-1,-0.5,1,2};
real [] c={0,1,1.5,-1,2};
real [] a1={1,2,3,4,4.7};

fs=dftesc(a,b,c);
xlimits(0,5);
ylimits(-1.2,2.2);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",Left,Ticks);
trace(fs,red);
tracesubd(a1,fs);
shipout("escalier3");
erase();




// Exemple 4, fonction en escalier, subdivision non adaptée


real [] a={1,2,3,4,4.7};
real [] b={-1,-1,1,2};
real [] c={0,1,1.5,-1,2};
real [] a1={1,3,4,4.7};

fs=dftesc(a,b,c);
xlimits(0,5);
ylimits(-1.2,2.2);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",Left,Ticks);
trace(fs,red);
tracesubd(a1,fs);
shipout("escalier4");
erase();



/// Exemple 5, fonction en escalier, subdivision adaptée

real [] a={1,2.5,3,3.4,4,4.7};
real [] b={-1,-0.5,-0.5,1,2};
real [] c={0,1,-0.5,1.5,-1,2};


fs=dftesc(a,b,c);
xlimits(0,5);
ylimits(-1.2,2.2);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",Left,Ticks);
tracei(fs);
trace(fs,red);
fctesc fm=fsubmin(fs);
 real [] sub=fm.sub;
  real [] vsub=fm.vsub;
//write(sub,vsub);
//tracesubd(a,fs);
shipout("escalier5");
erase();


// Exemple 6 et 7, encadrement d'un fonction par deux fonctions en escalier

real f(real x) {return (1+x^2);}

xlimits(-2,2);
ylimits(0,5.5);
real [] xpt;
xpt=sequence(-10,10)/5;
real [] vpt1,vpt2,vsub1,vsub2;
for (int i=0; i<11; ++i) {
  vpt1[i]=f(xpt[i]);
  vpt1[10+i]=f(xpt[10+i]);
  vpt2[i]=f(xpt[i+1]);
  vpt2[10+i]=f(xpt[10+i-1]);
}
vpt2[10]=1;
vpt1[10]=f(xpt[11]);
for (int i=0;i<10 ; ++i) {
  vsub1[i]=vpt1[i];
  vsub1[i+10]=vpt1[i+11];
  vsub2[i]=vpt2[i];
  vsub2[i+10]=vpt2[i+11];
}
fs1=dftesc(xpt,vsub1,vpt1);
fs2=dftesc(xpt,vsub2,vpt2);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",XZero,Ticks);
encadesc(fs1,fs2);
trace(fs1,red);
trace(fs2,green);
draw(graph(f,-2,2));
shipout("escalier6");
erase();
xlimits(-2,2);
ylimits(-0.5,1);
fctesc fs3=soust(fs1,fs2);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",XZero,Ticks);
tracei(fs3);
trace(fs3,blue);
//write(fs3.vpt);
//write(vfctesc(fs2,2-0.1));
//write (search(fs2.sub,2-0.1),fs2.sub.length);
shipout("escalier7");
erase();





// Exemple 6a et 7a, encadrement d'un fonction par deux fonctions en escalier

real f(real x) {return (1+x^2);}
int n=20;
xlimits(-2,2);
ylimits(0,5.5);
real [] xpt;
xpt=sequence(-2*n,2*n)/n;
real [] vpt1,vpt2,vsub1,vsub2;
for (int i=0; i<2*n+1; ++i) {
  vpt1[i]=f(xpt[i]);
  vpt1[2*n+i]=f(xpt[2*n+i]);
  vpt2[i]=f(xpt[i+1]);
  vpt2[2*n+i]=f(xpt[2*n+i-1]);
}
vpt2[2*n]=1;
vpt1[2*n]=f(xpt[2*n+1]);
for (int i=0;i<2*n ; ++i) {
  vsub1[i]=vpt1[i];
  vsub1[i+2*n]=vpt1[i+2*n+1];
  vsub2[i]=vpt2[i];
  vsub2[i+2*n]=vpt2[i+2*n+1];
}
fs1=dftesc(xpt,vsub1,vpt1);
fs2=dftesc(xpt,vsub2,vpt2);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",XZero,Ticks);
encadesc(fs1,fs2,blue,true);
trace(fs1,red);
trace(fs2,green);
draw(graph(f,-2,2));
shipout("escalier6a");
erase();
//size(300,0);
unitsize(1.5cm,6cm);
xlimits(-2,2);
ylimits(-0.1,.3);
fctesc fs3=soust(fs1,fs2);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",XZero,Ticks);
tracei(fs3);
trace(fs3,blue);
//write(fs3.vpt);
//write(vfctesc(fs2,2-0.1));
//write (search(fs2.sub,2-0.1),fs2.sub.length);
shipout("escalier7a");
erase();







// Exemple 8 fonction petite entraîne intégrale petite mais réciproque fausse
size(300,0);
unitsize(0,0);
size(300,0);
xlimits(0,5);
ylimits(-0.3,7);
real [] a={0,.5,1,1.5,2,2.5,3,3.5,4,4.5,4.8,5};
real [] b={.1,.05,.2,.3,.2,.3,.2,.3,.2,.1,7};
real [] c={1,1,1,1,1,1,1,1,1,1,1,1};
fs=dftesc(a,b,c);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",Left,Ticks);
tracei(fs);
trace(fs,blue);
shipout("escalier8");
erase();


// Exemple 9 fonction tordue mais intégrable
xlimits(0,1);
ylimits(-0.1,1.1);
xaxis("$x$", YZero,Ticks);
yaxis("$y$",Left,Ticks);
int n=10;
for (int i=1; i<n+1;++i ) {
  for (int j=1; j<i+1;++j ) {
    dot((j/i,1/i),blue);
  }
}
path g=(0,0)--(0,1/n)--(1,1/n)--(1,0)--cycle;
filldraw(g,pattern("hatch"));
string leg="$f(x)=\begin{cases} 0 & \text{ si $x=0$ ou
$x\in\mathbb{R}\setminus\mathbb{Q}$} \\ 1/q & \text{ si $x=p/q$ avec $p\wedge q=1$}\end{cases}$";
label(leg,(0.05,.9),E);
string leg2="$\displaystyle\int_0^1 f(x)dx=0$";
label(leg2,(0.3,.7),E);
shipout("escalier9");
erase();