 ## Functions Implementation: Square Root, Logarithm, Sine, Cosine, Arctangent. C Code and Octave Script

### 3.1. Square Root. Taylor/Maclaurin series. Chebyshev Approximation

##### Taylor series: ##### Using the case: ##### Then (4) will corresponds to the form of Maclaurin series: ##### Let’s choose f(x) as: ##### Using (6) and (7): ##### Let’s take only the first 4 terms of the series: ##### Using the Octave script FunctionsRealization.m (see the array sqrt_relative_mistake2) is calculated the relative error: Δsqrt(1+x)/sqrt (1+x) by the formula (9) in the range x = -0.3…0.4. See the Figure 3-1. In this case the maximum relative error reaches:  Figure 3-1: Relative mistake: Δsqrt(1+x)/sqrt(1+x) of the square root calculation with 4 elements of the Maclaurin Series

##### Let’s rewrite (9) in the form: ##### To find the coefficients, the three equations shall be made. We will use the Chebyshev nodes to reduce the Runge phenomenon (fluctuations are minimized): ##### In our case: ##### Using (13), (14), (15), we make up 3 equations for finding the coefficients: ##### Solving (16): ##### Then: #####  Figure 3-2: Relative mistake: Δsqrt(1+x)/sqrt(1+x) of the square root calculation with 4 elements of the Chebyshev Approximation

### 3.2. Square Root. Newton’s Method.

##### The x shall be found in (21): ##### We can use for it the iterative Newton formula: ##### The formula is easy to prove. Let there is an estimate of the square root: Suppose the error of this estimate is Δx. Then: ##### Opening brackets: ##### If Δx is small, then we can neglect the Δx^2. Then the Δx value: ##### Then: ##### Newton Method. Calculation example: ##### Let’s take the x initial value: ##### Iteration 1: ##### Iteration 2: ##### Iteration 3: ### 4. Logarithm

##### To calculate the logarithm, we use Maclaurin series (6): ##### Reasoning in the same way as in paragraph 3.1, we will leave the first 3 terms of the series: ##### Any positive number: a > 0 can be represented as: ##### Using the Octave script FunctionsRealization.m (see the ln_absolute_mistake) is calculated the absolute error: Δln (1+x) by (28) in the range x = -0.32…0.36. See the Figure 4-1. In this case the maximum absolute error reaches:  Figure 4-1: Absolute mistake: Δln(1+x) of the natural logarithm calculation with 3 elements of the Maclaurin Series

##### As in paragraph 3.1 to increase the accuracy of the calculations we will use the approximation Chebyshev nodes: ##### from (14). Let’s rewrite (28) in the form: ##### Using (32) and (14), we make up 3 equations for finding the coefficients: ##### Then: ##### Using the Octave script FunctionsRealization.m (see the ln_absolute_mistake1) is calculated the absolute error: Δln(1+x) by (35) in the range x = -0.32…0.36. See the Figure 4-2. In this case, the maximum absolute error reaches: #####  Figure 4-2: Absolute mistake: Δln(1+x) of the natural logarithm calculation with 3 elements of the Chebyshev Approximation

• ##### In the future we will need a decimal logarithm. To move from the natural logarithm is used the basic logarithmic formula: ### 5.1 Sine

##### To calculate the sine, we first use the Maclaurin series (6): ##### Reasoning similarly, as in paragraph 3.1, we will leave the first 4 terms of the series: ##### Using the Octave script FunctionsRealization.m (see the sin_absolute_mistake) is calculated the absolute error: Δsin(x) by (41) in the range x = -π…π. See the Figure 5-1. In this case the maximum absolute error reaches:  Figure 5-1: Absolute mistake: Δsin(x) of the sine calculation with 4 elements of the Maclaurin Series

##### As in paragraph 3.1 to increase the accuracy of calculations we will use the Chebyshev approximation. Rewrite the (41) in the form: ##### Taking advantage of the odd sine: sin (- x) = – sin(x) is chosen 9 = 4*2+1 Chebyshev nodes and using only the first four nodes to determine the coefficients a1, …, a4 : ##### Using (43) and (44), we make up 4 equations for finding the coefficients of a: ##### Solving (45): ##### Then: ##### Using Octave script FunctionsRealization.m (see the sin_absolute_mistake1) is calculated the absolute error: Δsin(x) by (47) in the range x = -π…π. See the Figure 5-2. In this case the maximum absolute error reaches: ##### To reduce the number of multiplications in (47), another entry is used:  Figure 5-2: Absolute mistake: Δsin(x) of the sine calculation with 4 elements of the Chebyshev Approximation

### 5.2 Cosine

##### To calculate the cosine is used the formula: ### 6. Arctangent

##### In this case β = 0 … π/4, which corresponds to (51). Then the tan(α) can be calculated: ##### Rewrite (53) in the form: ##### Using (52), (54) and (51) we can write for x > 1: ##### Reasoning similarly for x < -1: ##### Reasoning in the same way as in the previous paragraphs, we will leave the first 4 terms of the series (51): ##### Using Octave script FunctionsRealization. m ( see the sin_absolute_mistake) is calculated the absolute error: Δarctan(x) by (57) in the range x = -1…1. See the Figure 6-1. In this case the maximum absolute error reaches:  Figure 6-1: Absolute mistake: Δatan(x) of the sine calculation with 4 elements of the Maclaurin Series

##### To increase the accuracy of calculations we will use the Chebyshev approximation. Rewrite (57) in the form: ##### Taking advantage of the odd arctan: arctan (- x) = – arctan(x) is chosen 9 = 4*2+1 Chebyshev nodes and using only the first four nodes to determine the coefficients a1, …, a4 : ##### Using (59) and (60) we make up 4 equations for finding the coefficients: ##### Solving (61): ##### Then: ##### Using Octave script FunctionsRealization.m (see the atan_absolute_mistake1) is calculated the absolute error: Δatan(x) by (63) in the range x = -1…1. See the Figure 6-2. In this case, the maximum absolute error reaches: ##### Comparing (58) and (64) the error has decreased by more than 590 times! Figure 6-2: Absolute mistake: Δatan(x) of the sine calculation with 4 elements of the Chebyshev Approximation

##### To reduce the number of multiplications in (63), another entry is used: ### 8.1 square_root.c/h

##### square_root.c/h is supported square root functions

/*
Square Root
Float Arithmetic
Polynomial sqrt(1+x) = 1 + x (0.50016487 + x (-0.12822807 + 0.059227649 x))
0.7 < 1 + x <= 1.4
-0.3 < x <= 0.4
Input
float input – input value
Return
float – sqrt(input) Note: if input number is 0.0 or negative then return value is 0.0
*/
float sgrt_polynomial_float(float input);

/*
Square Root
Fixed Point Arithmetic
Polynomial sqrt(1+x) = 1 + x (0.50016487 + x (-0.12822807 + 0.059227649 x))
0.7 < 1 + x <= 1.4
-0.3 < x <= 0.4
Input
uint16 input – input value: 0…65535
Return
uint16 – sqrt(input) in format: b15*2^7 + b14*2^6 + … + b8*2^0 + b7*2^-1 + b6*2^-2 + … + b0*2^-8
Note: Result = ReturnValue/OUT_SQRT_FACTOR, where OUT_SQRT_FACTOR = 256
*/
uint16 sgrt_polynomial_int(uint16 input);

/*
Square Root
Float Arithmetic
Newton Method x(n+1) = (x(n)+ a/x(n))/2
Input
float input – input value
Return
float – sqrt(input) Note: if input number is 0.0 or negative then return value is 0.0
*/
float sgrt_Newton_float(float input);

/*
Square Root
Fixed Point Arithmetic
Newton Method x(n+1) = (x(n)+ a/x(n))/2
Input
uint16 input – input value: 0…65535
Return
uint16 – sqrt(input) in format: b15*2^7 + b14*2^6 + … + b8*2^0 + b7*2^-1 + b6*2^-2 + … + b0*2^-8
Note: Result = ReturnValue/OUT_SQRT_FACTOR, where OUT_SQRT_FACTOR = 256
*/
uint16 sgrt_Newton_int(uint16 input);

##### The square root functions were tested in main() using the test arrays: input_array_float[ ] and input_array_int[ ]. The test results are shown in Figure 8-1, Figure 8-2, Figure 8-3, Figure 8-4 Figure 8-1: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_polynomial_float() for x=1…100. Max relative mistake is 0.000163 Figure 8-2: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_polynomial_int() for x= 65535, 1, 2, … 100. Max relative mistake is 0.001302 Figure 8-3: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_Newton_float() for x=1…100. Max relative mistake is 0.0000000946 Figure 8-4: Relative mistake: Δsqrt(x)/sqrt(x) of the sgrt_Newton_int() for x= 65535, 1, 2, … 100. Max relative mistake is 0.000913

### 8.2 logarithm.c/h

##### logarithm.c/h is supported logarithm functions

/*
Square Root
Float Arithmetic
Polynomial log10(1+x) = log10(e) * ln(1+x) = 0.43429449 ( x (1.0004359 + x (-0.52194273 + 0.33578411 x)))
0.68 < 1 + x <= 1.36
-0.32 < x <= 0.36
Input
float input – input value
Return
float – log10(input) Note: if input number is 0.0 or negative then return value is 0.0
*/
float logarithm10_float(float input);

/*
Decimal Logarithm
Fixed Point Arithmetic
Polynomial log10(1+x) = log10(e) * ln(1+x) = 0.43429449 ( x (1.0004359 + x (- 0.52194273 + 0.33578411 x)))
0.68 < 1 + x <= 1.36
-0.32 < x <= 0.36
Input
uint16 input – input value: 1…65535 => log10(1)=0 … log10(65535)=4.8164735
Return
sint16 – log10(input) in format: b14*2^2 + b13*2^1 + b12*2^0 + b11*2^-1 + … + b0*2^-12
Note: Result = ReturnValue/OUT_LOGARITHM_FACTOR, where OUT_LOGARITHM_FACTOR = 4096
*/
sint16 logarithm10_int(uint16 input);

##### The logarithm10 functions were tested in main() using the test arrays: input_array_float[ ] and input_array_int[ ]. The test results are shown in Figure 8-5, Figure 8-6 Figure 8-5: Absolute mistake: Δlogarithm of the logarithm10_float() for x= 0.1, 1, 2, …, 99, 100. Max absolute mistake is 0.000365 Figure 8-6: Absolute mistake: Δlogarithm of the logarithm10_int() for x= 65535, 1, 2, …, 99, 100. Max absolute mistake is 0.000539

### 8.3 sine.c/h

##### Sine.c/h is supported sine and cosine functions

/*
Sine
Float Arithmetic
Polynomial sin(x) = x (0.9992497 + x^2 (-0.16564582 + x^2 (7.9537760e-3 – 1.4482903e-4 x^2)))
-pi <= x <= pi
Input
float input – input value
Return
float – sin(input) Note: if input > pi or input < -pi then return value is 0.0
*/
float sin_float(float input);

/*
Cosine
Float Arithmetic
cos(x) = sin(x + pi/2)
-pi <= x <= pi
Input
float input – input value
Return
float – cos(input) Note: if input > pi or input < -pi then return value is 0.0
*/
float cos_float(float input);

/*
Sine
Fixed Point Arithmetic
Polynomial sin(x) = x (0.9992497 + x^2 (-0.16564582 + x^2 (7.9537760e-3 – 1.4482903e-4 x^2)))
-pi <= x <= pi
Input
sint16 input – input value: -PI_INT_VALUE=-25736 … PI_INT_VALUE=25736
input format 3Q13: b14*2^1 + b13*2^0 + b12*2^-1 + b11*2^-2 + … + b0*2^-13 and b15-sign
Return
sint16 – sin(input) in format 1Q15: b14*2^-1 + b13*2^-2 + b12*2^-3 + b11*2^-4 + … + b0*2^-15 and b15-sign
Note: Result = ReturnValue/OUT_SINE_FACTOR, where OUT_SINE_FACTOR = 32768
*/
sint16 sin_int(sint16 input);

/*
Cosine
Fixed Point Arithmetic
cos(x) = sin(x + pi/2)
-pi <= x <= pi
Input
sint16 input – input value: -PI_INT_VALUE=-25736 … PI_INT_VALUE=25736
input format 3Q13: b14*2^1 + b13*2^0 + b12*2^-1 + b11*2^-2 + … + b0*2^-13 and b15-sign
Return
sint16 – cos(input) in format 1Q15: b14*2^-1 + b13*2^-2 + b12*2^-3 + b11*2^-4 + … + b0*2^-15 and b15-sign
Note: Result = ReturnValue/OUT_SINE_FACTOR, where OUT_SINE_FACTOR = 32768
*/
sint16 cos_int(sint16 input);

##### The sine / cosine functions were tested in main() using the test arrays: input_array_float[ ] and input_array_int[ ]. The test results are shown in Figure 8-7, Figure 8-8 Figure 8-7: Absolute mistake: Δsine of the sin_float() for x= -π … +π. Max absolute mistake is 2.587e-04 Figure 8-8: Absolute mistake: Δsine of the sin_int() for x= -π … +π. Max absolute mistake is 2.694e-02

### 8.3 atan.c/h

##### atan.c/h is supported arctangent function

/*
Arctangent
Float Arithmetic
Polynomial atan(x) = x (0.99904621 + x^2 (-0.3198728 + x^2 (0.14370242 – 0.037537854 x^2)))
-1 <= x <= 1
Polynomial atan(x) = +/- pi/2 – x^-1 (0.99904621 + x^-2 (-0.3198728 + x^-2 (0.14370242 – 0.037537854 x^-2)))
|x| > 1
Input
float input – input value
Return
float – atan(input)
*/
float atan_float(float input);

##### The atan function was tested in main() using the test arrays: input_array_float 1/2/3[ ]. The test results are shown in Figure 8-9, Figure 8-10, Figure 8-11 Figure 8-9: Absolute mistake: Δatan of the atan_float() for x = -100 … -1. Max absolute mistake is 1.036e-4 Figure 8-10: Absolute mistake: Δatan of the atan_float() for x = -1 … 1. Max absolute mistake is 1.032e-04 Figure 8-11: Absolute mistake: Δatan of the atan_float() for x = 1 … 100. Max absolute mistake is 1.036e-4