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

import math
import decimal
from decimal import Decimal as D

PI = D("3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679")

PHI = D("1.6180339887498948482045868343656381177203091798057628621354486227052604628189024497072072041893911374")

E = D("2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274")

SQRT2 = D("1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727")



def pi():
  # Compute Pi to the current precision  
  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)


def atanEuler(x):
  """
    Calculate arctan(1/x) using Euler's accelerated formula.
    Each term is made from the previous term by 
    multiplying by 2n/((2n+1)(1+x^2)) 
  """
  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()



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

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))))


