> with(plots,polygonplot3d,pointplot,display3d):
# CantorSequence constructs a list containing the sequence p^i, i=0, 1,
# ... N-1 where p<1/2
# CantorSequence(4,1/3) will produce the first four elements of the
# sequence used to
# construct the Cantor ternary set [1, 1/3, 1/9, 1/27].
> CantorSequence := proc(N:integer,p)
> local i,lst:
> lst := []:
> for i from 0 to N-1 do
> lst := [op(lst),p^i]:
> od:
> lst:
> end:
# Draw a solid 3D box in color farbe.  p1 and p2 are lists containing
# the upper-front-left and lower-back-right corners respectively.
> box3d := proc(p1,p2,farbe)
> local half1,half2:
> half1 := polygonplot3d([ [p1[1],p1[2],p1[3]], [p1[1],p1[2],p2[3]],
> [p1[1],p2[2],p2[3]], [p1[1],p2[2],p1[3]], [p1[1],p1[2],p1[3]],
> [p2[1],p1[2],p1[3]], [p2[1],p2[2],p1[3]], [p1[1],p2[2],p1[3]],
> [p1[1],p1[2],p1[3]], [p2[1],p1[2],p1[3]],
> [p2[1],p1[2],p2[3]],[p1[1],p1[2],p2[3]], [p1[1],p1[2],p1[3]]
> ],style=PATCH,color=farbe):
> half2 := polygonplot3d([ [p2[1],p2[2],p2[3]], [p2[1],p2[2],p1[3]],
> [p2[1],p1[2],p1[3]], [p2[1],p1[2],p2[3]], [p2[1],p2[2],p2[3]],
> [p1[1],p2[2],p2[3]], [p1[1],p1[2],p2[3]], [p2[1],p1[2],p2[3]],
> [p2[1],p2[2],p2[3]], [p1[1],p2[2],p2[3]],
> [p1[1],p2[2],p1[3]],[p2[1],p2[2],p1[3]], [p2[1],p2[2],p2[3]]
> ],style=PATCH,color=farbe):
> {half1, half2}:
> end:
# GraphPsiHat draws the graph of the Cantor function on a 3D Cantor set
# defined by the
# sequence cseq.  GraphPsiHat(CantorSeq(4,1/3)); will draw the graph of
# the Cantor ternary
# function.  NOTE: the graph is 4 dimensional.  The function value is
# represented by the color.
# Blue is the lowest value and red is the highest with gradations being
# in between.
> GraphPsiHat := proc(cseq:list)
> local n,k,k1,k2,k3,crds,crds2,these,box1,box2,box3,h,farbe:
> global all,plotme:
> all := []:
> crds := [0,1]:
> for n from 2 to nops(cseq) do
> crds2 := []:
> these := []:
> for k from 0 to nops(crds)/2-1 do
> crds2 :=
> [op(crds2),crds[2*k+1],crds[2*k+1]+cseq[n],crds[2*k+2]-cseq[n],crds[2*
> k+2]]:
> these := [op(these),
> [crds[2*k+1],crds[2*k+1]+cseq[n],crds[2*k+2]-cseq[n],crds[2*k+2]] ]:
> od:
> crds := crds2:
> for k1 from 1  to nops(these) do
> for k2 from 1 to nops(these) do
> for k3 from 1 to nops(these) do
> h :=  ((2*k1-1)/2^(n-1)+(2*k2-1)/2^(n-1)+(2*k3-1)/2^(n-1))/3:
> all := [op(all), [ these[k1],these[k2],these[k3], h ]]:
> od:
> od:
> od:
> od:
> plotme := {}: 
> for n from 1 to nops(all) do
> farbe := COLOUR(RGB,all[n][4],0,1-all[n][4]):
> box1 :=
> box3d([all[n][1][1],all[n][2][1],all[n][3][2]],[all[n][1][4],all[n][2]
> [4],all[n][3][3]],farbe):
> box2 :=
> box3d([all[n][1][2],all[n][2][1],all[n][3][1]],[all[n][1][3],all[n][2]
> [4],all[n][3][4]],farbe):
> box3 :=
> box3d([all[n][1][1],all[n][2][2],all[n][3][1]],[all[n][1][4],all[n][2]
> [3],all[n][3][4]],farbe):
> plotme := plotme union {op(box1),op(box2),op(box3)}:
> od:
> display3d(plotme):
> end:
> GraphPsiHat(CantorSequence(3,1/3));
>