Real Analysis

Written by Gianluca Gorni, University of Udine, Italy. Latest revision: March 31, 2022.

Description: This file contains a package that helps study some objects of interest in Advanced Real Analysis: the (ternary) Cantor set, Vitalis's function ("Devil's staircase"), Peano's curve, and Van der Waerden's continuous but nowhere differentiable function.

Related MathSource notebooks: The Cantor Set and Mathematica, The Peano Curve (I still have to understand if it is the same curve as the one I define here), Plotting the Monsters of Real Analysis.

Instructions: load the package by clicking on the button below, then evaluate the examples.

Button[StyleForm["Click here to load package!",FontColor->RGBColor[1,0,0],FontSize->14],ButtonFunction:>FrontEndExecute[FrontEndToken["EvaluateInitialization"]],ButtonEvaluator->None,Active->True]

Mathematica package

Preamble

BeginPackage["RealAnalysis`"];

ClearAll["`*"];

CantorSetQ::usage=
"CantorSetQ[x] decides for an exact number x whether it belongs to the (ternary) Cantor Set. When x is rational the procedure terminates, while for x irrational it may not terminate.";

CantorSet::usage=
"CantorSet is a domain, like Integers, Rationals etc. If you write Element[x,CantorSet] for an exact number x Mathematica will try to decide whether x belongs to the (ternary) Cantor Set. When x is rational the procedure terminates, while for x irrational it may not terminate.";

Options[CantorSetPlot]={PlotStyle->{},Heights->Automatic,Tree->False};

CantorSetPlot::usage=
"CantorSetPlot[{n1,n2,...},opts] draws a plot of the approximants to Cantor Set of orders n1,n2,.... The options are the same as for Graphics, plus PlotStyle (behaving as in Plot) and Heights (Heights->{h1,h2,...} decides the vertical positions of the approximants; the default is Automatic).";

Vitali::usage=
"Vitali[x] calculates the value of Vitali's continuous function that has zero derivative almost everywhere but is not constant. When x is a floating-point real number, the result is a floating point number with as much Precision as the input allows. When x is an exact real number the result is an exact rational, unless the recursion limit is exceeded.";

InfInverseVitali::usage=
"InfInverseVitali[x] calculates the minimum of the inverse image of x through Vitali's function. When x is a floating-point real number, the result is a floating point number with as much Precision as the input allows. When x is an exact rational number the result is an exact rational, belonging to the Cantor Set. When x is an exact irrational numbers you can calculate the value of InfInverseVitali[x] to arbitrary precision with N[InfInverseVitali[x],n].
See SupInverseVitali";

SupInverseVitali::usage=
"SupInverseVitali[x] calculates the maximum of the inverse image of x through Vitali's function. When x is a floating-point real number, the result is a floating point number with as much Precision as the input allows. When x is an exact rational number the result is an exact rational, belonging to the Cantor Set. When x is an exact irrational numbers you can calculate the value of SupInverseVitali[x] to arbitrary precision with N[SupInverseVitali[x],n].
See InfInverseVitali";

VitaliPlot::usage=
"VitaliPlot[x0,x1,opts] plots Vitali's continuous nowhere differentiable function, a.k.a. Devil's Staircase, on the interval between x0 and x1. The options are the same as for Graphics, plus PlotPoints (set to a high value) and PlotStyle. The internal routine uses a Line primitive, not Plot.";

myPeanoCurve::usage=
"myPeanoCurve[t] calculates the value of Peano's continuous curve that fills the two-dimensional unit square. When t is a floating-point real number, the result is a couple of numbers, usually one exact rational and one floating-point. When t is rational the result is a couple of exact rationals. When t is an exact irrational number the result is a couple of exact numbers too, unless the RecursionLimit is exceeded.";

InversePeano::usage=
"InversePeano[x,y] calculates a t in the Cantor Set such that myPeanoCurve[t]={x,y}, when the couple x,y is in the unit square. When either x or y is a floating-point real number, the result is floating point too, with as much precision as the imput allows. When x and y are both non-rational exact numbers InversePeano[x,y] is not defined, but you can still calculate the numerical value to arbitrary precision n with N[InversePeano[x,y],n].";

PeanoCurveSeparatePlots::usage=
"PeanoCurveSeparatePlots[n,x0,x1,opts] with n=1 or n=2 plots respectively the first and second component of Peano Curve on the interval between x0 and x1. The options are the same as for Graphics, plus PlotPoints (set to a high value) and PlotStyle. The internal routine uses a Line primitive, not Plot. PeanoCurveSeparatePlots[{1,2},x0,x1,opts] makes a GraphicsArray with the two components.";

PeanoCurvePlot::usage=
"PeanoCurvePlot[n,opts] plots the polygonal trajectory in the unit square of the n-th approximation to Peano's Curve. With the default option RoundCorners->True the corners of the polygonal are rounded, avoiding self-intersections.";

Options[PeanoCurvePlot]={HeadLength->Automatic,RoundCorners->True,RoundingAmount->Automatic,AspectRatio->Automatic,PlotRange->All,PlotStyle->{},Frame->True,Ticks->{{0,1/4,1/2,3/4,1},{0,1/4,1/2,3/4,1}},FrameTicks->{{0,1/4,1/2,3/4,1},{0,1/4,1/2,3/4,1}}};

Options[PeanoCurveMovie]={ShowPath->True,FramesPerSmallestInterval->Automatic};

PeanoCurveMovie::usage=
"PeanoCurveMovie[n] makes a movie of the point moving along the polygonal trajectory in the unit square of the n-th approximation to Peano's Curve. Options are ShowPath->True and FramesPerSmallestInterval->Automatic";

VanDerWaerden::usage=
"VanDerWaerden[x] calculates the value of Van der Waerden's continuous nowhere differentiable function. When x is a floating-point real number, the result is a real number with the same Precision. When x is rational the result is an exact rational. When x is an exact irrational number you can calculate the value of VanDerWaerden[x] to arbitrary precision with N[VanDerWaerden[x],n].";

VanDerWaerdenPlot::usage=
"VanDerWaerdenPlot[x0,x1,opts] plots Van der Waerden's continuous nowhere differentiable function on the interval between x0 and x1. The options are the same as for Graphics, plus PlotPoints (set to a high value) and PlotStyle. The internal routine uses a Line primitive, not Plot.";

Begin["`Private`"];

ClearAll["`*"];

Cantor Set tests

Basics:

CantorSetQ[0]=True;
CantorSetQ[1]=True;

Basic exclusions from the Cantor Set:

CantorSetQ[x_?NumericQ/;Or[Im[x]=!=0,x<0,1/3<x<2/3,x>1]]=False;

Self-similarity recursive formula for rationals with denominator divisible by 3:

CantorSetQ[x_Rational/;IntegerQ[Denominator[x]/3]]:=CantorSetQ[Mod[3x,1]];

Self-similarity recursive formula for rationals with denominator not divisible by 3:

CantorSetQ[x_Rational/;!IntegerQ[Denominator[x]/3]]:=!Or@@Map[1/3<#<2/3&,NestList[Mod[3#,1]&,x,EulerPhi[Denominator[x]]-1]];

Self-similarity recursive formula for irrationals:

CantorSetQ[x_?NumericQ/;(Precision[x]===Infinity&&Or[0<x<1/3,2/3<x<1])]:=CantorSetQ[Mod[3x,1]];

Definition of the domain CantorSet:

CantorSet/:Element[x_,CantorSet]:=With[{q=CantorSetQ[x]},q/;Or[q===True,q===False]];

Cantor Set plot

All possible sequences of n digits, each being either 0 or 2:

digitSequences02[0]={{0}};digitSequences02[n_]:=Flatten[Outer[List,Sequence@@Table[{0,2},{n}]],n-1]

Left and right extrema of the intervals remaining at the n-th step of the construction of the Cantor Set:

leftExtrema[n_]:=Map[FromDigits[#,3]&,digitSequences02[n]]/3^n

leftAndRightExtrema[n_]:=Union[Flatten[{#,#+1/3^n}&[leftExtrema[n]]]];

The intervals left at the n-th step of the construction of the Cantor Set:

intervals[n_]:=Transpose[({#,#+1/3^n})&[leftExtrema[n]]]

The ticks appropriate for Cantor Set plots in a given interval:

CantorTicks[a_,b_,k_:2]/;a<b:=
Module[{n,x2,x3,unitIntervalTicks},
n=k+Max[0,Ceiling[Log[3,1/(b-a)]]];
{x2,x3}=3^(-n){Floor[x0*3^n],Ceiling[x1*3^n]};
unitIntervalTicks=leftAndRightExtrema[n];
Select[Union[Flatten[Map[#+unitIntervalTicks&,Range[Floor[a],Ceiling[b]]]]],(a≤#≤b)&]];

Plotting the approximants to Cantor Set:

CantorSetPlot[{n__Integer?NonNegative},opts___Rule]:=Module[{howMany,lineStyles,hts,raisedIntervals,raisedLines,base},
howMany=Length[{n}];
{lineStyles,hts,treeQ}={PlotStyle,Heights,Tree}/.{opts}/.Options[CantorSetPlot];
Which[Head[lineStyles]=!=List,lineStyles={{lineStyles}},!MatchQ[lineStyles,{{___}..}],lineStyles={{}}];
If[!MatchQ[hts,{__?NumericQ}|Automatic],hts=Automatic];
Which[hts===Automatic,hts=Reverse[Range[howMany]]];
raisedIntervals[k_]:=Map[Flatten,Outer[List,intervals[{n}[[k]]],{hts[[Mod[k,Length[hts],1]]]}],{2}];
raisedLines[k_]:=Join[lineStyles[[Mod[k,Length[lineStyles],1]]],Map[Line,raisedIntervals[k]]];
Graphics[{base=Map[raisedLines,Range[howMany]],If[treeQ,base/.Line[{p1_,p2_}]:>
If[p1[[-1]]>1,{Arrow[{(p1+p2)/2,5p1/6+p2/6-{0,1}}],
Text[0,((p1+p2)/2+5p1/6+p2/6-{0,1})/2,{1,-1}],
Arrow[{(p1+p2)/2,p1/6+5p2/6-{0,1}}],
Text[2,((p1+p2)/2+p1/6+5p2/6-{0,1})/2,{-1,-1}]},{}],Nothing]},FilterRules[{opts},Options@Graphics],AspectRatio->Min[howMany/2,1],Axes->{True,False},AxesOrigin->{0,0},PlotRange->{{0,1},{Min[0,hts-1],Max[hts+1,0]}},Ticks->{CantorTicks[0,1],Automatic}]];

Show[Graphics[{base[[1]],base[[1]]/.Line[{p1_,p2_}]:>
If[p1[[-1]]>1,{Arrow[{(p1+p2)/2,5p1/6+p2/6-{0,1}}],
Text[0,((p1+p2)/2+5p1/6+p2/6-{0,1})/2,{1,-1}],
Arrow[{(p1+p2)/2,p1/6+5p2/6-{0,1}}],
Text[2,((p1+p2)/2+p1/6+5p2/6-{0,1})/2,{-1,-1}]},{}]},
base[[2]]]]

CantorSetPlot[n_Integer?NonNegative,opts___Rule]:=CantorSetPlot[{n},opts];

Calculation of Vitali's function

SetAttributes[Vitali,Listable];

Reduction to the unit interval:

Vitali[x_?NumericQ/;And[Element[x,Reals],Floor[x]=!=0]]:=Floor[x]+Vitali[Mod[x,1]];

Value at the origin:

Vitali[0]=0;

Value in the middle third:

Vitali[x_/;And[Precision[x]===Infinity,1/3≤x≤2/3]]=1/2;

Self-similarity formula for rationals whose denominator is divisible by 3:

Vitali[x_Rational/;And[0≤x<1/3,IntegerQ[Denominator[x]/3]]]:=Vitali[x]=Vitali[Mod[3x,1]]/2;

Vitali[x_Rational/;And[2/3<x<1,IntegerQ[Denominator[x]/3]]]:=Vitali[x]=(1+Vitali[Mod[3x,1]])/2;

Self-similarity formula for rational whose denominator is not divisible by 3:

Vitali[x_Rational/;And[Or[0≤x<1/3,2/3<x<1],!IntegerQ[Denominator[x]/3]]]:=Vitali[x]=Module[{g,rl},
rl=g[y_]:>Which[
0≤y<1/3,g[3y]/2,
1/3≤y≤2/3,1/2,
2/3≤y<1,(1+g[Mod[3y,1]])/2];
g[x]/.Solve[g[x]==Nest[Expand[#/.rl]&,g[x],EulerPhi[Denominator[x]]],g[x]][[1]]];

Self-similarity formula for exact real numbers:

Vitali[x_/;And[Head[x]=!=Rational,Precision[x]===Infinity,0≤x<1/3]]:=Vitali[Mod[3x,1]]/2;

Vitali[x_/;And[Head[x]=!=Rational,Precision[x]===Infinity,2/3<x<1]]:=(1+Vitali[Mod[3x,1]])/2;

Vitali's function for floating-point real numbers. The result is an exact number.

Vitali[x_Real/;0≤x≤1]:=Vitali[FromDigits[RealDigits[x,3],3]];

Inverse Vitali

SetAttributes[InfInverseVitali,Listable];

SetAttributes[SupInverseVitali,Listable];

Basic values:

InfInverseVitali[0]=0;
InfInverseVitali[1/2]=1/3;
InfInverseVitali[1]=1;

SupInverseVitali[0]=0;
SupInverseVitali[1/2]=2/3;
SupInverseVitali[1]=1;

Reduction to the unit interval:

InfInverseVitali[x_?NumericQ/;And[Im[x]===0,Floor[x]=!=0]]:=Floor[x]+InfInverseVitali[Mod[x,1]];

SupInverseVitali[x_?NumericQ/;And[Im[x]===0,Floor[x]=!=0]]:=Floor[x]+SupInverseVitali[Mod[x,1]];

Self-similarity formulas for rationals with even denominator:

InfInverseVitali[x_Rational/;And[0≤x<1/2,EvenQ[Denominator[x]]]]:=InfInverseVitali[2x]/3;

SupInverseVitali[x_Rational/;And[0≤x<1/2,EvenQ[Denominator[x]]]]:=SupInverseVitali[2x]/3;

InfInverseVitali[x_Rational/;And[1/2<x<1,EvenQ[Denominator[x]]]]:=(2+InfInverseVitali[2x-1])/3;

SupInverseVitali[x_Rational/;And[1/2<x<1,EvenQ[Denominator[x]]]]:=(2+SupInverseVitali[2x-1])/3;

Self-similarity formulas for rationals with odd denominator:

InfInverseVitali[x_Rational/;OddQ[Denominator[x]]]:=InfInverseVitali[x]=
Module[{rl,g},
rl=g[y_]:>Which[
0≤y<1/2,g[2y]/3,
y==1/2,1/3,
1/2≤y<1,(2+g[Mod[2y,1]])/3];
g[x]/.Solve[g[x]==Nest[Expand[#/.rl]&,g[x],EulerPhi[Denominator[x]]],g[x]][[1]]];

SupInverseVitali[x_Rational/;OddQ[Denominator[x]]]:=InfInverseVitali[x];

For floating-point real numbers:

InfInverseVitali[x_Real/;0≤x<1]:=Module[{bin,ter},
bin=RealDigits[x,2]//.{a___,1,c:((0)..)}:>{a,0,Sequence@@Map[1&,{c}]};
ter={2bin[[1]],bin[[2]]};
FromDigits[ter,3]+(1.`1)2/3^(Length[ter[[1]]]-ter[[2]])];

SupInverseVitali[x_Real/;0≤x<1]:=Module[{bin,ter},
bin=RealDigits[x,2]//.{a___,0,c:((1)..)}:>{a,1,Sequence@@Map[0&,{c}]};
ter={2bin[[1]],bin[[2]]};
FromDigits[ter,3]+(1.`1)2/3^(Length[ter[[1]]]-ter[[2]])];

Vitali's function plot

Options[VitaliPlot]={PlotPoints->2^5,PlotStyle->{},Axes->True,AspectRatio->Automatic};

VitaliPlot[x0_,x1_,opts___]/;x0<x1:=Module[{segmentNumber,lineOptions,n,x2,x3,pts,tks,verttks},
{segmentNumber,lineOptions}={PlotPoints,PlotStyle}/.{opts}/.Options[VitaliPlot];n=Max[1,Ceiling[Log[2,segmentNumber]-Log[3,x1-x0]]];
{x2,x3}=3^(-n){Floor[x0*3^n],Ceiling[x1*3^n]};pts=Transpose[{#,Vitali[#]}&[Union[Join[CantorTicks[x2,x3,n],{x2,x3}]]]];
tks=CantorTicks[x0,x1];
verttks=Vitali[tks];
Show[Graphics[Flatten[lineOptions]~Join~{Line[pts]}],FilterRules[{opts},Options@Graphics],PlotRange->{{x0,x1},All},Ticks->{tks,Vitali[tks]},FrameTicks->{tks,verttks,Map[{#,""}&,tks],Map[{#,""}&,verttks]},
FilterRules[Options[VitaliPlot],Options@Graphics]]];

Peano's curve

SetAttributes[myPeanoCurve,Listable];
SetAttributes[px,Listable];
SetAttributes[py,Listable];

End points:

px[0]=0;
py[0]=0;
px[1]=1;
py[1]=1;

Value in the middle third:

px[t_/;1/3≤t≤2/3]=1/2;
py[t_/;1/3≤t≤2/3]=2-3t;

Self-similarity for rationals with denominator divisible by 3:

px[t_Rational/;And[0≤t<1/3,IntegerQ[Denominator[t]/3]]]:=px[t]=py[3t]/2;

py[t_Rational/;And[0≤t<1/3,IntegerQ[Denominator[t]/3]]]:=py[t]=px[3t];

px[t_Rational/;And[2/3≤t<1,IntegerQ[Denominator[t]/3]]]:=px[t]=(1+py[Mod[3t,1]])/2;

py[t_Rational/;And[2/3≤t<1,IntegerQ[Denominator[t]/3]]]:=py[t]=px[Mod[3t,1]];

Self-similarity for rationals with denominator not divisible by 3:

px[t_Rational/;And[Or[0≤t<1/3,2/3<t<1],!IntegerQ[Denominator[t]/3]]]:=px[t]=Module[{rl,f,g,itn},
rl=
{f[y_]:>Which[
0≤y<1/3,g[3y]/2,
1/3≤y≤2/3,1/2,
2/3≤y<1,(1+g[Mod[3y,1]])/2],
g[y_]:>Which[0≤y<1/3,f[3y],
1/3≤y≤2/3,2-3y,
2/3≤y<1,f[Mod[3y,1]]]};
itn=If[EvenQ[#],#,2#]&[EulerPhi[Denominator[t]]];
f[t]/.Solve[f[t]==Nest[Expand[#/.rl]&,f[t],itn],f[t]][[1]]];

py[t_Rational/;And[Or[0≤t<1/3,2/3<t<1],!IntegerQ[Denominator[t]/3]]]:=py[t]=Module[{rl,f,g,itn},
rl=
{f[y_]:>Which[
0≤y<1/3,g[3y]/2,
1/3≤y≤2/3,1/2,
2/3≤y<1,(1+g[Mod[3y,1]])/2],
g[y_]:>Which[0≤y<1/3,f[3y],
1/3≤y≤2/3,2-3y,
2/3≤y<1,f[Mod[3y,1]]]};
itn=If[EvenQ[#],#,2#]&[EulerPhi[Denominator[t]]];
g[t]/.Solve[g[t]==Nest[Expand[#/.rl]&,g[t],itn],g[t]][[1]]];

For floating-point or exact non-rational real numbers:

px[t_?NumericQ/;And[Head[t]=!=Rational,0≤t≤1]]:=Module[{v,t0,t1},
v=Vitali[t];
{t0,t1}={InfInverseVitali[v],SupInverseVitali[v]};
Interpolation[{{t0,px[t0]},{t1,px[t1]}},InterpolationOrder->1][t]];

py[t_?NumericQ/;And[Head[t]=!=Rational,0≤t≤1]]:=Module[{v,t0,t1},
v=Vitali[t];
{t0,t1}={InfInverseVitali[v],SupInverseVitali[v]};
Interpolation[{{t0,py[t0]},{t1,py[t1]}},InterpolationOrder->1][t]];

The function to export:

myPeanoCurve[t_/;0≤t≤1]:={px[t],py[t]};

Inverse Peano

Inverse Peano Curve for couples of rationals, of which at least one has even denominator:

InversePeano[x_Rational/;EvenQ[Denominator[x]],y_?NumericQ]:=(2Floor[2x]+InversePeano[y,Mod[2x,1]])/3;

InversePeano[x_?NumericQ/;x<1,y_Rational/;EvenQ[Denominator[y]]]:=(6Floor[2x]+2Floor[2y]+InversePeano[Mod[2x,1],Mod[2y,1]])/9;

Inverse Peano Curve for couples of rationals, both with odd denominator:

InversePeano[x_?NumericQ/;Or[MemberQ[{0,1},x],And[Head[x]===Rational,OddQ[Denominator[x]]]],y_?NumericQ/;Or[MemberQ[{0,1},y],And[Head[y]===Rational,OddQ[Denominator[y]]]]]:=InversePeano[x,y]=Module[{rl,h,itn},
rl=h[a_,b_]:>(6If[a<1,Floor[2a],1]+2If[b<1,Floor[2b],1]+h[If[a<1,Mod[2a,1],1],If[b<1,Mod[2b,1],1]])/9;
itn=EulerPhi[Denominator[x]]EulerPhi[Denominator[y]];
h[x,y]/.Solve[h[x,y]==Nest[Expand[#/.rl]&,h[x,y],itn],h[x,y]][[1]]];

Inverse Peano Curve for floating-point real numbers:

InversePeano[x_/;0≤x≤1,y_/;0≤y≤1]/;Min[N@Precision[x],N@Precision[y]]<Infinity:=Module[{pr,a,b},
pr=Floor[Min[Map[Log[2,10]N@Precision[#]-Ceiling[Log[2,#]]&,{x,y}]]];{a,b}=Map[RealDigits[#,2,pr,-1][[1]]&,{x,y}];
{a,b}=Map[Take[#,pr]&,{a,b}];
FromDigits[{2Flatten[Transpose[{a,b}]],0},3]+(1.`1)/3^(2pr-2)];

Numeric Inverse Peano Cuve for couples of non-rational infinite-precision numbers:

N[InversePeano[x_?NumericQ/;And[0≤x≤1,Precision[x]==Infinity],y_?NumericQ/;And[0≤y≤1,Precision[y]==Infinity]/;Or[Head[x]=!=Rational,Head[y]=!=Rational]],n_:16]:=InversePeano[N[x,Ceiling[n*Log[3,2]]],N[y,Ceiling[n*Log[3,2]]]];

Peano curve plots and movie

Plotting separately the two components of Peano's Curve:

Options[PeanoCurveSeparatePlots]={PlotPoints->2^7,PlotStyle->{},Axes->True,AspectRatio->Automatic};

PeanoCurveSeparatePlots[k_/;MemberQ[{1,2},k],x0_,x1_,opts___]/;x0<x1:=Module[{segmentNumber,lineOptions,n,x2,x3,pts,tks,verttks},
{segmentNumber,lineOptions}={PlotPoints,PlotStyle}/.{opts}/.Options[PeanoCurveSeparatePlots];n=Max[1,Ceiling[Log[2,segmentNumber]-Log[3,x1-x0]]];
{x2,x3}=3^(-n){Floor[x0*3^n],Ceiling[x1*3^n]};pts=Transpose[{#,({px,py}[[k]])[#]}&[Union[Join[CantorTicks[x2,x3,n],{x2,x3}]]]];
tks=CantorTicks[x0,x1];
verttks=({px,py}[[k]])[CantorTicks[x0,x1,4]];
Show[Graphics[Flatten[lineOptions]~Join~{Line[pts]},FilterRules[{opts},Options@Graphics],PlotRange->{{x0,x1},All},Ticks->{tks,verttks},FrameTicks->{tks,verttks,Map[{#,""}&,tks],Map[{#,""}&,verttks]},
FilterRules[Options[PeanoCurveSeparatePlots],Options@Graphics]]]];

PeanoCurveSeparatePlots[{1,2},x0_,x1_,opts___]/;x0<x1:=Show[GraphicsArray[Table[{PeanoCurveSeparatePlots[k,x0,x1,opts]},{k,1,2}],
FilterRules[Options[PeanoCurveSeparatePlots],OptionsGraphics]]];

Plotting the trajectory of the n-th approximant to Peano's Curve:

PeanoCurvePlot[n_Integer?Positive,opts___Rule]:=Module[{pts=myPeanoCurve[leftAndRightExtrema[n]]},If[(RoundCorners/.{opts}/.Options[PeanoCurvePlot])===True,
peanoRounded[pts,opts],
peanoAngular[pts,opts]]];

A couple of utilities:

vectorNorm[x_List]:=Sqrt[x.x];

normalizeVector[x_List]:=x/vectorNorm[x];

The plot in case no corner-rounding is asked for:

peanoAngular[{pts__},opts___Rule]:=Module[{smallestSegment,hdl,pltsty,makeArw,segmts},
smallestSegment=vectorNorm[{pts}[[3]]-{pts}[[2]]];
{hdl,pltsty}={HeadLength,{PlotStyle}}/.{opts}/.Options[PeanoCurvePlot];
If[hdl===Automatic,hdl=Min[smallestSegment/6,.08]];
makeArw[{pt1_List,pt2_List}]:={Arrowheads[hdl],Arrow[{pt1,(pt1+pt2)/2}],Line[{(pt1+pt2)/2,pt2}]};
segmts=Partition[{pts},2,1];
Show[Graphics[Join[{PointSize[Min[.05,hdl/2]]},Flatten[pltsty],{Point[{0,0}],Point[{1,1}]},Flatten[Map[makeArw,segmts]]]],FilterRules[{opts},Options@Graphics]]];

The plot in case corner-rounding is asked for:

peanoRounded[{pts__},opts___Rule]:=Module[{smallestSegment,hdl,rndamnt,pltsty},
smallestSegment=Sqrt[#.#]&[{pts}[[3]]-{pts}[[2]]];
{hdl,rndamnt,pltsty}={HeadLength,RoundingAmount,{PlotStyle}}/.{opts}/.Options[PeanoCurvePlot];
If[hdl===Automatic,hdl=Min[smallestSegment/5,.08]];
If[rndamnt===Automatic,rndamnt=Min[.05,smallestSegment/8]];
Graphics[Join[{PointSize[Min[.05,hdl/2]]},Flatten[pltsty],{Point[{0,0}],Point[{1,1}],Arrowheads[hdl],Arrow[{{pts}[[1]],({pts}[[1]]+{pts}[[2]])/2}],Line[{({pts}[[-2]]+{pts}[[-1]])/2,{pts}[[-1]]}]},Flatten[Map[roundSingleCorner[#,rndamnt,HeadLength->hdl,opts]&,Partition[{pts},3,1]]]],FilterRules[{opts},Graphics]]];

The corner-rounding routine:

roundSingleCorner[{x_List,y_List,z_List},ϵ_?Positive,opts___]:=Module[{a,b,c,θ,r,d,angls},
a=normalizeVector[x-y]//N;
b=normalizeVector[z-y]//N;
c=normalizeVector[a+b]//N;
θ=ArcCos[c.b]//N;
r=ϵ*Tan[θ];
d=ϵ/Cos[θ];
angls=Apply[ArcTan,{(c.a)a-c,(c.b)b-c},1]//N;
angls=If[Abs[Subtract@@angls]<Pi,Sort[angls,Less],Sort[angls,Greater]+{0,2Pi}];
{Circle[y+d*c,r,angls],Line[{(x+y)/2,y+ϵ*a}],Arrow[{y+ϵ*b,(y+z)/2}]}];

A movie of the point moving along Peano's curve:

PeanoCurveMovie[n_Integer?Positive,opts___Rule]:=Module[{sp,tms,bckgnd,pts,trpts,smallestSegment,sty,f,g,frm},
{sp,fpsi}={ShowPath,FramesPerSmallestInterval}/.{opts}/.Options[PeanoCurveMovie];
Which[fpsi===Automatic,fpsi=Max[1,5-n],
fpsi>0,fpsi=Ceiling[fpsi],
True,Max[1,5-n]];
tms=leftAndRightExtrema[n];
bckgnd=If[sp===True,PeanoCurvePlot[n,DisplayFunction->Identity],Graphics[{}]];
pts=myPeanoCurve[tms];
trpts=Transpose[pts];
smallestSegment=Sqrt[#.#]&[pts[[3]]-pts[[2]]];
sty=PlotStyle/.{opts}/.{PlotStyle->{PointSize[Min[.05,smallestSegment/3]]}};
f=Interpolation[N[Transpose[{tms,trpts[[1]]}]],InterpolationOrder->1];
g=Interpolation[N[Transpose[{tms,trpts[[2]]}]],InterpolationOrder->1];
frm[t_]:=Show[{bckgnd,Graphics[Join[Flatten[sty],{Point[{f[t],g[t]}]}]]},FilterRules[{opts},Options@Graphics],AspectRatio->Automatic,PlotRange->{{-.1,1.1},{-.1,1.1}},Frame->True,FrameTicks->False];
Animate[frm[Global`t],{Global`t,0,1,1/(fpsi*3^(n+1))}]];

Calculation of Van der Waerden's function

SetAttributes[VanDerWaerden,Listable];

Sawtooth function:

sawtooth[x_?NumericQ]:=Abs[x-Round[x]];

VdW function vanishes on the integers:

VanDerWaerden[_Integer]=0;

VdW function is periodic with period 1:

VanDerWaerden[x_?NumericQ/;(Element[x,Reals]&&Floor[x]=!=0)]:=VanDerWaerden[Mod[x,1]]

Self-similarity recursive formula for rationals with even denominator:

VanDerWaerden[x_Rational/;EvenQ[Denominator[x]]]:=sawtooth[x]+VanDerWaerden[2x]/2;

Self-similarity formula for other rationals:

VanDerWaerden[x_Rational/;OddQ[Denominator[x]]]:=VanDerWaerden[x]=Module[{g,rl},
g[y_?NumericQ/;Floor[y]=!=0]:=g[Mod[y,1]];
rl=g[y_]:>sawtooth[y]+g[2y]/2;
g[x]/.Solve[g[x]==Nest[Expand[#/.rl]&,g[x],EulerPhi[Denominator[x]]],g[x]][[1]]];

Arbitrary precision floating-point VdW function for real numbers with infinite precision. This is based on the Hölder inequality |f(x)-f(y)|≤(Sqrt[8|x-y|]), and the estimate k x≤f(x)≤(k+2/3)x for 0<x≤1/2, with k=Ceiling[-Log[2,x]]-1, that hold for VdW function.

N[VanDerWaerden[x_/;And[Element[x,Reals],Head[x]=!=Rational,Precision[x]===Infinity]],n_:16]:=Module[{x0,k,y0,p,bin},
x0=sawtooth[x];
k=Ceiling[-Log[2,x0]]-1;
y0=k*x0;
p=Ceiling[Log[2,8*10^(2n)/y0^2]];
bin=FromDigits[RealDigits[x0,2,p,-1],2];
VanDerWaerden[bin]+(1.`1)Sqrt[8/2^p]];

Floating-point VdW function for floating-point real numbers. This is based on the inequality (Subscript[r, 0])≤f(x)≤(Subscript[r, 1])+(2/3)((Subscript[x, 1])-(Subscript[x, 0])), where (Subscript[r, 0]):=Min[f((Subscript[x, 1])),f((Subscript[x, 0]))],  (Subscript[r, 0]):=Min[f((Subscript[x, 1])),f((Subscript[x, 0]))], (Subscript[r, 1]):=Max[f((Subscript[x, 0])),f((Subscript[x, 1]))], (Subscript[x, 0])=(k/(2^n)), (Subscript[x, 1])=((k+1)/(2^n)), (Subscript[x, 0])≤x≤(Subscript[x, 1]).

VanDerWaerden[x_Real]:=Module[{dg,x0,ϵ,x1,f0,f1},
dg=RealDigits[Mod[x,1],2];
x0=FromDigits[dg,2];
ϵ=FromDigits[dg/.{{a___Integer,b_},c_}:>{Join[Map[0&,{a}],{1}],c},2];
x1=x0+ϵ;
{f0,f1}=VanDerWaerden[{x0,x1}];
Min[f0,f1]+Max[Abs[f1-f0],2(x1-x0)/3](1.`1)];

Plotting Van der Waerden's function

Options[VanDerWaerdenPlot]={PlotPoints->2^9,PlotStyle->{},Axes->True,AspectRatio->Automatic};

VanDerWaerdenPlot[x0_,x1_,opts___Rule]/;x0<x1:=Module[{segmentNumber,lineOptions,otherOptions,n,x2,x3},
{segmentNumber,lineOptions}={PlotPoints,PlotStyle}/.{opts}/.Options[VanDerWaerdenPlot];n=Floor[Log[2,segmentNumber/(x1-x0)]];
n0=Max[0,Floor[Log[2,8/(x1-x0)]]];
{x2,x3}=2^(-n){Floor[x0*2^n],Ceiling[x1*2^n]};
tks=Select[Range[Floor[x0],Ceiling[x1],1/2^n0],x2≤#≤x3&];
Show[Graphics[Flatten[lineOptions]~Join~{Line[Table[{x,VanDerWaerden[x]},{x,x2,x3,2^(-n)}]]}],FilterRules[{opts},Graphics],PlotRange->{{x0,x1},All},Ticks->{tks,Automatic},FrameTicks->{tks,Automatic,Map[{#,""}&,tks],Automatic},
FilterRules[Options[VanDerWaerdenPlot],Options@Graphics]]];

Tail end

End[];

EndPackage[];

Examples

Most calculations are recursive. If you work with a fraction with high denominator, or want very high precision with floating-point result, you may need to increase the $RecursionLimit up from its default of 256: for example

$RecursionLimit=500;

Cantor Set

Start from the unit interval [0,1], and remove the middle third [1/3,2/3]. Then from the remaining two intervals remove their middle thirds, and so on: when repeated infinitely many times what you are left with is the (ternary) Cantor Set. With this package you can decide whether an exact real number belongs to the Cantor Set. The calculation uses recursively the self-similarity formula

x∈CantorSet⟺{({{3x∈CantorSet, if 0≤x≤1/3}, {False, if 1/3<x<2/3}, {3x-2∈CantorSet, if 2/3≤x≤1}})

0∈CantorSet

(1/2)∈CantorSet

(1/4)∈CantorSet

(1/(Sqrt[2]))∈CantorSet

(132/177)∈CantorSet

Take the rational number whose expansion in base 3 is 0.02002002002...

FromDigits[{{{2,0,0}},0},3]

The digit 1 does not appear in the expansion. Hence the number should belong to the Cantor set. Let's check:

CantorSetQ[9/13]

You can also see the successive approximations to the Cantor Set:

base=CantorSetPlot[{0,1,2,3,4},PlotRange->{-.1,5.2}]

CantorSetPlot[{0,1,2,3,4},PlotRange->All,Tree->True]

Vitali's function ("Devil's staircase")

Vitali's function can be described on the unit interval this way: if x belongs to the Cantor Set, write its expansion in base 3, using the periodic 2222... alternative in case of ambiguity, turn all digits 2 into digits 1, and finally interpret the result as an expansion in base 2. When x is not in the Cantor Set, make a linear interpolation between the nearest points of the Cantor Set to the left and right of x. This function is interesting because it is continuous, it is differentiable almost everywhere, its derivative vanishes almost everywhere, but the function is not constant.
This package calculates exactly Vitali's function for (practically) all exact real numbers, and gives a floating-point value for floating-point real entries. For exact real numbers the computation uses recursively the self-similarity formula

f[x]={({{((f[3x])/2), if  0≤x≤1/3}, {1/2, if  1/3<x≤2/3}, {((1+f[3x-2])/2), if  2/3<x<1}})

Sample values of Vitali's function:

Vitali[0]

Vitali[1/2]

Vitali[1/4]

Vitali[1/Sqrt[11]]

Vitali[1/Pi]

As Vitali's function is onto, it admits right-inverses. InfInverseVitali and SupInverseVitali calculate the smallest and, respectively, the largest of the inverse images, which are always points in the Cantor Set. Evaluating InfInverseVitali[x] or InfInverseVitali[x] for rational x in the unit interval will produce all rationals in the Cantor Set. InfInverseVitali[Vitali[x]] and SupInverseVitali[Vitali[x]] are the elements in the Cantor Set closest to x respectively to the left and to the right, when x is in the unit interval.

InfInverseVitali[1/2]

SupInverseVitali[1/2]

InfInverseVitali[7/11]

CantorSetQ[InfInverseVitali[7/11]]

Here are the elements of the Cantor Set immediately to the left and right of 12/47:

{InfInverseVitali[Vitali[#]],SupInverseVitali[Vitali[#]]}&[12/47]

The same for a floating-point number:

{InfInverseVitali[Vitali[#]],SupInverseVitali[Vitali[#]]}&[.7975]

The graph of Vitali's function is sometimes called "Devil's staircase":

VitaliPlot[0,1]

An enlargement of a detail:

VitaliPlot[1/4,1/3,Frame->True,Axes->False,ImageSize->500]

Peano's curve

Peano's curve {x[t],y[t]} for t∈CantorSet is defined this way: take the expansion of t in base 3:

t=0,(Subscript[t, 1])(Subscript[t, 2])(Subscript[t, 3])(Subscript[t, 4])(Subscript[t, 5])(Subscript[t, 6])(Subscript[…, (base 3)])

In case there is ambiguity choose the periodic 2222… form rather than the 0000… form. All digits of t are 0 or 2. Now separate the digits in odd positions and the digits in the even positions, and divide each of them by 2. This way you get two sequences of 0-1 digits, that you interpret as numbers in base 2:

x[t]=0,((Subscript[t, 1])/2)((Subscript[t, 3])/2)((Subscript[t, 5])/2)(Subscript[…, (base 2)]),          y[t]=0,((Subscript[t, 2])/2)((Subscript[t, 4])/2)((Subscript[t, 6])/2)(Subscript[…, (base 2)])

When t∈[0,1] is outside the Cantor Set, take the linear interpolation between the values of Peano's curve at the two closest points in the Cantor Set to the left and to the right of t.

Peano's curve is famous because it is a continuous function from the interval [0,1] onto the unit square [0,1]×[0,1].

This package calculates the value of Peano's Curve for exact rational values of t using the self-similarity relations

x[t]={({{y[3t]/2, if 0≤x<1/3}, {2-3t, if 1/3≤t≤2/3}, {(1+y[3t-2])/2, if 2/3<x≤1}})        y[t]={({{x[3t], if 0≤x<1/3}, {1/2, if 1/3≤t≤2/3}, {x[3t-2], if 2/3<x≤1}})

Examples:

myPeanoCurve[1/2]

myPeanoCurve[1/3]

myPeanoCurve[1/5]

For exact non-rational values of t the result is usually exact, unless t is very close to the Cantor Set and the recursion limit is hit:

myPeanoCurve[1/Sqrt[2]]//Simplify

For floating-point values of t the result is half-floating-point and half exact:

myPeanoCurve[.2]

Since Peano's Curve is onto the unit square, it admits a right inverse. The package defines one such right inverse with values in the Cantor Set when the entry is a couple of rational numbers:

InversePeano[1/4,4/7]

myPeanoCurve[InversePeano[1/4,4/7]]

The calculation is based on the self-similarity relations

h[x,y]={({{(2(Subscript[x, 1])+h[y,x'])/3}, {(6(Subscript[x, 1])+2(Subscript[y, 1])+h[x',y'])/9}, {(6(Subscript[x, 1])+2(Subscript[y, 1])+h[x',y'])/9}})

where (Subscript[x, 1])={({{Floor[2x], if x<1}, {1, if x=1}}), (Subscript[y, 1])={({{Floor[2y], if y<1}, {1, if y=1}}), x'={({{Mod[2x,1], if x<1}, {1, if x=1}}), y'={({{Mod[2y,1], if y<1}, {1, if y=1}})

The (right) inverse of Peano's curve can be calculated for floating-point entries too:

InversePeano[.3,.9]

For exact non-rational entries you can calculate the inverse to arbitrary precision:

N[InversePeano[1/E,1/Pi],40]

N[myPeanoCurve[N[InversePeano[1/E,1/Pi],40]]-{1/E,1/Pi}]

You can plot the x-component of Peano's curve on the interval [0,1]:

PeanoCurveSeparatePlots[1,0,1]

Here is a plot of the y-component:

PeanoCurveSeparatePlots[2,0,1]

You can also plot in the unit square the trajectory of the n-th approximant to Peano's curve.  In the plots horizontal and the vertical segments do belong to the full Peano's curve, while the oblique ones are an artefact of the approximant.

PeanoCurvePlot[1]

PeanoCurvePlot[2]

PeanoCurvePlot[8]

In the plots above the corners are artificially rounded to avoid self-intersections, so that you can follow the path. Here is  what  would happen without the rounding:

PeanoCurvePlot[2,RoundCorners->False]

You can also make a movie of the point moving along (the n-th approximant to) Peano's Curve:

PeanoCurveMovie[3]

Van der Waerden's function

Bartel Leendert van der Waerden introduced in 1930 a specially simple example of a continuous function that is not differentiable anywhere. The original article is

B. L. van der Waerden. Ein einfaches Beispiel einer nichtdifferenzierbaren Funktion. Math. Z. 32 (1930).

The definition is

f[x_]:=(∑_(n=0)^∞) (Subscript[a, 0])[(2^n)x]/(2^n),    where (Subscript[a, 0])[x_]:=Abs[x-Round[x]]

This package calculates the exact value of f[x] when x is rational, using the following rules:

f[n]=0 when n is integer;

f[x]=f[Mod[x,1]] for all x (periodicity condition);

f[x]=(Subscript[a, 0])[x]+f[Mod[2x,1]]/2 for all x (self-similarity),

and the fact that the sequence {({{(Subscript[b, 0])=x}, {(Subscript[b, n+1])=Mod[(Subscript[b, n]),1]}}) is eventually cyclical when x is rational.

With the package you get the exact value of the Van Der Waerden function in any exact rational point:

VanDerWaerden[0]

VanDerWaerden[1/2]

VanDerWaerden[37/2^4]

VanDerWaerden[1/3]

VanDerWaerden[1/18]

On exact non-rational numbers the result is symbolic:

VanDerWaerden[E]

VanDerWaerden[Pi]

You can calculate VdW function on exact non-rational real numbers with arbitrary precision. This is based on Hölder's inequality |f(x)-f(y)|≤(Sqrt[8|x-y|]), and the estimate f(x)≥(Ceiling[-Log[2,(Subscript[x, 1])]]-1) (Subscript[x, 1]) with (Subscript[x, 1])=|x-Round[x]|, that hold for VdW function.

N[VanDerWaerden[E]]

N[VanDerWaerden[Pi],31]

VdW function can take floating-point input. The output is still floating-point, with as much precision as allowed by the input.  This is based on the inequality (Subscript[r, 0])≤f(x)≤(Subscript[r, 1])+(2/3)((Subscript[x, 1])-(Subscript[x, 0])), where (Subscript[r, 0]):=Min[f((Subscript[x, 1])),f((Subscript[x, 0]))],  (Subscript[r, 1]):=Max[f((Subscript[x, 0])),f((Subscript[x, 1]))], (Subscript[x, 0])=(k/(2^n)), (Subscript[x, 1])=((k+1)/(2^n)), k∈ℤ (Subscript[x, 0])≤x≤(Subscript[x, 1]).

VanDerWaerden[.71]

You can plot Van der Waerden's function on the unit interval this way:

VanDerWaerdenPlot[0,1]

This is made with a Line graphic primitive with specially calculated subdivision. Of course you could as well use Plot[VanDerWaerden[x],{x,0,1}],  but the adaptive procedure of Plot is not well-suited to a nonsmooth function like ours: in different regions you may get either too fine or too raw a subdivision, and have symmetry messed up.

An enlargement of a detail:

VanDerWaerdenPlot[.2,.3,Frame->True]