Free cookie consent management tool by TermsFeed Policy Generator

source: stable/HeuristicLab.ExtLibs/HeuristicLab.NativeInterpreter/0.1/NativeInterpreter-0.1/lib/vdt/atan2.h @ 18091

Last change on this file since 18091 was 16269, checked in by bburlacu, 6 years ago

#2958: Add C++ source code

File size: 3.6 KB
Line 
1/*
2 * atan2.h
3 * The basic idea is to exploit Pade polynomials.
4 * A lot of ideas were inspired by the cephes math library (by Stephen L. Moshier
5 * moshier@na-net.ornl.gov) as well as actual code.
6 * The Cephes library can be found here:  http://www.netlib.org/cephes/
7 *
8 *  Created on: Sept 20, 2012
9 *      Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente
10 */
11
12/*
13 * VDT is free software: you can redistribute it and/or modify
14 * it under the terms of the GNU Lesser Public License as published by
15 * the Free Software Foundation, either version 3 of the License, or
16 * (at your option) any later version.
17 *
18 * This program is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
21 * GNU Lesser Public License for more details.
22 *
23 * You should have received a copy of the GNU Lesser Public License
24 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
25 */
26
27#ifndef ATAN2_H_
28#define ATAN2_H_
29
30#include "vdtcore_common.h"
31#include "atan.h"
32
33namespace vdt{
34
35inline double fast_atan2( double y, double x ) {
36    // move in first octant
37    double xx = std::fabs(x);
38    double yy = std::fabs(y);
39    double tmp (0.0);
40    if (yy>xx) {
41        tmp = yy;
42        yy=xx; xx=tmp;
43        tmp=1.;
44    }
45
46    // To avoid the fpe, we protect against /0.
47    const double oneIfXXZero = (xx==0.);
48    double t=yy/(xx+oneIfXXZero);
49    double z=t;
50
51    double s = details::PIO4;
52    double factor = details::MOREBITSO2;
53    t = (t-1.0) / (t+1.0);
54
55    if( z > details::T3PO8 ) {
56        s = details::PIO2;
57        factor = details::MOREBITS;
58        t = -1.0 / z ;
59    }
60    if ( z <= 0.66 ) {
61        s = 0.;
62        factor = 0.;
63        t = z;
64    }
65
66    const double t2 = t * t;
67
68    const double px = details::get_atan_px(t2);
69    const double qx = details::get_atan_qx(t2);
70
71    //double res = y +x * x2 * px / qx + x +factor;
72
73    const double poq=px / qx;
74
75    double ret = t * t2 * poq + t;
76    ret+=s;
77
78    ret+=factor;
79
80    // Here we put the result to 0 if xx was 0, if not nothing happens!
81    ret*= (1.-oneIfXXZero);
82
83    // move back in place
84    if (y==0) ret=0.0;
85    if (tmp!=0) ret = details::PIO2 - ret;
86    if (x<0) ret = details::PI - ret;
87    if (y<0) ret = -ret;
88
89
90    return ret;
91}
92
93inline float fast_atan2f( float y, float x ) {
94
95    // move in first octant
96    float xx = std::fabs(x);
97    float yy = std::fabs(y);
98    float tmp (0.0f);
99    if (yy>xx) {
100        tmp = yy;
101        yy=xx; xx=tmp;
102        tmp =1.f;
103    }
104
105    // To avoid the fpe, we protect against /0.
106    const float oneIfXXZero = (xx==0.f);
107   
108    float t=yy/(xx/*+oneIfXXZero*/);
109    float z=t;
110    if( t > 0.4142135623730950f ) // * tan pi/8
111            z = (t-1.0f)/(t+1.0f);
112
113    //printf("%e %e %e %e\n",yy,xx,t,z);
114    float z2 = z * z;
115   
116    float ret =(((( 8.05374449538e-2f * z2
117                    - 1.38776856032E-1f) * z2
118                    + 1.99777106478E-1f) * z2
119                    - 3.33329491539E-1f) * z2 * z
120                    + z );
121
122    // Here we put the result to 0 if xx was 0, if not nothing happens!
123    ret*= (1.f - oneIfXXZero);
124   
125    // move back in place
126    if (y==0.f) ret=0.f;
127    if( t > 0.4142135623730950f ) ret += details::PIO4F;
128    if (tmp!=0) ret = details::PIO2F - ret;
129    if (x<0.f) ret = details::PIF - ret;
130    if (y<0.f) ret = -ret;
131
132    return ret;
133
134}
135
136
137
138//------------------------------------------------------------------------------
139
140} // end namespace vdt
141
142
143#endif
Note: See TracBrowser for help on using the repository browser.