#########################################################
#   ACCA7022 - CW1 Advanced Statistical Programming
#   Student Name: Anand
#   File: pell.r
#   Description:
#     This program builds a data frame called pelldata
#     for the first 200 terms of three number sequences:
#       1. Fibonacci numbers
#       2. Pell numbers
#       3. Rational approximations to √2
#     It then calculates how close each approximation is
#     to the true √2 value and how many decimal places match.
#     Finally, it saves the output as an Excel file.
#########################################################

# Load the library needed to write Excel files
# (Install first if not already installed: install.packages("openxlsx"))
library(openxlsx)

# Let's start by deciding how many terms we want.
N <- 200

# Store the real square root of 2 (used for comparison).
true_sqrt2 <- sqrt(2)

# Create empty vectors to fill later.
n_values <- integer(N)        # This will store index numbers (0,1,2,…)
fib_values <- numeric(N)      # Fibonacci sequence
pell_values <- numeric(N)     # Pell sequence
approx_values <- numeric(N)   # Approximations to √2
diff_values <- numeric(N)     # Difference between approx and real √2
accuracy_values <- integer(N) # How many decimals are correct

# -----------------------------
# 1. Create the index values
# -----------------------------
# We’ll simply fill 0 to 199 into our n_values.
for (i in 1:N) {
  n_values[i] <- i - 1
}

# -----------------------------
# 2. Generate Fibonacci numbers
# -----------------------------
# Rules: F0 = 0, F1 = 1, and F(n) = F(n-1) + F(n-2)
fib_values[1] <- 0
fib_values[2] <- 1

# Start the loop from 3rd element
for (i in 3:N) {
  fib_values[i] <- fib_values[i - 1] + fib_values[i - 2]
}

# -----------------------------
# 3. Generate Pell numbers
# -----------------------------
# Rules: P0 = 0, P1 = 1, and P(n) = 2*P(n-1) + P(n-2)
pell_values[1] <- 0
pell_values[2] <- 1

for (i in 3:N) {
  pell_values[i] <- 2 * pell_values[i - 1] + pell_values[i - 2]
}

# -----------------------------
# 4. Approximation of √2
# -----------------------------
# We’ll build rational approximations (called convergents)
# using the continued fraction for √2: [1; 2, 2, 2, 2, ...]
# Formula:
#   p[k] = a[k]*p[k-1] + p[k-2]
#   q[k] = a[k]*q[k-1] + q[k-2]
#   a0 = 1, all others a[k] = 2
# Approximation = p[k] / q[k]

p_prev2 <- 0  # p[-2]
p_prev1 <- 1  # p[-1]
q_prev2 <- 1  # q[-2]
q_prev1 <- 0  # q[-1]

for (k in 0:(N - 1)) {
  # Decide the coefficient a[k]
  if (k == 0) {
    a_k <- 1
  } else {
    a_k <- 2
  }
  
  # Compute numerator and denominator
  p_k <- a_k * p_prev1 + p_prev2
  q_k <- a_k * q_prev1 + q_prev2
  
  # Store the approximation value
  approx_values[k + 1] <- p_k / q_k
  
  # Shift for next iteration
  p_prev2 <- p_prev1
  p_prev1 <- p_k
  q_prev2 <- q_prev1
  q_prev1 <- q_k
}

# -----------------------------
# 5. Find difference and accuracy
# -----------------------------
# For each approximation, we’ll:
#   - Find the difference vs real √2
#   - Round it to 5 decimal places
#   - Count how many decimal digits are correct

for (i in 1:N) {
  # (a) Difference from real sqrt(2)
  diff_values[i] <- round(approx_values[i] - true_sqrt2, digits = 5)
  
  # (b) Compare digit by digit to count matching decimals
  real_str <- formatC(true_sqrt2, format = "f", digits = 12)
  approx_str <- formatC(approx_values[i], format = "f", digits = 12)
  
  real_parts <- strsplit(real_str, "\\.")[[1]]
  approx_parts <- strsplit(approx_str, "\\.")[[1]]
  
  real_int <- real_parts[1]
  real_dec <- strsplit(real_parts[2], "")[[1]]
  approx_int <- approx_parts[1]
  approx_dec <- strsplit(approx_parts[2], "")[[1]]
  
  if (real_int != approx_int) {
    accuracy_values[i] <- 0
  } else {
    count <- 0
    for (j in 1:length(real_dec)) {
      if (real_dec[j] == approx_dec[j]) {
        count <- count + 1
      } else {
        break
      }
    }
    accuracy_values[i] <- count
  }
}

# -----------------------------
# 6. Combine all into one table
# -----------------------------
pelldata <- data.frame(
  n = n_values,
  Fn = fib_values,
  Pn = pell_values,
  Sqrt2_Approx = approx_values,
  Diff_vs_sqrt2_5dp = diff_values,
  Accuracy_DP = accuracy_values
)

# -----------------------------
# 7. Show sample results
# -----------------------------
cat("Here are the first 10 rows:\n")
print(pelldata[1:10, ])

cat("\nAnd the last 5 rows (just to check the end):\n")
print(pelldata[(N - 4):N, ])

# -----------------------------
# 8. Save to Excel file
# -----------------------------
# The Excel file will be created in your working directory.
write.xlsx(pelldata, "pelldata.xlsx", rowNames = FALSE)

cat("\nExcel file 'pelldata.xlsx' has been created successfully in your working directory.\n")

#########################################################
# End of Script
#########################################################