Computing unit 1

22. Any cubic equation can be put into the form x3 + ax2 + bx + c = 0. Assuming the coefficients are real, such an equation has at least one real root, which can be computed in closed form as follows. Make the substitution y = x + a/3. The equation then reduces to y3 + py + q = 0, p = 3b − a2 3 , q = 2a3 27 − ab 3 + c. If the discriminant D = (p/3)3 + (q/2)2 is positive, then there is only one real root, and it is given by x = − a 3 + 3 √ − q 2 + √D + 3 √ − q 2 − √D. Write a function which uses this method (in real arithmetic) to compute the unique real root of a cubic equation.

# input coefficients of the cubic equation f12<-function(a,b,c){ p<-(3*b-a^2)/3 q<-2*a^3/27-a*b/3+c D<-(p/3)^3+(q/2)^2 if (D>0) { root<- -a/3+sign(-q/2+sqrt(D))*abs(-q/2+sqrt(D))^(1/3)+sign(-q/2-sqrt(D))*abs(-q/2-sqrt(D))^(1/3) root } else { writeLines("The polynomial does not have a unique root") } }

1. Calculate the remainder after dividing 31 079 into 170 166 719.

170166719 %% 31079

Using seq() and rep() as needed, create the vectors 0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 and 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

a<-rep(seq(1:5)-1, each =5) b<-rep(seq(1:5),5) print(a) print(b)

7. Using seq() and rep() as needed, create the vector 1 2 3 4 5 2 3 4 5 6 3 4 5 6 7 4 5 6 7 8 5 6 7 8 9

c<-rep(seq(1,5),5)+rep(seq(1:5)-1, each =5) print(c)

19. Write a function which computes the sum of the digits of a given natural number. Hint: you can use %% and %/% to extract and eliminate the last digit, respectively.

checksum_func <- function(x) { sum(x %/% 10^(0:(nchar(x))) %% 10) }

20. Write a function which computes the total number of digits (in the decimal representation) in a given sequence of natural numbers

f10<- function(x){ leng<-nchar(x) sum(leng) }

21. Write a function which takes three numbers as arguments and returns the sum of the squares of the two larger numbers.

f11<-function(x,y,z){ all<-c(x,y,z) all<-all[!(all==min(all))] sum(all^2) }

10. Write a function which saves the 26 upper-case letters of the Roman alphabet to a given file.

f1<-function(file_name) save(LETTERS, file = file_name)

12. Write a function which outputs whether a given number is positive or negative.

f2<-function(x){ if (x>0) {writeLines("The number is positive")} if (x==0) {writeLines("The number is zero")} if (x<0) {writeLines("The number is negative")} }

13. Write a function which outputs whether a given number exceeds π in absolute value

f3<-function(x){ if (abs(x)>pi) {writeLines("The number exceeds Pi in absolute value")} else {writeLines("The number does not exceed Pi in absolute value")} }

16. Write a function which outputs a given character string once for each character in the string. (Hint: ?"for".)

f6<-function(x){ for (i in 1:nchar(x)) cat(x, "\n") }

17. Write a function which for a given natural number n outputs all natural numbers ≤ n in reverse order, each on a separate line ("countdown").

f7<-function(x){ for (i in x:1) cat(i, "\n") }

18. For two natural numbers x and y, x %/% y and x %% y give the result and remainder, respectively, of the integer division of x by y. Write a function which outputs whether a given x is divisible by a given y (i.e., whether the remainder is zero, or equivalently, whether x is a multiple of y).

f8<-function(x,y){ if (x%%y==0) writeLines("The first input is divisible by the second input") else writeLines("The first input is not divisible by the second input") }

8. Write a function which, for a given natural number n, returns a sequence where each i ≤ n is repeated i times, in ascending order. E.g., for n = 4 the function should return c(1, 2, 2, 3, 3, 3, 4, 4, 4, 4)

f<-function(n){ sequence<-c() for (k in 1:n) sequence<-c(sequence,rep(k,k)) sequence }

23. Write a function factorial to compute n! = ∏n i=1 i using prod(). (a) Compute 10!, 50!, 100!, and 1000!. (b) Use factorial to write a function which computes the binomial coefficient ( n m ). Com- pute (4 2 ), (50 20 ), and (5000 2000 ). (c) The sum() function can used to sum up all elements of its arguments, while the log() and exp() functions take the natural logarithm and exponential of all elements in their arguments. Use these functions to create an improved function to compute the binomial coefficients. Compute (4 2 ), (50 20 ), and (5000 2000 ) using the improved version.

factorial<-function(x) prod(1:x) # a factorial(10) factorial(50) factorial(100) factorial(1000) # b binom<-function(n,m){ factorial(n)/(factorial(m)*factorial(n-m)) } binom(4,2) binom(50,20) binom(5000,2000) # c improved<-function(n,m){ exp(sum(log(1:n))-sum(log(1:m))-sum(log(1:(n-m)))) } improved(4,2) improved(50,20) improved(5000,2000) # The improved approach scales the numbers down and sums # them instead of multiplying them, while the regular approach struggles with # multiplying infinities.

14. Write a function which returns the geometric (or the harmonic) mean of two given numbers.

mean_func <- function(kind, x, y) { if (kind == "geometric") { return(sqrt(x*y)) } else if (kind == "harmonic") { return(2/(1/x+1/y)) } else { return("Enter a valid mean") } } mean_func("geometric", 10, 20) mean_func("harmonic", 10, 20) mean_func("arithmetic", 10, 20)

Calculate the sum sn = ∑n i=1 i and compare with n(n + 1)/2 for n = 100, 200, 400, 800. Use the quick formula to compute the sum for all values of n between 1 and 100, and store the results in a vector.

n = c(100, 200 ,400, 800) for (i in n) { long <- sum(1:i) short <- i*((i+1)/2) print(c(long, short)) print(long == short) } # part 2 results <- (1:100)*((1:100)+1)/2

Calculate the sum sn = ∑n i=1 i2 and compare with n(n + 1)(2n + 1)/6 for n = 200, 400, 600, 800. Use the quick formula to compute the sum for all values of n between 1 and 100, and store the results in a vector.

n = c(200, 400, 600, 800) for (i in n){ long <- sum((1:i)^2) short <- i * (i + 1) * (2 * i + 1)/6 print(c(long, short)) print(long == short) } # part 2 n <- (1:100) n * (n + 1) * (2 * n + 1)/6

9. Create a vector called numbers which contains 3 5 8 10 12 Dump numbers to a file called 'numbers.R' and delete numbers using rm(). Using ls(), confirm that numbers has been deleted. Now, use source() to retrieve the vector numbers

numbers<-c(3,5,8,10,12) dump("numbers", file = "numbers.R") rm("numbers") ls() source("numbers.R")

11. Plot the graph of the function f (x) = { 3x + 2, x ≤ 3 2x − 0.5x2, x > 3 on the interval [0, 6]. (Do not worry about the discontinuity at x = 3.)

plot(NULL, xlim=c(0,6), ylim=c(-10,15)) curve( 3*x + 2, from=0, to=3,add=TRUE) curve(2*x-0.5*x^2,from=3,to=6,add=TRUE)

Calculate the sum sn = ∑n i=1 ri for r = 1.08 and compare it to (rn+1 − 1)/(r − 1) − 1 for n = 10, 20, 30, 40. Use the quick formula to compute the sum for all values of n between 1 and 100, and store the results in a vector.

r <- 1.08 n <- c(10,20,30,40) for (i in n) { x <- 1:i long <- sum(r^x) short <- ((r^(i+1) - 1) / (r-1)) - 1 print(c(long, short)) print(long == short) } # part 2 (r^((1:100) +1) - 1) / (r - 1) - 1

24. Write a function to compute ρn = Γ((n − 1)/2) Γ(1/2)Γ((n − 2)/2) . 2 Note that a straightforward implementation using gamma to compute the Gamma function will fail for large n. Why? (Try n = 2000 for example). Instead, compute ρn using the logarithm of the Gamma function (provided by lgamma). Can you guess the limit of ρn/√n for n → ∞?

rhoofntest<-function(n){gamma((n-1)/2)/(gamma(1/2)*gamma((n-2)/2))} rhoofntest(2000) rhoofn<-function(n){exp(lgamma((n-1)/2)-(lgamma(1/2)+lgamma((n-2)/2)))} rhoofn(2000) # The implementation using gamma() does not work, because of the fast rising # nature of the gamma function leads to an overflow # resulting in the expression Inf/Inf rhoofn(1000)/sqrt(1000) rhoofn(10000)/sqrt(10000) rhoofn(100000)/sqrt(100000) rhoofn(1000000)/sqrt(1000000) # the limit of this expression tends as n -> Inf is ~0.4 # This is our guess for the limit : 1/sqrt(2*pi) # For this we used the Stirling approximation for n!

Calculate the sum sn = ∑n i=1 1/i, and compare with log(n) for n = 500, 1000, 2000, 4000, 8000. Do you know the limit of sn − log(n) for n → ∞?

seq_sum3<-c() logarithm<-c() n3<-c(500, 1000, 2000, 4000, 8000) for (i in 1:5){ seq_sum3[i]<-sum(1/(1:n3[i])) logarithm[i]<-log(n3[i])} print(logarithm) print(seq_sum3) sum(1/(1:1000000))-log(1000000) #the difference between harmonic series and log(n) converges to #Euler's constant gamma

15. Write a function which outputs which of two given character strings is shorter. (Hint: ?"nchar".)

string_length_func <- function(char1, char2) { if (nchar(char1) > nchar(char2)){ return(paste0(char1, " is longer than ", char2)) } else if (nchar(char1) == nchar(char2)) { return(paste0(char1, " is the same length as ", char2)) } else { return(paste0(char1, " is shorter than ", char2)) } } string_length_func("abc","abc")