Solutions

Exercise 3

1. Factorial:

   Using While loop

   fact1 <- function(n){
     if(n<0) return(NA)
     result <- 1
     while(n > 0){
       result <- result * n
       n <- n - 1
     }
         result
   }
   fact1(5)

   [1] 120

   Using repeat

   fact2 <- function(n){
     if(n<0) return(NA)
     result <- 1
     repeat{
       if(n < 2) break
       result <- result * n
       n <- n - 1
     }
     result
   }
   fact2(6)

   [1] 720

   Using For-loop

   fact3 <- function(n){
     if(n<0) return(NA)
     result <- 1
     for(i in seq_len(n)){
       result <- result * i 
     }
     result
   }
   fact3(4)

   [1] 24

2. Displaying the Multiplication Table

   {
     for(i in seq(9)) cat("\t", i)
     cat("\n")
     for(i in seq_len(9)){
       cat(i)
       for(j in seq_len(9)){
         cat("\t", i*j)
       }
       cat("\n")
     }
   }

        1   2   3   4   5   6   7   8   9
   1    1   2   3   4   5   6   7   8   9
   2    2   4   6   8   10  12  14  16  18
   3    3   6   9   12  15  18  21  24  27
   4    4   8   12  16  20  24  28  32  36
   5    5   10  15  20  25  30  35  40  45
   6    6   12  18  24  30  36  42  48  54
   7    7   14  21  28  35  42  49  56  63
   8    8   16  24  32  40  48  56  64  72
   9    9   18  27  36  45  54  63  72  81

3. grade(x)

   grade <- function(x){
     if(x >= 90) 'A'
     else if(x >= 80) 'B'
     else if(x >= 70) 'C'
     else if(x >= 60) 'D'
     else if(x >= 50) 'E'
     else 'F'
   }
   grade(93)

   [1] "A"

4. grades(x)

   grades <- function(x){
     for(i in x) cat(grade(i), " ")
   }
   grades(c(73, 92, 80, 49))

   C  A  B  F  

5. \(\pi\) estimation:

   result <- 0
   for(i in 0:200000){
     result <- result + 4*(-1)^i/(i*2+1)
   }
   sprintf("%.30f", result)

   [1] "3.141597653564761838396179882693"

   sprintf("%.30f", pi)

   [1] "3.141592653589793115997963468544"

6. \(\pi\) estimation:

   const <- 2*sqrt(2)/9801
   res <- 0
   for(k in 0:2){
     res <- res + factorial(4*k)*(1103 + 26390*k)/factorial(k)^4/396^(4*k)
   }
   1/const/res

   [1] 3.141593

   sprintf("%.30f", pi)

   [1] "3.141592653589793115997963468544"

   sprintf("%.30f", 1/const/res)

   [1] "3.141592653589792671908753618482"

7. sqrt_avg

   sqrt_avg <- function(x){
     begin = 0
     end = x
     repeat{
       mid <- (begin + end)/2
       err <- x - mid * mid
       if(abs(err) < 1e-8) break
       if(err > 0) begin <- mid
       else end <- mid
     }
     mid
   }
   sqrt_avg(2)

   [1] 1.414214

8. sqrt_heron

   sqrt_heron <- function(x){
     x0 <- 1
     repeat{
       x_new <- (x0 + x/x0)/2
       if(abs(x_new - x0) < 1e-8) break
       x0 <- x_new
     }
     x0
   }
   sqrt_heron(2)

   [1] 1.414214

Exercise 4

1. Fibonacci:

   fib1 <- function(x){
     a = 0
     b = 1
     while(x>1){
       temp <- a
       a <- b 
       b <- b + temp
       x <- x - 1
     }
     b
   }
   fib1(10)

   [1] 55

2. GCD

   gcd <- function(a, b){
     if(b == 0) a
     else gcd(b, a%%b)
   }

3. LCM

   lcm <- function(a, b){
     a * b /gcd(a, b)
   }

4. Taylor Series:

   exp(x) : Will correctly approximate values from\(-11< x < 55\). Notice that for numbers outside this interval, we need more iterations.

   my_exp <- function(x){
     result <- fact <- i <- 1 #taken the first term as result
     while(i <= 100){
       result <- result + x^i/fact
       i <- i + 1
       fact <- fact * i
     }
     result
   }

   log(x) : Used the first 50 terms. Only approximates \(0.2\leq x\leq2\) . The error is huge for numbers outside this interval

   my_log_restricted <- function(x){
     if(x>2 | x<0.2) return(NA)
     i <- 1
     result <- 0
     while(i <= 100){
       result <- result + (-1)^(i-1) * (x-1)^i/i
       i <- i + 1
     }
     result
   }

   sin(x) : \(-\infty <x<\infty\)

   my_sin <- function(x){
     x <- x %% (2*pi)
     fact <- i <- j <- 1
     result <- 0
     while(i <= 50){
       result <- result + (-1)^(i-1)*x^(2*i-1)/fact
       i <- i + 1
       fact <- fact * (2*i - 1 - 1) * (2*i - 1)
     }
     result
   }

5. Significant Digits:

   signif_digits <- function(x, base){
     x <- abs(x)
     if(round(x,6) >= 1 & x < base) 0
     else {
       pow <- (-1)^(x < 1)
       pow + signif_digits(x/base^pow, base)
     }
   }

   signif_digits_1 <- function(x, base){
     x <- abs(x)
     result <- 0
     pow <- if(x < 1) -1 else 1
     while(round(x, 6) < 1 | x >= base){
       result <- result + pow
       x <- x / base^pow
     }
     result  
   }

6. Digit Sum

   digit_sum <- function(n){
     n <- abs(n)
     if(n < 10) return(n)
     n %% 10 + Recall(n%/%10)
   }

7. Digital Root

   digit_root <- function(n){
     if(abs(n)>10) digit_root(digit_sum(n))
     else n
   }

8. Logarithms:

   my_log_using10 <- function(x){
      const <- 2.302585092994045901094
      b <- signif_digits(x, 10)
      a <- x/10^b
      if(a>2) {
        a <- a/10
        b <- b + 1
      }
      my_log_restricted(a) + b*const
   }

   my_log_using2 <- function(x){
      const <- 0.6931471805599452862268
      b <- signif_digits(x, 2)
      a <- x/2^b
      my_log_restricted(a) + b*const
   }

   my_log2 <-function(x){
     if(x==1) return(0)
     w <- 1
     if(x<1) {
       x <- 1/x
       w <- -1
     }
     res <- 0
     n <- 0
     for(i in 1:20){
       m <- 0
       while(x<2){
         m <- m+1
         x <- x**2
       }
       x <- x/2
       n <- m + n
       res <- res + 2^-n
     }
     w * res
   }

9. Optimization

   x <- 1
   repeat{
     f <- x^2*exp(3*x) - 10
     f_prime <- 2*x*exp(3*x) + 3*x^2*exp(6*x)
     x_new <- x - f/f_prime
     if(abs(x_new - x)<1e-8) break
     x <- x_new
   }
   x

   [1] 0.8645552

   x^2*exp(3*x)

   [1] 10

   x <- 1
   repeat{
     f <- x^3 + 8 
     f_prime <- 3*x^2
     x_new <- x - f/f_prime
     if(abs(x_new - x)<1e-8) break
     x <- x_new
   }
   x

   [1] -2

   x^3

   [1] -8

   my_sqrt <- function(x){
     y <- 1
     repeat{
       f <- y^2 - x 
       f_prime <- 2*y
       y_new <- y - f/f_prime
       if(abs(y_new - y)<1e-8) break
       y <- y_new
     }
     y
   }
   my_sqrt(49)

   [1] 7

   my_sqrt(2)

   [1] 1.414214

10. my_log_newton <- function(x){
      if(x==1)return(0)
      y <- 1
      repeat{
        f <- my_exp(y) - x 
        f_prime <- my_exp(y)
        y_new <- y - f/f_prime
        if(abs(y_new - y)<1e-16) break
        y <- y_new
      }
      y
    }

11. my_optim <- function(x0, f, fprime){
      repeat{
        x <- x0 - f(x0)/fprime(x0)
        if(abs(x - x0)<1e-8) break
        x0 <- x
      }
      x0
    }
    
    my_optim(1, function(z) z^3 + 8, function(x)3*x^2)

    [1] -2

12. Numerical Differentiation/Derivatives

    my_derive <- function(x, f){
       h <- 1e-8
      (f(x+h) - f(x))/h
    }

    my_optim2 <- function(f){
      x <- 1
       repeat{
        x_new <- x - f(x)/my_derive(x, f)
        if(abs(x_new - x)<1e-8) break
        x <- x_new
       }
      x
    }
    
    my_optim2(function(x)x^2*exp(3*x)-10)

    [1] 0.8645552

13. Better Fibonacci Recursion

    my_fib <- function(x, start = 0, end = 1){
      if(x == 1)end
      else Recall(x - 1, end, start + end)
    }
    my_fib(30) # Compare the times with the previous fib function

    [1] 832040

    my_fib(200)# Do not run fib(200)

    [1] 2.805712e+41

14. Tower of Hanoi

    tower_of_hanoi <- function(height, from, to, via){
      if(height == 1) 
        cat("Move disk",height,"from", from, "to", to, "\n")
      else {
        tower_of_hanoi(height-1, from, via, to)
        cat("Move disk",height,"from", from, "to", to, "\n")
        tower_of_hanoi(height-1, via, to, from)
      }
    }
    tower_of_hanoi(4, "A", "B", "C")

    Move disk 1 from A to C 
    Move disk 2 from A to B 
    Move disk 1 from C to B 
    Move disk 3 from A to C 
    Move disk 1 from B to A 
    Move disk 2 from B to C 
    Move disk 1 from A to C 
    Move disk 4 from A to B 
    Move disk 1 from C to B 
    Move disk 2 from C to A 
    Move disk 1 from B to A 
    Move disk 3 from C to B 
    Move disk 1 from A to C 
    Move disk 2 from A to B 
    Move disk 1 from C to B