Talk:CORDIC

Text of introduction on main page

FYI the introduction only talks about the tan function but the problem description is to implement sin and cos. We should probably put the problem description there in the page, but could mention that CORDIC can also be used to implement other functions such as tan or log_10 (e.g. the implementation at https://github.com/bikerglen/cordic/blob/master/log10.cpp )

Validity

The task description currently says that the algorithm is valid for values of alpha up to π/2 but it seems that that limit should be π/4. The accuracy of the algorithm at values slightly smaller than π/4 isn't very good (for example: at 0.784 radians), but it's *much* better than the accuracy of the algorithm at values slightly larger than π/4 (for example 0.786 radians). --Rdm (talk) 13:31, 10 July 2023 (UTC)

Using the Algol 68 version, I get the following:

   angle  cordic cos       cos     difference  cordic sin       sin     difference  cordic tan       tan     difference
  0.7810:   0.710185  0.710185  2.3725799e-10 |  0.704015  0.704015  2.3933777e-10 |  0.991311  0.991311  6.6818373e-10
  0.7819:   0.709571  0.709571  3.73174525e-9 |  0.704634  0.704634  3.75788911e-9 |  0.993043  0.993043  1.05185740e-8
  0.7828:   0.708956  0.708956  1.7983937e-10 |  0.705253  0.705253  1.8078339e-10 |  0.994778  0.994778  5.0734306e-10
  0.7837:   0.708340  0.708340  3.44724094e-9 |  0.705872  0.705872  3.45929518e-9 |  0.996515  0.996515  9.73335601e-9
  0.7845:   0.707724  0.707724  6.95373681e-9 |  0.706489  0.706489  6.96588398e-9 |  0.998256  0.998256  1.96510264e-8
  0.7854:   0.707107  0.707107  1.1102230e-16 |  0.707107  0.707107  0.00000000e+0 |  1.000000  1.000000  1.1102230e-16
  0.7863:   0.706489  0.706489  3.27581606e-9 |  0.707724  0.707724  3.27010385e-9 |  1.001747  1.001747  9.27351729e-9
  0.7871:   0.705872  0.705872  6.79132328e-9 |  0.708340  0.708340  6.76765832e-9 |  1.003497  1.003497  1.92424936e-8
  0.7880:   0.705253  0.705253  3.22337856e-9 |  0.708956  0.708956  3.20654492e-9 |  1.005250  1.005250  9.14117892e-9
  0.7889:   0.704634  0.704634  6.74470002e-9 |  0.709571  0.709571  6.69777700e-9 |  1.007006  1.007006  1.91443004e-8
  0.7898:   0.704015  0.704015  3.17049509e-9 |  0.710185  0.710185  3.14294724e-9 |  1.008765  1.008765  9.00724273e-9

--Tigerofdarkness (talk) 14:18, 10 July 2023 (UTC)

Among the scientific articles I have read on the topic, a few mention a [0,π/4] validity but the others stick to [0,π/2]. For me, the above issue comes from accumulation of rounding errors: the algorithm starts by trying to withdraw π/4 from the angle: if it can't, it will then substract arctan(0.1) ≈ 3π/100 a number of times, generating a bigger rounding error. To improve accuracy, it probably makes sense to find the appropriate 'CORDIC window' by calling a trigonometric identity, such as tan(2a) = 2tan(a)/(1-tan²(a)), before starting the calculation. Successors of the HP-35, the first pocket calculator with CORDIC inside, compute mathematical functions with 'extended real numbers' (mantissa with 3 extra digits) to avoid these errors.

--Aerobar (talk) 20:22, 10 July 2023 (UTC)

I see I should have said a bit more... Doesn't the accuracy depend on how many values you have in the arctan table and the value of epsilon ?

It seemed OK to me and not particularly less accurate around pi/4. For some reason, I assumed we should try to minimise the number of iterations so I made epsilon 10^-9 in the Algol 68 sample - with smaller epsilons, more accurate values can be obtained.

--Tigerofdarkness (talk) 21:40, 10 July 2023 (UTC)

The size of the arctan table is completely defined by the internal floating point precision, but the key for accuracy shall lie in the appropriate selection of epsilon - there even might be different values for it, depending on the number to compute. One can also decide to exit the main loop after a predefined number of iterations, considering that overcomputing destroys precision. The expertise of early computer scientists was in such tweaking!

--Aerobar (talk) 20:53, 25 July 2023 (UTC)

OK, I hadn't realised that in the pseudo code:

constant θ[n] = arctan 10^(-n) // or simply 10^(-n) depending on floating point precision,

the comment meant that the table should be extended as necessary with smaller values. After studying the RPL sample, I see it starts with a 6 element arctan table but this grows to include values down to 1e-12.

So for IEEE 64-bit floats which can represent values with about 16 decimal digits, epsilon should be 1e-16 ?

Using epsilon = 1e-16 and an arctan table extending down to 1e-16, I get:

    angle cordic cos       cos     difference  cordic sin       sin     difference  cordic tan       tan     difference
 -9.0000:  -0.911130 -0.911130  1.1102230e-16 | -0.412118 -0.412118  0.00000000e+0 |  0.452316  0.452316  0.00000000e+0
  0.0000:   1.000000  1.000000  0.00000000e+0 |  0.000000  0.000000  0.00000000e+0 |  0.000000  0.000000  0.00000000e+0
  1.5000:   0.070737  0.070737  7.4940054e-16 |  0.997495  0.997495  0.00000000e+0 | 14.101420 14.101420  1.5099033e-13
  6.0000:   0.960170  0.960170  1.1102230e-16 | -0.279415 -0.279415  5.5511151e-16 | -0.291006 -0.291006  6.1062266e-16
    angle cordic cos       cos     difference  cordic sin       sin     difference  cordic tan       tan     difference
  0.7810:   0.710210  0.710210  3.3306691e-16 |  0.703990  0.703990  2.2204460e-16 |  0.991242  0.991242  6.6613381e-16
  0.7819:   0.709576  0.709576  2.2204460e-16 |  0.704629  0.704629  1.1102230e-16 |  0.993028  0.993028  4.4408921e-16
  0.7828:   0.708942  0.708942  3.3306691e-16 |  0.705267  0.705267  4.4408921e-16 |  0.994817  0.994817  9.9920072e-16
  0.7837:   0.708307  0.708307  6.6613381e-16 |  0.705905  0.705905  5.5511151e-16 |  0.996609  0.996609  1.6653345e-15
  0.7846:   0.707671  0.707671  5.5511151e-16 |  0.706542  0.706542  3.3306691e-16 |  0.998405  0.998405  1.3322676e-15
  0.7855:   0.707035  0.707035  1.1102230e-16 |  0.707179  0.707179  1.1102230e-16 |  1.000204  1.000204  2.2204460e-16
  0.7864:   0.706398  0.706398  1.1102230e-16 |  0.707815  0.707815  2.2204460e-16 |  1.002006  1.002006  4.4408921e-16
  0.7873:   0.705761  0.705761  1.1102230e-16 |  0.708450  0.708450  1.1102230e-16 |  1.003811  1.003811  2.2204460e-16
  0.7882:   0.705123  0.705123  3.3306691e-16 |  0.709085  0.709085  5.5511151e-16 |  1.005619  1.005619  1.3322676e-15
  0.7891:   0.704484  0.704484  5.5511151e-16 |  0.709720  0.709720  4.4408921e-16 |  1.007431  1.007431  1.3322676e-15

--Tigerofdarkness (talk) 11:57, 26 July 2023 (UTC)