> ?polygon
# This Maple sheet allows you to generate pictures of the system of
# triangles made along the way by a policy of alternately converting the
# middle triangle of all red triangles to blue, and then converting the
# middle of all blue triangles to red. Note that it doesn't settle down
# to a uniform beige. There are, in the end, two kinds of triangles, one
# mostly red, the other mostly blue, and they are mirror images of each
# other. The connection to induction is that you can prove that there is
# a limiting percentage of the area that ends up red; it's 4/7. The are
# calculation is this: In the blue fairy's move, the pair [r,b] of
# market shares (adding to 1), is replaced with the new pair
# [r,b]->[3*r/4,b+r/4]. The next step carries the new
# [r,b]->[r+b/4,3*b/4].
> blushare:=u->[3*u[1]/4,u[2]+u[1]/4];

             blushare := u -> [3/4 u[1], u[2] + 1/4 u[1]]

> 
> redshare:=u->[u[1]+u[2]/4,3*u[2]/4];

             redshare := u -> [u[1] + 1/4 u[2], 3/4 u[2]]

> bothsteps:=redshare@blushare;

                    bothsteps := redshare@blushare

> bothsteps([1,0]);

                               13
                              [--, 3/16]
                               16

> bothsteps(");

                               181  75
                              [---, ---]
                               256  256

> 

> 
> 
> 
> evalf(");

                      [.7070312500, .2929687500]

> bothsteps([r,b]);

                     13
                    [-- r + 1/4 b, 3/4 b + 3/16 r]
                     16

> solve({r*13/16+b/4=r,r+b=1},{r,b});

                          {r = 4/7, b = 3/7}

# To get at the question of whether this steady state is reached and if
# so how fast, introduce the relevant quantity: r-4/7=rs (for rsurplus).
# Then r=4/7+rs,b=3/7-rs. 
# Thus 13/16*r+b/4-4/7=13/16*(rs+4/7)+(3/7-rs)/4-4/7=9/16*rs
> 13/16*(rs+4/7)+(3/7-rs)/4-4/7;

                               9/16 rs

# And this tells you that the red surplus, over the value it will in the
# limit have, drops by a factor of *9/16 over each two-move cycle. 
# 
# The rest of this sheet is Maple code which actually does the steps and
# produces pictures. If you're working with a black and white monitor or
# printer, use this definition of
# makepix2:=subs({red=white,blue=black},****the rest of the code for
# makepix2*******)end;
> with(plottools);

[arc, arrow, circle, cone, cuboid, curve, cutin, cutout, cylinder,

    disk, dodecahedron, ellipse, ellipticArc, hemisphere, hexahedron,

    hyperbola, icosahedron, line, octahedron, pieslice, point,

    polygon, rectangle, rotate, scale, semitorus, sphere, stellate,

    tetrahedron, torus, transform, translate]

> 
> 
> with(plots);

[animate, animate3d, changecoords, complexplot, complexplot3d,

    conformal, contourplot, contourplot3d, coordplot, coordplot3d,

    cylinderplot, densityplot, display, display3d, fieldplot,

    fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d,

    inequal, listcontplot, listcontplot3d, listdensityplot, listplot,

    listplot3d, loglogplot, logplot, matrixplot, odeplot, pareto,

    pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d,

    polyhedraplot, replot, rootlocus, semilogplot, setoptions,

    setoptions3d, spacecurve, sparsematrixplot, sphereplot, surfdata,

    textplot, textplot3d, tubeplot]


> midpt:=proc(u,v)[(u[1]+v[1])/2,(u[2]+v[2])/2]end;

  midpt := proc(u, v) [1/2*u[1] + 1/2*v[1], 1/2*u[2] + 1/2*v[2]] end

> 
> 
> bluefairy:=proc(triangle)local
> tcolor,poly,m1,m2,m3,p1,p2,p3,t1,t2,t3,t4;
> tcolor:=triangle[2];if tcolor=red then
> poly:=triangle[1];p1:=poly[1];p2:=poly[2];p3:=poly[3];m1:=midpt(p1,p2)
> ;m2:=midpt(p2,p3);m3:=midpt(p3,p1);t1:=[[p1,m1,m3],red];t2:=[[m1,p2,m2
> ],red];t3:=[[m2,p3,m3],red];t4:=[[m1,m2,m3],blue];[t1,t2,t3,t4];RETURN
> (t1,t2,t3,t4)else RETURN(triangle)fi end;

bluefairy := proc(triangle)
local tcolor, poly, m1, m2, m3, p1, p2, p3, t1, t2, t3, t4;
    tcolor := triangle[2];
    if tcolor = red then
        poly := triangle[1];
        p1 := poly[1];
        p2 := poly[2];
        p3 := poly[3];
        m1 := midpt(p1, p2);
        m2 := midpt(p2, p3);
        m3 := midpt(p3, p1);
        t1 := [[p1, m1, m3], red];
        t2 := [[m1, p2, m2], red];
        t3 := [[m2, p3, m3], red];
        t4 := [[m1, m2, m3], blue];
        [t1, t2, t3, t4];
        RETURN(t1, t2, t3, t4)
    else RETURN(triangle)
    fi
end

> bl:=bluefairy([[[0,0],[1,0],[1/2,evalf(sqrt(3)/2)]],red]);

bl := [[[0, 0], [1/2, 0], [1/4, .4330127020]], red],

    [[[1/2, 0], [1, 0], [3/4, .4330127020]], red], [

    [[3/4, .4330127020], [1/2, .8660254040], [1/4, .4330127020]], red

    ], [[[1/2, 0], [3/4, .4330127020], [1/4, .4330127020]], blue]

> makepix2:=proc(u)polygon(u[1],color=u[2])end;

         makepix2 := proc(u) polygon(u[1], color = u[2]) end

> redfairy:=proc(triangle)local
> tcolor,poly,m1,m2,m3,p1,p2,p3,t1,t2,t3,t4;
> tcolor:=triangle[2];if tcolor=blue then
> poly:=triangle[1];p1:=poly[1];p2:=poly[2];p3:=poly[3];m1:=midpt(p1,p2)
> ;m2:=midpt(p2,p3);m3:=midpt(p3,p1);t1:=[[p1,m1,m3],blue];t2:=[[m1,p2,m
> 2],blue];t3:=[[m2,p3,m3],blue];t4:=[[m1,m2,m3],red];[t1,t2,t3,t4];RETU
> RN(t1,t2,t3,t4)else RETURN(triangle)fi end;

redfairy := proc(triangle)
local tcolor, poly, m1, m2, m3, p1, p2, p3, t1, t2, t3, t4;
    tcolor := triangle[2];
    if tcolor = blue then
        poly := triangle[1];
        p1 := poly[1];
        p2 := poly[2];
        p3 := poly[3];
        m1 := midpt(p1, p2);
        m2 := midpt(p2, p3);
        m3 := midpt(p3, p1);
        t1 := [[p1, m1, m3], blue];
        t2 := [[m1, p2, m2], blue];
        t3 := [[m2, p3, m3], blue];
        t4 := [[m1, m2, m3], red];
        [t1, t2, t3, t4];
        RETURN(t1, t2, t3, t4)
    else RETURN(triangle)
    fi
end

> rbl:=map(redfairy,[bl]);

rbl := [[[[0, 0], [1/2, 0], [1/4, .4330127020]], red],

    [[[1/2, 0], [1, 0], [3/4, .4330127020]], red], [

    [[3/4, .4330127020], [1/2, .8660254040], [1/4, .4330127020]], red

    ], [[[1/2, 0], [5/8, .2165063510], [3/8, .2165063510]], blue], [

    [[5/8, .2165063510], [3/4, .4330127020], [1/2, .4330127020]],

    blue], [

    [[1/2, .4330127020], [1/4, .4330127020], [3/8, .2165063510]],

    blue], [

    [[5/8, .2165063510], [1/2, .4330127020], [3/8, .2165063510]], red

    ]]

> nops([bl]);

                                  4

> bl1:=[bl];

bl1 := [[[[0, 0], [1/2, 0], [1/4, .4330127020]], red],

    [[[1/2, 0], [1, 0], [3/4, .4330127020]], red], [

    [[3/4, .4330127020], [1/2, .8660254040], [1/4, .4330127020]], red

    ], [[[1/2, 0], [3/4, .4330127020], [1/4, .4330127020]], blue]]

> recursionstep:=proc(rb)local
> k,rbr,rbrb;rbr:=[seq(bluefairy(rb[k]),k=1..nops(rb))];rbrb:=[seq(redfa
> iry(rbr[k]),k=1..nops(rbr))];rbrb end;

recursionstep := proc(rb)
local k, rbr, rbrb;
    rbr := [seq(bluefairy(rb[k]), k = 1 .. nops(rb))];
    rbrb := [seq(redfairy(rbr[k]), k = 1 .. nops(rbr))];
    rbrb
end

> r2:=recursionstep(rbl):
> nops(r2);

                                  40

> pix2:=map(makepix2,r2):display(pix2);

> r3:=recursionstep(r2):pix3:=map(makepix2,r3):display(pix3);

> r4:=recursionstep(r3):pix4:=map(makepix2,r4):display(pix4);


>