HP-35 style trig functions

[Druzyek]

I've been looking at HP-35 CORDIC routines for calculating tangent using decimal numbers (here and here) so I can use them for a calculator I'm building. The site shows that once you find tangent, you can use square root and multiplies to find sine and cosine. Is there another algorithm you can use to find just sine or just cosine that doesn't require square root and multiply? Unlike the HP-35 designers, I have plenty of space for extra algorithms and tables I did try implementing the standard powers of 2 CORDIC, but since I'm working with BCD numbers, it turned out extremely slow.

[pier4r]

only as info, eric with his great works saved the 35 site as well.

https://archived.hpcalc.org/laporte/Trigonometry.htm

[Thomas Klemm]

(01-13-2019 08:32 PM)Druzyek Wrote:  Is there another algorithm you can use to find just sine or just cosine that doesn't require square root and multiply?

You may have a look at the Python examples in the attachment of Exploring the CORDIC algorithm with the WP-34S.

You can rotate by an angle and pre-calculate the radius K:

Code:

P, d = 1.0, 1.0
for i in range(N):
    P *= (1 + d*d)
    d /= 2

K = sqrt(P)

This value is a constant depending only on N.

Thus you can use (1/K, 0) as a starting value to calculate the values of cos(arg) and sin(arg) directly:

Code:

# Trigonometrics
arg = 0.4
cordic(0, 1, 1/K, 0, arg, 1.0, ATAN)
# (0.9210609940028851, 0.3894183423086506, -9.860820082244643e-20)
# cos(arg) = 0.9210609940028851
# sin(arg) = 0.3894183423086506

Kind regards
Thomas

[Druzyek]

Hi Thomas,
Thanks for the link! I checked out the Python script and it looks like you are using powers of 2 there but powers of 10 in your original post. I can get tangent with powers of 10 like you show in the post but I don't see how to get sine and cosine from that without the square root identity or powers of 2.

[Dan]

I used BCD and CORDIC algorithms in the firmware I wrote for the AriCalculator. There are CORDIC algorithms for sine and cosine but I didn't use them, instead I used a CORDIC algorithm to find tan and then the identities to find sine and cosine, as you mentioned.

Keep in mind that if you are using BCD your multiplication and division routines will be slowed by the conversion to and from binary if you use shift and add algorithms or the processor's hardware instructions for multiplication and division.

What hardware and software are you using for your project?

[Druzyek]

Hi Dan, I have admired your AriCalculator for a while!

Quote:There are CORDIC algorithms for sine and cosine but I didn't use them

Right, but do any of them use a table of atan(10^(-i))? That lets you shift X and Y to divide by 10, which is easy in BCD. Everything else I've seen (other than the links I posted) uses a table of atan(2^(-i)), which means you have to halve BCD values and that is turning out to be extremely slow.

I'm keeping everything in BCD without ever converting to binary. My numbers are 16 bytes of packed BCD. What I have so far can do multiplies in 16k cycles, division in 25k cycles, and tan in about 1.1 million cycles. 90% of the cycles in tan are for halving X and Y down before adding them. I have a lookup table for 8 bit shifts (ie 0-0x99 BCD divided by 256) and finish with 1 bit shifts but with the last CORDIC cycle requiring 106 divides by 2, it requires a huge amount of shifting. Using the method I posted above with divides by 10 only takes 60k-80k cycles! Calculating sine and cosine from that might take 2-4 times longer with the square root and multiplies, which is still very fast, but it would be nice to know how to do it without that (I'm starting to think it is not possible).

I'm using an MSP430 programmed in assembly, so no hardware multiply or divide. I used the same chip programmed in C for my first calculator: RPN Scientific Calculator. The chip has 16 bit registers and a BCD add instruction, so it is not bad to work with in assembly. Using the slow divide by 2 CORDIC I have the four functions plus the 6 trig and inverse functions in only 4k of firmware.

[Thomas Klemm]

(01-14-2019 02:34 AM)Druzyek Wrote:  Using the method I posted above with divides by 10 only takes 60k-80k cycles! Calculating sine and cosine from that might take 2-4 times longer with the square root and multiplies, which is still very fast, but it would be nice to know how to do it without that (I'm starting to think it is not possible).

Using powers of 10 you can't tell a priori how often you have to rotate by the angle \(\arctan 10^{-k}\).
The problem is that each rotation is also a dilation by the factor \(\sqrt{1+10^{-2k}}\).
And somehow you have to take that into account.

You could keep these values calculated in a table but you'd still have to multiply them to get the final radius.
I'd estimate that the average amount of multiplications are about 5 times the number of digits.

Otherwise you end up with \(x=r \cos \phi\) and \(y=r \sin \phi\) for an unknown radius \(r\).
Now you can either calculate \(\tan \phi = \frac{y}{x}\) or calculate the radius \(r=\sqrt{x^2+y^2}\) and then from this get \(\cos \phi = \frac{x}{r}\) and \(\sin \phi = \frac{y}{r}\).
This means two multiplication, one square root and two divisions.

Cheers
Thomas

[Thomas Klemm]

(01-14-2019 02:34 AM)Druzyek Wrote:  I'm starting to think it is not possible

Since you have plenty of space for lookup-tables you could pre-compute COS and SIN:

Code:

from math import radians, sin, cos

COS = [ [ cos(radians(10**i * j)) for j in range(10)] for i in range(1, -18, -1)]
SIN = [ [ sin(radians(10**i * j)) for j in range(10)] for i in range(1, -18, -1)]

arg = radians(17.188733853924695)
digits = [1, 7, 1, 8, 8, 7, 3, 3, 8, 5, 3, 9, 2, 4, 6, 9, 5]

x, y = 1.0, 0.0
for i, j in enumerate(digits):
    u, v = COS[i][j], SIN[i][j]
    x, y = x * u - y * v, x * v + y * u

print(x, y)
print(cos(arg), sin(arg))

This results in:

Code:

# 0.9553364891256062 0.2955202066613395
# 0.955336489125606 0.29552020666133955

However now 4 multiplications are needed for each digit.
Here you can use degrees for the angle.
But of course you could use radians instead.

Cheers
Thomas

[Dan]

(01-14-2019 02:34 AM)Druzyek Wrote:  Do any of them use a table of atan(10^(-i))? That lets you shift X and Y to divide by 10, which is easy in BCD. Everything else I've seen (other than the links I posted) uses a table of atan(2^(-i)), which means you have to halve BCD values and that is turning out to be extremely slow.

The algorithm presented on page 1 - 9 of this paper uses base 2 but it can be modified to use base 10. However doing so places a limit on the size of the angle since atan(10^0) + atan(10^-1) + atan(10^-2) +...< 52 degrees, whereas in the original base 2 algorithm atan(2^0) + atan(2^-1) + atan(2^-2) +... < 100 degrees, meaning the first quadrant is covered. But the author states on page 18 that

"The preceding approach could be reworked using powers of 1/10 instead of powers of 1/2; it is only necessary (as was done for the hyperbolic examples) to work with lists that contain many duplicate entries".

I recall other authors mentioning that base 10 can be used but I've never implemented CORDIC algorithms for sine and cosine so I don't know the details.

(01-14-2019 02:34 AM)Druzyek Wrote:  I'm keeping everything in BCD without ever converting to binary. My numbers are 16 bytes of packed BCD.

So numbers are 32 digits long? If you are not converting to binary then what algorithms are you using for multiplication and division? Are you using CORDIC algorithms for these as well?
EDIT: I've just seen that you have an article on your webpage on BCD multiplication, which I look forward to reading!

(01-14-2019 02:34 AM)Druzyek Wrote:  I'm using an MSP430 programmed in assembly, so no hardware multiply or divide. I used the same chip programmed in C for my first calculator: RPN Scientific Calculator. The chip has 16 bit registers and a BCD add instruction, so it is not bad to work with in assembly. Using the slow divide by 2 CORDIC I have the four functions plus the 6 trig and inverse functions in only 4k of firmware.

Great work Druzyek, I look forward to checking out the link!

[Thomas Klemm]

(01-15-2019 03:43 AM)Dan Wrote:  "The preceding approach could be reworked using powers of 1/10 instead of powers of 1/2; it is only necessary (as was done for the hyperbolic examples) to work with lists that contain many duplicate entries".

You only have to add an additional loop in the cordic function and adjust K:

Code:

from math import sqrt, sin, cos, atan

N = 16
ATAN  = [atan(10**(-i)) for i in range(N)]

P, d = 1.0, 1.0
for i in range(N):
    P *= (1 + d*d)
    d /= 10

K = sqrt(P)**9

def cordic(x, y, z, d, A):
    for a in A:
        for _ in range(9):
            s = cmp(z, 0)
            z -= s * a
            x, y = x - y * s * d, y + x * s * d
        d /= 10
    return x, y, z

This results in:

arg = 0.4
cordic(1/K, 0, arg, 1.0, ATAN)
(0.9210609940028851, 0.38941834230865136, -2.986136694229353e-16)

Compare this to:

>>> cos(arg), sin(arg)
(0.9210609940028851, 0.3894183423086505)


The cost is a possible unnecessary flipping between two states. But you can avoid the calculation of the square root.

Cheers
Thomas

[Druzyek]

Thomas,
Thanks for the examples. I had not seen the digit by digit method. The last example looks like exactly what I need. I think I can eliminate a lot of cycles in the inner loop if sign changes on the first iteration or two, since I can skip iterations that would cancel out anyway.

Quote:If you are not converting to binary then what algorithms are you using for multiplication and division? Are you using CORDIC algorithms for these as well?

Multiplication is the one exception where I convert the bytes of one argument to binary and store them so they can be halved with the Russian peasant algorithm. I leave the bytes of the other argument in BCD, so doubling is just adding the byte to itself and halving is just a shift since the other byte is in binary. This is a lot faster than the quarter squares method. Division is just subtract and shift. I get a big speed up shifting by nibbles rather than bytes.

[Druzyek]

Reviving an old thread, but it's on the same topic:

Base 10 CORDIC works fine for sine, cosine, and arctan. I'd like to use it for arcsin as well so I don't have to use identities that would require extra calculations. I found an example that calculates arcsin by comparing Y to the argument to arcsin every iteration instead of comparing Y to 0 as with arctan: link

It works with the base 2 algorithm but I can't get it to work with the changes made for base 10. Is there a way to do it?

[Thomas Klemm]

I had to slightly adjust the algorithm:

Code:

from math import sqrt, atan, asin

N = 16
ATAN  = [atan(10**(-i)) for i in range(1, N)]

P, d = 1.0, 0.1
for i in range(1, N):
    P *= (1 + d*d)
    d /= 10

K = sqrt(P)**9

def cmp(a, b):
    return (a > b) - (a < b)

def arcsin(z, d, A):
    x, y = 1/K, 0
    arg = 0
    for a in A:
        for _ in range(9):
            s = cmp(z, y)
            arg += s * a
            x, y = x - y * s * d, y + x * s * d
        d /= 10
    return x, y, arg

To calculate \(\sin^{-1}(0.4)\) use:

Code:

z = 0.4
arcsin(z, 0.1, ATAN)

(0.9165151389911695, 0.3999999999999998, 0.4115168460674864)

Compare it to the proper value:

Code:

asin(z)

0.41151684606748806

Also note that \(x = \cos(\arg)\) is calculated as well.
But you can ignore this unless you need it.

---

Just a few notes:

• I had to skip the first value in ATAN. Otherwise \(1/K\) is too small and we'd spiral once around the origin before hitting \(z\).
• As a consequence the input range is limited to \(0 \le z < 0.8\).
• I haven't figured out the proper value yet. It might be slightly above that.
• They dropped the cmp function in Python 3 so I had to provide that.