# ============================================================================
# PELL NUMBERS CALCULATOR
# Script to calculate Pell numbers and their approximations to sqrt(2)
# ============================================================================

# Clear the environment
rm(list = ls())

# ============================================================================
# TASK 1: Calculate first 200 Pell numbers and Fibonacci numbers
# ============================================================================

# Initialize the number of terms to calculate
n_terms <- 200

# Create vectors to store the sequences
index <- 1:n_terms
fibonacci <- numeric(n_terms)
pell <- numeric(n_terms)

# Set initial values for Fibonacci sequence
# F(1) = 0, F(2) = 1
fibonacci[1] <- 0
fibonacci[2] <- 1

# Set initial values for Pell sequence
# P(1) = 0, P(2) = 1
pell[1] <- 0
pell[2] <- 1

# Calculate Fibonacci numbers using a loop
# Formula: F(n) = F(n-1) + F(n-2)
for (i in 3:n_terms) {
  fibonacci[i] <- fibonacci[i-1] + fibonacci[i-2]
}

# Calculate Pell numbers using a loop
# Formula: P(n) = 2*P(n-1) + P(n-2)
for (i in 3:n_terms) {
  pell[i] <- 2 * pell[i-1] + pell[i-2]
}

# Create the dataframe with three columns
pelldata <- data.frame(
  n = index,
  Fibonacci = fibonacci,
  Pell = pell
)

# Display first few rows to verify
cat("\n=== First 10 rows of pelldata ===\n")
print(head(pelldata, 10))

# ============================================================================
# TASK 2: Calculate approximation to sqrt(2)
# ============================================================================

# Initialize vector for approximations
sqrt2_approx <- numeric(n_terms)

# Calculate approximation using loop
# Formula: sqrt(2) ≈ (P(n-1) + P(n)) / P(n)
# For n=1, we'll handle specially to avoid division issues
sqrt2_approx[1] <- NA  # First term is undefined

for (i in 2:n_terms) {
  # Approximation formula: (P(n-1) + P(n)) / P(n)
  sqrt2_approx[i] <- (pell[i-1] + pell[i]) / pell[i]
}

# Add approximation column to dataframe
pelldata$sqrt2_approximation <- sqrt2_approx

# Display results
cat("\n=== Dataframe with sqrt(2) approximations (rows 2-11) ===\n")
print(pelldata[2:11, ])

# ============================================================================
# TASK 3: Calculate difference from actual sqrt(2)
# ============================================================================

# Calculate the actual value of sqrt(2) using R's sqrt() function
actual_sqrt2 <- sqrt(2)

# Initialize vector for differences
difference <- numeric(n_terms)

# Calculate the difference for each approximation using a loop
for (i in 1:n_terms) {
  if (is.na(sqrt2_approx[i])) {
    difference[i] <- NA
  } else {
    # Calculate absolute difference and round to 5 decimal places
    difference[i] <- round(abs(sqrt2_approx[i] - actual_sqrt2), 5)
  }
}

# Add difference column to dataframe
pelldata$difference_from_sqrt2 <- difference

# Display results
cat("\n=== Dataframe with differences (rows 2-11) ===\n")
print(pelldata[2:11, ])

# ============================================================================
# TASK 4: Calculate accuracy (number of correct decimal places)
# ============================================================================

# Initialize vector for accuracy
accuracy <- numeric(n_terms)

# Calculate accuracy using a loop
for (i in 1:n_terms) {
  if (is.na(sqrt2_approx[i])) {
    # First term is undefined
    accuracy[i] <- NA
  } else {
    # Calculate the absolute difference
    diff <- abs(sqrt2_approx[i] - actual_sqrt2)
    
    # Determine number of correct decimal places
    # A decimal place is correct if the difference is less than 0.5 * 10^(-decimal_places)
    # We'll count how many decimal places match
    
    if (diff == 0) {
      # Perfect match (unlikely with floating point)
      accuracy[i] <- 15  # Maximum precision in R
    } else {
      # Calculate number of correct decimal places
      # Using -log10(diff) gives us approximately the number of correct places
      # We subtract 0.3 to be conservative (ensures we count only fully correct places)
      correct_places <- floor(-log10(diff) - 0.3)
      
      # Ensure non-negative
      if (correct_places < 0) {
        correct_places <- 0
      }
      
      accuracy[i] <- correct_places
    }
  }
}

# Add accuracy column to dataframe
pelldata$decimal_accuracy <- accuracy

# ============================================================================
# FINAL OUTPUT
# ============================================================================

# Display final dataframe structure
cat("\n=== Final dataframe structure ===\n")
str(pelldata)

# Display first 20 rows of complete dataframe
cat("\n=== First 20 rows of complete pelldata ===\n")
print(pelldata[1:20, ])

# Display last 10 rows to show convergence
cat("\n=== Last 10 rows showing high accuracy ===\n")
print(pelldata[191:200, ])

# Summary statistics
cat("\n=== Summary Statistics ===\n")
cat("Actual sqrt(2):", actual_sqrt2, "\n")
cat("Approximation at n=200:", sqrt2_approx[200], "\n")
cat("Difference at n=200:", difference[200], "\n")
cat("Decimal accuracy at n=200:", accuracy[200], "places\n")

# Create a summary table for key values
cat("\n=== Key Approximations ===\n")
key_indices <- c(2, 5, 10, 20, 50, 100, 150, 200)
summary_table <- pelldata[key_indices, ]
print(summary_table)

# Optional: Save the dataframe to a CSV file
# Uncomment the line below if you want to save the results
# write.csv(pelldata, "pell_numbers_results.csv", row.names = FALSE)

cat("\n=== Script completed successfully ===\n")