1 /*---------------------------------------------------------------------------*
2 Project: Horizon
3 File: math_Arithmetic.cpp
4
5 Copyright (C)2009 Nintendo Co., Ltd. All rights reserved.
6
7 These coded instructions, statements, and computer programs contain
8 proprietary information of Nintendo of America Inc. and/or Nintendo
9 Company Ltd., and are protected by Federal copyright law. They may
10 not be disclosed to third parties or copied or duplicated in any form,
11 in whole or in part, without the prior written consent of Nintendo.
12
13 $Rev: 24106 $
14 *---------------------------------------------------------------------------*/
15
16 #include <nn/math.h>
17
18 #include <nn/math/math_Arithmetic.h>
19 #include <cstdlib>
20
21 namespace nn { namespace math {
22
23 namespace
24 {
25 //---- 指数関数のテーブル
26 struct ExpTbl
27 {
28 f32 exp_val;
29 f32 exp_delta;
30 };
31 struct ExpTbl sExpTbl[32+1] =
32 {
33 { 0.500000000f, 0.022136891f }, // Exp(-0.693147181) = Exp(ln2 * (-16)/16)
34 { 0.522136891f, 0.023116975f }, // Exp(-0.649825482) = Exp(ln2 * (-15)/16)
35 { 0.545253866f, 0.024140451f }, // Exp(-0.606503783)
36 { 0.569394317f, 0.025209240f }, // Exp(-0.563182084)
37 { 0.594603558f, 0.026325349f }, // Exp(-0.519860385)
38 { 0.620928906f, 0.027490871f }, // Exp(-0.476538687)
39 { 0.648419777f, 0.028707996f }, // Exp(-0.433216988)
40 { 0.677127773f, 0.029979008f }, // Exp(-0.389895289)
41 { 0.707106781f, 0.031306292f }, // Exp(-0.346573590)
42 { 0.738413073f, 0.032692340f }, // Exp(-0.303251891)
43 { 0.771105413f, 0.034139753f }, // Exp(-0.259930193)
44 { 0.805245166f, 0.035651249f }, // Exp(-0.216608494)
45 { 0.840896415f, 0.037229665f }, // Exp(-0.173286795)
46 { 0.878126080f, 0.038877963f }, // Exp(-0.129965096)
47 { 0.917004043f, 0.040599237f }, // Exp(-0.086643398)
48 { 0.957603281f, 0.042396719f }, // Exp(-0.043321699) = Exp(ln2 * (-1)/16)
49 { 1.000000000f, 0.044273782f }, // Exp(0.000000000) = Exp(ln2 * (0)/16)
50 { 1.044273782f, 0.046233950f }, // Exp(0.043321699) = Exp(ln2 * (1)/16)
51 { 1.090507733f, 0.048280902f }, // Exp(0.086643398)
52 { 1.138788635f, 0.050418480f }, // Exp(0.129965096)
53 { 1.189207115f, 0.052650697f }, // Exp(0.173286795)
54 { 1.241857812f, 0.054981743f }, // Exp(0.216608494)
55 { 1.296839555f, 0.057415992f }, // Exp(0.259930193)
56 { 1.354255547f, 0.059958015f }, // Exp(0.303251891)
57 { 1.414213562f, 0.062612584f }, // Exp(0.346573590)
58 { 1.476826146f, 0.065384679f }, // Exp(0.389895289)
59 { 1.542210825f, 0.068279507f }, // Exp(0.433216988)
60 { 1.610490332f, 0.071302499f }, // Exp(0.476538687)
61 { 1.681792831f, 0.074459330f }, // Exp(0.519860385)
62 { 1.756252160f, 0.077755926f }, // Exp(0.563182084)
63 { 1.834008086f, 0.081198475f }, // Exp(0.606503783)
64 { 1.915206561f, 0.084793439f }, // Exp(0.649825482) = Exp(ln2 * 15/16)
65 { 2.000000000f, 0.088547565f }, // Exp(0.693147181) = Exp(ln2 * 16/16)
66 };
67
68
69
70 //---- 対数関数のテーブル
71 struct LogTbl
72 {
73 f32 log_val;
74 f32 log_delta;
75 };
76 struct LogTbl sLogTbl[256+1] =
77 {
78 { 0.000000000f, 0.003898640f }, // Log(1.00000000)
79 { 0.003898640f, 0.003883500f }, // Log(1.00390625)
80 { 0.007782140f, 0.003868477f }, // Log(1.00781250)
81 { 0.011650617f, 0.003853569f }, // Log(1.01171875)
82 { 0.015504187f, 0.003838776f }, // Log(1.01562500)
83 { 0.019342963f, 0.003824096f }, // Log(1.01953125)
84 { 0.023167059f, 0.003809528f }, // Log(1.02343750)
85 { 0.026976588f, 0.003795071f }, // Log(1.02734375)
86 { 0.030771659f, 0.003780723f }, // Log(1.03125000)
87 { 0.034552382f, 0.003766483f }, // Log(1.03515625)
88 { 0.038318864f, 0.003752350f }, // Log(1.03906250)
89 { 0.042071214f, 0.003738322f }, // Log(1.04296875)
90 { 0.045809536f, 0.003724399f }, // Log(1.04687500)
91 { 0.049533935f, 0.003710579f }, // Log(1.05078125)
92 { 0.053244515f, 0.003696862f }, // Log(1.05468750)
93 { 0.056941376f, 0.003683245f }, // Log(1.05859375)
94 { 0.060624622f, 0.003669729f }, // Log(1.06250000)
95 { 0.064294351f, 0.003656311f }, // Log(1.06640625)
96 { 0.067950662f, 0.003642991f }, // Log(1.07031250)
97 { 0.071593653f, 0.003629768f }, // Log(1.07421875)
98 { 0.075223421f, 0.003616640f }, // Log(1.07812500)
99 { 0.078840062f, 0.003603608f }, // Log(1.08203125)
100 { 0.082443669f, 0.003590668f }, // Log(1.08593750)
101 { 0.086034337f, 0.003577821f }, // Log(1.08984375)
102 { 0.089612159f, 0.003565066f }, // Log(1.09375000)
103 { 0.093177225f, 0.003552402f }, // Log(1.09765625)
104 { 0.096729626f, 0.003539827f }, // Log(1.10156250)
105 { 0.100269453f, 0.003527341f }, // Log(1.10546875)
106 { 0.103796794f, 0.003514942f }, // Log(1.10937500)
107 { 0.107311736f, 0.003502631f }, // Log(1.11328125)
108 { 0.110814366f, 0.003490405f }, // Log(1.11718750)
109 { 0.114304771f, 0.003478264f }, // Log(1.12109375)
110 { 0.117783036f, 0.003466208f }, // Log(1.12500000)
111 { 0.121249244f, 0.003454235f }, // Log(1.12890625)
112 { 0.124703479f, 0.003442344f }, // Log(1.13281250)
113 { 0.128145823f, 0.003430535f }, // Log(1.13671875)
114 { 0.131576358f, 0.003418807f }, // Log(1.14062500)
115 { 0.134995165f, 0.003407158f }, // Log(1.14453125)
116 { 0.138402323f, 0.003395589f }, // Log(1.14843750)
117 { 0.141797912f, 0.003384098f }, // Log(1.15234375)
118 { 0.145182010f, 0.003372684f }, // Log(1.15625000)
119 { 0.148554694f, 0.003361348f }, // Log(1.16015625)
120 { 0.151916042f, 0.003350087f }, // Log(1.16406250)
121 { 0.155266129f, 0.003338901f }, // Log(1.16796875)
122 { 0.158605030f, 0.003327790f }, // Log(1.17187500)
123 { 0.161932820f, 0.003316753f }, // Log(1.17578125)
124 { 0.165249573f, 0.003305788f }, // Log(1.17968750)
125 { 0.168555361f, 0.003294896f }, // Log(1.18359375)
126 { 0.171850257f, 0.003284075f }, // Log(1.18750000)
127 { 0.175134332f, 0.003273325f }, // Log(1.19140625)
128 { 0.178407657f, 0.003262646f }, // Log(1.19531250)
129 { 0.181670303f, 0.003252035f }, // Log(1.19921875)
130 { 0.184922338f, 0.003241494f }, // Log(1.20312500)
131 { 0.188163832f, 0.003231021f }, // Log(1.20703125)
132 { 0.191394853f, 0.003220615f }, // Log(1.21093750)
133 { 0.194615468f, 0.003210276f }, // Log(1.21484375)
134 { 0.197825743f, 0.003200003f }, // Log(1.21875000)
135 { 0.201025746f, 0.003189795f }, // Log(1.22265625)
136 { 0.204215541f, 0.003179653f }, // Log(1.22656250)
137 { 0.207395194f, 0.003169575f }, // Log(1.23046875)
138 { 0.210564769f, 0.003159560f }, // Log(1.23437500)
139 { 0.213724329f, 0.003149609f }, // Log(1.23828125)
140 { 0.216873938f, 0.003139720f }, // Log(1.24218750)
141 { 0.220013658f, 0.003129893f }, // Log(1.24609375)
142 { 0.223143551f, 0.003120127f }, // Log(1.25000000)
143 { 0.226263679f, 0.003110422f }, // Log(1.25390625)
144 { 0.229374101f, 0.003100778f }, // Log(1.25781250)
145 { 0.232474879f, 0.003091193f }, // Log(1.26171875)
146 { 0.235566071f, 0.003081667f }, // Log(1.26562500)
147 { 0.238647738f, 0.003072199f }, // Log(1.26953125)
148 { 0.241719937f, 0.003062790f }, // Log(1.27343750)
149 { 0.244782726f, 0.003053437f }, // Log(1.27734375)
150 { 0.247836164f, 0.003044142f }, // Log(1.28125000)
151 { 0.250880306f, 0.003034904f }, // Log(1.28515625)
152 { 0.253915210f, 0.003025721f }, // Log(1.28906250)
153 { 0.256940931f, 0.003016594f }, // Log(1.29296875)
154 { 0.259957524f, 0.003007521f }, // Log(1.29687500)
155 { 0.262965046f, 0.002998503f }, // Log(1.30078125)
156 { 0.265963548f, 0.002989539f }, // Log(1.30468750)
157 { 0.268953087f, 0.002980628f }, // Log(1.30859375)
158 { 0.271933715f, 0.002971770f }, // Log(1.31250000)
159 { 0.274905486f, 0.002962965f }, // Log(1.31640625)
160 { 0.277868451f, 0.002954212f }, // Log(1.32031250)
161 { 0.280822663f, 0.002945510f }, // Log(1.32421875)
162 { 0.283768173f, 0.002936860f }, // Log(1.32812500)
163 { 0.286705033f, 0.002928260f }, // Log(1.33203125)
164 { 0.289633293f, 0.002919710f }, // Log(1.33593750)
165 { 0.292553003f, 0.002911210f }, // Log(1.33984375)
166 { 0.295464213f, 0.002902760f }, // Log(1.34375000)
167 { 0.298366973f, 0.002894358f }, // Log(1.34765625)
168 { 0.301261331f, 0.002886005f }, // Log(1.35156250)
169 { 0.304147335f, 0.002877700f }, // Log(1.35546875)
170 { 0.307025035f, 0.002869442f }, // Log(1.35937500)
171 { 0.309894478f, 0.002861232f }, // Log(1.36328125)
172 { 0.312755710f, 0.002853069f }, // Log(1.36718750)
173 { 0.315608779f, 0.002844952f }, // Log(1.37109375)
174 { 0.318453731f, 0.002836881f }, // Log(1.37500000)
175 { 0.321290612f, 0.002828856f }, // Log(1.37890625)
176 { 0.324119469f, 0.002820876f }, // Log(1.38281250)
177 { 0.326940345f, 0.002812941f }, // Log(1.38671875)
178 { 0.329753286f, 0.002805051f }, // Log(1.39062500)
179 { 0.332558337f, 0.002797205f }, // Log(1.39453125)
180 { 0.335355542f, 0.002789402f }, // Log(1.39843750)
181 { 0.338144944f, 0.002781643f }, // Log(1.40234375)
182 { 0.340926587f, 0.002773927f }, // Log(1.40625000)
183 { 0.343700514f, 0.002766253f }, // Log(1.41015625)
184 { 0.346466767f, 0.002758622f }, // Log(1.41406250)
185 { 0.349225390f, 0.002751033f }, // Log(1.41796875)
186 { 0.351976423f, 0.002743486f }, // Log(1.42187500)
187 { 0.354719909f, 0.002735980f }, // Log(1.42578125)
188 { 0.357455889f, 0.002728515f }, // Log(1.42968750)
189 { 0.360184404f, 0.002721090f }, // Log(1.43359375)
190 { 0.362905494f, 0.002713706f }, // Log(1.43750000)
191 { 0.365619200f, 0.002706362f }, // Log(1.44140625)
192 { 0.368325561f, 0.002699057f }, // Log(1.44531250)
193 { 0.371024618f, 0.002691792f }, // Log(1.44921875)
194 { 0.373716410f, 0.002684565f }, // Log(1.45312500)
195 { 0.376400975f, 0.002677378f }, // Log(1.45703125)
196 { 0.379078353f, 0.002670229f }, // Log(1.46093750)
197 { 0.381748581f, 0.002663117f }, // Log(1.46484375)
198 { 0.384411699f, 0.002656044f }, // Log(1.46875000)
199 { 0.387067743f, 0.002649008f }, // Log(1.47265625)
200 { 0.389716751f, 0.002642009f }, // Log(1.47656250)
201 { 0.392358761f, 0.002635048f }, // Log(1.48046875)
202 { 0.394993808f, 0.002628122f }, // Log(1.48437500)
203 { 0.397621931f, 0.002621233f }, // Log(1.48828125)
204 { 0.400243164f, 0.002614381f }, // Log(1.49218750)
205 { 0.402857545f, 0.002607563f }, // Log(1.49609375)
206 { 0.405465108f, 0.002600782f }, // Log(1.50000000)
207 { 0.408065890f, 0.002594035f }, // Log(1.50390625)
208 { 0.410659925f, 0.002587324f }, // Log(1.50781250)
209 { 0.413247249f, 0.002580647f }, // Log(1.51171875)
210 { 0.415827895f, 0.002574004f }, // Log(1.51562500)
211 { 0.418401899f, 0.002567396f }, // Log(1.51953125)
212 { 0.420969295f, 0.002560821f }, // Log(1.52343750)
213 { 0.423530116f, 0.002554280f }, // Log(1.52734375)
214 { 0.426084395f, 0.002547772f }, // Log(1.53125000)
215 { 0.428632167f, 0.002541297f }, // Log(1.53515625)
216 { 0.431173465f, 0.002534856f }, // Log(1.53906250)
217 { 0.433708320f, 0.002528446f }, // Log(1.54296875)
218 { 0.436236767f, 0.002522069f }, // Log(1.54687500)
219 { 0.438758836f, 0.002515725f }, // Log(1.55078125)
220 { 0.441274561f, 0.002509412f }, // Log(1.55468750)
221 { 0.443783972f, 0.002503130f }, // Log(1.55859375)
222 { 0.446287103f, 0.002496880f }, // Log(1.56250000)
223 { 0.448783983f, 0.002490661f }, // Log(1.56640625)
224 { 0.451274644f, 0.002484473f }, // Log(1.57031250)
225 { 0.453759117f, 0.002478316f }, // Log(1.57421875)
226 { 0.456237433f, 0.002472189f }, // Log(1.57812500)
227 { 0.458709623f, 0.002466092f }, // Log(1.58203125)
228 { 0.461175715f, 0.002460026f }, // Log(1.58593750)
229 { 0.463635741f, 0.002453989f }, // Log(1.58984375)
230 { 0.466089730f, 0.002447982f }, // Log(1.59375000)
231 { 0.468537712f, 0.002442004f }, // Log(1.59765625)
232 { 0.470979715f, 0.002436055f }, // Log(1.60156250)
233 { 0.473415770f, 0.002430135f }, // Log(1.60546875)
234 { 0.475845905f, 0.002424244f }, // Log(1.60937500)
235 { 0.478270148f, 0.002418381f }, // Log(1.61328125)
236 { 0.480688529f, 0.002412546f }, // Log(1.61718750)
237 { 0.483101076f, 0.002406740f }, // Log(1.62109375)
238 { 0.485507816f, 0.002400962f }, // Log(1.62500000)
239 { 0.487908777f, 0.002395211f }, // Log(1.62890625)
240 { 0.490303988f, 0.002389487f }, // Log(1.63281250)
241 { 0.492693475f, 0.002383791f }, // Log(1.63671875)
242 { 0.495077267f, 0.002378122f }, // Log(1.64062500)
243 { 0.497455389f, 0.002372480f }, // Log(1.64453125)
244 { 0.499827870f, 0.002366865f }, // Log(1.64843750)
245 { 0.502194735f, 0.002361276f }, // Log(1.65234375)
246 { 0.504556011f, 0.002355714f }, // Log(1.65625000)
247 { 0.506911724f, 0.002350177f }, // Log(1.66015625)
248 { 0.509261902f, 0.002344667f }, // Log(1.66406250)
249 { 0.511606569f, 0.002339182f }, // Log(1.66796875)
250 { 0.513945751f, 0.002333723f }, // Log(1.67187500)
251 { 0.516279474f, 0.002328290f }, // Log(1.67578125)
252 { 0.518607764f, 0.002322881f }, // Log(1.67968750)
253 { 0.520930646f, 0.002317498f }, // Log(1.68359375)
254 { 0.523248144f, 0.002312140f }, // Log(1.68750000)
255 { 0.525560284f, 0.002306806f }, // Log(1.69140625)
256 { 0.527867090f, 0.002301497f }, // Log(1.69531250)
257 { 0.530168587f, 0.002296212f }, // Log(1.69921875)
258 { 0.532464799f, 0.002290952f }, // Log(1.70312500)
259 { 0.534755751f, 0.002285715f }, // Log(1.70703125)
260 { 0.537041466f, 0.002280503f }, // Log(1.71093750)
261 { 0.539321969f, 0.002275314f }, // Log(1.71484375)
262 { 0.541597282f, 0.002270149f }, // Log(1.71875000)
263 { 0.543867431f, 0.002265007f }, // Log(1.72265625)
264 { 0.546132438f, 0.002259888f }, // Log(1.72656250)
265 { 0.548392326f, 0.002254792f }, // Log(1.73046875)
266 { 0.550647118f, 0.002249720f }, // Log(1.73437500)
267 { 0.552896838f, 0.002244670f }, // Log(1.73828125)
268 { 0.555141508f, 0.002239643f }, // Log(1.74218750)
269 { 0.557381150f, 0.002234638f }, // Log(1.74609375)
270 { 0.559615788f, 0.002229655f }, // Log(1.75000000)
271 { 0.561845443f, 0.002224695f }, // Log(1.75390625)
272 { 0.564070138f, 0.002219757f }, // Log(1.75781250)
273 { 0.566289895f, 0.002214840f }, // Log(1.76171875)
274 { 0.568504735f, 0.002209946f }, // Log(1.76562500)
275 { 0.570714681f, 0.002205073f }, // Log(1.76953125)
276 { 0.572919754f, 0.002200221f }, // Log(1.77343750)
277 { 0.575119974f, 0.002195391f }, // Log(1.77734375)
278 { 0.577315365f, 0.002190581f }, // Log(1.78125000)
279 { 0.579505946f, 0.002185793f }, // Log(1.78515625)
280 { 0.581691740f, 0.002181026f }, // Log(1.78906250)
281 { 0.583872766f, 0.002176279f }, // Log(1.79296875)
282 { 0.586049045f, 0.002171554f }, // Log(1.79687500)
283 { 0.588220599f, 0.002166848f }, // Log(1.80078125)
284 { 0.590387447f, 0.002162163f }, // Log(1.80468750)
285 { 0.592549610f, 0.002157498f }, // Log(1.80859375)
286 { 0.594707108f, 0.002152853f }, // Log(1.81250000)
287 { 0.596859961f, 0.002148229f }, // Log(1.81640625)
288 { 0.599008190f, 0.002143624f }, // Log(1.82031250)
289 { 0.601151813f, 0.002139038f }, // Log(1.82421875)
290 { 0.603290851f, 0.002134473f }, // Log(1.82812500)
291 { 0.605425324f, 0.002129926f }, // Log(1.83203125)
292 { 0.607555250f, 0.002125399f }, // Log(1.83593750)
293 { 0.609680650f, 0.002120892f }, // Log(1.83984375)
294 { 0.611801541f, 0.002116403f }, // Log(1.84375000)
295 { 0.613917944f, 0.002111933f }, // Log(1.84765625)
296 { 0.616029877f, 0.002107482f }, // Log(1.85156250)
297 { 0.618137360f, 0.002103050f }, // Log(1.85546875)
298 { 0.620240410f, 0.002098637f }, // Log(1.85937500)
299 { 0.622339046f, 0.002094242f }, // Log(1.86328125)
300 { 0.624433288f, 0.002089865f }, // Log(1.86718750)
301 { 0.626523153f, 0.002085506f }, // Log(1.87109375)
302 { 0.628608659f, 0.002081166f }, // Log(1.87500000)
303 { 0.630689826f, 0.002076844f }, // Log(1.87890625)
304 { 0.632766670f, 0.002072540f }, // Log(1.88281250)
305 { 0.634839209f, 0.002068253f }, // Log(1.88671875)
306 { 0.636907462f, 0.002063984f }, // Log(1.89062500)
307 { 0.638971446f, 0.002059733f }, // Log(1.89453125)
308 { 0.641031179f, 0.002055499f }, // Log(1.89843750)
309 { 0.643086679f, 0.002051283f }, // Log(1.90234375)
310 { 0.645137961f, 0.002047084f }, // Log(1.90625000)
311 { 0.647185045f, 0.002042902f }, // Log(1.91015625)
312 { 0.649227947f, 0.002038737f }, // Log(1.91406250)
313 { 0.651266683f, 0.002034589f }, // Log(1.91796875)
314 { 0.653301272f, 0.002030458f }, // Log(1.92187500)
315 { 0.655331730f, 0.002026343f }, // Log(1.92578125)
316 { 0.657358073f, 0.002022245f }, // Log(1.92968750)
317 { 0.659380318f, 0.002018164f }, // Log(1.93359375)
318 { 0.661398482f, 0.002014099f }, // Log(1.93750000)
319 { 0.663412582f, 0.002010051f }, // Log(1.94140625)
320 { 0.665422633f, 0.002006019f }, // Log(1.94531250)
321 { 0.667428651f, 0.002002003f }, // Log(1.94921875)
322 { 0.669430654f, 0.001998003f }, // Log(1.95312500)
323 { 0.671428657f, 0.001994019f }, // Log(1.95703125)
324 { 0.673422675f, 0.001990050f }, // Log(1.96093750)
325 { 0.675412726f, 0.001986098f }, // Log(1.96484375)
326 { 0.677398824f, 0.001982161f }, // Log(1.96875000)
327 { 0.679380985f, 0.001978240f }, // Log(1.97265625)
328 { 0.681359225f, 0.001974334f }, // Log(1.97656250)
329 { 0.683333559f, 0.001970444f }, // Log(1.98046875)
330 { 0.685304003f, 0.001966569f }, // Log(1.98437500)
331 { 0.687270572f, 0.001962709f }, // Log(1.98828125)
332 { 0.689233281f, 0.001958864f }, // Log(1.99218750)
333 { 0.691192146f, 0.001955035f }, // Log(1.99609375)
334 { 0.693147181f, 0.001951220f }, // Log(2.00000000)
335 };
336
337
338
339 /*!--------------------------------------------------------------------------*
340 Name: FExpLn2
341
342 @brief 指数関数の値をテーブル引きで求めます。
343 e^x を求めます。
344
345 @param[in] x 指数の値。
346 - loge(2) < x < loge(2) でなければなりません。
347
348 @return e^x の近似値。
349 *---------------------------------------------------------------------------*/
350 NN_MATH_INLINE f32
FExpLn2(f32 x)351 FExpLn2(f32 x)
352 {
353 f32 fidx;
354 u16 idx;
355 f32 r;
356
357 // |x| < ln2 でなくてはならない
358 // map [-ln2, ln2) -> [0, 32)
359 fidx = (x + F_LN2) * (32 / (2 * F_LN2));
360 idx = F32ToU16(fidx);
361 r = fidx - U16ToF32(idx);
362
363 return sExpTbl[idx].exp_val + r * sExpTbl[idx].exp_delta;
364 }
365
366
367
368 /*!--------------------------------------------------------------------------*
369 Name: FLog1_2
370
371 @brief 対数関数の値をテーブル引きで求めます。
372
373 @param[in] x 対数を求める値。
374 1 <= x < 2 でなければなりません。
375
376 @return loge(x) の近似値。
377 *---------------------------------------------------------------------------*/
378 NN_MATH_INLINE f32
FLog1_2(f32 x)379 FLog1_2(f32 x)
380 {
381 f32 fidx;
382 u16 idx;
383 f32 r;
384
385 // 1 <= x < 2 でなくてなならない
386 // map [1, 2) -> [0, 256)
387 fidx = (x - 1) * 256.f; // 0 <= fidx < 256
388 idx = F32ToU16(fidx);
389 r = fidx - U16ToF32(idx);
390
391 return sLogTbl[idx].log_val + r * sLogTbl[idx].log_delta;
392 }
393
394
395 } // end of anonymous namespace
396
397
398 namespace internal
399 {
400
401
402 /*!--------------------------------------------------------------------------*
403 Name: FExp
404
405 @brief 指数関数の値をテーブル引きで求めます。
406 e^x を求めます。
407
408 @param[in] x 指数の値。
409
410 @return e^x の近似値。
411 *---------------------------------------------------------------------------*/
412 // -10 から 10 までを両端を含んで2^16分割し、各サンプル点において誤差を計測
413 // sExpTbl[32]の場合 std に対する最大相対誤差 0.0235 % 平均相対誤差 0.0156 %
414 // sExpTbl[256]の場合 std に対する最大相対誤差 0.000427 % 平均相対誤差 0.000245 %
415 f32
FExp(f32 x)416 FExp(f32 x)
417 {
418 /*
419 exp(x) = exp(xn + k * ln2)
420 = exp(xn) * exp(k * ln2)
421 = exp(xn) * 2^k
422
423 |xn| < ln2
424 */
425
426 s16 k;
427 f32 kf;
428 f32 xn;
429 f32 expxn;
430
431 // x を xn + k * ln2 の形に変形します
432 // ただし |xn| < ln2, k は整数 です
433 // -128 <= k <= 127 でなければなりません
434 k = F32ToS16(F_INVLN2 * x);
435 kf = S16ToF32(k);
436 xn = x - kf * F_LN2;
437
438 // exp(xn) をテーブル引きで求めます
439 expxn = FExpLn2(xn);
440
441 // exp(x) = exp(xn) * 2^k
442 // 直接指数部に足しこみます
443 // Expの結果は負にならないので桁上がり分はマスクして消します
444 u32 expx = (F32AsU32(expxn) + (k << 23)) & 0x7FFFFFFF;
445 return U32AsF32(expx);
446 }
447
448
449
450
451
452 /*!--------------------------------------------------------------------------*
453 Name: FLog
454
455 @brief 対数関数の値をテーブル引きで求めます。
456
457 @param[in] x 対数を求める値。
458
459 @return loge(x) の近似値。
460 *---------------------------------------------------------------------------*/
461 // 0.001 から 10 までを両端を含んで2^16分割し、各サンプル点において誤差を計測
462 // sLogTbl[32]の場合 std に対する最大相対誤差 1.528 % 平均相対誤差 0.00857 %
463 // sLogTbl[128]の場合 std に対する最大相対誤差 0.386 % 平均相対誤差 0.000623 %
464 // sLogTbl[256]の場合 std に対する最大相対誤差 0.192 % 平均相対誤差 0.000167 %
465 // 相対誤差はx=1の前後で最大
466 f32
FLog(f32 x)467 FLog(f32 x)
468 {
469 /*
470 log(x) = log(xn * 2^k)
471 = log(xn) + log(2^k)
472 = log(xn) + k * ln2
473
474 1 < xn < 2
475 */
476
477 s32 k;
478 f32 kf;
479 f32 xn;
480 f32 logxn;
481
482 // x を xn * 2^k の形に変形します
483 k = FGetExpPart(x);
484 xn = FGetMantPart(x);
485 kf = S16ToF32(static_cast<s16>(k));
486
487 // log(xn) をテーブル引きで求めます
488 logxn = FLog1_2(xn);
489
490 // log(x) = log(xn) + k * ln2
491 return logxn + kf * F_LN2;
492 }
493
494 } // end of namespace internal
495
496
497
498 /*!--------------------------------------------------------------------------*
499 Name: Bezier
500
501 @brief ベジェ補間を行います
502
503 @param[in] p1 点1での値
504 @param[in] p2 点1での制御値
505 @param[in] p3 点2での制御値
506 @param[in] p4 点2での値
507 @param[in] s 補間対象位置。(点1:0.0～1.0:点2)
508
509 @return 補間結果の値。
510 *---------------------------------------------------------------------------*/
Bezier(f32 p1,f32 p2,f32 p3,f32 p4,f32 s)511 f32 Bezier(f32 p1, f32 p2, f32 p3, f32 p4, f32 s)
512 {
513 // もう少しマシな計算方法にすること
514 f32 t = 1.f - s;
515 f32 tt = t * t;
516 f32 ss = s * s;
517
518 f32 a1 = tt * t;
519 f32 a2 = tt * s * 3.f;
520 f32 a3 = ss * t * 3.f;
521 f32 a4 = ss * s;
522
523 return a1*p1 + a2*p2 + a3*p3 + a4*p4;
524 }
525
526
527 /*!--------------------------------------------------------------------------*
528 Name: CatmullRom
529
530 @brief Catmull-Rom 補間 を行います
531
532 @param[in] p0 点1での制御値
533 @param[in] p1 点1での値
534 @param[in] p2 点2での値
535 @param[in] p3 点2での制御値
536 @param[in] s 補間対象位置。(点1:0.0～1.0:点2)
537
538 @return 補間結果の値。
539 *---------------------------------------------------------------------------*/
CatmullRom(f32 p0,f32 p1,f32 p2,f32 p3,f32 s)540 f32 CatmullRom(f32 p0, f32 p1, f32 p2, f32 p3, f32 s)
541 {
542 // もう少しマシな計算方法にすること
543 return Hermite(p1, 0.5f * (p0 + p2),
544 p2, 0.5f * (p1 + p3),
545 s);
546 }
547
548
549
550 /*!--------------------------------------------------------------------------*
551 Name: CntBit1
552
553 @brief ビット列の中で1となっているビットの数を数えます。
554
555 @param[in] x 走査するビット列。
556
557 @return x 中で1となっているビットの数を返します。
558 *---------------------------------------------------------------------------*/
CntBit1(u32 x)559 u32 CntBit1(u32 x)
560 {
561 x = x - ((x >> 1) & 0x55555555);
562 x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
563 x = (x + (x >> 4)) & 0x0f0f0f0f;
564 x = x + (x >> 8);
565 x = x + (x >> 16);
566 return x & 0x0000003f;
567 }
568
569
570
571 /*!--------------------------------------------------------------------------*
572 Name: CntBit1
573
574 @brief ビット列の中で 1 となっているビットの数を数えます。
575
576 @param[in] first ビット列の先頭ワードへのポインタ。
577 @param[in] last ビット列の最終ワードの次のワードへのポインタ。
578
579 @return ビット列中で1となっているビットの数を返します。
580 *---------------------------------------------------------------------------*/
CntBit1(const u32 * first,const u32 * last)581 u32 CntBit1(const u32* first, const u32* last)
582 {
583 const u32 n = u32(last - first);
584
585 u32 i, j, lim;
586 unsigned int s = 0;
587 unsigned int ss, x;
588
589 for (i = 0; i < n; i += 31)
590 {
591 lim = (n < i + 31) ? n : (i + 31);
592 ss = 0;
593 for (j = i; j < lim; ++j)
594 {
595 x = *(first + j);
596 x -= ((x >> 1) & 0x55555555);
597 x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
598 x = (x + (x >> 4)) & 0x0f0f0f0f;
599 ss += x;
600 }
601 x = (ss & 0x00ff00ff) + ((ss >> 8) & 0x00ff00ff);
602 x = (x & 0x0000ffff) + (x >> 16);
603 s += x;
604 }
605 return s;
606 }
607
608
609 /*!--------------------------------------------------------------------------*
610 Name: DistBit
611
612 @brief ビット列間の距離を求めます。
613 ビット列間の距離はビット列同士をビットごとに比較した場合の
614 等しいビットの数です。
615
616 @param[in] first1 ビット列 1 の先頭ワードへのポインタ。
617 @param[in] last1 ビット列 1 の最終ワードの次のワードへのポインタ。
618 @param[in] first2 ビット列 2 の先頭ワードへのポインタ。
619
620 @return ビット列 1 とビット列 2 の距離を返します。
621 *---------------------------------------------------------------------------*/
DistBit(const u32 * first1,const u32 * last1,const u32 * first2)622 u32 DistBit(const u32* first1, const u32* last1, const u32* first2)
623 {
624 unsigned int n = 0;
625 while(first1 != last1)
626 {
627 n += DistBit(*first1, *first2);
628 ++first1;
629 ++first2;
630 }
631 return n;
632 }
633
634
635 /*!--------------------------------------------------------------------------*
636 Name: RevBit
637
638 @brief ビット列の順序を反転します。
639
640 @param[in] x 反転するビット列
641
642 @return x のビット順序を逆順にしたビット列を返します。
643 *---------------------------------------------------------------------------*/
RevBit(u32 x)644 u32 RevBit(u32 x)
645 {
646 x = ((x & 0x55555555) << 1) | ((x >> 1) & 0x55555555);
647 x = ((x & 0x33333333) << 2) | ((x >> 2) & 0x33333333);
648 x = ((x & 0x0f0f0f0f) << 4) | ((x >> 4) & 0x0f0f0f0f);
649 x = (x << 24) | ((x & 0xff00) << 8) | ((x >> 8) & 0xff00) | (x >> 24);
650
651 return x;
652 }
653
654
655 /*!--------------------------------------------------------------------------*
656 Name: IExp
657
658 @brief 整数の非負整数乗を求めます。
659
660 @param[in] x 底となる値
661 @param[in] n 指数の値。非負でなければなりません。
662
663 @return x^n を返します
664 *---------------------------------------------------------------------------*/
665 int
IExp(int x,u32 n)666 IExp(int x, u32 n)
667 {
668 int p, y;
669
670 y = 1;
671 p = x;
672 for (;;)
673 {
674 if (n & 1) { y *= p; }
675 n >>= 1;
676 if (n == 0) { return y; }
677 p *= p;
678 }
679 }
680
681
682 /*!--------------------------------------------------------------------------*
683 Name: ILog10
684
685 @brief 整数の常用対数を計算し、結果を整数で返します。
686
687 @param[in] x 常用対数を求める値。
688
689 @return x の常用対数を整数で返します。
690 *---------------------------------------------------------------------------*/
ILog10(u32 x)691 u32 ILog10(u32 x)
692 {
693 u32 y;
694 static u32 table2[11] = {
695 0, 9, 99, 999, 9999, 99999, 999999, 9999999, 99999999, 999999999, 0xFFFFFFFF
696 };
697
698 y = (19 * (31 - CntLz(x))) >> 6;
699 y += ((table2[y + 1] - x) >> 31);
700
701 return y;
702 }
703
704
705 namespace internal
706 {
707
708 /*!--------------------------------------------------------------------------*
709 Name: CntLz_
710
711 @brief ビット列の左側からの連続する 0 のビットの数を数えます。
712 専用の命令を持たない PC 用です。
713
714 @param[in] x ビット列。
715
716 @return x 中で左側からの 0 のビットが連続する数を返します。
717 *---------------------------------------------------------------------------*/
718 u32
CntLz_(u32 x)719 CntLz_(u32 x)
720 {
721 unsigned int y;
722 int n, c;
723
724 n = 32;
725 c = 16;
726 do {
727 y = x >> c;
728 if (y != 0)
729 {
730 n -= c; x = y;
731 }
732 c >>= 1;
733 } while(c != 0);
734 return n - x;
735 }
736
737 } // internal
738
739
740 }} // nw::math
741