BLI_math_rotation: properly name the quaternion power function.
[blender.git] / source / blender / blenlib / intern / math_base.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * The Original Code is: some of this file.
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  * */
25
26 /** \file blender/blenlib/intern/math_base.c
27  *  \ingroup bli
28  */
29
30 #include "BLI_math.h"
31
32 #include "BLI_strict_flags.h"
33
34 int pow_i(int base, int exp)
35 {
36         int result = 1;
37         BLI_assert(exp >= 0);
38         while (exp) {
39                 if (exp & 1) {
40                         result *= base;
41                 }
42                 exp >>= 1;
43                 base *= base;
44         }
45
46         return result;
47 }
48
49 /* from python 3.1 floatobject.c
50  * ndigits must be between 0 and 21 */
51 double double_round(double x, int ndigits)
52 {
53         double pow1, pow2, y, z;
54         if (ndigits >= 0) {
55                 pow1 = pow(10.0, (double)ndigits);
56                 pow2 = 1.0;
57                 y = (x * pow1) * pow2;
58                 /* if y overflows, then rounded value is exactly x */
59                 if (!isfinite(y))
60                         return x;
61         }
62         else {
63                 pow1 = pow(10.0, (double)-ndigits);
64                 pow2 = 1.0; /* unused; silences a gcc compiler warning */
65                 y = x / pow1;
66         }
67
68         z = round(y);
69         if (fabs(y - z) == 0.5)
70                 /* halfway between two integers; use round-half-even */
71                 z = 2.0 * round(y / 2.0);
72
73         if (ndigits >= 0)
74                 z = (z / pow2) / pow1;
75         else
76                 z *= pow1;
77
78         /* if computation resulted in overflow, raise OverflowError */
79         return z;
80 }