# decTrig.py
# https://docs.python.org/3/library/decimal.html#recipes

import math
import decimal
from decimal import Decimal as D


def pi():
  # Compute Pi to the current precision  
  # https://math.stackexchange.com/questions/3028868/what-is-the-formula-for-pi-used-in-the-python-decimal-library
  decimal.getcontext().prec += 2  
          # extra digits for intermediate steps
  t, prevS, s = D(3), 0, 3
  n, na, d, da = 1, 0, 0, 24
  while s != prevS:
    prevS = s
    n, na = n+na, na+8
    d, da = d+da, da+32
    t = (t * n) / d
    s += t
  decimal.getcontext().prec -= 2
  return +s         # unary plus applies the new precision


def cos(x):
  # Return the cosine of x as measured in radians.
  # cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + ...
  decimal.getcontext().prec += 2
  i, prevS, s = 0, 0, 1 
  fact, num, sign = 1, 1, 1
  while s != prevS:
    prevS = s  
    i += 2
    fact *= i * (i-1)
    num *= x * x
    sign *= -1
    s += num / fact * sign 
  decimal.getcontext().prec -= 2    
  return +s


def sin(x):
  # Return the sine of x as measured in radians
  # sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...
  decimal.getcontext().prec += 2
  i, prevS, s = 1, 0, x
  fact, num, sign = 1, x, 1
  while s != prevS:
    prevS = s  
    i += 2
    fact *= i * (i-1)
    num *= x * x
    sign *= -1
    s += num / fact * sign 
  decimal.getcontext().prec -= 2    
  return +s


def atan(x):
  if x == 0:
    return D(0)
  elif x > 0:
    return atanPos(x)
  else:
    return D(-1)*atanPos(-x)

def atanPos(x):
  if x > D(1):
    return pi()/D(2) - atanEuler(x)
  else:
    return atanEuler(D(1)/x)

"""
  Calculate arctan(x) using Euler's accelerated formula.
  https://en.wikipedia.org/wiki/Inverse_trigonometric_functions
  Each term is made from the previous term by 
  multiplying by 2n/((2n+1)(1+x^2)) 
"""
def atanEuler(x):
  decimal.getcontext().prec += 2

  x2p1 = D(1) + (x * x)
  term = x / x2p1
  total = term
  n2 = D(2)
  while True:
    term *= n2 / ((n2+1) * x2p1)
    if term == D(0):
      break
    total += term
    n2 += D(2)
  decimal.getcontext().prec -= 2
  return +total


def radians(degrees):
  return degrees * pi()/D(180)

def degrees(radians):
  return radians * D(180)/pi()


def printGrouped(num, spaces=4):
  # print the digits in groups of 5
  digs = list(spaces*" "+str(num))
  for i in range(len(digs)):
    print(digs[i], end='')
    if i % 5 == 0:
      print(" ", end ='')
    if i % 60 == 0:
      print()
  print()


def countDigits(digs):
  counts = [0]*10
  for i in range(len(digs)):
    if digs[i].isdigit():
      counts[int(digs[i])] += 1
  return counts


def findFact(d):
  f = 1
  n = 2
  while math.log10(f) < d:
    f = f*n
    n += 1
  print("log10(fact("+ str(n-1) + ")) >=", d)

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

if __name__ == "__main__":
  decimal.getcontext().prec = 100

  degree = float(input("degree? "))
  angle = math.radians(degree)

  print('\n   Math sin:', math.sin(angle))
  print('Decimal sin:', sin(D(angle)))

  print('\n   Math cos:', math.cos(angle))
  print('Decimal cos:', cos(D(angle)))

  print("\n   Math pi:", math.pi)
  print("Decimal pi:",  pi())

  x = float(input("x? "))

  print("\n   Math atan:", math.degrees(math.atan(x)))
  print("Decimal atan:",  degrees(atan(D(x))))