# File: pell.r
# Module: ACCA7022 – Advanced Statistical Programming
# Assignment: Data Manipulation (Task 1 – Code)
# Author: Karthik
# Date: Oct 22, 2025
# Purpose:
#   1) Generate the first 200 Pell numbers P_n and Fibonacci numbers F_n using LOOPS.
#   2) Build a dataframe 'pelldata' with columns:
#        n, Fn (Fibonacci), Pn (Pell), Sqrt2_Approx, Diff_vs_sqrt2_5dp, Accuracy_DP
#   3) The √2 approximation uses the classical Pell-based convergents:
#        H_n / P_n  where  H_n are the half-companion Pell numbers.
#      Using the identity H_n = P_n + P_{n-1}, we can compute this with loops only.
#   4) For each n, compute:
#        - Signed difference (Approx − sqrt(2)), displayed to five decimals.
#        - Accuracy as the number of correct decimal places compared to sqrt(2).
###############################################################################

# Housekeeping
# Ensure reproducible numeric display in any quick checks (optional)
options(scipen = 999)

# Store the true value of sqrt(2) once (to be reused)
true_sqrt2 <- sqrt(2)

# Configuration
N <- 200                     # Number of terms required
# Pre-allocate vectors
# I keep indexing aligned to mathematical convention:
#   n = 0, 1, 2, ..., 199  (total 200 terms)
n_idx    <- 0:(N - 1)        # Index column n

F <- numeric(N)              # Fibonacci numbers  F_0=0, F_1=1, F_n = F_{n-1} + F_{n-2}
P <- numeric(N)              # Pell numbers       P_0=0, P_1=1, P_n = 2*P_{n-1} + P_{n-2}

approx_s2        <- rep(NA_real_, N)  # √2 approximation column
diff_signed_5dp  <- rep(NA_real_, N)  # Signed difference to 5 decimals (Approx − sqrt(2))
accuracy_dp      <- rep(NA_integer_, N) # Correct decimal places (integer)

# Build Fibonacci by LOOP
# Base cases
F[1] <- 0  # F_0
if (N >= 2) F[2] <- 1  # F_1

# Loop for remaining terms: for n = 2,3,... (in R's 1-based indexing: positions 3..N)
if (N >= 3) {
  for (i in 3:N) {
    # i corresponds to mathematical index n = i-1
    F[i] <- F[i - 1] + F[i - 2]
  }
}

# Build Pell by LOOP
# Pell base cases
P[1] <- 0  # P_0
if (N >= 2) P[2] <- 1  # P_1

# Loop for remaining terms: P_n = 2*P_{n-1} + P_{n-2}
if (N >= 3) {
  for (i in 3:N) {
    P[i] <- 2 * P[i - 1] + P[i - 2]
  }
}

# Compute √2 approximations, signed differences, and accuracy
# I use the classical convergents to √2:
#   H_n / P_n  with H_n = P_n + P_{n-1}.
# Note:
#   - For n = 0 (i == 1), P_0 = 0 so the fraction is undefined; we leave NA.
#   - For n >= 1 (i >= 2), Approx_n = (P_n + P_{n-1}) / P_n.
# Accuracy metric:
#   Number of correct decimal places d such that:
#       Approx - sqrt(2)| < 0.5 * 10^{-d}
#   We compute d = floor( -log10( 2 * |error| ) ), with bounds [0, 15].
#   (15 is a practical cap due to double precision.)
max_cap_dp <- 15L

for (i in 1:N) {
  if (i == 1) {
    # n = 0: undefined approximation because P_0 = 0
    approx_s2[i]       <- NA_real_
    diff_signed_5dp[i] <- NA_real_
    accuracy_dp[i]     <- NA_integer_
  } else {
    # n = i-1 >= 1
    # H_n = P_n + P_{n-1}
    Hn <- P[i] + P[i - 1]
    # Approximation to sqrt(2)
    approx_val <- Hn / P[i]
    approx_s2[i] <- approx_val
    
    # Signed difference (Approx − sqrt(2)), shown to five decimals
    diff_signed_5dp[i] <- round(approx_val - true_sqrt2, 5)
    
    # Absolute error for accuracy calculation
    abs_err <- abs(approx_val - true_sqrt2)
    
    # Handle the exact-equality edge case (theoretically never happens for irrational √2)
    if (abs_err == 0) {
      accuracy_dp[i] <- max_cap_dp
    } else {
      # d = floor( -log10( 2 * |error| ) )
      # If error is large, this can go negative; enforce lower bound at 0.
      raw_dp <- floor(-log10(2 * abs_err))
      if (is.na(raw_dp) || raw_dp < 0) raw_dp <- 0L
      if (raw_dp > max_cap_dp) raw_dp <- max_cap_dp
      accuracy_dp[i] <- as.integer(raw_dp)
    }
  }
}

# Assemble the dataframe 'pelldata'
pelldata <- data.frame(
  n                 = n_idx,
  Fn                = F,
  Pn                = P,
  Sqrt2_Approx      = approx_s2,
  Diff_vs_sqrt2_5dp = diff_signed_5dp,
  Accuracy_DP       = accuracy_dp,
  # Ensure integers display as such (optional; keeps numerics for arithmetic)
  check.names = FALSE
)

#  Light-touch verification / preview
# Print top rows to visually confirm structure and early convergents:
# Expected early convergents for √2: 1/1=1.00000, 3/2=1.50000, 7/5=1.40000, 17/12≈1.41667, ...
# which correspond to n = 1,2,3,4 with (H_n, P_n) pairs.
cat("\nPreview (first 10 rows):\n")
print(head(pelldata, 10))

# Show the best (highest accuracy) within the first 200 terms:
best_row <- NA_integer_
best_acc <- -1L
for (i in 1:N) {
  # Skip NA rows
  if (!is.na(accuracy_dp[i]) && accuracy_dp[i] > best_acc) {
    best_acc <- accuracy_dp[i]
    best_row <- i
  }
}
cat("\nBest within first", N, "terms:\n")
if (!is.na(best_row)) {
  print(pelldata[best_row, , drop = FALSE])
} else {
  cat("No valid approximation rows found.\n")
}

# ---------------------------
# Save to CSV for documentation/testing convenience
# ---------------------------
write.csv(pelldata, file = "/Users/mac/Desktop/ACCA7002 Advanced Statistical Programming/Karthik/pelldata.csv", row.names = FALSE)