# dirty.py # Andrew Davison, ad@coe.psu.ac.th, June 2026 ''' Functions implemented here: * sqrt() * log10() * power10() --> exp, ln, pow --> sinh(), asinh() * atan() --> asin(), acos() * sin() --> cos() * tan() The program includes an extensive test-rig at the end. The original BASIC versions of these functions come from: 'Quick and Dirty Functions for BASIC' Dennis Allison Dr. Dobb's Journal Volume 3, Number 9, 1978 and Dr. Dobb's Journal - Vol 3, 1978 https://archive.org/details/dr_dobbs_journal_vol_03/ ''' import math PI = 3.141593 HALF_PI = PI/2 QUARTER_PI = PI/4 TWO_PI = 2*PI SQRT10 = 3.162278 LOG10_E = 0.4342945 # log_10(e) LN10 = 2.302585 # ln(10) def sqrtD(x): if x < 0: raise ValueError("Math domain error: sqrt of negative number") if x == 0: return 0 # Scale x down by powers of 100 t = 1.0 while x > 1: x = x / 100 t *= 10 # A linear guess at sqrt(x) y = 0.115442 + 1.15442*x # y = 1 # 0.5 # Newton-Raphson iteration until convergence count = 0 while True: oldY = y y = 0.5 * (y + (x / y)) count += 1 if abs(y - oldY) < 1e-6: break # print("No. iterations:", count) # ranges from 3 to 6 return y * t # ---------- log, power, exponent-------- def log10D(x): if x <= 0: raise ValueError("Math domain error: log of non-positive number") sign = 1 if x < 1: sign = -1 x = 1/x # range reduction to [0.1, 1.0] n = 0 while x >= 1: x = x / 10 n += 1 # count powers of 10 # Transformation to range [1/sqrt(10, sqrt(10)] s = SQRT10 * x z = (s-1)/(s+1) z2 = z*z # Coefficients for the rational approximation # Hart p.110; precision = 8.66, N = M = 2; table LOG10 2323 p = 3.187822 + z2*(-2.655809 + z2*0.2668633) q = 3.670116 + z2*(-4.280973 + z2) y = (n + z*p/q - 0.5) * sign return y def pow10D(x): """ 10^x Separates the integer and fractional parts. """ isInvert = (x < 0) if isInvert: x = -x # range reduction to [0, 0.5] n = int(x) # Center the fractional part around 0.5 xf = x - n - 0.5 xf2 = xf * xf # Rational approximation coefficients # Hart: p.104 (EXPN = N^x) # Precision = 8.87, N = 1, M = 2 # EXPDC 1642 p = 69.55697 + xf2*8.751758 q = 60.41641 + xf2*(34.29514 + xf2) z = (2.0*xf*p)/(q - xf*p) # Multiply by sqrt(10) to correct the -0.5 shift # and by 10^n to shift back tn = 1 for i in range(n): tn *= 10 y = (z+1) * SQRT10 * tn if isInvert: y = 1.0 / y return y def expD(x): return pow10D(x * LOG10_E) def lnD(x): return log10D(x) * LN10 def powD(x, y): return expD(y * lnD(x)) # -------- hyperbolics ---------- def sinhD(x): # sinh(x) = (e^x − e^{−x​})/2 ex = expD(x) return (ex - (1/ex))/2 def asinhD(x): return lnD(x + sqrtD(x*x + 1.0)) # ---------- inverse trig ------- def atanD(x): sign = 1 if x < 0: sign = -1 x = -x isInvert = (x > 1) if isInvert: x = 1/x # range reduction to [0, tan pi/4 = 1] x2 = x * x # Rational approximation coefficients # Hart p.130; ARCTN 5092 p = 80.78998 + x2 * (72.58144 + x2 * 11.11774) q = 80.78999 + x2 * (99.51128 + x2 * (28.13286 + x2)) y = x * p/q if isInvert: y = HALF_PI - y # atan(x) = pi/2 - atan(1/x) return y * sign def asinD(x): ''' Could use arcsin(x) = arctan(x/sqrt(1−x^2)​) but a more stable alternative is: arcsin⁡(x) = 2*arctan⁡(x/(1+sqrt(1−x^2)) ) ''' if x < -1 or x > 1: raise ValueError("asin domain error") if x == 1: return HALF_PI if x == -1: return -HALF_PI # return atanD(x/sqrtD(1 - x*x)) return 2 * atanD(x/(1 + sqrtD(1 - x*x)) ) def acosD(x): if x < -1 or x > 1: raise ValueError("acos domain error") return HALF_PI - asinD(x) # ---------------- trig ------------------- def sinD(x): # Shift phase to use cosine approximation x = HALF_PI - x if abs(x) > 1000: raise ValueError("Input too large") if x < 0: x = -x # Range reduction to [0, 2*pi) x = x % TWO_PI # Determine quadrant (0 to 3; ccw) and sign quadrant = int(x / HALF_PI) sign = 1 if quadrant == 1 or quadrant == 2: sign = -1 x = PI - x elif quadrant == 3: x = x - TWO_PI x2 = x * x # Polynomial approximation for cosine # Hart p.118, COS 3502 y = 0.9999999 + x2*(-0.4999991 + x2*(0.04166359 + \ x2*(-0.001385370 + x2*(0.00002315393)))) return y * sign def tanD(x): # Reduce to [-pi/2, pi/2] x = x % PI if x > HALF_PI: x -= PI if abs(abs(x) - HALF_PI) < 1e-12: raise OverflowError( "tan() undefined near odd multiples of pi/2" ) sign = 1 if x < 0: sign = -1 x = -x if x <= QUARTER_PI: # Direct evaluation y = 4*x/PI return sign * tan_82s(y) else: # Use tan(x) = 1 / tan(pi/2 - x) xr = HALF_PI - x y = 4*xr/PI return sign/tan_82s(y) def tan_82s(x): '''Computes tan(pi*x/4) for x in [0, 1] 'A Guide to Approximations' Jack G. Ganssle Accurate to about 8.2 decimal digits over the range [0, pi/4]. . ''' c1= 211.849369664121 c2= -12.5288887278448 c3= 269.7350131214121 c4= -71.4145309347748 x2 = x*x return x*(c1 + c2*x2)/(c3 + x2*(c4+x2)) # --------------------------------------------- # Compare custom funcs with Python's math funcs if __name__ == "__main__": tests = { "sqrt": (sqrtD, math.sqrt, [0.5, 2, 25, 144, 0.001]), "log10": (log10D, math.log10, [0.1, 1, 10, 50, 1000]), "pow10": (pow10D, lambda x: 10**x, [0, 0.5, 1, 2, -1, -0.5]), "atan": (atanD, math.atan, [0, 0.5, 1, 5, -1]), "sin": (sinD, math.sin, [0, math.pi/4, math.pi/2, math.pi, 2*math.pi]), "tan": (tanD, math.tan, [0, math.pi/6, math.pi/4, math.pi/2 - 1e-6, -math.pi/3, 1000]) } tests2 = { "exp": (expD, math.exp, [-5, -1, 0, 1, 5]), "ln": (lnD, math.log, [0.1, 0.5, 1, math.e, 10]), "pow": ( lambda xy: powD(xy[0], xy[1]), lambda xy: math.pow(xy[0], xy[1]), [(2, 3), (9, 0.5), (10, -1), (5, 2.5)] ), "sinh": (sinhD, math.sinh, [-3, -1, 0, 1, 3]), "asinh": (asinhD, math.asinh, [-10, -1, 0, 1, 10]), "asin": (asinD, math.asin, [-1, -0.5, 0, 0.5, 1]), "acos": (acosD, math.acos, [-1, -0.5, 0, 0.5, 1]) } allTests = {} allTests.update(tests) allTests.update(tests2) print(f"{'Func':<6} | {'Input':<10} | {'Dirty':<12} | " f"{'Python':<12} | {'Error'}") print("-" * 63) for name, (dirtyF, pyF, vals) in allTests.items(): for val in vals: try: resD = dirtyF(val) resPy = pyF(val) err = abs(resD - resPy) if isinstance(val, tuple): valStr = str(val) elif isinstance(val, int): valStr = f"{val: d}" else: valStr = f"{val: .6f}" print(f"{name:<6} | {str(valStr):<10} | " f"{resD:< 12.6f} | {resPy:< 12.6f} | {err:.2e}") except Exception as e: print(f"{name:<6} | {str(val):<10} | ERROR: {e}") print()