import decimal
import numpy as np

from decimal import *
getcontext().prec = 200

# de functie np.log werkt niet in combinatie met de module decimal, dus alles moet met de functie np.log10
# de module math is geen optie, want die is veel te onnauwkeurig
# de constanten e en pi moet je zelf nauwkeurig meegeven, want binnen decimal zijn die onnauwkeurig

e = Decimal ('2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901')
pi = Decimal ('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196')

# 10^0

x = Decimal ('1')
getalexact = Decimal ('1')
exponentexact = Decimal ('0')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))

# 10^1

x = Decimal ('10')
getalexact = Decimal ('3.6288')
exponentexact = Decimal ('6')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))

# 10^2

x = Decimal ('100')
getalexact = Decimal ('9.3326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864')
exponentexact = Decimal ('157')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))

# 10^3

x = Decimal ('1000')
getalexact = Decimal ('4.0238726007709377354370243392300398571937486421071463254379991042993851239862902059204420848696940480047998861019719605863166687299480855890132382966994459099742450408707375991882362772718873252013122')
exponentexact = Decimal ('2567')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))

# 10^4

x = Decimal ('10000')
getalexact = Decimal ('2.8462596809170545189064132121198688901480514017027992307941799942744113400037644437729907867577847758158840621423175288300423399401535187390524211613827161748198241998275924182892597878981242531043265')
exponentexact = Decimal ('35659')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))

# 10^5

x = Decimal ('100000')
getalexact = Decimal ('2.8242294079603478742934215780245355184774949260912248505789180865429779509010630178725517714138311636107136117373619629514749961831239180227260734090938324220055569688667840380377379444961268442101980')
exponentexact = Decimal ('456573')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))

# 10^6

x = Decimal ('1000000')
getalexact = Decimal ('8.2639316883312400623766461031726662911353479789638730451677758855633796110356450844465305113114639733516068042108785885414647469506478361823012109754232995901156417462491737988838926919341416003974215')
exponentexact = Decimal ('5565708')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))

# 10^7

x = Decimal ('10000000')
getalexact = Decimal ('1.2024234005159034561401534879443075697676801824947563081172508508669676920037337277075969393672212361228664996104650925368418125437973581763283743859957027075176261489825713568024383286731249165156653')
exponentexact = Decimal ('65657059')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))

# 10^8

x = Decimal ('100000000')
getalexact = Decimal ('1.6172037949214623863387731856128040432923745306487350797898134626667371544075043865805008769291907857306593069977095263790720953162951796541671790735130090066341293377470615910986852970969058342604960')
exponentexact = Decimal ('756570556')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))

# 10^9

x = Decimal ('1000000000')
getalexact = Decimal ('9.9046265792229937372808211050657043217172213900293119277105585378984692554496119536646786298014599266365248886088755466982273949436921553304917026919138198142995066054469923125356087065662932655248803')
exponentexact = Decimal ('8565705522')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))

# 10^10

x = Decimal ('10000000000')
getalexact = Decimal ('2.3257962056730833651049447199498788053980140675288356663530707441440170169965578757085252986087824677345004147127681644274203880836311510117131391958249885618171595910465617288605909670517066642034499')
exponentexact = Decimal ('95657055186')
logexact = Decimal (exponentexact) + Decimal (np.log10 (Decimal (getalexact)))
lnstirling = Decimal (x) * Decimal (np.log10 (Decimal (x))) / Decimal (np.log10 (Decimal (e))) \
 - Decimal (x) \
 + Decimal (np.log10 (Decimal (np.sqrt (Decimal ('2') * Decimal (pi) * Decimal (x))))) / Decimal (np.log10 (Decimal (e))) \
 + Decimal ('1') / (Decimal ('12') * Decimal (x)) \
 - Decimal ('1') / (Decimal ('360') * Decimal (x) ** Decimal ('3'))
logstirling = Decimal (lnstirling) * Decimal (np.log10 (Decimal (e)))
exponentstirling = Decimal (np.floor (Decimal (logstirling)))
getalstirling = Decimal ('10') ** (Decimal (logstirling) - Decimal (exponentstirling))
ratio = Decimal ('10') ** (Decimal (logstirling) - Decimal (logexact))
relerror = Decimal (np.log10 (Decimal (np.abs (Decimal (ratio) - Decimal ('1')))))

print (x)
print ('{:42.40f}'.format(getalstirling))
print (exponentstirling)
print ('{:42.40f}'.format(ratio))
print ('{:4.2f}'.format(relerror))