/* Program to accompany the article Wicklin, R. (2016) "Visualize the Cantor function" published 29JUN2016 on The DO Loop blog: http://blogs.sas.com/content/iml/2016/06/29/viz-cantor-function-in-sas.htm */ proc iml; /* rows of B contain all 8-digit combinations of 0s and 2s */ B = expandgrid({0 2}, {0 2}, {0 2}, {0 2}, {0 2}, {0 2}, {0 2}, {0 2}); B = B[2:nrow(B),]; /* remove first row of zeros */ k = ncol(B); /* k = 8 */ v = 3##(-(1:k)); /* vector of powers 3^{-i} */ t = B * v`; /* x values: right-hand endpts of middle-third intervals */ u = 2##(-(1:k)); /* vector of powers 2^{-i} */ f = B/2 * u`; /* y values: Cantor function on Cantor set */ call series(t, f); /* graph the Cantor function */ quit; /* More elaborate program that adds tick marks to the graph */ proc iml; B = expandgrid({0 2}, {0 2}, {0 2}, {0 2}, {0 2}, {0 2}, {0 2}, {0 2}); *B = expandgrid({0 2}, {0 2}, {0 2}, {0 2}); B = B[2:nrow(B),]; /* remove first row of zeros */ k = ncol(B); v = 3##(-(1:k)); /* vector of powers 3^{-i} */ t = B * v`; /* some points in the Cantor set */ u = 2##(-(1:k)); /* vector of powers 2^{-i} */ f = B/2 * u`; /* values of Cantor function on Cantor set */ title "Cantor Function"; call series(t, f) yvalues=do(0,1.01,0.125); /* For each row of a matrix A, find the rightmost column that contains a positive number. For example: A = {1 0 1 1 0 0 0, 0 1 0 1 1 1 0, 1 0 0 0 0 0 0}; p = RightmostPosCol(A); * = {4, 6, 1} ; The reason: - In the first row, column 4 is the last nonpositive value in row - In the second row, column 6 is the last nonpositive value in row - In the third row, column 1 is the last nonpositive value in row */ start RightmostPosCol(A); cols = j(nrow(A), 1, 0); do j = 1 to ncol(A); idx = loc(A[,j] > 0); if ncol(idx)>0 then cols[idx] = j; end; return (cols); finish; /* remove first row of zeros */ depth = RightmostPosCol(B); z = t - 3##(-depth); /* left hand endpts of middle-third intervals */ x = 0 // colvec(z || t) // 1; /* points in Cantor set */ y = 0 // colvec(f || f) // 1; /* function values */ /* extract labels for the first few intervals */ d = 1 // colvec(depth || depth) // 1; /* step at which interval was added */ idx = loc(d <= 2); xvalues = x[idx]; xlabels = '"' + putn(xvalues, "fract5.") + '"'; idx = loc(d <= 3); yvalues = unique(y[idx]); ylabels = '"' + putn(yvalues, "fract5.") + '"'; *print xvalues xlabels yvalues ylabels; create Cantor var {x y}; append; close; submit xvalues xlabels yvalues ylabels; proc sgplot data=Cantor2; series x=x y=y; xaxis values=(&xvalues) valuesdisplay=(&xlabels) grid; yaxis values=(&yvalues) valuesdisplay=(&ylabels) grid; run; endsubmit; quit;