Talk:CORDIC

From Rosetta Code
Latest comment: 4 months ago by Tigerofdarkness in topic Validity

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)Reply[reply]

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)Reply[reply]

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)Reply[reply]

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)Reply[reply]
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)Reply[reply]
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)Reply[reply]