Author Topic: (Library) Arbitrary precision ArcTan, ArcTan2 and ArcSin for Arduino  (Read 16643 times)

Bob101

  • Jr. Member
  • **
  • Posts: 71
    • View Profile
UPDATE 17/10/2012: RECODED LIBRARY, PLEASE SEE LATER POST FOR INSTRUCTIONS
UPDATE 17/07/2012: ADDED FINISH() TO LIBRARY

== Introduction ==
Thanks to the great work of Nick Gammon in porting the big number library to the Arduino and Gabriel in implementing arbitrary precision Sin and Cos functions using big numbers it was possible to perform most of the sun tracking calculations to arbitrary precision on the Arduino.

Unfortunately, the much needed ArcTan, ArcTan2 and ArcSin functions did not work well as implemented (leading the Arduino to crash in most cases) and it was necessary to resort to using the built in low accuracy (typically 2 digit) atan, atan2 and asin functions. Though these functions are accurate enough for nearly all applications they are in principle a hard limit of the accuracy of the calculations.

I believe I have identified the reason for the failure of ArcTan to return and I have altered the implementation to avoid this. I have wrapped all of Gabriel's BigNumberMathFunctions.ino into a library BigNumberMath and have added the newly implemented ArcTan, ArcTan2 and ArcSin along with other functions required to evaluate them efficiently. The new library is attached.

== Instructions ==
To Install (on Arduino 1.0.1):
1) Ensure you have the BigNumber library installed (http://arduino.cc/forum/index.php/topic,85692.0.html)
2) Unzip the BigNumberMath folder into the libraries directory in your Arduino IDE

To Use (on Arduino 1.0.1):
1) Include the BigNumber and BigNumberMath libraries, at the top of your sketch write the lines:
Code: [Select]
#include <BigNumber.h>
#include <BigNumberMath.h>
2) In the void setup() function of your sketch initialise the BigNumberMath library, this automatically initialises the BigNumber library with the correct scale, you are strongly advised to not initialise BigNumber separately. Example:
Code: [Select]
BigNumberMath::begin(30);
initialises the BigNumberMath library and the BigNumber library with a scale of 30. This scale seems to provide a good balance of accuracy and speed on the Arduino Uno. Increase this number to work to higher precision (at the cost of speed), note that if this number is too large the Arduino will not be able to fit the required big numbers in memory and will crash without warning or give incorrect answers.
3) Use the static constants and functions of the library by referencing them directly. Example:
Code: [Select]
BigNumberMath::arctangent(BigNumber("5"),10);
would return atan(5) to the accuracy obtained by summing the first 10 terms of the expansion. This number of terms typically provides good accuracy and speed on the Arduino Uno.

Example of the library in use:
Code: [Select]
#include <BigNumber.h>
#include <BigNumberMath.h>

void setup() {
  Serial.begin (9600);
  BigNumberMath::begin(30);
  printBigNum(BigNumberMath::arctangent(BigNumber("5"),10));
  BigNumberMath::finish(); // v2 addition
}

void loop() {
  // null
}

// function to display a big number and free it afterwards
void printBigNum (BigNumber n)
{
  char * s = n.toString ();
  Serial.println (s);
  free (s);
}

This will output over serial:
Code: [Select]
1.373400766945015874925235233734

For comparison the actual value of ArcTan[5]=1.373400766945015860861271926445.

== List of features ==
CONSTANTS:
Code: [Select]
BigNumber pi;
BigNumber zero;
BigNumber one;
BigNumber two;
BigNumber oneEighty;
BigNumber recipRootThree;

TRIGONOMETRIC FUNCTIONS
Code: [Select]
BigNumber cosine(BigNumber z, int p);
BigNumber sine(BigNumber z, int p);
BigNumber arcsine(BigNumber z, int p); // new
BigNumber arctan2(BigNumber y, BigNumber x, int p); // new
BigNumber arctangent(BigNumber z, int nmax); // new
BigNumber within2pi(BigNumber angle);
BigNumber to_BigRad(BigNumber angle);
BigNumber to_BigDeg(BigNumber angle);

COMBINATORIAL FUNCTIONS
Code: [Select]
BigNumber fallingFactorial2(int x, int n); // new - definition: x!!/n!!
BigNumber factorial2(int x); // new - definition: x!!
BigNumber factorial(int x);

BASIC FUNCTIONS
Code: [Select]
BigNumber raiseToPower(BigNumber base, long exponent);

I have also included BigNumberMathTest an Arduino 1.0.1 compatible sketch that prints out tables for Cos, Sin, ArcSin, ArcTan plus special checks (including of ArcTan2). Note that the "-100" displayed in "F[-100]" etc... is not the argument with which the function has been evaluated, it was used to automate checking of the results in Mathematica. All of the functions have been checked against Mathematica for a small range (201 points that span a "reasonable" part of the domain), these Mathematica checking files have been included in the BigNumberMathTest/checks folder for those who are interested.

== Discussion ==

This library is designed to be a drop-in replacement for Gabriel's BigNumberMathFunctions.ino. My hope is that this library can be included in Gabriel's Sun Tracking program to allow full arbitrary precision sun tracking.  All that is required is to replace every occurence of function with BigNumberMath::function and ensure that a sensible number of terms in the expansion is chosen. To begin with, I would recommend 20 terms for Cos and Sin and 10 terms for ArcTan, ArcTan2 and ArcSin on the Arduino Uno.

On the Arduino Uno for the few cases tested these new implementations typically return >7 digits of precision for a big number scale = 30 with 10 terms of the expansion, this takes less than 3 seconds to evaluate in most cases.

== Disclaimer ==
I am in no way claiming that this implementation is optimal or even accurate and I am very happy for any improvements to be made to this library. This is the first revision of my first attempt at arbitrary precision calculations and my first c++ and Arduino library! There are still lots of things to be fixed. Please feel free to get involved, there are lots of people out there far more experienced than me. I do hope that any changes you do make will be published back here so that everybody can benefit from them. Please report any bugs/issues and possible improvements!

I would like to make clear that this work is heavily influenced by the implementation by Gabriel, many of the functions are directly copied. Euler's approximation is still used in the approximation of ArcTan, however, some numerical tricks, implemented by me, have been used to express this in a more computationally friendly way. I will post more on this shortly. I think I found a new (probably not) efficient way to compute the coefficients of the Euler expansion using lots of numerically easy operations (floor, power), this trick is used in the ArcTan function implemented in this library.

Thanks to Gabriel for his dedication to this cause! I am building a solar tracker currently and hope to be more involved in this community. Looking forward to hearing from you guys!
« Last Edit: October 17, 2012, 02:23:52 PM by Bob101 »


Bob101

  • Jr. Member
  • **
  • Posts: 71
    • View Profile
Re: (Library) Arbitrary precision ArcTan, ArcTan2 and ArcSin for Arduino
« Reply #1 on: July 16, 2012, 12:42:04 PM »
To help understand precisely why the old implementation of ArcTan often crashed the Arduino and to explain some of the mathematics behind the new ArcTan function I have attached a short (2 page) summary of what changes were implemented and why.

Hope this helps!

« Last Edit: July 17, 2012, 01:55:45 AM by Bob101 »


Gabriel

  • Administrator
  • Hero Member
  • *****
  • Posts: 655
    • View Profile
Re: (Library) Arbitrary precision ArcTan, ArcTan2 and ArcSin for Arduino
« Reply #2 on: July 17, 2012, 07:09:34 AM »
This looks awesome!
Thanks for the arctan fix/addition. I spent a lot of time trying to implement it, but I never was able to figure out a way around the memory overflow issue inherent to the large factorials in Euler's Approximation.

I will have to add these improved trig functions to the Arduino Mega version of the Sun Tracking / Heliostat program once I finish updating some other things first. If it works out, we should have the full accuracy of the algorithm, which gives the sun's position to within 0.01 degrees if I recall correctly. Not bad for an 8-bit microcontroller!

Thanks for your help! I hope to hear more from you. It's always nice to have another programmer around to watch my back.

Bob101

  • Jr. Member
  • **
  • Posts: 71
    • View Profile
Re: (Library) Arbitrary precision ArcTan, ArcTan2 and ArcSin for Arduino
« Reply #3 on: July 17, 2012, 08:04:33 AM »
Thanks!

Don't worry, I plan to stick around. As I said I am involved in building a hobby solar tracker and fully intend to make the best use of your considerable effort in implementing the tracking algorithms! Your program is really looking great.

One Note (before you start implementing):
- I have added a new function "finish()" which explicitly destructs the BigNumberMath object freeing all used memory. I will upload this version later.

I have checked the BigNumberMath library for memory leaks on return from each of the functions (using the Memory free library) and everything is looking great, no memory appears to be leaking in my test cases.

Bob101

  • Jr. Member
  • **
  • Posts: 71
    • View Profile
Re: (Library) Arbitrary precision ArcTan, ArcTan2 and ArcSin for Arduino
« Reply #4 on: July 17, 2012, 12:39:59 PM »
Updated library now attached to first post.


Bob101

  • Jr. Member
  • **
  • Posts: 71
    • View Profile
Re: (Library) Arbitrary precision ArcTan, ArcTan2 and ArcSin for Arduino
« Reply #5 on: September 18, 2012, 02:29:28 PM »
Not sure if this will help with the use of this library but I can see a few ways to improve it. This will slightly increase the memory footprint but should save the library repeatedly declaring the same BigNumbers it also allows us to ensure they are properly destructed.

Specifically:
- Add declarations for BigNumber(-1), BigNumber(6) and destruct them properly.
- Replace lines like BigNumber sum=0 -> BigNumber sum=zero and if(y.isZero() && x<BigNumber(0)) -> if(y.isZero() && x<zero) and sum += pi/BigNumber(6) -> sum += pi/six etc...

The library will then not declare any magic numbers, it will only declare BigNumbers inside functions if they are to be returned or when they must be generated at runtime eg... output=output*BigNumber(i)

I will implement these changes soon but I seriously doubt they will fix the issues we are experiencing using this library in the main program.

Bob101

  • Jr. Member
  • **
  • Posts: 71
    • View Profile
Re: (Library) Arbitrary precision ArcTan, ArcTan2 and ArcSin for Arduino
« Reply #6 on: October 17, 2012, 02:22:01 PM »
== Introduction ==
To improve performance and ease usage of the library it has now been completely recoded from scratch.

This version of the library is based on libmath.b from bc-1.06 by Philip A. Nelson. All implemented functions use a range reduction and then a taylor expansion. Functions are expanded until the next term in the taylor expansion is smaller than the current BigNumber scale meaning that in most cases numbers returned can be trusted in all except the last digit. During evaluation the BigNumber scale is dynamically adjusted to maintain both speed and accuracy.

This port required largely trivial translation of the b code for cosine, sine and arctangent. I have also used these functions as before to provide arcsine and arctan2. In addition to this I altered the range reduction for arctan to guarantee range reduction without a loop, this gives a performance increase.

== Instructions ==
To Install (on Arduino 1.0.1):
1) Ensure you have the BigNumber library installed (http://arduino.cc/forum/index.php/topic,85692.0.html) (NOTE: a new version is available since the last release of this library, it fixes a memory leak spotted in sqrt)
2) Unzip the BigNumberMath folder into the libraries directory in your Arduino IDE

To Use (on Arduino 1.0.1):
1) Include the BigNumber and BigNumberMath libraries, at the top of your sketch write the lines:
Code: [Select]
#include <BigNumber.h>
#include <BigNumberMath.h>
2) In the void setup() function of your sketch initialise the BigNumberMath library, this automatically initialises the BigNumber library with the correct scale, you are strongly advised to not initialise BigNumber separately. Example:
Code: [Select]
BigNumberMath::begin(10);initialises the BigNumberMath library and the BigNumber library with a scale of 10. This scale seems to provide a good balance of accuracy and speed on the Arduino Uno. Increase this number to work to higher precision (at the cost of speed), note that if this number is too large the Arduino will not be able to fit the required big numbers in memory and will crash without warning or give incorrect answers. The scale should not be increased above 30 due to the limited accuracy of hardcoded constants (eg... pi,pi/2,pi/4) included in the library.
3) Use the static constants and functions of the library by referencing them directly. Example:
Code: [Select]
BigNumberMath::arctangent(BigNumber("5"));would return atan(5) to the highest accuracy possible at the current BigNumber scale (the last digit should not be trusted due to the rounding convention of BigNumbers).
NOTE: If you use the non-integer constants you should divide them by one to reduce them to the current scale.
Example:
Code: [Select]
BigNumber myPI=BigNumberMath::pi/BigNumberMath::one; // reduce constant pi to the current BigNumber scale
Example of the library in use:
Code: [Select]
#include <BigNumber.h>
#include <BigNumberMath.h>

void setup() {
  Serial.begin (9600);
  BigNumberMath::begin(30);
  printBigNum(BigNumberMath::arctangent(BigNumber("5"))); // v3 change
  BigNumberMath::finish(); // v2 addition
}

void loop() {
  // null
}

// function to display a big number and free it afterwards
void printBigNum (BigNumber n)
{
  char * s = n.toString ();
  Serial.println (s);
  free (s);
}

This will output over serial:
Code: [Select]
1.373400766945015860861271926444
For comparison the actual value of ArcTan[5]=1.373400766945015860861271926445.

== List of features ==
CONSTANTS:
Code: [Select]
BigNumber pi;
BigNumber piOverOneEighty;
BigNumber piOverTwo;
BigNumber piOverFour;
BigNumber piOverSix;
BigNumber rootThree;
BigNumber twoPlusRootThree;
BigNumber twoMinusRootThree;
BigNumber recipRootThree;
BigNumber minusOne;
BigNumber zero;
BigNumber one;
BigNumber two;
BigNumber three;

TRIGONOMETRIC FUNCTIONS:
Code: [Select]
BigNumber cosine(BigNumber z);
BigNumber sine(BigNumber z);
BigNumber arcsine(BigNumber z);
BigNumber arctan2(BigNumber y, BigNumber x);
BigNumber arctangent(BigNumber z);
BigNumber to_BigRad(BigNumber angle);
BigNumber to_BigDeg(BigNumber angle);

BASIC FUNCTIONS:
Code: [Select]
BigNumber raiseToPower(BigNumber base, long exponent);
LIBRARY FUNCTIONS:
Code: [Select]
void setScale(const int scale = 0); // dynamically alters the current BigNumber scale
int getScale(); // returns the current BigNumber scale

I have also included BigNumberMathTest an Arduino 1.0.1 compatible sketch that prints out tables for Cos, Sin, ArcSin, ArcTan plus special checks (including of ArcTan2). Note that the "-100" displayed in "F[-100]" etc... is not the argument with which the function has been evaluated, it was used to automate checking of the results in Mathematica (actual argument is this 1/100th of this value). All of the functions have been checked against Mathematica for a small range (201 points that span a "reasonable" part of the domain), these Mathematica checking files have been included in the BigNumberMathTest/checks folder for those who are interested. NOTE: All checks now completely agree with Mathematica as there is no early truncation of the Taylor series in the new implemenation.

== Discussion ==
The most important changes in this library are:
1) A precision should no longer be specified.
2) The functions are much more acccurate so a smaller scale may be used (eg... 10 is now suitable where 30 was used previously).

Library attached to first post of this thread.

Thanks!

Bob101

  • Jr. Member
  • **
  • Posts: 71
    • View Profile
Re: (Library) Arbitrary precision ArcTan, ArcTan2 and ArcSin for Arduino
« Reply #7 on: October 17, 2012, 02:48:25 PM »
Since the biggest change I made was in the range reduction... here is some explanation.

In libmath.b the range reduction was implemented as:
Code: [Select]
f=0;
while (x > .2) {
    f += 1;
    x = (x-.2) / (1+x*.2);
  }

This works, but, for example, if we arbitrarily pick x=10 the loop will run as follows:
Code: [Select]
f=0; x=10
f=1; x=3.26667
f=2; x=1.85484
f=3; x=1.20706
f=4; x=0.811221
f=5; x=0.525897
f=6; x=0.294881
f=7; x=0.0895974
We have to range reduce 7 times before we can use our Taylor expansion! This means 7 multiplications, 14 additions, 7 divisions, 7 subtractions and then we can finally use our Taylor expansion.

I studied repeated application of this range reduction and saw no clear way to generate the "nth" application of the reduction without just calculating it. Maybe somebody else can show me how? NOTE: The next post examines this issue and shows how to do this in general.

Instead I considered using the range reduction function f(x)=(x-1/Sqrt(3))/(1+x*1/Sqrt(3)). On the surface this looks very similar. We have just replaced .2 with 1/Sqrt(3), not a big deal right? Wrong. Now we can quickly see a pattern emerging from repeated application of this function. In fact, if we apply this function 6 times we will get back to where we started! Exercise for the reader: show this, explain why this happens.

This leads us to the question: can we utilise this to perform an efficient range reduction? The answer is: Yes! Since there are only 5 possible results no matter how many times we apply this range reduction we can simply graph them and find in which domain we should use them.

Eye Candy:


The above graph shows the modulus of the 6 functions generated by repeated application of f(x). We can see that if we pick the right function we are guaranteed to reduce the modulus of our argument to < 2-Sqrt(3). Our Taylor expansion is valid here... perfect! Working out which function we should pick for each argument we arrive at the code:
Code: [Select]
    if (z > twoPlusRootThree) {
        f = three;
    z = minusOne/z;
    }
    else if (z > one) {
        f = two;
    z = (z-rootThree)/(one+rootThree*z);
    }
    else if (z > twoMinusRootThree) {
        f = one;
    z = (three*z-rootThree)/(three+rootThree*z);
    }
    else { // else no range reduction necessary
        f = zero;
    }
Now in the worst case we have to do 2 multiplications, 1 addition, 1 division, 1 subtraction. In most cases we can be lazy and just do 1 division.

Harrah!
« Last Edit: October 19, 2012, 10:43:41 AM by Bob101 »

Bob101

  • Jr. Member
  • **
  • Posts: 71
    • View Profile
Re: (Library) Arbitrary precision ArcTan, ArcTan2 and ArcSin for Arduino
« Reply #8 on: October 18, 2012, 12:13:11 PM »
I spoke to a friend about the range reduction algorithm and he showed me a very neat way to study repeated application of these functions in general. Take a deep breath... things are about to get real*!

*complex

Note that the original range reduction function:
f(x)=(x-1/5)/(1+x/5) = (5x-1)/(5+x),
has the form of a Möbius transformation:
m(x)=(ax+b)/(cx+d),
where ad-bc ≠ 0.

This can be represented as a projective matrix:
A = (a  b)
      (c  d),
where the coefficients of the nth application of the transformation correspond to the values of A^n. That is, the composition of the function can be mapped to matrix multiplication.

Unfortunately, it is not trivial to compute the nth power of A. To try to rectify this we can diagonalize A. For the particular function being studied it has eigenvalues (5+i) and (5-i) with corresponding eigenvectors (1,-i)^T and (1,i)^T (where "^T" means transpose). Thus it can be diagonalized by the matrix:
P=(1  1)
    (-i   i)

A' = P^(-1) A P = (5+i    0  )
                            (  0   5-i)

And now A^n can be computed easily from: A^n=P (A')^n P^(-1), where (A')^n has only diagonal elements (5+i)^n and (5-i)^n.

Multiplying out P (A')^n P^(-1) and expressing (5+i) and (5-i) in polar form we arrive at the following closed form solution:
A^n = ( Sqrt[26]^n*Cos(n*ArcTan[1/5])     -Sqrt[26]^n*Sin(n*ArcTan[1/5]) )
           ( Sqrt[26]^n*Sin(n*ArcTan[1/5])       Sqrt[26]^n*Cos(n*ArcTan[1/5]) )

Which means, the range reduction could instead read:
f^n(x) = (ax+b)/(cx+d)
where a=d=Sqrt[26]^n*Cos(n*ArcTan[1/5]) and b=-c=Sqrt[26]^n*Sin(n*ArcTan[1/5]).

I doubt this would be that efficient as it requires computing both Cos and Sin (or Sin and then using a trick to get Cos) just to range reduce ArcTan! Still... I thought this was very interesting and actually I was not expecting it!
« Last Edit: October 18, 2012, 12:44:58 PM by Bob101 »