# dPade.py
# Andrew Davison, ad@coe.psu.ac.th, June 2026
"""
Obtain a Pade approximation for a function f(x), 
using its Taylor series coefficients. 
The series' expansion point is stored in the x0s
list.

A list of functions is presented, and the chosen one is 
given a Pade approximation based on the specified (n,m) 
order for that function.

The approx is evaluated against a series of values, and 
the results are compared to the actual function. 

The functions listed here are the ones implemented using
rational approximations in dirty.py, a Python
version of the code in 'Quick and Dirty Functions for BASIC'
by Dennis Allison.
"""

import math
from fractions import Fraction
import faUtils
import sympy as sp



def errReport(fName, func, p, q, testVals, x0=0):
  print(f"\n{fName} Comparisons:")
  print("      x        Actual         Pade"
        "         Abs Err       Rel Err" )
  print( "-"*70)

  tot = 0.0
  validCount = 0
  for xVal in testVals:
    try:
      actual = func(xVal)
      approx = float( faUtils.evalPade(xVal, p, q, x0))
    except (ValueError, ZeroDivisionError, TypeError):
      print(f"{xVal:8.3f}   (domain error, skipped)")
      continue

    absErr = abs(approx - actual)
    relErr = abs((approx - actual) / actual) \
                           if actual != 0 else 0.0
    tot += relErr
    validCount += 1
    print(f"{xVal:8.3f}   {actual:12.8f}   {approx:12.8f}   "
            f"{absErr:12.8f}  {relErr:12.8f}")

  if validCount == 0:
    print("\nNo valid test points.")
    return None
  avg = tot / validCount
  print(f"\nAverage Relative Error = {avg:.8e}")
  return avg


# -----------------------

x = sp.symbols('x')

fNames = ["log10(x)", "pow(10,x)", "atan(x)"]   # was log10(x)

spFns = [lambda x: sp.log(x,10), lambda x: sp.Pow(10,x), 
         sp.atan ]
fns = [lambda x: math.log10(x), lambda x: math.pow(10,x), 
       math.atan ]

dataLists = [ [0.2, 0.4, 0.6],  [0.1, 0.3, 0.5],
              [0.1, 0.5, 0.9] ]
orders = [ (2,2), (1,2), (2,3)]   #(n,m)
x0s = [0.5, 0.25, 0.5]


print("Functions:", fNames)
i = int(input("Enter index of function name: "))
n, m = orders[i]
print(f"\nApproximating {fNames[i]} for order {orders[i]}")
p, q, _ = faUtils.padeApprox(spFns[i], n, m, x0=x0s[i])
errReport(fNames[i], fns[i], p, q, dataLists[i], x0=x0s[i])