```---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
R.  3. Other examples of use  (more)

Click → here to see "Examples of use" this is an extension of.

From A to P there are developments of the parts that bear the same name.  Click the
letter to see them.
From S01 to S20 there are new parts.
Click here to see an index of this paper, here for one on its own page, here for
an index of subjects.

USE as CALCULATOR

# A-B You must first write to R (only once):

source("http://macosa.dima.unige.it/r.R")

# source executes the commands that are stored in the indicated external file. We can
# use source to load commands, exercises, … from files in my computer or on Internet.
# See here for an example. See s18 for another one. Here for an example of "worksheet".

# If you expect to be unable to use the network, you can save a copy of
# http://macosa.dima.unige.it/r.R  on your computer and then use "source" from the

# C  Other functions:

# 3       2         1        -1          0       2       -0.1

#  [ Although (-8)^(1/3) is -2, software calculates x^(1/3) only when x ≥ 0; so
#    For the calculation of derivatives, use sqrt(k^2) instead of abs(k)            ]

# INI_N(x) calculates the initial decimal place of the number x
INI_N(0.25); INI_N(37*1000); INI_N(-57/10000000)
#    2            3               5

# If I divide M by N, M%/%N is the integer quotient, M%%N is the remainder
750 %% 360; 750 %/% 360
# 30          2

factorial(4);  choose(10,3)             # factorial, binomial coeff.     (see combn)
#  24  (= 4*3*2*1)  120  (= 10/3*9/2*8/1)

factorial(170); factorial(171);
# 7.257416e+306   Inf
# 171! is too large. If you type .Machine you have that 1.797693e+308 is the maximum
# usable number.

# To get a help for other functions, type:   (see here)
?Math

# Other π (pi) another mathematical constant identified with a specific name, e, is the
# Napier's constant; the function exp(x) is equivalente to e^x
e; exp(1); more(e)                # see "more", down
# 2.718282  2.718282  2.71828182845905

# Besides round:

trunc(2.6);  floor(2.6)           # truncate (to integer) instead of round
#  2           2
trunc(-7.4); floor(-7.4)          # with differences on negative numbers (see)
#  -7         -8
Trunc(26.47,1); Floor(26.47,1)    # these also truncate to non-integers
#   26.4             26.4
Trunc(-26.47,-1); Floor(-26.47,-1)
#   -20             -30
print(20/3, 2); print(20/3, 15)   # round to 2, 15 digits; see "more", down.
#   6.7       6.66666666666667
signif(20/3, 2); signif(20/3, 7)  # round to 2, 7 digits (7 significant digits at most)
#   6.7       6.666667
# more visualizes 15 digits, last() is the value of the last calculation
(5/13+1/7)/17; fraction(last()); more((5/13+1/7)/17)
#   0.0310278      48/1547       0.0310277957336781
# [ if A is an integer of N digits more(A) displays A correctly even if it has 16
#   digits; if A has 17 digits more(A) displays all the digits but the last one can
#   be wrong:  a=12345678901234567; more(a)  gives 12345678901234568  Why? see here]

# Some uses of circular functions

180*degrees; 180*degrees/pi; sin(30*degrees); cos(30*degrees); tan(45*degrees)
# 3.141593       1        0.5            0.8660254      1
cos(30*degrees); fraction(cos(30*degrees)/sqrt(3))
# 0.8660254    1/2   (I understand that cos(30°) = √3/2)

# I can also use WolframAlpha by introducing 0.8660254 or cos(30°) or cos(30*degrees);
# I obtain √3/2  ( "degrees" is defined so:  degrees = pi/180 )

# asin, from [-1,1] to [-π/2,π/2], acos, from [-1,1] to [0,π], atan, from (-∞,∞) to
# (-π/2,π/2), are the inverse functions; for the output in degrees add "/degrees".

# D-E2 str compactly display the internal structure of an R object. Example:

str( divisors(123456) )
# num [1:28] 1 2 3 4 6 8 12 16 24 32 ...  (they are 28: 1 2 3 4 6 ...)
str(sin); str(choose)
# function (x)  function (n,k)    (a function of 1 variable, a function of 2 variables)

# How to print a table of data. Various methods.

z = c(1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6,7,8,9,0,1,2,3,4); length(z)
# 24   (the data are 24; I put then in 8 columns, separated by ", "
write(z, file="",sep=', ',ncolumns=8)
#  1, 2, 3, 4, 5, 6, 7, 8
#  9, 0, 1, 2, 3, 4, 5, 6
#  7, 8, 9, 0, 1, 2, 3, 4
write(z, file="",sep='\t',ncolumns=8)  # I separe the data by tabs
#  1       2       3       4       5       6       7       8
#  9       0       1       2       3       4       5       6
#  7       8       9       0       1       2       3       4
# [ I can "write" data to file; in these cases I must put file=""; if I want to write
#  to the file "XX.txt" in the working directory I use:  write(z, file="XX.txt", … )
#  otherwise I can specify the path:                   write(z, file="…/XX.txt", … ) ]

# cat, that concatenates data output (and with '\n' put them in new lines), produces:
cat('',z[1:12],'\n',z[13:24],'\n')
#  1 2 3 4 5 6 7 8 9 0 1 2
#  3 4 5 6 7 8 9 0 1 2 3 4
#       ( \n cab be used also with "type" and other writing commands )

# If the data have different precisions, the tabs may give rise to bad effects:
cat(20/3, 72, 15/4,'\n', 3, 1/9, 7,'\n')
# 6.666667 72 3.75
#  3 0.1111111 7

# I can make up width a tab and a rounding:
k = function(x) round(x,3)
cat('',k(20/3),'\t',k(72),'\t',k(15/4),'\n',k(3),'\t',k(1/9),'\t',k(7),'\n')
#  6.667   72      3.75
#  3       0.111   7
# or using \b to introduce a "backspace":
cat('',k(20/3),'\t\b',k(72),'\t\b\b',k(15/4),'\n',k(3),'\t\b',k(1/9),'\t\b\b',k(7),'\n')
#  6.667  72     3.75
#  3      0.111  7
# or in this way:
k = function(x) signif(x,3)
cat('',k(20/3),'\t',k(72),'\t',k(15/4),'\n',k(3),'\t',k(1/9),'\t',k(7),'\n')
#  6.67    72      3.75
#  3       0.111   7
# or opportunely by more tabs:
cat('',20/3,'\t',72,'\t\t',15/4,'\n',3,'\t\t',1/9,'\t',7,'\n')
#  6.666667        72              3.75
#  3               0.1111111       7

# cat is also convenient for printing tables. Let's see, for example, how I can print
# some binomial coefficients  (we use sep="" otherwise " " is used as a separator):

# [ for(var in seq) expr  repeats "expr" for all the values of "var" listed in "seq" ]

n = 5; for(k in 0:n) cat("C(", n, ",", k, ") = ", choose(n,k), sep="", '\n')
# C(5,0) = 1
# C(5,1) = 5
# C(5,2) = 10
# C(5,3) = 10
# C(5,4) = 5
# C(5,5) = 1
# See here for Pascal's triangle

# How to use cat to write (the beginning of) an irrational number ...

N=function(m) {cat("0."); for(i in 0:m) {for(j in 1:2^i) cat("0"); cat("1")}; cat("\n") }
N(0); N(1); N(2); N(5)
# 0.01
# 0.01001
# 0.0100100001
# 0.010010000100000000100000000000000001000000000000000000000000000000001

# CAT(X,N) print the sequence X of numbers or strings in rows of N columns

a=c(123,54321,356,123456,654321,13579,246810,"End"); CAT(a,4)

# 123     54321   356     123456
# 654321  13579   246810  End

x=NULL; N=4; for(i in 1:N) for(j in 1:N) x=c(x,i*j); CAT(x,N)

# 1       2       3       4
# 2       4       6       8
# 3       6       9       12
# 4       8       12      16

# Cat(x,N) produces a string that begins with the string x and ends with a sequence of
# " " so to get one of at least N characters [Cat("be",3) -> "be "; Cat("be",0) -> "be"].

m=17; n=8; p=0
cat(sep="",Cat("Mario Bianchi",m),Cat("Genova",n),Cat("via Napoli, 7",p), "\n",
Cat("Veronica Rossi",m),  Cat("Milano",n),Cat("via Roma, 21",p),  "\n",
Cat("Luisa Neri",m),      Cat("Roma",n),  Cat("via Genova, 13",p),"\n" )

#  Mario Bianchi    Genova  via Napoli, 7
#  Veronica Rossi   Milano  via Roma, 21
#  Luisa Neri       Roma    via Genova, 13
#( 1234567890123456712345678                  m=17; n=8)

# caT(x,N) produces a string that begins with a sequence of " " and ends with the
# string x so to get one of at least N characters.

m=14; n=9; p=17
cat(sep="",caT("Mario Bianchi",m),caT("Genova",n),caT("via Napoli, 7",p), "\n",
caT("Veronica Rossi",m),  caT("Milano",n),caT("via Roma, 21",p),  "\n",
caT("Luisa Neri",m),      caT("Roma",n),  caT("via Genova, 13",p),"\n" )

#   Mario Bianchi   Genova    via Napoli, 7
#  Veronica Rossi   Milano     via Roma, 21
#      Luisa Neri     Roma   via Genova, 13
#( 1234567890123412345678912345678901234567   m=14; n=9; p=17)

Other CALCULI

# F-G1 I can use an indexed variable (or vector), for example I can write XX[3]=8,
# only if I have already used XX. A simplwe way is to make this assignment: XX = NULL.
# Example:

for(i in 1:9) W[i] = i*2; W                    # If I have never used W:
#  Error in W[i] = i * 2 : object "W" not found
W = NULL; for(i in 1:9) W[i] = i*2; W
#  [1]  2  4  6  8 10 12 14 16 18

# A more general way is to define the length of the vector XX:
Z = vector( length=9) ; for(i in 1:9) Z[i] = i*2; Z
#  [1]  2  4  6  8 10 12 14 16 18

# prod (like sum) returns the product of all the values present in its arguments:

x = c(2, 3, 7, 100); prod(x)
#  4200

# reverse returns the arguments in reverse order
x = c(2, 3, 7, 100); y = reverse(x); y
# [1]   100   7   3   2

# combn(x,N) generate all combinations of the elements of x taken N at a time
x = c(2, 3, 7, 100); combn(x,2)
#      [,1] [,2] [,3] [,4] [,5] [,6]
# [1,]    2    2    2    3    3    7
# [2,]    3    7  100    7  100  100      They are 6 = choose(4,2) (see)
combn(x,3)
#      [,1] [,2] [,3] [,4]
# [1,]    2    2    2    3
# [2,]    3    3    7    7
# [3,]    7  100  100  100                They are 4 = choose(4,3)

# Interval computation.   +  -  *  /  ^:

# The area of the rectangle with one side between 30.5 and 31.5 and the other between
# 1.3 and 1.4 cm   [put as operand a certain number or the pair of extremes of the
# interval in which it is]:
approx( c(30.5, 31.5), c(1.3, 1.4), "*")
#   [1] min  39.65    [2] max  44.10
approx2( c(30.5, 31.5), c(1.3, 1.4), "*")
#   [1] center  41.875    [2] radius  2.225
# It is in [39.65 cm², 44.10 cm²], or:  41.9 ± 2.3 cm²

# The double of the first side:
approx( c(30.5, 31.5), 2, "*")
#    [1] min  61    [2] max  63     it is 62 ± 1 cm
# Idem for   +, -, /, ^.

# Another example:
x = 2.1 ± 0.1,  y = 2.9 ± 0.1,  1/z = 1/x + 1/y.  How much is z?
x = c(2.0, 2.2); y = c(2.8, 3.0)
rx = approx(1,x, "/"); rx           # 0.4545455  0.5000000
ry = approx(1,y, "/"); ry           # 0.3333333  0.3571429
rz = approx(rx,ry, "+"); rz         # 0.7878788  0.8571429
z = approx(1,rz, "/"); z
#   [1] min  1.166667   [2] max  1.269231
# or:
approx2(1,rz, "/")
#   [1] center  1.21794872   [2] radius  0.05128205
# or alternatively:  1.22 ± 0.06

# Other example (see here):

#     catheti: 22 ≤ a ≤ 23,  38 ≤ b ≤ 39  →  hypotenuse: 43.9 ≤ c ≤ 45.3

# In other situations (if there are other functions or a variable appears more than
# once) I can use EPS(r) that generates random values from -r to r. An example.

# A marble is launched in the air with initial speed v=4±0.1 m/s.  After t=0.6±0.06 s
# it reaches the height h = v*t-g*t^2/2, where g = 9.80±0.001 m/s^2. How tall is h?
t = 0.6+EPS(0.06); v = 4+EPS(0.1); g = 9.8+EPS(0.001)
h = function(t,v,g) v*t-g*t^2/2; m=min( h(t,v,g) ); M=max( h(t,v,g) ); m; M
#  0.4394657  0.7852536     h is between these value, or:
(m+M)/2; (M-m)/2
#  0.6123597  0.1728939     h = 0.612 +/- 0.173, or   h = 0.6 +/- 0.2 (m)
# (the solutions with more figures do not make sense!)

FIGURES

# H-I If I put ASP = k before PLANE (monometric plane) I can change the aspect ratio
# between the units of y-axis and x-axis. If I move the window ASP doesn't change.

BF=2.6; HF=2.6
ASP = 1/2; PLANE(-2,2, -4,4); circle(0,0, 2, "blue")
ASP = 3;   PLANE(-6,6, -2,2); circle(0,0, 2, "brown")

# You can save the current plot and to replay it (after you've opened a new graphics
# device). If I type Name = recordPlot() (Name is any name) the plot is saved; if I
# open a new grapics window later and type replayPlot(Name) the old object reappears
# and I can work it again [if you want delete the latest device you can use dev.off();
# dev.new() open a new device; the two commands can be useful to realize animations; I
# use dev.new(width=b,height=h) - or dev.new(width=BF,height=HF) - if I want a device
# like that obtained with BF=b; HF=h].

# With  ARC(xCenter,yCenter, radius, from,to, col)  I draw arcs of a circle:
HF=2.5; BF=2.5; PLANE(-4,4, -4,4); circle(0,0, 4, "violet")
ARC(-3,0, 1, 0,180, "blue"); ARC(1,0, 3, 180,360, "red")

# Another example (a right-angle triangle):
PLANE(-1,10, -1,10)
segm(9,1, 3,8, "red"); inclination(9,1, 3,8)
# [1] -49.39871    (inclination of red segment)
# I draw a 5 long segment form 3,8 towards the direction by 90° decreased
Direction(3,8, inclination(9,1, 3,8)-90, 5, "blue")
# I draw a segment from the arrival point to 9,1
segm(Directionx,Directiony, 9,1, "violet")
# I draw the arc (with center 3,8 and radius 2) that goes from side to side
ARC(3,8, 2, inclination(9,1, 3,8)-90, inclination(9,1, 3,8), "orange")
type(3,5,"90°")

# The third figure uses the command ARCd, which draws an arrow to give a direction
# to the arc.
PLANE(-1,10, -1,10)
ARCd(5,5, 4, 45,180, "blue"); ARCd(5,5, 4, 45,-180, "red")
Directio(5,5, 180,4, "seagreen"); Directio(5,5, 45,4, "seagreen")
CLEAN(4,6, 7,8); text(5,7.5, "135",font=4)
CLEAN(4,6, 2,3); text(5,2.5,"-225",font=4)

# circl and circl2, circl3p arc arcd draw thinner lines than circle, circle3p, ARC, ARCd.
# halfl, l2p and p2p draw thinner lines than halfline, line2p and perp2p.
# polylin and polyl draw thinner or much thinner lines than polyline; polyli makes use
# of dashed lines.

See here how to make this image of a building made with the Meccano.
# arcC(a,b,R, u,v,C) draws the circular sector (with center a,b and radius R) from u°
# to v° in color C
PLANEww(-5,5,-5,5); B="brown"
arcC(0,0,5, -30,50,"green"); arcC(0,0,5, 50,130,"blue")
arcC(0,0,5, 130,210,"yellow")
ARC(0,0,5, -30,210, B)
Direction(0,0,-30,5,B); Direction(0,0,50,5,B)
Direction(0,0,130,5,B); Direction(0,0,210,5,B)

# "Black" can be specified by 1. Other colors can be specified by giving an index.
# With Ncolor() I can quote the numbers associated with these colors:
Ncolor()
#1 black, 2 red, 3 green, 4 blue, 5 cyan, 6 violet, 7 yellow, 8 grey
# (then: 9 is black, 10 is red, ...)
HF=1; BF=4; PLANE(0,21,-1,1); for (c in 1:20) POINT(c,0, c)

# Alternatively, colors can be specified directly with a string of the form rgb(a,b,c)
# where each of numbers a,b,c has a value in the range 0 to 1 that states the quantity
# of red,green,blue: POINT(2,3, rgb(1,0,0) )  draw (2,3) red, POINT(4,5, rgb(1,0.5,0) )
# draw (4,5) orange, POINT(3,2, rgb(0.7,0.7,0.7) )  draw (3,2) grey.

# I can also use colors() to get more names than those obtainable with Colors().

# Other pictures here

# J-J6 If I type Quadrilateral(a,b,c,d)  width 4 directions (numbers ranging from 0
# to 360) instead of a,b,c,d  I obtain:

#     [ the first figure has been obtained with  Quadrilateral(0,80,195.67,286) ]
# You can easily build similar "problems".

# P() allows you to click on a point and read its coordinates. Example:

HF=3; BF=3; PLANE(0,10, 0,10)
P()        # I click (more or less) the point 5,9 and I obtain
# xP = 4.989773 ;  yP = 9.005114      or something like that
# See here.

# PPP() store the coordinates of the clicked point by the names xP and yP
#  (we have seen here how it can be usefully employed in the basic school)

# NW() prepares a new window of the same dimensions as the last one.

# polyS, polyS2 act like polyC but they emphasize the polygon by shading lines (± dense)
# with a given angle (rather than with a color); see the left figure.

BF=3;HF=3; PLANE(-1,5, -1,5)
x = c(-1,1,5,2); y = c(0,-1,2,4); polyS(x,y, 45, "red")
x1 = c(0,3,5);  y1 = c(5,-1,1); polyS2(x1,y1, -45, "seagreen")

# See here for the construction of the center and left figures, and for the use of the
# commands:  polySN(x,y, angle, N, col) for other shading lines,
# text(x,y,text, font=N) and the options to change alignment, size, color, font, string,
# rotation, family,
# dart2, arrow2, Dart2, Arrow2 that connect the top of the preceding arrow with the top
# of the new one.  With arroW2 and ArroW2 I have an intermediate thickness between
# "dart" and "arrow".# Some examples of writing (on a coloured background):
PLANE(0,6, 0,6)
x=c(1,6,6,1); y=c(1,1,6,6); polyC(x,y, "grey")
type(2,2,"A\nB"); typeC(3,2,"grey","C\nD")
typeC(4,2, 0,"E\nF"); typeR(2,4, 30,"G\nH")
typeRC(3,4, 30, 0,"I\nL"); text(4.5,5, font=4,srt=45,"M\nN")
text(5.3,2.5, font=4, cex=1.2, srt=180,"O\nP")
# If you want superimpose any white rectangle before writing see S18
# With type and text, I can use bquote(2*pi) and bquote(sqrt(3)) to obtain 2π and √3.
# For other formulas that I can obtain with bquote, see here. Examples:

# A lot of geometrical problems can be easily solved.
See here for an example.
# Direction(x,y, d,r, col) with the color 0 doesn't draw the "radius", but in Directionx
# and Directiony there are the coordinates of the arrival point. An example: (see here)

# Here are some other commands that are useful for drawing figures.
# I can shift with a translation, rotate and multiply the points by:
# MulTraRot(x,y, mx,my, tx,ty, ang) and RotTraMul(x,y, ang, tx,ty, mx,my).   The
# first multiply the coordinates by mx and my, then shift by tx,ty and finally rotates
# the point about (0,0) by ang°.  The second produces first a rotation about (0,0),
# then a translation and finally a multiplication of the coordinates. If I add [1] or
# [2] to the command I have the 1st or the 2nd coordinate. To operate on curves see S01.
# We have seen here symm(x,y, A,B, C,D) and symm2(x,y, A,B, d) to obtain the symmetric
# of x,y with respect to the line through (A,B) and (C,D) or through (A,B) with the
# direction of d degrees.
# RotP(x,y, ang, x0,y0) rotates (x,y) about (x0,y0) by ang°.
# RotPx(x,y, ang, x0,y0), RotPy(x,y, ang, x0,y0) give his two coordinates.
# MulRotP(x,y, m, ang, x0,y0) multiply by m the distance of (x,y) from (x0,y0) and
# rotates it about (x0,y0) by ang°.
# See here for an example, that generates the following immage.

# The immage on the right recalls some activities with the pentominoes described here.

# After drawing a line or an arrow (with segm, line, arrow, dash, ...), the end point
# coordinates are stored in xend and yend.  This can be useful for building various
# types of figures. An initial point can be traced with Dot: its coordinates are put
# in xend, yend. See here for an example, that produce the figure above, on the right.
# Dot2 traced a point large 3 pixel.

GRAPHS

# K Other ways of putting inscriptions along the axes

# On the left:
year = c(1880,1890,1900,1910,1920,1930,1940,1950,1960,1970,1980)
fruit = c(19,  21,  23,  28,  31,  27,  23,  37,  65,  88,  79)
BF=4; HF=2.5
Grid(1880,1980, 0,90)
polyline(year,fruit, "blue"); POINT(year,fruit, "brown")
underx("year"); underX(1880,1880); underX(1980,1980)
# underx(n) writes the name n in the center, underX(n, p) in the position p
undery("fruit"); aboveY(90,90); aboveY(0,0); aboveY(50,50)

# With AXES(a,b,"Col") I can trace two axes of color Col through (a,b). In the picture:
AXES(1880,0,"brown")
# AXES2 tracks them thicker.   underx2 writes in the line below.
underx2("kg of fruit consumed annually by an Italian")

# On the right:
# ...
UnderY("fruit",50); AboveY(90,90); AboveY(0,0); AboveY(50,50)
# or:  n=c(0,50,90); for(i in 1:3) AboveY(n[i],n[i])
# Other commands are  underY, aboveX, AboveX, UnderX  (in addition to  abovex, abovey).

# If in x and y there are the coordinates of some points, PlaneP(x,y) automatically
# chooses a "plane" that contains all of them (it is simliar to graphF for functions);
# PLANEP(x,y) chooses it monometric. For the sake of brevity we use the preceding data:
PlaneP(year,fruit); polyline(year,fruit, "blue")
PLANEP(year,fruit); polyline(year,fruit, "blue")

# range (see) gives min/max of some data, RANGE estimate min/max of a function between
# 2 input (if the function had a bump we could use maxmin - see):
range(fruit)
#19 88
f = function(x) (x+1/3)^2+1; RANGE(f, -5,10)
#"~ min F(min)  Max F(Max)   mfMf[1:4] = "
# -0.3333303   1.0000000  10.0000000 107.7777778
RANGE(f,-0.5,0)
#"~ min F(min)  Max F(Max)   mfMf[1:4] = "
# -0.3333333  1.0000000  0.0000000  1.1111111
mfMf[1]; mfMf[2]
# -0.3333333   1.0000000
RANGE(sin, -5,10)
#"~ min F(min)  Max F(Max)   mfMf[1:4] = "
# 4.71239 -1.00000 -4.71239  1.00000

Two of the various min/max __/

# RANGE0 is less precise but faster than RANGE
RANGE0(f, -0.5,10)
#"~ min F(min)   Max F(Max)   mfMf[1:4] = "
# -0.3330333   1.0000001  10.0000000 107.7777778
#
# Using these commands I can easily get any type of chart; it can be useful for teachers:
# two examples referring to typical activities of the first years of school (see):

# K1-K2 In addition to graphs of 1 variable functions, I can plot curves described as
# functions of 2 variables using the command CURVE(K, col). Examples:

PLANE(-5,5, -5,5)
f = function(x,y) y-(x^2-3);  CURVE(f, "red")       # parabola
g = function(x,y) x-(y^2-3);  CURVE(g, "blue")      # parabola
h = function(x,y) x^2+y^2-16; CURVE(h, "violet")    # circle
k = function(x,y) x*y-1;      CUR(k, "black")       # hyperbola
t = function(x,y) dist_taxi(x,y, -2.5,-3)-2
CURVE(t,"brown"); POINT(-2.5,-3, "brown")

# where dist_taxi(x1,y1,x2,y2) is the taxicab (or Manhattan) distance: |Δx|+|Δy|

# For other ways to describe curves, see here [the conics] and the S01 section.

# To solve eq. of type F(x,y) = 0 with respect to x (fixed y=k) or y (fixed x=k) see
# here

# CUR tracks thin-stroke curves; CURV uses an intermediate thickness.
# An example [see:
# here]:
#
# 2^(x*y) = x+y
# ( and y = x )

# Another:
f=function(x,y) x^y - y^x
BF=2.5; HF=2.5
PLANE(0,8, 0,8)
CURV(f, "brown")
# [ when x>0, y>0;  but
#   we have also, ie:
#   (-2)^(-4) = 1/2^4 =
#   0.0625 = (-4)^(-2)  ]
# (e,e) is a double point

# An egg:

A=c(-4,0); B=c(3,0)
F = function(x,y) 2*point_point(A[1],A[2],x,y)+3*point_point(B[1],B[2],x,y)-24
BF=2.5; HF=2; PLANE(-5,5, -4,4); CURVE(F, "brown")
# How to color the interior (see here):
for(i in 1:10) Diseq1(F,0, "pink")
CURVE(F, "brown"); POINT(A[1],A[2],"blue"); POINT(B[1],B[2],"blue")

# Likewise, graph1F and graph1 work in the same way as graphF and graph but with
# thin-stroke curves. graph2F and graph2 make them with an intermediate thickness.

Plane(-5,5, -5,15)
q = function(x) x^2-4; graph(q,-5,5, "blue")
s = function(x) x^2-2; graph1(s,-5,0, "red"); graph2(s,0,5, "seagreen")
#      f, g, h, k                       q, s

# Not all curves have a curved shape:
F = function(x,y) sqrt(x^2+y^2)-2*sqrt(x*y); PLANE(-3,3, -3,3); CURV(F,"seagreen")
G = function(x,y) x^2*y-x*y^2+x*y*sqrt(x*y); PLANE(-3,3, -3,3); CURV(G,"brown")

# NOTE. graph and graphF take enough time to plot the graphs, but they make them very
# precisely. See here to pay attention to the use of CURVE and plot (that we have not
# discussed so far). The following figure shows the graphs of some functions that have
# jumps, both those obtained with graph and those, incorrect, obtained with CURVE and
# plot.
# Here you see the use of GRAPH, GRAPH1, GRAPH2, faster than graph, graph1, graph2.

# The graph of a function does not display the possible holes: points where it is not
# defined (or where it is non continuos). Here's how to proceed in these cases.

HF=2.5; BF=3; F=function(x) (x^5+x^4-3*x^2-2*x+1)/(x^4-3*x+1); graphF(F, -2,4, "brown")
# I get (A). But F is non defined where x^4-3*x+1=0. Where?

graphF(F, -2,4, "brown"); G = function(x) x^4-3*x+1; graph1(G,-4,6, "red")
A=solution(G,0, 0,1); B=solution(G,0, 1,1.5); A; B; POINT(c(A,B),c(0,0),"black")
#   0.3376668    1.307486
E=1e-10; POINT(A+E,F(A+E),"yellow"); POINT(B+E,F(B+E),"yellow")
# For painting the holes I drew (in yellow color) points on the graph with abscissa
# close to that of the holes. I get (B). The final graph, (C):
graph2F(F,-2,4, "brown"); pointO(A+E,F(A+E),"blue"); pointO(B+E,F(B+E),"blue")

# We have seen how to calculate the distance of a point from another one or a straight
# line. distPF(x,y, F, a,b) calculates the distance between a point x,y and the graph
# of a 1-input-1-output function F in the domain a,b.
# Let's see, for example, how to calculate the distances of (2,0) - blue point - from
# both K and G (I want to understand which of the two curves is closest to the point)

PLANE(-1,2, -1,2)
K = function(x) 2*x-0.5; graph(K,-1,2, "red")
G = function(x) x^2; graph(G,-1,2, "brown")
POINT(2,-1, "blue"); distPF(2,-1, K, -1,2)
#dist. & nearest point (punto pił vicino):
#  2.012461    0.200000 -0.100000
POINT(0.2,-0.1, "black"); line(0.2,-0.1, 2,-1, "blue")
# in this case I would also find the distance with the point_line command
point_line(2,-1,  0,-1/2, 1,1.5)
# 2.012461
# For G, whose chart is not a straight line, I can not do it.
distPF(2,-1, G, -1,2)
#dist. & nearest point (punto pił vicino):
# 1.9490881  0.5535740 0.3064442
# I can improve accuracy starting from a smaller interval.
distPF(2,-1, G, 0.5,0.6)
# 1.9490881   0.5535738 0.3064440
POINT(0.5535740, 0.3064442, 1); line(0.5535740, 0.3064442, 2,-1, "blue")
# The curve G is closer
# Let us calculate the distance of the red point, (0,1), from G:
POINT(0,1, "red"); distPF(0,1, G, -1,2)
#  0.8660254   -0.7071070  0.5000003       I can improve precision:
distPF(0,1, G, -0.8,-0.6)
#  0.8660254   -0.7071068  0.5000000
# These numbers look like all square roots of fractions between integers. Check:
fraction(distPF(0,1, G, -0.8,-0.6)^2)
#  0.8660254   -0.7071068  0.5000000
#  3/4  1/2 1/4      The distance is √3/2
# I plot (0.707…, 0.5)
POINT(sqrt(1/2), 1/2, "seagreen"); line(0,1, sqrt(1/2),1/2, "blue")

# NOTE:  CURVE(h,…) to build the graph of h(x,y)=0 looks where h changes the sign; if
# h(x,y)≥0 or h(x,y)≤0 everywhere, the graph is not traced.  Then try to operate, as
# illustrated below, with the CURVEN (or CURN) command (see here for more information).

PLANE(-4,4, -4,4)
g = function(x,y) x^2+4*y^2-4*x*y; CURVE(g,"green")
# I get nothing ( though in this case I know the curve is
# x-2*y = 0  because  x^2+4*y^2-4*x*y = (x-2*y)^2 )
CURVEN(g,1, "blue")                    # CURN(g,1, "brown")
#[ y = F(x) +/- 1e-04 ]

# I can plot grids and write coordinates in non-decimal format. This can be useful,
# for example, in the case of trigonometric functions, but not only.  I can indicate
# the distances between the sides of the grid with the TICKx and TICKy commands and
# use the Plane2 or PLANE2 commands to trace the coordinate plane. See here for graphs
# like the following:

# K3 You can solve an equation F(x)=G(x) by defining H = F-G and using the command
# solution. But you can also use solution2(F,G, a,b), if the solution is between a and
# b.
f = function(x) x^2-2; g = function(x) x^3+1
graphF(f, -3,3, "blue"); graph(g, -3,3, "brown")
x = solution2(f,g, -2,0); x; more(x); POINT(x,f(x), "red")
#   -1.174559   -1.17455941029298

f = function(x) sin(x); g = function(x) 1-x
graphF(f, -2,2, "brown"); graph(g, -3,3, "red")
z = solution2(f,g, 0,1); more(z); POINT(z,f(z), "blue")
# 0.510973429388569

# You can also refer to WolframAlpha, but the outputs are not always easy to interpret:
# see.

# You can solve the polynomial equations where the only variable that appears is the
# unknown:   solution of   -78 + 73*x - 57*x^2 + 73*x^3 + 21*x^4 = 0
#                   [ the coefficients must be set by increasing degrees ]

q = c(-78, 73, -57, 73, 21); solpol(q)     [q -> reverse(q) if the order is decreasing]
# 0.857142857142857  -4.33333333333333
#           [from the polynomial remainder theorem, we know that a univariate
#               polynomial equation of degree N has at most N solutions]
# The outputs are temporarily put into solut (solut[1],solut[2],...)
X = solut; more(X); fraction(X)
#   0.857142857142857   -4.33333333333333
#   6/7                 -13/3
# Graphic control (see below to the left):
f = function(x) -78+73*x-57*x^2+73*x^3+21*x^4
graphF(f, -5,2, "blue")
# We can use the command polfun(a,x) to create the function a[1]+a[2]*x+a[3]*x^2+...
g = function(x) polfun(q,x); graphF(g, -5,2, "blue")

# Here's another example:   solution of  -2/9 - 2/9*x + x^3 + x^4 = 0
# There is no second degree term. I put 0 in place of it.

u = c(-2/9, -2/9, 0, 1, 1); solpol(u)
#   0.60570686427738    -1
f = function(x) -2/9 - 2/9*x + x^3 + x^4
# or:
f = function(x) polfun(u,x)
graphF(f, -1.5,1, "blue")
POINT(solut,c(0,0), "red")

# To get even the complex solutions (see here for complex numbers) instead of solpol we
# use:

u = c(-2/9, -2/9, 0, 1, 1); polyroot(u)
#   0.6057069-0i  -0.3028534+0.5245575i  -0.3028534-0.5245575i  -1+0i

# For more complex polynomial calculations, use Wolframalpha (see). Just that I put:

# quotient and remainder of  (x^5 + x^3 - 1) / (x^2 + x - 5)
# to obtain:  quotient:  x^3 - x^2 + 7*x - 12   remainder:  47*x - 61
# or   (x^5 + x^3 - 1) * (x^2 + x - 5)   to have:
# x^7 + x^6 - 4*x^5 + x^4 - 5*x^3 - x^2 - x + 5
# or   simplify  (x^2 + 1) / (x^5 - x)   to have:   1/(x^3 - x)
# or   lcm x^5 - x, x^2 + 1   to have:  (x^2 + 1)*(x^3 - x)
# or   gcd x^5 - x, x^2 + 1   to have:   x^2 + 1

# [ Here I consider polynomials where a single variable appears.  If I consider
#  a*x+a*y-y-x  like a polynomial of a single indeterminate I must specify what the
#  indeterminate is (the concept of "polynomial without indeterminates" does not make
#  sense). Here's what happens with WolframAlpha using the Collect command:
#  Collect[ a*x+a*y-y-x, a] -> a*(x+y)-x-y, Collect[ a*x+a*y-y-x, x] -> (a-1)*x+a*y-y,
#  Collect[ a*x+a*y-y-x, y] -> (a-1)*y+a*x-x.   Collect[ u*x^2+u*x+x+u,u] ->
#  u*(x^2+x+1)+x,  Collect[ u*x^2+u*x+x+u,x] ->  u*x^2+(u+1)*x+u.    See here, also
#  to factorize or expand terms:     factor a^2-b^2 -> (a-b)*(a+b),
#                                    expand (a-b)*(a+b) -> a^2-b^2    ]

# You can still add and multiply polynomials with sump and prodp.   To work on P =
# x^5 + x^3 - 1 and Q = x^2 + x - 5 you must give the coefficients as input (you need
# to put them by decreasing degree):
P = c(1,0,1,0,0,-1); Q = c(1,1,-5)
# and do this:
sump(P,Q); prodp(P,Q); sump(P,-2*Q)    # I have also subtracted from P the double of Q
#  1  0  1  1  1 -6
#  1  1 -4  1 -5 -1 -1  5
#  1  0  1 -2 -2  9
# You obtain:   x^5+x^3+x^2+x-6   x^7+x^6-4x^5+x^4-5x^3-x^2-x+5   x^5+x^3-2x^2-2x+9

# If you want use solpol you must reverse the coefficients. The solution of P-2*Q=0:
solpol(reverse(sump(P,-2*Q)))
# -1.3914711989105  the (only) solution of x^5+x^3-2x^2-2x+9=0; ie solpol(c(9,-2,-2,1,0,1))

# To find  (x+1)^7  multiply  1  by  x+1  7 times:
A = 1; B = c(1,1); for(i in 1:7) A = prodp(A,B); A              # "for" -> D-E2
#  1  7 21 35 35 21  7  1
# That is:  1 + 7x + 21x^2 + 35x^3 + 35x^4 + 21x^5 + 7x^6 + x^7

# You can multiply 2, 3, … 8 polynomials with  prodp2,…, prodp8.  Example:
# (x^5 + x^3 - 1) * (x^2 + x - 5) * (2*x - 1/4)
prodp3(c(1,0,1,0,0,-1),c(1,1,-5),c(2,-1/4))
#   2.00   1.75  -8.25   3.00 -10.25  -0.75  -1.75  10.25  -1.25
fraction(prodp3(c(1,0,1,0,0,-1),c(1,1,-5),c(2,-1/4)))
#   2   7/4 -33/4     3 -41/4  -3/4  -7/4  41/4  -5/4
# ie:  2*x^8 + 7/4*x^7 - 33/4*x^6 + 3*x^5 - 41/4*x^4 - 3/4*x^3 - 7/4*x^2 +  41/4*x - 5/4
#
# With divp you can divide a polynomial by a polynomial. See here.
#
# To repeat a sequence of commands or to modify the execution depending on the value of
# certain variables, use commands: if, ifelse, for, while. It is also often convenient
# to use connectives. Here are some simple examples:

for(i in -2:2) if(i > 0) print(i)                        #  1  2
for(i in -2:2) if(i >= 0) print(i)                       #  0  1  2
for(i in -2:2) if(i == 0) print(i)                       #  0
for(i in -2:2) if(i != 0) print(i)                       #  -2  -1  1  2
for(i in -2:2) if( !(i==0) ) print(i)                    #  -2  -1  1  2
for(i in -2:2) if(i > 0 | i <0) print(i)                 #  -2  -1  1  2
for(i in -2:2) if(i >= 0 & i <= 0) print(i)              #  0
for(i in -2:2) if(i == 0) print(i) else print (1/i)      #  -0.5  -1  0  1  0.5
for(i in -2:2) print( ifelse(i==0, i, 1/i) )             #  -0.5  -1  0  1  0.5
for(k in c(1,5,10,50)) print(1/k)                        #  1  0.2  0.1  0.02
n = 5; k = n; p = 1; while(k > 0) {print(k); p = p*k; k = k-1}; p
#  5  4  3  2  1  120
p = 1; n = 1; while(TRUE) {if(p > 10000) break else {print(p); n=n+1; p = p*n} }
#  1  2  6  24  120  720  5040
p = 1; n = 1; while(! FALSE) {if(p > 10000) break else {print(p); n=n+1; p = p*n} }
#  1  2  6  24  120  720  5040

# { and } allow you to group instructions. Try to understand by yourself what the
# various commands do; then use  help("if")   help("ifelse")   help("&")   help("TRUE")
# You can also open from the Help menu the "Html Guide" and then "Packages" (in Windows
# you can open the PDF file "R Reference" that contains a detailed help).
# break allows you to exit from the inside of a "cycle" (or even of the definition of a
# "function").
# In order to signal the end of a long cycle you can introduce the command alarm() which
# gives an audible (or visual, in some platforms) signal to the user.

# Loops (sequences of statements which may be carried out several times, like the ones
# seen above) are the basis of the computer's operation (but they are also the basis of
# every human mathematical activity). An example: the description of an algorithm to
# calculate the reciprocal using only sums, sign change and multiplication, eg 1/7:

A = 7; Y = 0.1         # Y is taken in such a way that A*Y<2
Y = Y*(2-A*Y); more(Y)      # 0.13
Y = Y*(2-A*Y); more(Y)      # 0.1417
Y = Y*(2-A*Y); more(Y)      # 0.14284777
Y = Y*(2-A*Y); more(Y)      # 0.14285714224219
Y = Y*(2-A*Y); more(Y)      # 0.142857142857143
Y = Y*(2-A*Y); more(Y)      # 0.142857142857143

# Another command: switch(n, a,b,c,...) is the n-th value in a,b,c,...
switch(4, 1/1, 1/2, 1/3, 1/4, 1/5); switch(4, "H","O","U","S","E")
#      0.25                                "S"
# switch(word, p1=X, p2=Y, ...) is X or Y or ... if word is p1 or p2 or ...
who = function(colo) switch(colo, blond="Bea", red="Lia", dark="Rosa"); who("red")
#    "Lia"
age = function(name) switch(name, Bea=19,Lia=21,Rosa=22); age("Lia")
#    21
K = c("A","B","C"); K1 = NULL; for(i in 1:3) K1[i]=switch(K[i],"A"=1,"B"=2,"C"=3); K1
#    1 2 3
K2 = NULL; for(i in 1:3) K2[i]=switch(K1[i],"1"="A","2"="B","3"="C"); K2
#    "A" "B" "C"

# subset allow you to extract "subsets" from a collection of objects specifying the
# condition that characterizes them.
data = c(118, 135, 127, 119, 121, 124, 130, 132, 128, 122, 121, 123, 127)
data1 = subset(data, data > 120); data2 = subset(data, data-100 > 20 & data-100 < 30)
data1; data2
# 135 127 121 124 130 132 128 122 121 123 127
# 127 121 124 128 122 121 123 127

# An example: the first integer smaller than 10000 with the maximum number of divisors
n=0; for(i in 1:9999) {m=length(divisors(i)); if(m>n) {n=m; j=i}}; c(j,n)
#  7560   64     7560 has 64 divisors; you can see them with divisors(7560)

# ifelse is often used to define functions to be computed in different ways in
# different parts of the domain. A simple example:
H = function(x) ifelse(x<5, -2, ifelse(x>10, x*2, x) )
# -2 if x < 5, 2x if x > 10, x elsewher, that is in [5,10]
BF=2.5; HF=2; graphF(H, -5, 15, "blue")

# Here (on the right) another example (a trick to restrict the domain):
K = function(x) ifelse(x >= -1 & x <= 2,abs(-x^2+1),x/0); V = function(x) K(x+1/2)+1
BF=2.5; HF=2; Plane(-2,3, 0,5); graph(K, -2,3, "brown"); graph(V, -2,3, "seagreen")

# The other graphic is obtained with the function DOM, that is 0 if the input belongs
# to the domain of the function FDOM, is not defined otherwise
Plane(-2,3, 0,5); FDOM=K; graph(DOM,-3,3, "brown")
FDOM=V; DOM1 = function(x) DOM(x)+2; graph(DOM1,-2,3, "seagreen")
# DOM is useful to study functions whose domain and range are not known.
# With RANGE (see) we can find the image of a function when input is between 2 values:
RANGE0(K, -1,2)
# "~ min F(min)   Max F(Max)   mFMF[1:4] = "
#    -1    0       2    3                         from 0 (in -1) to 3 (in 2)

# I can use the connectives to represent inequalities too (see):
A = function(x,y) 1 <= x^2+y^2
B = function(x,y) (1 <= x^2+y^2 & x^2+y^2 <= 2) | x^2+y^2 <= 0.5
PLANE(-2,2, -2,2); Diseq2(A,0,"seagreen"); POINT(1/2,0, "red")
PLANE(-2,2, -2,2); for(i in 1:10) Diseq2(B,0,"brown"); POINT(1/2,0, "yellow")
A(1/2, 0); B(1/2, 0)
#  FALSE      TRUE

# K4 You can find the minimum degree polynomial whose graph goes through some points
# with the commands:  xy_2 (2 points), #&133;, xy_7 (7), xy_fr (fractional form). See.
# From the points:    (-2,4) (-1,5) (2,3) (3,-1) (4,-3) (7,4)    to the polynomial:
#     211/30 + 341/360*x - 3/2*x^2 - 71/288*x^3 + 37/240*x^4 - 19/1440*x^5

# You can automatically plot the graph of a 2nd-degree polynomial function and find its
# vertex and focus, with the Parabola and ParabolaM commands: See
#        y = -2/3*x^2 + 3*x + 5/3      V: 9/4 121/24      F: 9/4 14/3

STATISTICS (and Probability)

# L-N3 By using BAR2 I can build bar charts that do not start at 0: see here.

# Note: it is meaningless (although used to fool the reader) to use inclined pie charts.
# See here.

# Besides what we have seen in (L) and (N1), there are many ways to get simply
# bar and circular charts:  see here.

# See here for the construction of tables like the following:

# Let's go deep into what we discussed in (M),(N),(N2),(N3) (and here).

# I can compare boxplots related to a given phenomenon at different times using the
# same scale, if boxAB is assigned the interval in which to plot the boxplot. Example:

# heights of 20-year-old males in 1880 in the intervals of extremes 145,150,…,185 (cm)
#              (below 145 and above 185 the percentages are negligible)

freq1 = c(4.0, 6.8, 20.4, 28.4, 23.5, 11.7, 3.7, 1.5)
interv1 = seq(145,185, 5)
BF=4; HF=3; histoclas(interv1,freq1)                 # see below, on the left

BF=4; HF=1.4; boxAB = c(145,195); morestat()        # see below, the upper boxplot

# heights of 20-year-old males in 1992 in the intervals of extremes 155,160,…,195 (cm)

freq2 = c(1.7, 6.8, 18.6, 29.1, 25.2, 13.2, 4.3, 1.1)
interv2 = seq(155,195, 5)
BF=4; HF=3; histoclas(interv2,freq2)                 # see the histogram on the right
BF=4; HF=1.4; boxAB = c(145,195); morestat()
#                               see the other boxplot (some writings have been added)

# When you have individual data that can be classified into different intervals it may
# be convenient to represent them with rectangles whose area is the fraction of the
# data falling in the intervals, so that the histogram has a total area of one, as we
# have seen for data already classified in intervals (see here). To do this we use the
# hist command. Let's see some examples regarding the beans considered in the "part 2".

# dev.new() opens a new device, abline adds grid lines through the current plot.

# [ abline(h=…)/abline(v=…) draw hor./ver. lines; abline(a,b,…) draws the line y=a+b*x
#  in the commands abline and plot, lty=3 causes the lines to be dashed, col=… colored,
#  lwd=… determines the thickness, ...]

# The number of classes is automatically selected or can be specified, roughly, by
# nclass (or n). With cex.axis I can size the characters. right=FALSE specifies that
# the intervals where the data are classified are of the form […,…).

# The essential thing is:  hist(beans, right=FALSE, probability=TRUE).

# "Density" indicates that (on the y axis) the density is represented, that is the
# unitary percentage frequencies (instead of probability=TRUE I can use freq=FALSE).
# [ see here for the unitary percentage frequencies ]

# [ with probability=FALSE or freq=TRUE I get the histogram of the absolute frequencies ]

dev.new()
hist(beans, n=5, right=FALSE, probability=TRUE, cex.axis=0.8, col="grey")
abline(h=axTicks(2),lty=3)

# and what I get with  n=20  and with  n=10  (or without specifying  nclass)

# If you want to add a histogram to an open window, append ,add=TRUE to the options.
# If you do not want to color the inside do not put ,col=…. If you want to color the
# board differently, add ,border=…. To not have the title of the window, add ,main="".
# For more, see the "details".

# To know the height of a histogram (maximum density), for example in the 1st case:

max( hist(beans, nclass=5, right=FALSE, probability=TRUE, plot=FALSE)\$density )
#   2.461538

# To know all the absolute or relative frequencies you can use hist(…)\$counts or
# hist(…)\$density.  An example with other data:

da = c(1,2,2,3,3,3,4,4,4,4,5,5,5,5,5,6,7,8,9)   #  I classify in [0,3), [3,6), [6,10)
A=hist(da,c(0,3,6,10), plot=FALSE,right=FALSE)\$counts;  A;    A/length(da)*100
#                                                3  12  4     15.78947 63.15789 21.05263

# \$ is used to extract outputs from a "command". See here and here for other examples.

# We recall that mean and percentiles for non-classified data are meaningful only if
# the numbers are rounded. For example, if I have numbers truncated to integers I must
# first convert them by adding 1/2 (data = data+1/2); if they are truncated to tenths I
# must add 1/2 tenth (data = data+ 0.1 / 2); if they are truncated to the hundreds I have
# to add 1/2 hundred (data = data + 100/2).

# To simulate (with runif and sample, and set.seed) phenomena and analyze probabilistic
# problems (use of Var, Sd, SdM, study of gaussian and other distributions), and to
# encode/decode messages, see the developments (by the end).

#         Developments HERE

# For an example of random number generation let's see the commands  Die(n), DieB(n)
# (that hide the function runif, discussed above). Die(1) and DieB(1) produce the ouput
# of 1 die. Die(N) and DieB(N) with N ≥ 2 produce N outputs and display the distribution
# histogram. One of the two dice is made of wood, the other of thin cardboard, as shown
# below. What are the two dice histograms?

# Test with N≤20000, and evaluate the probability that the output is less than or equal
# to 3 in one case and the other. Which of Die(1) and DieB(1) would you use for a non-
# make-up game?
# Die2(N) and DieB2(N) (with N>1) dispalay the distribution of N outputs of the throw
# of 2 dice.  DIE display how to build a die.

# By textAnalysis(text) I can analyze the texts statistically.  This command does not
# parse the accented letters. If you are in Windows you can use a command with the same
# name (textAnalysis) that also analyzes these characters by loading it with:

source("http://macosa.dima.unige.it/R/textanalysis.R")

# See here for using the command (and maxIsto, maxIstoN, usable to compare histograms
# of different phenomena in the same scale). Here are some of the outputs you get by
# comparing the same piece of "Alice in Wonderland" in Italiano and in English.

MORE

# O By using rm(list=ls()) you can remove every definition. Then you must run
# again  source("http://macosa.dima.unige.it/r.R")
# rm(v) and rm(v1,v2,...) remove only the variables v, v1, v2, ...

# The warning messages are ignored: if I introduce   5/sqrt(-1)
# I obtain only   NaN   (it is not a number)
# To print them (as they occur) you must type:  options(warn=0). In this case I obtain:
#     Warning message:   sqrt(-1)   NaN
# To deactivate you must type:  options(warn=-1)

# It is easy to define sequences by recursion. Three examples (see here for "while"):

# First example:   x(1)=1, x(n+1) = ( x(n) + A/x(n) ) / 2   (→ √A)
A=49/9; N=4; x = 1; for (i in 1:N) x = (x+A/x)/2; x
# 2.333335
A=49/9; N=10; x = 1; for (i in 1:N) x = (x+A/x)/2; x
# 2.333333       (L = (L+A/L)/2  ->  L^2 = A  ->  L = √A)

# Second example:   u(1) = 1, u(2) = 2, u(n+2) = ( 3*u(n+1)-u(n) )/2
U1 = 1; U2 = 2; N=10; for(i in 3:N) {U = (3*U2-U1)/2; U1 = U2; U2 = U}; U
# 2.996094
U1=1; U2=2; N=20; for(i in 3:N) {U=(3*U2-U1)/2; U1=U2; U2=U}; U
# 2.999996
U1=1; U2=2; N=30; for(i in 3:N) {U=(3*U2-U1)/2; U1=U2; U2=U}; U
# 3         ( If U1=A and U2=A+1 the limit is A+2. Why? )

# Third example:   v(1) = 1, v(n+1) = 1+1/v(n)
N=10; x = 1; for (i in 1:N) x = 1+1/x; more(x)
# 1.61797752808989
N=100; x = 1; for (i in 1:N) x = 1+1/x; more(x)
# 1.61803398874989
N=1000; x = 1; for (i in 1:N) x = 1+1/x; more(x)
# 1.61803398874989
# I can find the same value by solving x=1+1/x, ie x^2-x-1=0, for x>0
f=function(x) x^2-x-1; more(solution(f,0, 0,3))
# 1.61803398874989     It is the golden ratio (1+sqrt(5))/2
# Or:
N=100; x = 1; for (i in 1:N) x = sqrt(1+x); more(x)
# 1.61803398874989

# Another example, where I do not proceed by iteration:
exchange = function(x,y) {q=x; x=y; y=q}      # I exchange x and y
m=56; p=49; while(p>0){if(p>m) exchange(p,m); k=m; m=p; p=k%%p}; m
#  ( I put in order m,n; I replace m with m-n; I replace n with the rest of m/n;
#    I get the greatest common divisor of m and n )
# 7
m=4321; p=1234; while(p>0){if(p>m) exchange(p,m); exchange(p,m%%p)}; m
# 1
# I verify using the GCD command:
GCD( c(56,49) ); GCD( c(4321,1234) )
#  7   1

# For a more complex recursive function see here.

# P-P3  mmpaper allows to create millimeter paper (more freely than mmpaper1/3). Using
# rect for drawing rectangles I can trace histograms too. Examples:

mmpaper(40,25); x = c(1,23,38); y = c(23,9,17); polyline(x,y,"red") # the 1st picture
mmpaper(40,25)
rect(5,0, 10,17, col="red");  rect(10,0, 15,22, col="green")  # 5,0 and 10,17 are two
rect(15,0, 20,9, col="blue"); rect(20,0, 25,15, col="cyan")   #     opposite vertices
rect(25,0, 30,19,col="violet"); rect(30,0, 35,12, col="yellow")       # the 2nd picture

# mmPaper traces the millimeter paper on the same window. The third figure has been
# obtained by these further lines (lwd is the "line width"):

mmPaper(40,25)
rect(5,0, 10,17,border="brown",lwd=3);  rect(10,0, 15,22,border="brown",lwd=3)
rect(15,0, 20,9,border="brown",lwd=3);  rect(20,0, 25,15,border="brown",lwd=3)
rect(25,0, 30,19,border="brown",lwd=3); rect(30,0, 35,12,border="brown",lwd=3)

# I can add the names of the coordinates using the text command:
mmpaper(100,60)
text(100,-1.7,"100",cex=0.8)
text(0,-1.7,"0",cex=0.8)
text(-2.5,0,"0",cex=0.8,srt=90)
text(-2.5,20,"20",cex=0.8,srt=90)
text(-2.5,60,"60",cex=0.8,srt=90)
Point(0,0, "blue"); Point(100,60, "blue")
line(0,0, 100,60, "blue")
arrow(20*100/60,20, 20*10/6,0, "blue")
arrow(0,20, 20*100/60, 20, "blue")
text(20*10/6,-1.7,"~33",cex=0.8)

# I can represent circles like the ones below, on the same scale. For details see here
# (the following pictures are reduced).

New Sections

S 01 Parametric and Polar Equations. Geometric Transformations

Examples:

[commands:  param  polar  distPC  rottramul  multra  multiply ... ZZZ; see also #J]

See HERE for the commands (and more).

S 02 Inequalities and Systems of Inequalities (and plane regions)

Examples:

[commands:   diseq   diseq2   FIGURE ]

See HERE for the commands (and more).

S 03 Vector Calculus

[commands: prodv (or crossp) prods (or dotp) modulus (or magnitude) versor; see also J-2, #J]

See HERE for the commands (and more).

S 04 Lenghts and Areas

1                      1.3660254                  1.5707963

3.75              3                    1.333333          -1.333333

3.141593                     length = 24area = 6 π

12.8

[commands: areaF areaPar areaPolar lengFun lengPar lengPolar lengPar3 startP x/y Click]

See HERE for the commands (and more).

S 05 Equation systems

8x + 7y + 6z = 7.4
7x +10y +12z = 8.65
x + 2y + 4z = 1.75

x=0.55 y=0.30 z=0.15

x=15/13, y=75/13
x=3,     y=3
[commands:  eqSystem   solut]

See HERE for the commands (and more).

S 06 To estimate the best fit straight line

x:  12±1  17±1  24±1  33±1.5
y:  18±2  27±2  31±2.5  45±3

[commands:  pointDiff  pointDiff0  pointDiff2  colRect  rep]

See HERE for the commands (and more).

S 07 Centroid (and Barycenter)

C = centrePol(X,Y);  C;  POINT(C[1],C[2],"red")

[Centroid coordinates are not the average of X and Y]

[commands:  centrePol  xTab  yTab  Line  LIne]

See HERE for the commands (and more).

S 08 Integrals and Derivatives (and Limits). Taylor polynomial. Volumes

[commands:  deriv,deriv2,…,deriv6  derivxy  integral  Gintegral  eval …]

See HERE for the commands (and more).

S 09 Differential equations (1st and 2nd order). Direction fields

y'(x) = x - y(x)

[commands:   soledif    diredif]

See HERE for the commands (and more).

S 10 Calculations with "big" numbers. Strings <-> Numbers. Bases ≠ ten

x*(10+x)+5*x^3 for x = 123456789012345 -> 9408381861768148891412848380186556157340600

171! -> 124101807021766782342484052410310399261660557750169318538895180361199607522169
175299275197812048758557646495950167038705280988985869071076733124203221848436
431047357788996854827829075454156196485215346831804429323959817369689965723590
3947616152278558180061176365108428800000000000000000000000000000000000000000

as.roman(1944)  ->  MCMXLIV

MNb(4,100, 2,55)   (4/100 in base 2, 55 digits)
0. 0 0 0 0 1 0 1 0 0 0 1 1 1 1 0 1 0 1 1 1 0 0 0 0 1 0 1 0 0
0 1 1 1 1 0 1 0 1 1 1 0 0 0 0 1 0 1 0 0 0 1 1 1 1   remainder 36

## 105!
#                                    1
#                                   081
#                                  39675
#                                 8240290
#                                900504101
#                               30580032964
#                              9720646107774
#                             902579144176636
#                            57322653190990515
#                           3326984536526808240
#                          339776398934872029657
#                         99387290781343681609728
#                        0000000000000000000000000

[commands:  add pro dif app  strtoi nchar substr toString NUMBER paste0  as.roman MNb]

See HERE for the commands (and more).

S 11 Perspective. 3D

[commands:   casa()   GRAF3D   BOH()   seba(0)   seba(1) ]

See HERE for the commands (and more).

S 12 Regression. Splines. Curve Fitting. Moving average (or mean).

[commands:   regression    regression1/2/3     spline  splineC   FRUN…]

See HERE for the commands (and more).

S 13 Bezier curves

[commands:   bezier ]

See HERE for the commands (and more).

S 14 Escher

[commands:   escher ]

See HERE for the commands (and more).

S 15 Trees and other Graphs

[commands:   PQ  Ptext  ptext  PQtext  Ptext2  polyD  CLEAN ]

See HERE for the commands (and more).

S 16 Windows without the margins

[commands:    BOXW   boxW]

See HERE for the commands (and more).

S 17 Multiple graphs (in the same window)

[commands:   rowcol  colrow    Boxm_  Box_    BoxW  BoxmW    stepline  stepl]

See HERE for the commands (and more).

S 18 Animations

w*sin(x*u)+k*sin(x*v)

[commands:  Hour Minute Second Sec wait  BOXw BOXg CLEAN  fun3P curve4P polar3P  slope]

See HERE for the commands (and more).

S 19 Complex numbers. Matrices. Functions of two variables. Series

See HERE for the commands (and more).

S 20 Tree structure (and other graphs)

tree( B+C/D*E )

`+`
B
`*`
`/`
C
D
E

[commands:  tree   (and other commands) ]

See HERE for the commands (and more).

S 21 Random generator of exercises

# An example of how to use runif to randomly generate exercises

x = 3+runif(1)*5+runif(10000)^(runif(1)+0.7)*(runif(1)+5)
BF=4; HF=1.5; statistic(x)
x = 3+runif(1)*5+runif(10000)^(runif(1)+0.7)*(runif(1)+5)
BF=4; HF=1.5; statistic(x)

# and to use set.seed to generate specific exercises
set.seed(153); x = 3+runif(1)*5+runif(10000)^(runif(1)+0.7)*(runif(1)+5)
BF=4; HF=1.5; statistic(x)
set.seed(153); x = 3+runif(1)*5+runif(10000)^(runif(1)+0.7)*(runif(1)+5)
BF=4; HF=1.5; statistic(x)            # the same exercise, specified by the "seed" 153

See HERE for the commands (and more).

Command INDEX
Introduction        A
>                         source("http://macosa.dima.unige.it/r.R")
Use as calculator   B - E
* /                       ;                       .
div                       divisors
Developments    E1 - E2
Trees
Here
%/%                       quotient                remainder
factorial                 choose                  Inf  .Machine  ?Math
round  floor              e  exp                  print
signif                    more                    last
degrees                   str                     write
\n \t \b                  cat                     for
Other calculi       F - G
=                         c(…)                    sort
[ ]                       min                     max
length                    sum
Developments    G1
->  <-
Here
prod                      reverse                 combn
approx approx2            EPS
Figures             H - J
PLANE                     BF HF                   circle
POINT Point               circle3p                center3p
halfline                  line2p                  line
point_point               inclination             angle
mid2p                     perp2p                  polyline
leng                      […:…]                   polyC
Developments    J1 - J6
areaPol                   xy                      arrow
dart                      dirArrow                Direction
line_line                 solut                   point_inclina
line_line2                noPrint                 circleA
xrot  yrot                circle3pA               type
circle_circle             polylineR               polyRC
triLLL                    triLAL                  triALA
incentre                  P_tan
Here
ASP                       record/replayPlot       dev.off
ARC  ARCd                 circl circl3p arc       halfl l2p p2p
polyli                    Ncolor                  rgb
NW                        polyS polyS2            text
dart2 arrow2 D… A…        bquote                  Direction
MulTraRot                 RotTraMul               xend yend  Dot
Graphs              K
Plane                     abovex                  abovey
function                  graph                   solution
Developments    K1 - K3
graphF                    coldash                 NpF NPF NpsF
maxmin
Here
underx …X   undery U…Y    aboveY A…               RANGE
PlaneP  PLANEP            CURVE cur               graph1F graph1 …2…
plot                      distPF                  CURVEN
TICKx …y                  Plane2 PLANE2           solution2
solpol solut              polyroot                sump prodp
if ifelse for while       == != | &     { }       TRUE  break  subset
Statistics          L - N
Bar                       Pie                     Strip
ABOVE UNDER above         distribution            min
max                       length                  stem
histogram                 noClass                 median
mean                      Histogram               morestat
Developments    N1 - N4
BAR BarH PIE PIEh         rep                     Histogram
stem/scale                histoclas               pointO/D/V
Here (also Probability)
BAR2 BARM                 PIEH PieC BARC Pie_ …   boxAB
hist                      dev.new                 abline
nclass                    cex.axis                right=FALSE
runif                     sample                  set.seed
Var Sd SdM                Die DieB                textAnalysis
More                O - P
clock                     GCD                     primes
LCM                       Pie(0)                  gonio(0)
GONIO(0)                  Strip(0)                mmpaper1/2/3
SQUARES                   INPUT                   OUTPUT
r/RATIO
Developments    P1 - P4
seq  :                    PLANE/W/w/ww            Plane/W/w/ww
circleC                   GRID/w Grid/w           gridH/V GridH/V
Gonio                     CROSS CIRCLE LCROSS COUNT THERMO/2
affine
Here
mmpaper mmPaper           CiRcLe dNmM …

New Sections

S01  Parametric and Polar Equations. Geometric Transformations
param  polar  distCC distPC  rottramul  multra  multiply … ZZZ; see also #J
S02  Inequalities and Systems of Inequalities
diseq   diseq2
S03  Vector Calculus
prodv (or crossp) prods (or dotp) modulus (or magnitude) versor; see also J-2, #J
S04  Lenghts and Areas
areaF areaPar areaPolar lengFun lengPar lengPolar lengPar3 startP  x/yClick
S05  Equation systems
eqSystem   solut
S06  To estimate the best fit straight line
pointDiff  pointDiff0  pointDiff2  colRect  rep
S07  Centroid (and Barycenter)
centrePol  xTab  yTab  Line  LIne
S08  Integrals and Derivatives (and Limits). Taylor polynomial. Volumes
deriv,deriv2,…,deriv6  derivxy  integral  Gintegral  eval …
S09  Differential equations (1st and 2nd order). Direction fields
soledif    diredif
S10  Calculations with "big" numbers. Strings <-> Numbers. Bases ≠ ten
add pro dif app  strtoi nchar substr toString NUMBER paste0  as.roman MNb
S11  Perspective
casa()   GRAF3D   BOH()   seba(0)   seba(1)
S12  Regression. Splines. Curve Fitting. Moving Mean.
regression    regression1/2/3     spline
S13  Bezier curves
bezier
S14  Escher
escher
S15  Trees and other Graphs
PQ  Ptext  ptext  PQtext  Ptext2  polyD  CLEAN
S16  Windows without the margins
BOXW   boxW
S17  Multiple graphs (in the same window)
rowcol  colrow    Boxm_  Box_    BoxW  BoxmW    stepline  stepl
S18  Animations
Hour Minute Second Sec wait  BOXw BOXg CLEAN  fun3P curve4P polar3P  slope
S19  Complex numbers. Matrices. Functions of two variables. Series