Visually Exploring Computational Pseudo-Randomness: Complexity, Fractals, and Chaos

Dainel Carvalho, Wolfram, International Business Development Executive Brazil

#WolframTechConf

Computer Engineer
MSc Electrical Engineering

Started using Mathematica at version 4 at undergrad college LAB on calculus classes a long time ago!
Computing math on PC, instead of hand calculation ever since.

ABSTRACT

Graphical exploration of the intricate interplay between computational randomness, pseudo-random algorithms, and complexity. View how pseudo-randomness emerges from deterministic computational processes, revealing fractal patterns and chaotic behaviors that challenge traditional notions of randomness. By examining the computational complexity of these systems, we highlight how simple rules can generate intricate, unpredictable outputs akin to natural phenomena. Through practical demonstrations in the Wolfram Language, showcase tools for modeling and analyzing pseudo-random systems, with their applications in cryptography, simulation, and complex systems research, AI and data science. Bridges theoretical insights from NKS with practical computational techniques, offering a fresh perspective on randomness and complexity in modern computing.

ImageSynthesize["Create an black and white kind of pixel image to illustrate the abstract and this notebook text"]

Comparing Real Mathematical Constants and Random Numbers

With both visual and numerical visualization, it seams that some deterministic math real functions can be a source of pseudo random numbers:

References:
https://demonstrations.wolfram.com/ComparingRealMathematicalConstantsAndRandomNumbers/

Histograms for the Digits of Four Mathematical Constants

As more we have digits of real math numbers, better it converges to a uniform distribution, having a kind of pseudo random behaviour:

Interactive visualization of random integer sequences, it tends towards a uniform distribution:

SeedRandom[92177,Method-> "Rule30CA"]
Manipulate[
random=RandomInteger[9,size];
GraphicsGrid[{{
Histogram[random,ColorFunction->"SolarColors",PlotLabel-> "histogram of\nrandom integers"],
SmoothHistogram[random,ColorFunction->"SolarColors",PlotLabel-> "smooth histogram\nof random integers"],
PieChart[Last/@#,ChartLabels->First/@#,ColorFunction->"SolarColors",PlotLabel-> "pie chart\nof random integers"]&[Sort[Tally[random]]]}}],
{size,20,100000,1}]

The digits of well know deterministic mathematical constants tends towards to a uniform distribution at the histogram as having more digits.

Reference:
 https://demonstrations.wolfram.com/HistogramsForTheDigitsOfFourMathematicalConstants/

Exploring White Noise

Controlling the distribution, we can create white noise kind of pseudo random numbers:

𝒫=WhiteNoiseProcess[UniformDistribution[{-1,1}]];
sequenceWhiteNoise=RandomFunction[𝒫,{0,2000}];

White noise as source of randomness has a uniform distribution too at the histogram, for a graphical analysis:

Grid[{ArrayPlot[Partition[Values[#],Floor[Sqrt[Length[Values[#]]]]],ColorFunction->"DarkRainbow", Frame->None,PlotLabel-> Row[{Keys[#]," random\nsequence 2D lattice"}]],
MatrixPlot[Log[Abs[Fourier[Partition[Values[#],Floor[Sqrt[Length[Values[#]]]]]]]],ColorFunction->"SolarColors", PlotLabel-> Row[{Keys[#]," matrix\nabs fourier"}]],
ListPlot[Abs[Fourier[Values[#]]],ColorFunction->"DarkRainbow",GridLines->Automatic, PlotLabel-> Row[{Keys[#],"\nabs fourier list"}]],
ListLogLogPlot[Abs[Fourier[Values[#]]],ColorFunction->"DarkRainbow",GridLines->Automatic,PlotLabel-> Row[{Keys[#],"\nlog log abs fourier"}]],
Histogram[Values[#],100,ColorFunction->"SolarColors", PlotLabel-> Row[{Keys[#]," histogram"}]]}&/@
{"White Noise"->(Last /@First[Normal[sequenceWhiteNoise]])} ]
({{, , , , }})

Visualization of a white noise pseudo random numbers:

ListPlot[sequenceWhiteNoise,Filling->Axis,ColorFunction->"SolarColors"]

Fit visualization of the white noise pseudo random numbers:

ListFitPlot[sequenceWhiteNoise,ColorFunction->"SolarColors"]

Function to visualize statistics properties of the list:

descriptiveStatistics[serie_]:=
Grid[Partition[Labeled[AngularGauge[Last[#],{-1,1},GaugeLabels->Last[#],ImageSize->Tiny],First[#]]&/@
({
"Standard Deviation"->StandardDeviation[#],
"Variance"->Variance[#],
"Covariance"->Covariance[#],
"Kurstosis"->Kurtosis[#],
"Mean"->Mean[#],
"Mode"->Commonest[#,3],
"Median"->Median[#],
"Max"->Max[#],
"Min"->Min[#]} & [serie]),UpTo[3]]];

Visualization of descriptive statistical properties of the pseudo random list, white noise generated:

descriptiveStatistics[sequenceWhiteNoise]

Exploring Noise for Pseudo Randomness

Exploring a set of different noises interactively:

Manipulate[Grid[{{ArrayPlot[Partition[Values[#],Floor[Sqrt[Length[Values[#]]]]],ColorFunction->"DarkRainbow", Frame->None,PlotLabel-> Row[{Keys[#]," random\nsequence 2D lattice"}]],
MatrixPlot[Log[Abs[Fourier[Partition[Values[#],Floor[Sqrt[Length[Values[#]]]]]]]],ColorFunction->"SolarColors", PlotLabel-> Row[{Keys[#]," matrix\nabs fourier"}]]},{
ListPlot[Abs[Fourier[Values[#]]],ColorFunction->"DarkRainbow",GridLines->Automatic, PlotLabel-> Row[{Keys[#]," abs fourier list"}]],
ListLogLogPlot[Abs[Fourier[Values[#]]],ColorFunction->"DarkRainbow",GridLines->Automatic,PlotLabel-> Row[{Keys[#]," log log abs fourier"}]]},
{Histogram[Values[#],100,ColorFunction->"SolarColors", PlotLabel-> Row[{Keys[#]," histogram"}]],SpanFromLeft} }]&[
(noisesType->ResourceFunction["ColoredNoise"][size,noisesType] )],
{size,100,3000,1},
{{noisesType,"White"},{"Violet","Purple","Blue","White","Pink","Flicker","Red","Brown"}},
SaveDefinitions->True]

Exploring Noise for Pseudo Randomness
Side by Side View

Visual analysis side by side of the pseudo real random numbers matrices:

noisesType={"Violet","Purple","Blue","White","Pink","Flicker","Red","Brown"};

{ArrayPlot[Partition[Values[#],Floor[Sqrt[Length[Values[#]]]]],ColorFunction->"DarkRainbow", Frame->None,PlotLabel-> Row[{Keys[#]," random\nsequence 2D lattice"}]],
MatrixPlot[Log[Abs[Fourier[Partition[Values[#],Floor[Sqrt[Length[Values[#]]]]]]]],ColorFunction->"SolarColors", PlotLabel-> Row[{Keys[#]," matrix\nabs fourier"}]],
ListPlot[Abs[Fourier[Values[#]]],ColorFunction->"DarkRainbow",GridLines->Automatic, PlotLabel-> Row[{Keys[#]," abs fourier list"}]],
ListLogLogPlot[Abs[Fourier[Values[#]]],ColorFunction->"DarkRainbow",GridLines->Automatic,PlotLabel-> Row[{Keys[#]," log log abs fourier"}]],
Histogram[Values[#],100,ColorFunction->"SolarColors", PlotLabel-> Row[{Keys[#]," histogram"}]]}&/@
(#->ResourceFunction["ColoredNoise"][2000,#]&/@ noisesType) // Grid

Randomness Tests on Binary Sequences

Creating a binary random  numbers sequence with ECA rule 30 central cell, starting with one central cell as initial condition:

ArrayPlot[CellularAutomaton[30,{{1},0},100]]

Getting the binary value of the ECA 30 rule, 10K steps:

eca30central=#[[All,Ceiling[Length[#[[1]]]/2]]]&[CellularAutomaton[30,{{1},0},10000]];

Notes that it starts with one, the initial condition, 1 at the central cell.

Take[eca30central,10]

View as a binary image:

ArrayPlot[Partition[eca30central,Floor[Sqrt[Length[eca30central]]]],ColorFunction->"SolarColors",Frame->None]

Very close counts of zeroes and ones:

BarChart[Sort[Tally[eca30central]],ColorFunction->"SolarColors"]

Computing the randomness tests with the same binary random sequence:

Grid[Table[
{seq=ResourceFunction[randomTest][eca30central,"PValue"],
Labeled[AngularGauge[N@seq,{0,1},GaugeLabels->N@seq],randomTest<>"\nPValue"],
seqt=ResourceFunction[randomTest][eca30central,"TestStatistic"],
Labeled[AngularGauge[N@seqt,{0,1},GaugeLabels->N@seqt],randomTest<>"\nTestStatistic"]},
{randomTest,{"SerialRandomnessTest","ArcsineLawRandomnessTest",(* "RunCountRandomnessTest", *) "ChiSquareRandomnessTest"}}],Frame->All]

Randomness Tests on Real Sequences
Interactive

Run Count Randomness Test, pseudo real random matrix:

Manipulate[
SeedRandom[921779,Method->"Rule30CA"];
randomReal= RandomReal[{0,1},size];

Grid[
{{ArrayPlot[Partition[randomReal,Floor[Sqrt[Length[randomReal]]]],ColorFunction->"DarkRainbow",Frame->None],seq=ResourceFunction["RunCountRandomnessTest"][randomReal,"PValue"],
Labeled[AngularGauge[N@seq,{0,1},GaugeLabels->N@seq],"RunCountRandomnessTest\nPValue"]},

{SpanFromAbove,seqt=ResourceFunction["RunCountRandomnessTest"][randomReal,"TestStatistic"],
Labeled[AngularGauge[N@seqt,{0,1},GaugeLabels->N@seqt],"RunCountRandomnessTest\nTestStatistic"]}},
Frame->All],
{size,10,1000,1},SaveDefinitions->True]

P-Values higher, tending to one are good for pseudo random usage
P-Values lower, tending to zero or negative value is BAD (inappropriate) for usage as pseudo random source

More Experiments Measuring Pseudo Randomness with Noise Sequences

Another interactive example, with entropy and P-Value for standard randomness test, for measurement:

rangeAdjustment[list_]:=Module[{minVal,maxVal},
    minVal=Min[list];
    maxVal=Max[list];
    (list-minVal)/(maxVal-minVal)
]

rangeAdjustment[Range[-10,10]] // N

Manipulate[RandomSeed[921779];
sequenceNoise=ResourceFunction["ColoredNoise"][{size,size},noisesType];
Grid[{{plot2D[sequenceNoise,ColorFunction->"DarkRainbow",Frame->None,ImageSize->Medium],
Labeled[AngularGauge[N@#,{0,1},GaugeLabels->N@#],"RunCountRandomnessTest\nPValue"]&[
ResourceFunction["RunCountRandomnessTest"][rangeAdjustment[Flatten[sequenceNoise]],"PValue"]],
Labeled[AngularGauge[N@#,{0,1},GaugeLabels->N@#],"RunCountRandomnessTest\nTestStatistic"]&[
ResourceFunction["RunCountRandomnessTest"][rangeAdjustment[Flatten[sequenceNoise]],"TestStatistic"]],
Labeled[ThermometerGauge[N@#,{0,14},GaugeLabels->N@#],"Entropy"]&[
N[Entropy[Flatten[sequenceNoise]]]]}}],
{size,100,1000,1},
{{plot2D,ArrayPlot},{ArrayPlot, MatrixPlot}},
{{noisesType,"White"},{"Violet","Purple","Blue","White","Pink","Flicker","Red","Brown"}},
SaveDefinitions->True]

For a simple control test, this BAD repetitive random number sequence who’d have a very low P-Value:

Manipulate[RandomSeed[921779];
sequenceNoise=Table[#,size]&[RandomReal[9,size]];
Grid[{{plot2D[sequenceNoise,ColorFunction->"DarkRainbow",Frame->None,ImageSize->Medium],
Labeled[AngularGauge[N@#,{0,1},GaugeLabels->N@#],"RunCountRandomnessTest\nPValue"]&[
ResourceFunction["RunCountRandomnessTest"][rangeAdjustment[Flatten[sequenceNoise]],"PValue"]],
Labeled[AngularGauge[N@#,{0,1},GaugeLabels->N@#],"RunCountRandomnessTest\nTestStatistic"]&[
ResourceFunction["RunCountRandomnessTest"][rangeAdjustment[Flatten[sequenceNoise]],"TestStatistic"]],
Labeled[ThermometerGauge[N@#,{0,14},GaugeLabels->N@#],"Entropy"]&[
N[Entropy[Flatten[sequenceNoise]]]]}}],
{size,100,1000,1},
{{plot2D,ArrayPlot},{ArrayPlot, MatrixPlot}},
SaveDefinitions->True]

Fractals as Sources for Pseudo Random Sequences
Fractal Mandelbrot Set

Just for fun!! :-)

At first look it seams to be a good way to go, but it can get really risk quickly.
Navigating in the fractal space we can converge to just one color, or a limited number of colors and patterns, getting in a very low entropy set of values.

There are both sophisticated randomness and repetitive patterns at fractals complexity:

MandelbrotSetPlot[{-0.65+0.47 I,-0.4+0.72 I},PlotLegends->Automatic,ColorFunction->Hue]

We can find intricate sophisticated patterns, and boring patterns too navigating on the complex and real planes:

Manipulate[c=-real+complex I;
GraphicsGrid[{{MandelbrotSetPlot[Epilog->{Red,PointSize[Medium],Point[{Re[c],Im[c]}]},ColorFunction->Hue],JuliaSetPlot[c,ColorFunction->Hue]}}],
{real,0,1,.1},
{complex,0,1,.1},
SaveDefinitions->True]

To measure, lets try to get a standard initial condition 2D image data RGB:

ms=MandelbrotSetPlot[]
ArrayPlot[msData=First[First[ms]]]

Check image size:

Dimensions[msData]

View a histogram of all values that compose the 2D image:

Histogram[Flatten[msData],ColorFunction->"SolarColors"]

And the smooth histogram:

SmoothHistogram[Flatten[msData],ColorFunction->"SolarColors"]

A sample of frequence view with abs of Fourier:

ListPlot[Abs[Fourier[Flatten[RandomSample[Flatten[msData],3000]]]],ColorFunction->"SolarColors"]

Measure the randomness of the fractal 2D matrix data, it is to low for pseudo random usage:

ResourceFunction["RunCountRandomnessTest"][Flatten[msData],"PValue"]

ResourceFunction["ChiSquareRandomnessTest"][Flatten[msData],"PValue"]

Fractals as Sources for Pseudo Random Sequences
Fractal Julia Set

There are both sophisticated randomness and repetitive patterns at fractals complexity:

JuliaSetPlot[0.365-0.37 I,PlotLegends->Automatic,ColorFunction->Hue]

We can find intricate sophisticated patterns, and boring patterns too navigating interactively on the complex and real planes:

Manipulate[c=real-complex I;
GraphicsGrid[{{JuliaSetPlot[-1,Epilog->{Red,PointSize[Medium],Point[{Re[c],Im[c]}]},ColorFunction->Hue],JuliaSetPlot[c,ColorFunction->Hue]}}],
{real,0,1,.1},
{complex,0,1,.1},
SaveDefinitions->True]

Get the image data from Julia Set Fractal, the 2d matrix RGB values:

js=JuliaSetPlot[-1]
ArrayPlot[jsData=First[First[js]]]

Check image size:

Dimensions[jsData]

View a histogram of all values that compose the 2D image:

Histogram[Flatten[jsData],ColorFunction->"SolarColors"]

And the smooth histogram:

SmoothHistogram[Flatten[jsData],ColorFunction->"SolarColors"]

A sample of frequence view with abs of Fourier:

ListPlot[Abs[Fourier[RandomSample[Flatten[jsData],5000]]],ColorFunction->"SolarColors"]

Measure the randomness of the fractal 2D matrix data, it is to low for pseudo random usage:

ResourceFunction["RunCountRandomnessTest"][Flatten[jsData],"PValue"]

ResourceFunction["ChiSquareRandomnessTest"][Flatten[jsData],"PValue"]

As fractals can go to a space with simple patterns and converges to some or one color with lower entropy, to use Fractals as source of pseudo random, some strong boundaries should be set to minimize risks.

Elementary Cellular Automata Patterns as Source of Pseudo Randomness

Exploring statistical tests to assess the randomness of binary sequences for Cellular Automata patterns on a lattice. The focus is on understanding and applying different methods such as the CUSUM Max Randomness Test, and Run Length Randomness Test. Each of these tests provides unique insights into the properties of randomness in a sequence of zeros and ones, with applications ranging from cryptographic analysis to random number generation. We will delve into the mechanics of these tests, interpret their results, and discuss and experiment interactively their assumptions and limitations.

The CUSUM Max Randomness Test is a statistical test used to assess the randomness of a binary sequence, specifically sequences of zeros and ones. It computes a test statistic based on the maximum value achieved by a cumulative sums random walk derived from the sequence. The test then returns a two-tailed p-value to indicate the likelihood that the sequence is random. It is important to note that this test is valid only for sequences of length 100 or more.

We can measure the randomness level of zeroes and ones lists with "CUSUMMaxRandomnessTest":

There is a very low values for BAD randomness, and so, for non random binary lists:

Table[1,200]
ResourceFunction["CUSUMMaxRandomnessTest"][%,"PValue"]

We have good higher values for randomness, the P-Value tending towards 1:

RandomInteger[1,200]
ResourceFunction["CUSUMMaxRandomnessTest"][%,"PValue"]

Repetitive or low entropy data, has low values, here an binary identity matrix:

Flatten[IdentityMatrix[Floor[Sqrt[200]]]]
ResourceFunction["CUSUMMaxRandomnessTest"][%,"PValue"]

Repeating a short random least, looking for an average P-Value this way:

Flatten[Table[RandomInteger[1,10],20]]
ResourceFunction["CUSUMMaxRandomnessTest"][%,"PValue"]

The p-value is a statistical measure used to determine the significance of an observed result. In the context of randomness tests, like the CUSUM Max Randomness Test or the Run Length Randomness Test, the p-value helps assess whether a sequence is likely to be random or not.

Here’s how the p-value is interpreted: 
• Low p-value: This suggests that the observed pattern in the data is unlikely to have occurred by random chance. Therefore, the sequence may not be random.  
• High p-value: This indicates that the observed pattern in the data is consistent with what would be expected by random chance. Therefore, the sequence is likely random.

In summary, the p-value quantifies the probability of obtaining a result at least as extreme as the one observed, assuming that the null hypothesis (e.g., the sequence is random) is true. The lower the p-value, the stronger the evidence against the null hypothesis.

Setting the seed to random...

SeedRandom[92177,Method->"Rule30CA"]

Function to compute the compress ration of a matrix, list of lists:

compressionRatio[matrix_]:=Module[{originalSize,compressedSize},
  originalSize=ByteCount[matrix];
  compressedSize=ByteCount[Compress[matrix]];
  originalSize/compressedSize
]

Using 3 measurements for randomness level of a ECA patters at the lattice:

The entropy that came from physics to measure the order and chaos, and can be applied to information theory to measure randomness on data, interesting data is something in between order (repetitive boring patters) and chaos (random information)

Compression ratio, ordered data compression is high, more chaotic data has a low compression ratio

Cum Sum Max P-Value is low for repetitive patterns and higher rations for chaotic information

Putting it all together, we can with this demonstration, investigate visually and interactively the ECA patters and the level of randomness...

init=RandomInteger[1,1000];
Manipulate[ca=CellularAutomaton[rule,Take[init,size],size];
Row[{ArrayPlot[ca,Frame->None,ImageSize->300,PlotLabel->"eca rule "<> ToString[rule]],
VerticalGauge[N[Entropy[Flatten[ca]]],{0,1},GaugeLabels->"Value",PlotLabel->"entropy"],
VerticalGauge[compressionRatio[ca],{0,100},GaugeLabels->"Value",PlotLabel->"compression\nratio"],
Labeled[AngularGauge[ResourceFunction["CUSUMMaxRandomnessTest"][Flatten[ca],"PValue"],{0,1},GaugeLabels->"Value"],"Cum Sum\nP-Value",Top]}],
{{rule,100}, 0,255,1},
{rule,{30,75,110,151}},
{{size,100}, 10,1000,1},
SaveDefinitions->True]

Totalistic Cellular Automata Patterns as Source of Pseudo Randomness

The Run Length Randomness Test evaluates randomness by examining the lengths of consecutive identical elements (runs) within a sequence of zeros and ones, or real number between 0 and 1. By analyzing the frequency and distribution of these runs, the test determines whether the sequence deviates significantly from what would be expected in a random sequence. Both tests are useful in fields like cryptography and data analysis to ensure data sequences are genuinely random.

Test run length randomness "RunLengthRandomnessTest":

Random kinds of values:

ResourceFunction["RunLengthRandomnessTest"][(RandomInteger[2,1000]/.{0->.1,1->.5,2->.9}),"PValue"]

ResourceFunction["RunLengthRandomnessTest"][RandomReal[{0,1},1000],"PValue"]

Repetitive kinds of values:

ResourceFunction["RunLengthRandomnessTest"][Table[.2,1000],"PValue"]

ResourceFunction["RunLengthRandomnessTest"][(Flatten[IdentityMatrix[{100,100}]]/.{1-> .1,0->.9}),"PValue"]

Another view with 3 colors totalistic CA, r=1, k=3:

initTCA=RandomInteger[2,1000];
Manipulate[tca=CellularAutomaton[{rule,{3,1}},Take[initTCA,size],size];
Row[{ArrayPlot[tca,Frame->None,ImageSize->300,PlotLabel->"tca rule "<> ToString[rule],ColorRules->If[!color,{0->White,1->Gray,0->Black},{0->Yellow,1->Orange,2->Red}]],
VerticalGauge[N[Entropy[Flatten[tca]]],{0,1},GaugeLabels->"Value",PlotLabel->"entropy"],
VerticalGauge[compressionRatio[tca],{0,100},GaugeLabels->"Value",PlotLabel->"compression\nratio"],
Labeled[AngularGauge[ResourceFunction["RunLengthRandomnessTest"][(Flatten[tca]/.{0->.1,1->.5,2->.9}),"PValue"],{0,.1},GaugeLabels->"Value"],"run length\nP-Value",Top] }],
{{rule,100}, 0,2186,1},
{rule,{1765,1766,1769,1770}},
{{size,100}, 10,1000,1},
{{color,False},{False, True}},
SaveDefinitions->True]

Shallow[tca]

Test for the number of colors:

Tally[Sort[Flatten[tca]]]
PieChart[Last/@ %,ChartLabels->First/@ %]

Through the application of the CUSUM Max Randomness Test, and Run Length Randomness Test, this interactive notebook has provided a comprehensive overview of randomness testing for binary sequences matrices with NCA and TCA lattice patterns. These tests serve as powerful tools in identifying non-random patterns and ensuring the integrity of random sequences used in various fields. Understanding the strengths and limitations of each test allows for more informed decisions in selecting the appropriate method for different applications. The insights gained from this exploration emphasize the importance of rigorous randomness testing in maintaining the robustness and security of systems relying on random sequences for chaotic or more "boring" patterns.

Further Work

Randomness is an intriguing area for computation, has been studied by many computer scientists, engineer and mathematicians.
It goes from a theoretical field to real world practical applications.
We have a perceptions of randomness all over.
But it still hard to find randomness and measure it.
There is a lot of possibilities to study and advance in this field.

Those visualization can help to explore the fascinating interplay between computational pseudo-randomness, fractals, chaos, and complexity. Through a series of visualizations and statistical tests, we’ve delved into how deterministic processes can yield seemingly random patterns, challenging our traditional notions of randomness. From analyzing the statistical properties of real mathematical constants and noise sequences to examining the chaotic beauty of fractals and cellular automata, we have highlighted the subtleties and complexities inherent in generating and measuring randomness. The insights gained underscore the importance of rigorous analysis in applications like cryptography, data science, and complex systems research, where the robustness and unpredictability of random sequences are crucial. This exploration not only bridges theoretical concepts from New Kind of Science (NKS) with practical computational techniques but also opens avenues for further research into the nature of randomness and its applications in modern computing. [Gen AI - Wolfram Notebook Assistant - LLM]

Books and references links:

https://www.wolframscience.com/nks/p25--how-do-simple-programs-behave/

https://mathworld.wolfram.com/Entropy.html

https://mathworld.wolfram.com/ElementaryCellularAutomaton.html

https://reference.wolfram.com/language/ref/Compress.html

https://resources.wolframcloud.com/FunctionRepository/resources/CUSUMMaxRandomnessTest

https://resources.wolframcloud.com/FunctionRepository/resources/RunCountRandomnessTest/

https://brianmcfee.net/dstbook-site/content/intro.html

https://random.tastemaker.design/

Donald E. Knuth’s seminal series, "The Art of Computer Programming", Volume 2, Seminumerical Algorithms

https://www.wolfram.com/language/elementary-introduction/3rd-ed/

Márcio Rosa, Laboratórios Virtuais – Cálculo, Álgebra e Geometria, UNICAMP, 2020.