• Welcome to the new Internet Infidels Discussion Board, formerly Talk Freethought.

Calculating Ad Approximating Math Functions

S

steve_bank

Diabetic retinopathy and poor eyesight. Typos ...
Joined
Nov 9, 2017
Messages
13,769
Location
seattle
Basic Beliefs
secular-skeptic
Working solutions to c;calculating math functions.
 

Calculating sin with a Taylor series. The problem becomes the size of the factorials. A table of inverse factorials is created for speed.

As x s raised to a power the range of x would hee to be limited to 0-2*pi to avoid overflow. Sin is periodic. sn(pi/4) = 1. Sni() of any integer multiple of pi/4 is 1. For x > 2*pi x would have to be scaled down to 0-2*pi.

Code: 
#Python sin Taylor series
import math as ma
import array as ar
import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt
import warnings

# table 1/n! n = 1,3,5.....
factab = [
       1.0,
        0.16666666666666666,
        0.008333333333333333,
        0.0001984126984126984,
        2.7557319223985893e-06,
        2.505210838544172e-08,
        1.6059043836821613e-10,
        7.647163731819816e-13,
        2.8114572543455206e-15,
        8.22063524662433e-18,
        1.9572941063391263e-20,
        3.8681701706306835e-23,
        6.446950284384474e-26,
        9.183689863795546e-29,
        1.1309962886447718e-31,
        1.2161250415535181e-34,
        1.151633562077195e-37,
        9.67759295863189e-41,
        7.265460179153071e-44,
        4.902469756513544e-47,
        2.9893108271424046e-50,
        1.6552108677421951e-53,
        8.359650847182804e-57,
        3.866628513960594e-60,
        1.643974708316579e-63,
        6.446959640457174e-67,
        2.3392451525606576e-70,
        7.876246304918039e-74,
        2.4674957095607893e-77,
        7.210682961895936e-81,
        1.9701319568021682e-84,
        5.043860616493007e-88,
        1.2124664943492803e-91,
        2.74189618803546e-95,
        5.843768516699617e-99,
        1.1758085546679308e-102,
        2.2370786808750587e-106,
        4.030772397973078e-110,
        6.887854405285506e-114,
        1.117795262136564e-117,
        1.724992688482352e-121,
        2.53451761457883e-125,
        3.549744558233656e-129,
        4.7443792545223946e-133,
        6.057685462873333e-137,
        7.396441346609687e-141,
        8.64474210683694e-145,
        9.680562269694223e-149,
        1.039579281539328e-152,
        1.0715102881254669e-156]

def sin_taylor(x):
    n = 50
    k = 1
    s = 0
    for i in range(n):
        # i even add term else subtract term
        if (i% 2) == 0:
             s += pow(x,k)*factab[i]
        else:
            s -= pow(x,k)*factab[i]
        k += 2
    return s


n = 20
xmax = ma.pi*2  # sin() 360 degrees
x = np.ndarray(shape=(n,1),dtype = np.float64)
y = np.ndarray(shape=(n),dtype = np.float64)  # interpolated sin
ys = np.ndarray(shape=(n),dtype = np.float64) # sin() library function
x = np.linspace(0,xmax,n)  # 2*pi radians

print("sin taylor")
for i in range(n):
    y[i] = sin_taylor(x[i])
    ys[i] = ma.sin(x[i])
    print("%2.8f  %+2.15f   %+2.15f" %(x[i],ys[i],y[i]))


[fig, ax] = plt.subplots()
ax.plot(x, y, linewidth=2.0,color='r')
ax.plot(x, ys, linewidth=2.0,color='k')
ax.grid(color='k', linestyle='-', linewidth=2)
plt.show()


A double has around 15 decimal digits depending on the number. 16 digits shown.


Code: 
x radians       library sin()         Taylor series
0.00000000  +0.0000000000000000   +0.0000000000000000
0.33069396  +0.3246994692046835   +0.3246994692046835
0.66138793  +0.6142127126896678   +0.6142127126896678
0.99208189  +0.8371664782625285   +0.8371664782625283
1.32277585  +0.9694002659393304   +0.9694002659393303
1.65346982  +0.9965844930066698   +0.9965844930066700
1.98416378  +0.9157733266550575   +0.9157733266550577
2.31485774  +0.7357239106731317   +0.7357239106731318
2.64555171  +0.4759473930370737   +0.4759473930370737
2.97624567  +0.1645945902807340   +0.1645945902807341
3.30693964  -0.1645945902807338   -0.1645945902807332
3.63763360  -0.4759473930370735   -0.4759473930370736
3.96832756  -0.7357239106731313   -0.7357239106731308
4.29902153  -0.9157733266550573   -0.9157733266550566
4.62971549  -0.9965844930066698   -0.9965844930066703
4.96040945  -0.9694002659393305   -0.9694002659393290
5.29110342  -0.8371664782625288   -0.8371664782625288
5.62179738  -0.6142127126896680   -0.6142127126896658
5.95249134  -0.3246994692046837   -0.3246994692046831
6.28318531  -0.0000000000000002   +0.0000000000000099
 
Last edited:
en.wikipedia.org

CORDIC - Wikipedia

en.wikipedia.org en.wikipedia.org

CORDIC (for "coordinate rotation digital computer"), also known as Volder's algorithm, or: Digit-by-digit method Circular CORDIC (Jack E. Volder),[1][2] Linear CORDIC, Hyperbolic CORDIC (John Stephen Walther),[3][4] and Generalized Hyperbolic CORDIC (GH CORDIC) (Yuanyong Luo et al.),[5][6] is a simple and efficient algorithm to calculate trigonometric functions, hyperbolic functions, square roots, multiplications, divisions, and exponentials and logarithms with arbitrary base, typically converging with one digit (or bit) per iteration. CORDIC is therefore also an example of digit-by-digit algorithms. CORDIC and closely related methods known as pseudo-multiplication and pseudo-division or factor combining are commonly used when no hardware multiplier is available (e.g. in simple microcontrollers and FPGAs), as the only operations it requires are additions, subtractions, bitshift and lookup tables. As such, they all belong to the class of shift-and-add algorithms. In computer science, CORDIC is often used to implement floating-point arithmetic when the target platform lacks hardware multiply for cost or space reasons.
Click to expand...
Code from the link. The code works only in the 1st quadrant 0-pi/2 or 90 degrees. The input x would gave to be processed to find the quadrant for x>pi/2 ro make it work 0-2*pi.


Code: 
from math import atan2, sqrt, sin, cos, radians
import math as ma
import array as ar
import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt
import warnings


ITERS = 100
theta_table = [atan2(1, 2**i) for i in range(ITERS)]

def compute_K(n):
    """
    Compute K(n) for n = ITERS. This could also be
    stored as an explicit constant if ITERS above is fixed.
    """
    k = 1.0
    for i in range(n):
        k *= 1 / sqrt(1 + 2 ** (-2 * i))
    return k

def CORDIC(alpha, n):
    K_n = compute_K(n)
    theta = 0.0
    x = 1.0
    y = 0.0
    P2i = 1  # This will be 2**(-i) in the loop below
    for arc_tangent in theta_table:
        sigma = +1 if theta < alpha else -1
        theta += sigma * arc_tangent
        x, y = x - sigma * y * P2i, sigma * P2i * x + y
        P2i /= 2
    return x * K_n, y * K_n

n = 20
xmax = ma.pi/2  # sin() 360 degrees
x = np.ndarray(shape=(n,1),dtype = np.float64)
y = np.ndarray(shape=(n),dtype = np.float64)  # interpolated sin
ys = np.ndarray(shape=(n),dtype = np.float64) # sin() libray function
yc = np.ndarray(shape=(n),dtype = np.float64) # sin() libray function
x = np.linspace(0,xmax,n)  # 2*pi raduans


for i in range(n):
     yc[i], y[i] = CORDIC(x[i], ITERS)
     #y[i] = sinx
     ys[i] = ma.sin(x[i])
     print("%2.8f  %+2.16f   %+2.16f" %(x[i],ys[i],y[i]))


[fig, ax] = plt.subplots()
ax.plot(x, y, linewidth=2.0,color='r')
ax.plot(x, ys, linewidth=2.0,color='k')
ax.plot(x, yc, linewidth=2.0,color='k')
ax.grid(color='k', linestyle='-', linewidth=2)
plt.show()

for 100 iterations

Code: 
x, sin(), CORDIC
0.00000000  +0.0000000000000000   +0.0000000000000000
0.08267349  +0.0825793454723323   +0.0825793454723323
0.16534698  +0.1645945902807339   +0.1645945902807339
0.24802047  +0.2454854871407991   +0.2454854871407992
0.33069396  +0.3246994692046835   +0.3246994692046835
0.41336745  +0.4016954246529694   +0.4016954246529693
0.49604095  +0.4759473930370735   +0.4759473930370735
0.57871444  +0.5469481581224268   +0.5469481581224267
0.66138793  +0.6142127126896678   +0.6142127126896681
0.74406142  +0.6772815716257410   +0.6772815716257410
0.82673491  +0.7357239106731316   +0.7357239106731317
0.90940840  +0.7891405093963936   +0.7891405093963936
0.99208189  +0.8371664782625285   +0.8371664782625287
1.07475538  +0.8794737512064891   +0.8794737512064885
1.15742887  +0.9157733266550574   +0.9157733266550580
1.24010236  +0.9458172417006346   +0.9458172417006353
1.32277585  +0.9694002659393304   +0.9694002659393305
1.40544935  +0.9863613034027223   +0.9863613034027227
1.48812284  +0.9965844930066698   +0.9965844930066698
1.57079633  +1.0000000000000000   +0.9999999999999999
 
The usual method of calculating square roots and the like is Newton's method:
\(x^n = a\)

\( \displaystyle{ x \gets \frac{1}{n} \left( (n - 1) x + \frac{a}{x^{n-1}} \right) } \)
  • Square root (n = 2) -- \( x \gets \frac{1}{2} (x + a/x) \)
  • Inverse square root (n = -2) -- \( x \gets \frac{1}{2} x (3 - a x^2) \)
  • Reciprocal (n = -1) -- \( x \gets x (2 - a x) \)
I'm mentioning reciprocals because that's what one uses to do division in hardware that lacks that capability. In practice, it's all done behind the scenes, either in a chip itself as microcode or in some run-time library.

For fast convergence in Newton's method, one needs a good initial guess. Here is a notable one for the inverse square root. It directly uses the floating-point representation of a number, taking minus one half of the number's exponent. One can do similar speedups for the other functions.

 Fast inverse square root and algorithm - John Carmack's Unusual Fast Inverse Square Root (Quake III) - Stack Overflow and Fast Inverse Square Root — A Quake III Algorithm - YouTube
 
To calculate inverse trigonometric functions, one can use what I call Archimedes's method and François Viète's method.

Archimedes (~287 BCE - ~212 BCE): repeatedly use a half-angle formula:

\( \displaystyle{ \tan \frac{x}{2} = \frac{\tan x}{1 + \sqrt{1 + \tan^2 x} } } \)

and then after n iterations do

\( \displaystyle{ x \gets 2^n \tan \frac{x}{2^n} } \)

Viète (1540 - 1603) or Vieta: repeatedly use another half-angle formula: \( \sin x = 2 \sin (x/2) \cos (x/2) \) After n iterations:

\( \displaystyle{ x \gets 2^n \sin \frac{x}{2^n} \ = \frac{\sin x}{ \prod_{k=1}^n \cos \frac{x}{2^k} } } \)

where \( \cos x = \sqrt{1 - \sin^2 x} \) and \( \displaystyle{ \cos \frac{x}{2} = \sqrt{ \frac{1 + \cos x}{2} } } \)

Both methods were used to calculate pi by their inventors.
 
From the links the CORDIC algorithm was driven by a need to imminent trig functions on a machine without a ha4dware multiplier. Shit and add, add, and subtract integers. Circa early 1950s. There us a block diagram in one of the links.

With the speed of modern micro-controllers a Taylor series can easily be used.

I call it successive approximation. The table is arc tangents at successively smaller values. An angle is summed adding and subtracting values as the value approaches the target angle. At each step the x and y values on a unit circle are calculated. On a unit circle the y coordinate is the sine.

An example f a real situation. A sensor s attached t a shaft that rotates 360 degrees. The sensor outputs a voltage from 0-12 volts. An 8 bit analog to digital converter is simulated to show the effects of quantization errors.

The algorithm works in the first quadrant 0-90 degrees, so for angles > 90 the quadrant has to be found.

Code: 
from math import atan2, sqrt, sin, cos, radians
import math as ma
import array as ar
import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt
import warnings

#--------------------------------------------------------------
def shaft_angle(volts, n):
    # sensor roration 0-360 degrees  0-12 volts
    alpha = 0.
    ang = 2*ma.pi*volts/12.  #angle in radiabs
    sgn = 1
    if(ang > 0. and ang <= ma.pi/2):   #1st quadrant
        sgn = 1.
        alpha = ang
    if(ang > ma.pi/2. and ang <= 2*ma.pi):  #2nd quadrant
        sgn = 1.
        alpha = ma.pi - ang ;
    if(ang > ma.pi and ang <= 1.5*ma.pi): #3rd quadrant
        sgn = -1.
        alpha = ang - ma.pi;
    if(ang > 1.5*ma.pi and ang <= 2*ma.pi): #4th quadrabt
        sgn = -1.
        alpha = 2*ma.pi - ang

    k_n = 1.0
    for i in range(n): k_n *= 1 / sqrt(1 + 2 ** (-2 * i))
    atan_table = [atan2(1, 2**i) for i in range(n)]
    theta = 0.0
    x = 1.0
    y = 0.0
    p2i = 1
    
    for i in range(n):
        sigma = +1 if theta < alpha else -1
        theta += sigma *atan_table[i]
        x, y = x - sigma * y * p2i, sigma * p2i * x + y
        p2i /= 2

    return sgn*y * k_n
#---------------------------------------------------------------------
def adc(nbits,vin):
    lsb = 12./pow(2,nbits)
    if vin < lsb: return 0
    vout = 0.
    n = pow(2,nbits)-1
    for i in range(n):
        vout += lsb
        if vout >= vin: break
    return vout
#----------------------------------------------------------------------
niterations = 50
ndata = 20
y = np.ndarray(shape=(ndata),dtype = np.float64)  # interpolated sin
volts = np.ndarray(shape=(ndata),dtype = np.float64)
volts = np.linspace(0,12,ndata)  #0-12 volts
vadc = np.ndarray(shape=(ndata),dtype = np.float64)

print("continous volage, digitized voltage, sine of shaft angle")
for i in range(ndata):
    vadc[i] = adc(8,volts[i])
    #print(vadc)
    y[i] = shaft_angle(vadc[i], niterations)
    print("% 2.8f   %2.8f  %+2.16f" %(volts[i],vadc[i],y[i]))


[fig, ax] = plt.subplots()
ax.plot(vadc, y, linewidth=2.0,color='k')
ax.grid(color='k', linestyle='-', linewidth=2)
plt.show()



Code: 
continous volage, digitized voltage, sine of shaft angle
 0.00000000   0.00000000  -0.0000000000000002
 0.63157895   0.65625000  +0.3368898533922198
 1.26315789   1.26562500  +0.6152315905806262
 1.89473684   1.92187500  +0.8448535652497073
 2.52631579   2.53125000  +0.9700312531945439
 3.15789474   3.18750000  +0.9951847266721973
 3.78947368   3.79687500  +0.9142097557035308
 4.42105263   4.45312500  +0.7242470829514663
 5.05263158   5.06250000  +0.4713967368259979
 5.68421053   5.71875000  +0.1467304744553612
 6.31578947   6.32812500  -0.1709618887602995
 6.94736842   6.98437500  -0.4928981922297844
 7.57894737   7.59375000  -0.7409511253549603
 8.21052632   8.25000000  -0.9238795325112880
 8.84210526   8.85937500  -0.9972904566786909
 9.47368421   9.51562500  -0.9637760657954396
 10.10526316   10.12500000  -0.8314696123025450
 10.73684211   10.78125000  -0.5956993044924332
 11.36842105   11.39062500  -0.3136817403988908
 12.00000000   11.95312500  -0.0245412285229123
 
You must log in or register to reply here.
Back
Top Bottom