/usr/include/mapnik/fastmath.hpp is in libmapnik-dev 2.2.0+ds1-6build2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 | /*****************************************************************************
*
* This file is part of Mapnik (c++ mapping toolkit)
*
* Copyright (C) 2011 Artem Pavlenko
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
*****************************************************************************/
#ifndef MAPNIK_FASTMATH_HPP
#define MAPNIK_FASTMATH_HPP
/* Timings:
* fast_sin(not inlined) 8.95s
* fast_sin(inlined) 6.64s
* sin 28.7s
* => 4.3x speedup
* worst case accuracy abs(fast_sin(x)/sin(x) - 1) = 0.000018 (at 0.664, M_PI-0.664, M_PI+0.664; 2*M_PI-0.664)
*/
static inline double fast_sin(double x)
{
bool negative = false;
double result;
while (x > 2*M_PI) x -= 2*M_PI;
while (x < 0) x += 2*M_PI;
if (x > M_PI) {
x -= M_PI;
negative = true;
}
if (x < 0.664 || x > M_PI-0.664) {
//series expansion at x=0: x-x^3/6+x^5/120-x^7/5040+...
if (x > M_PI-0.664) x = M_PI - x;
result = x*(1. + x*x*(-1/6. + x*x/120.));
} else {
//series expansion at x=pi/2
//1-x^2/2+x^4/24-x^6/720+x^8/40320+...
x -= M_PI/2;
result = 1. + x*x*(-1/2.+x*x*(1/24. + x*x*(-1/720.)));
}
return negative?-result:result;
}
static inline double fast_cos(double x)
{
return fast_sin(x + M_PI/2);
}
static inline double atan_helper(double x)
{
//Series expansion at x=0:
// x-x^3/3+x^5/5-x^7/7+...
if (x < 0.30) {
return x * (1 + x*x*(-1/3. + x*x*(1/5.) + x*x*(-1/7. + x*x*(1/9.))));
}
else if (x < 0.71) {
//Series expansion at x=0.5
//atan(1/2)+(4 x)/5-(8 x^2)/25-(16 x^3)/375+(96 x^4)/625 +...
x -= 0.5;
return 0.463647609000806116 /*atan(0.5) */ + x *(4./5. + x *(-8./25. + (-16./375.*x)));
} else {
//series expansion at x=1:
//pi/4+x/2-x^2/4+x^3/12-x^5/40+...
x -= 1;
return (M_PI/4.) + x * (1/2. + x*(-1/4. +x*(1/12. + x * (-1/40.))));
}
}
/*
* fast_atan(not inlined) 6.74s
* fast_atan(everything inlined) 6.78s
* fast_atan(only helper inlined) 6.75
* atan 27.5s
* => 4x speedup
* worst case accuracy abs(fast_atan(x)/atan(x) - 1) = 0.000271 (at 1.411)
*/
double inline fast_atan(double x)
{
double negative = false;
double result;
if (x < 0) {
x = -x;
negative = true;
}
if (x <= 1) result = atan_helper(x); else result = M_PI/2 - atan_helper(1/x);
return negative?-result:result;
}
static inline double fast_atan2(double y, double x)
{
double result = M_PI/2;
if (x == 0 && y == 0) return 0;
if (x != 0) result = fast_atan(y/x);
if (x < 0 && y >= 0) return result + M_PI;
if (x <= 0 && y < 0) return -M_PI + result;
return result;
}
#endif // MAPNIK_FASTMATH_HPP
|