From 03486b80e0b8bac1af77234dc7f0325113dde8eb Mon Sep 17 00:00:00 2001 From: Theodor Chikin Date: Tue, 7 Oct 2025 19:22:47 +0300 Subject: [PATCH] created and debugged precise fixed-pint sin and cos --- naive_DFT.py | 111 +++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 103 insertions(+), 8 deletions(-) diff --git a/naive_DFT.py b/naive_DFT.py index 2788976..bcf45ce 100755 --- a/naive_DFT.py +++ b/naive_DFT.py @@ -3,7 +3,7 @@ from math import sin, cos, pi, sqrt import plotly.graph_objs as go from plotly.subplots import make_subplots -FP_acc = 1e2 +FP_acc = 1e3 INP_L = 1000 @@ -20,6 +20,10 @@ def abs_f(re, im): def abs_FP(re, im): return int(sqrt(re*re + im*im)) + +def sqrt_FP(val): + print(val) + return int(sqrt(val)) def DFT_naive(inp, out): for f in range(len(out)): @@ -39,9 +43,90 @@ def abs_FP(re, im): # return sqrt(re*re + im*im) return int(sqrt(re*re + im*im)) +trigon_debug = 0 def sin_FP(phi_fp): - return sin(pi*(phi_fp/FP_acc)/(pi_FP/FP_acc)) + if (trigon_debug): + print("sin_FP========") + print("phi:", phi_fp) + if phi_fp < 0: + if (trigon_debug): + print("phi < 0. recursive inversion...") + return -1 *sin_FP(-1*phi_fp) + + while phi_fp >= 2*pi_FP: + if (trigon_debug): + print("phi is bigger than 2Pi. Decreasing...") + phi_fp -= 2*pi_FP + if (trigon_debug): + print("phi:", phi_fp) + + if phi_fp >= pi_FP: + if (trigon_debug): + print("phi > pi_FP. recursive inversion...") + print(phi_fp, pi_FP) + return -1*sin_FP(phi_fp - pi_FP) + + if phi_fp == pi_FP/2: + return 1*FP_acc + if phi_fp == 0: + return 0 + + if phi_fp > pi_FP/2: + if (trigon_debug): + print("phi > pi_FP/2. recursive inversion...") + return sin_FP(pi_FP - phi_fp) + + #now phi should be inside [0, Pi/2). checking... + if phi_fp < 0: + raise ValueError('error in sin_FP. after all checks phi < 0') + if phi_fp >= pi_FP/2: + raise ValueError('error in sin_FP. after all checks phi > pi_FP/2') + + #now phi is inside [0, Pi/2). So, cos(phi) > 1 always + if (trigon_debug): + print("phi:", phi_fp) + return sin_FP_constrained(phi_fp) + + +sin_05_debug = 1 + +def sin_FP_constrained(phi_fp): + phi_trh = pi_FP/16 + if (trigon_debug): + print("sin_FP_constrained===========") + print("phi:", phi_fp) + print("check is phi inside [0, Pi/2)") + if phi_fp < 0: + raise ValueError('error in sin_FP. after all checks phi < 0') + if phi_fp >= pi_FP/2: + raise ValueError('error in sin_FP. after all checks phi > pi_FP/2') + if (trigon_debug): + print("Ok") + if (phi_fp > phi_trh): + if (sin_05_debug): + print("phi_fp:", phi_fp," >",phi_trh,"... recusion") + print("phi_fp/2:", phi_fp/2) + sin_phi_05 = sin_FP_constrained(phi_fp/2) + sin_phi = (2*sin_phi_05 * sqrt_FP(FP_acc**2 - sin_phi_05*sin_phi_05))/FP_acc + if (sin_05_debug): + print("sin_phi_05:", sin_phi_05) + print("sin_phi:",sin_phi) + return sin_phi + else: + res = int((phi_fp * 3.141592653589793238462643383279502884197169399375105820974944 * FP_acc)/FP_acc - ((phi_fp * 3.141592653589793238462643383279502884197169399375105820974944 * FP_acc)/FP_acc)**3/(6*FP_acc**2)) +# res = int((phi_fp * 1.0 * FP_acc)/FP_acc - ((phi_fp * 1.0 * FP_acc)/FP_acc)**3/(6*FP_acc**2)) + if (trigon_debug): + print("calculating sin(x) as x:",phi_fp, res) + + return res + + + + + +# return sin(pi*(phi_fp/FP_acc)/(pi_FP/FP_acc)) + def cos_FP(phi_fp): return sin_FP(phi_fp - pi_FP/2) @@ -54,8 +139,11 @@ def DFT_naive_FP(inp_float, out): val_im = 0 for n in range(len(inp)): phi = 2*pi_FP*f*n/INP_L - val_re += inp[n] * sin_FP(phi) /INP_L - val_im += inp[n] * cos_FP(phi) /INP_L + phi_sin = sin_FP(phi) + phi_cos = cos_FP(phi) + print(phi, phi_sin, phi_cos) + val_re += inp[n] * phi_sin /INP_L + val_im += inp[n] * phi_cos /INP_L val_abs = abs_FP(val_re, val_im) print("F, val_abs:",f, val_abs) @@ -82,16 +170,23 @@ def main(): def sin_tester(): - N = 400 - angs = [(N/2 - i)/100 for i in range(N)] + N = 4000 + angs = [(i - N/2)/1000 for i in range(N)] res_f = [] res_FP = [] + res_cos_FP = [] for phi in angs: res_f.append(sin(phi*pi)) - res_FP.append(sin_FP(phi*FP_acc*pi_FP)) +# print(phi, phi*FP_acc*pi_FP/FP_acc) + val_fp = sin_FP(phi*FP_acc*pi_FP/FP_acc) + val_cos_fp = cos_FP(phi*FP_acc*pi_FP/FP_acc) + print("angle, sin, cos:",phi, val_fp, val_cos_fp) + res_FP.append(val_fp) + res_cos_FP.append(val_cos_fp) chart = go.Figure() chart.add_trace(go.Scatter(x = angs, y=res_f, name="sin_float", mode="markers+lines")) - chart.add_trace(go.Scatter(x = angs, y=res_f, name="sin_FP", mode="markers+lines")) + chart.add_trace(go.Scatter(x = angs, y=[val/FP_acc for val in res_FP], name="sin_FP", mode="markers+lines")) + chart.add_trace(go.Scatter(x = angs, y=[val/FP_acc for val in res_cos_FP], name="cos_FP", mode="markers+lines")) chart.show()