Improve fixed-point tangent calculation for division-less operations and accuracy

When calculating the twin_tan(), division is used, specifically $\frac{\sin x}{\cos x}$

// twin_tan()
return ((s << 15) / c) << 1;

To minimize the use of the division operator, you can refer to the division implementation mentioned in https://hackmd.io/@sysprog/linux-intdiv, which describes the implementation in jemalloc to avoid direct use of the division operator.

Currently, testing and experimental data are needed to confirm that using this division method in twin_tan within mado can achieve faster calculation times.

As for the other issue, it's whether twin_tan can be improved further, including enhancements in precision and calculation speed.

Perhaps coordinate rotation digital computer(CORDIC) cloud be a solution.

The note here provide a quick overview of the mathematical principles behind conventional CORDIC and scaling-free CORDIC.

Reference:

An array architecture for fast computation of discrete Hartley transform by A. S. Dhar and S. Banerjee, 1991, mentions the use of small-angle approximations for sine and cosine functions to compensate for the need for a scaling factor $K$.

In the method, a magic number $(n−2.852)/3$ is referenced, where this magic number imposes the restriction that the number of bits shifted $j$ in the iteration must be greater than $\left\lceil \frac{n-2.585}{3} \right\rceil$.

This magic number is derived from the error estimation obtained from the remainder terms of the Taylor series approximation for sine and cosine functions. If the initial number of shifted bits $j$ does not exceed this magic number, the result of the approximation will not be accurate.

The detailed derivation process is recorded in the note: Scaling-free method of CORDIC.

By using the hybrid CORDIC, also known as the coarse-fine CORDIC algorithm [1], the lookup table size, which originally requires 15 twin_angle_t elements, can be improved to only require 5 twin_angle_t elements.

static const twin_angle_t atan_table[] = {
    0x1000, /* arctan(2^0)  = 45°    -> 4096  */
    0x0972, /* arctan(2^-1) = 26.565° -> 2418 */
    0x04fe, /* arctan(2^-2) = 14.036° -> 1278 */
    0x0289, /* arctan(2^-3) = 7.125°  -> 649  */
    0x0146, /* arctan(2^-4) = 3.576°  -> 326  */
-    0x00a3, /* arctan(2^-5) = 1.790°  -> 163  */
-    0x0051, /* arctan(2^-6) = 0.895°  -> 81   */
-    0x0029, /* arctan(2^-7) = 0.448°  -> 41   */
-    0x0014, /* arctan(2^-8) = 0.224°  -> 20   */
-    0x000a, /* arctan(2^-9) = 0.112°  -> 10   */
-    0x0005, /* arctan(2^-10) = 0.056° -> 5    */
-    0x0003, /* arctan(2^-11) = 0.028° -> 3    */
-    0x0001, /* arctan(2^-12) = 0.014° -> 1    */
-    0x0001, /* arctan(2^-13) = 0.007° -> 1    */
-    0x0000, /* arctan(2^-14) = 0.0035° -> 0   */
};

    /* Coarse rotate */
+    for (int i = 0; i < 5; i++) {
        twin_fixed_t temp_x = x;
        if (y > 0) {
            x += (y >> i);
            y -= (temp_x >> i);
            angle += atan_table[i];
        } else {
            x -= (y >> i);
            y += (temp_x >> i);
            angle -= atan_table[i];
        }
    }
    /* Fine rotate */
+    for (int i = 5; i < 15; i++) {
+        twin_fixed_t temp_x = x;
+        if (y > 0) {
+            x += (y >> i);
+            y -= (temp_x >> i);
+            angle += (0x1000 >> i);
+        } else {
+            x -= (y >> i);
+            y += (temp_x >> i);
+            angle -= (0x1000 >> i);
+        }
+    }

Since arctan can be approximated by a first-order Taylor expansion, the fine rotation angle, relative to the coarse rotation angle, can be represented by $2^{-j}$ for $i=j,\dots,n-1$ without the need for a lookup table. This allows the fine rotation part to calculate the angle using simple bit-shift operations.

Refs:

@ndsl7109256

However, adopting this method causes the arc ellipse calculation to deviate from expectations, resulting in "gaps" in the generated graphics. More analysis of the effective bit precision in the arc ellipse calculation is required.

Another interesting finding worth exploring is that adding 1 to the return value of twin_acos as

+    return twin_atan2_first_quadrant(y, x) + 1;

fills the gaps between the purple and red arc ellipses in the flowchart, regardless of whether the original method or hybrid CORDIC is used.

Left side are the arc ellipse by original CORDIC, right side are the arc ellipse by original by hybrid CORDIC. The black solid-line box indicates the location of the gap.

@ndsl7109256 I’d love to hear your thoughts or any observations you might have on this. Let me know if you’ve noticed anything interesting!

I think increasing the number of lookup table values could help. Below is a list of errors comparing the use of lookup table values to a first-order Taylor expansion, where i represents the number of lookup table values used.

According to Wikipedia [1], human vision can resolve a visual angle of one arcminute (approximately 0.017°). Based on this, we might consider making a trade-off using the provided data.

i	5	6	7
arctan(2^{-i})	1.790°	0.895°	0.448°
(0x1000<<i) * 45° / 4096	1.406°	0.703°	0.352°
error	0.384°	0.192	0.096°

And here is the comparison of using different i .

As for adding 1 to the return value of twin_acos, I didn't quite grasp the intuition behind this. Could you explain it in more detail?
Ref:
[1] https://en.wikipedia.org/wiki/Minute_and_second_of_arc#Human_vision

Consider the changes below:

--- a/src/trig.c
+++ b/src/trig.c
@@ -36,7 +36,26 @@ twin_fixed_t twin_tan(twin_angle_t a)
     if (s == 0)
         return 0;
 
-    return ((s << 15) / c) << 1;
+    /* Avoid division by using CORDIC-like approach */
+    /* For better accuracy, we use extended precision */
+    int64_t sin_ext = (int64_t) s;
+    int64_t cos_ext = (int64_t) c;
+
+    /* Handle negative values */
+    bool negative = (sin_ext < 0) ^ (cos_ext < 0);
+    if (sin_ext < 0)
+        sin_ext = -sin_ext;
+    if (cos_ext < 0)
+        cos_ext = -cos_ext;
+
+    /* Scale up for precision, then perform division */
+    int64_t result = (sin_ext << 16) / cos_ext;
+
+    /* Clamp to valid range */
+    if (result > TWIN_FIXED_MAX)
+        result = TWIN_FIXED_MAX;
+
+    return negative ? -(twin_fixed_t) result : (twin_fixed_t) result;
 }
 
 static inline twin_fixed_t sin_poly(twin_angle_t x)

It enhances twin_tan with extended precision division

Replace ((s << 15) / c) << 1 pattern that could cause overflow
Use int64_t intermediate calculations for better precision
Add XOR-based sign handling for negative values
Implement range clamping to prevent extreme value instability
Scale with 16-bit shift instead of 15+1 for improved accuracy

Preliminary evaluations:

Accuracy: <0.005 error for normal angles, 2.9 max for extreme (89°)
Memory: Zero additional overhead, no lookup tables required

twin_tan() is gone: #163 removed it together with twin_atan2(), since nothing in the tree called them, so the division in its body no longer exists. The CORDIC accuracy discussion here still applies to twin_acos() (src/trig.c), which arc flattening uses in src/path.c; that deserves its own issue if someone picks it up.

Navigation Menu

Repository navigation

Improve fixed-point tangent calculation for division-less operations and accuracy #33

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Improve fixed-point tangent calculation for division-less operations and accuracy #33

Description

Activity

jouae commented on Sep 7, 2024

jouae commented on Mar 9, 2025

jouae commented on Mar 11, 2025

ndsl7109256 commented on Mar 12, 2025

jserv commented on Jul 30, 2025

jserv commented on Oct 19, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions