Cycles: Code cleanup, spaces around keywords
[blender-staging.git] / intern / cycles / util / util_math.h
1 /*
2  * Copyright 2011-2013 Blender Foundation
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16
17 #ifndef __UTIL_MATH_H__
18 #define __UTIL_MATH_H__
19
20 /* Math
21  *
22  * Basic math functions on scalar and vector types. This header is used by
23  * both the kernel code when compiled as C++, and other C++ non-kernel code. */
24
25 #ifndef __KERNEL_OPENCL__
26
27 #ifdef _MSC_VER
28 #  define _USE_MATH_DEFINES
29 #endif
30
31 #include <float.h>
32 #include <math.h>
33 #include <stdio.h>
34
35 #endif
36
37 #include "util_types.h"
38
39 CCL_NAMESPACE_BEGIN
40
41 /* Float Pi variations */
42
43 /* Division */
44 #ifndef M_PI_F
45 #define M_PI_F          ((float)3.14159265358979323846264338327950288)          /* pi */
46 #endif
47 #ifndef M_PI_2_F
48 #define M_PI_2_F        ((float)1.57079632679489661923132169163975144)          /* pi/2 */
49 #endif
50 #ifndef M_PI_4_F
51 #define M_PI_4_F        ((float)0.785398163397448309615660845819875721)         /* pi/4 */
52 #endif
53 #ifndef M_1_PI_F
54 #define M_1_PI_F        ((float)0.318309886183790671537767526745028724)         /* 1/pi */
55 #endif
56 #ifndef M_2_PI_F
57 #define M_2_PI_F        ((float)0.636619772367581343075535053490057448)         /* 2/pi */
58 #endif
59
60 /* Multiplication */
61 #ifndef M_2PI_F
62 #define M_2PI_F         ((float)6.283185307179586476925286766559005768)         /* 2*pi */
63 #endif
64 #ifndef M_4PI_F
65 #define M_4PI_F         ((float)12.56637061435917295385057353311801153)         /* 4*pi */
66 #endif
67
68 /* Float sqrt variations */
69
70 #ifndef M_SQRT2_F
71 #define M_SQRT2_F       ((float)1.41421356237309504880)                                         /* sqrt(2) */
72 #endif
73
74 #ifndef M_LN2_F
75 #define M_LN2_F      ((float)0.6931471805599453)        /* ln(2) */
76 #endif
77
78 #ifndef M_LN10_F
79 #define M_LN10_F     ((float)2.3025850929940457)        /* ln(10) */
80 #endif
81
82 /* Scalar */
83
84 #ifdef _WIN32
85
86 #ifndef __KERNEL_OPENCL__
87
88 ccl_device_inline float fmaxf(float a, float b)
89 {
90         return (a > b)? a: b;
91 }
92
93 ccl_device_inline float fminf(float a, float b)
94 {
95         return (a < b)? a: b;
96 }
97
98 #endif
99
100 #endif
101
102 #ifndef __KERNEL_GPU__
103
104 ccl_device_inline int max(int a, int b)
105 {
106         return (a > b)? a: b;
107 }
108
109 ccl_device_inline int min(int a, int b)
110 {
111         return (a < b)? a: b;
112 }
113
114 ccl_device_inline float max(float a, float b)
115 {
116         return (a > b)? a: b;
117 }
118
119 ccl_device_inline float min(float a, float b)
120 {
121         return (a < b)? a: b;
122 }
123
124 ccl_device_inline double max(double a, double b)
125 {
126         return (a > b)? a: b;
127 }
128
129 ccl_device_inline double min(double a, double b)
130 {
131         return (a < b)? a: b;
132 }
133
134 /* These 2 guys are templated for usage with registers data.
135  *
136  * NOTE: Since this is CPU-only functions it is ok to use references here.
137  * But for other devices we'll need to be careful about this.
138  */
139
140 template<typename T>
141 ccl_device_inline T min4(const T& a, const T& b, const T& c, const T& d)
142 {
143         return min(min(a,b),min(c,d));
144 }
145
146 template<typename T>
147 ccl_device_inline T max4(const T& a, const T& b, const T& c, const T& d)
148 {
149         return max(max(a,b),max(c,d));
150 }
151
152 #endif
153
154 ccl_device_inline float min4(float a, float b, float c, float d)
155 {
156         return min(min(a, b), min(c, d));
157 }
158
159 ccl_device_inline float max4(float a, float b, float c, float d)
160 {
161         return max(max(a, b), max(c, d));
162 }
163
164 #ifndef __KERNEL_OPENCL__
165
166 ccl_device_inline int clamp(int a, int mn, int mx)
167 {
168         return min(max(a, mn), mx);
169 }
170
171 ccl_device_inline float clamp(float a, float mn, float mx)
172 {
173         return min(max(a, mn), mx);
174 }
175
176 #endif
177
178 ccl_device_inline int float_to_int(float f)
179 {
180         return (int)f;
181 }
182
183 ccl_device_inline int floor_to_int(float f)
184 {
185         return float_to_int(floorf(f));
186 }
187
188 ccl_device_inline int ceil_to_int(float f)
189 {
190         return float_to_int(ceilf(f));
191 }
192
193 ccl_device_inline float signf(float f)
194 {
195         return (f < 0.0f)? -1.0f: 1.0f;
196 }
197
198 ccl_device_inline float nonzerof(float f, float eps)
199 {
200         if(fabsf(f) < eps)
201                 return signf(f)*eps;
202         else
203                 return f;
204 }
205
206 ccl_device_inline float smoothstepf(float f)
207 {
208         float ff = f*f;
209         return (3.0f*ff - 2.0f*ff*f);
210 }
211
212 /* Float2 Vector */
213
214 #ifndef __KERNEL_OPENCL__
215
216 ccl_device_inline bool is_zero(const float2 a)
217 {
218         return (a.x == 0.0f && a.y == 0.0f);
219 }
220
221 #endif
222
223 #ifndef __KERNEL_OPENCL__
224
225 ccl_device_inline float average(const float2 a)
226 {
227         return (a.x + a.y)*(1.0f/2.0f);
228 }
229
230 #endif
231
232 #ifndef __KERNEL_OPENCL__
233
234 ccl_device_inline float2 operator-(const float2 a)
235 {
236         return make_float2(-a.x, -a.y);
237 }
238
239 ccl_device_inline float2 operator*(const float2 a, const float2 b)
240 {
241         return make_float2(a.x*b.x, a.y*b.y);
242 }
243
244 ccl_device_inline float2 operator*(const float2 a, float f)
245 {
246         return make_float2(a.x*f, a.y*f);
247 }
248
249 ccl_device_inline float2 operator*(float f, const float2 a)
250 {
251         return make_float2(a.x*f, a.y*f);
252 }
253
254 ccl_device_inline float2 operator/(float f, const float2 a)
255 {
256         return make_float2(f/a.x, f/a.y);
257 }
258
259 ccl_device_inline float2 operator/(const float2 a, float f)
260 {
261         float invf = 1.0f/f;
262         return make_float2(a.x*invf, a.y*invf);
263 }
264
265 ccl_device_inline float2 operator/(const float2 a, const float2 b)
266 {
267         return make_float2(a.x/b.x, a.y/b.y);
268 }
269
270 ccl_device_inline float2 operator+(const float2 a, const float2 b)
271 {
272         return make_float2(a.x+b.x, a.y+b.y);
273 }
274
275 ccl_device_inline float2 operator-(const float2 a, const float2 b)
276 {
277         return make_float2(a.x-b.x, a.y-b.y);
278 }
279
280 ccl_device_inline float2 operator+=(float2& a, const float2 b)
281 {
282         return a = a + b;
283 }
284
285 ccl_device_inline float2 operator*=(float2& a, const float2 b)
286 {
287         return a = a * b;
288 }
289
290 ccl_device_inline float2 operator*=(float2& a, float f)
291 {
292         return a = a * f;
293 }
294
295 ccl_device_inline float2 operator/=(float2& a, const float2 b)
296 {
297         return a = a / b;
298 }
299
300 ccl_device_inline float2 operator/=(float2& a, float f)
301 {
302         float invf = 1.0f/f;
303         return a = a * invf;
304 }
305
306
307 ccl_device_inline float dot(const float2 a, const float2 b)
308 {
309         return a.x*b.x + a.y*b.y;
310 }
311
312 ccl_device_inline float cross(const float2 a, const float2 b)
313 {
314         return (a.x*b.y - a.y*b.x);
315 }
316
317 #endif
318
319 #ifndef __KERNEL_OPENCL__
320
321 ccl_device_inline bool operator==(const int2 a, const int2 b)
322 {
323         return (a.x == b.x && a.y == b.y);
324 }
325
326 ccl_device_inline float len(const float2 a)
327 {
328         return sqrtf(dot(a, a));
329 }
330
331 ccl_device_inline float2 normalize(const float2 a)
332 {
333         return a/len(a);
334 }
335
336 ccl_device_inline float2 normalize_len(const float2 a, float *t)
337 {
338         *t = len(a);
339         return a/(*t);
340 }
341
342 ccl_device_inline float2 safe_normalize(const float2 a)
343 {
344         float t = len(a);
345         return (t)? a/t: a;
346 }
347
348 ccl_device_inline bool operator==(const float2 a, const float2 b)
349 {
350         return (a.x == b.x && a.y == b.y);
351 }
352
353 ccl_device_inline bool operator!=(const float2 a, const float2 b)
354 {
355         return !(a == b);
356 }
357
358 ccl_device_inline float2 min(float2 a, float2 b)
359 {
360         return make_float2(min(a.x, b.x), min(a.y, b.y));
361 }
362
363 ccl_device_inline float2 max(float2 a, float2 b)
364 {
365         return make_float2(max(a.x, b.x), max(a.y, b.y));
366 }
367
368 ccl_device_inline float2 clamp(float2 a, float2 mn, float2 mx)
369 {
370         return min(max(a, mn), mx);
371 }
372
373 ccl_device_inline float2 fabs(float2 a)
374 {
375         return make_float2(fabsf(a.x), fabsf(a.y));
376 }
377
378 ccl_device_inline float2 as_float2(const float4 a)
379 {
380         return make_float2(a.x, a.y);
381 }
382
383 #endif
384
385 #ifndef __KERNEL_GPU__
386
387 ccl_device_inline void print_float2(const char *label, const float2& a)
388 {
389         printf("%s: %.8f %.8f\n", label, (double)a.x, (double)a.y);
390 }
391
392 #endif
393
394 #ifndef __KERNEL_OPENCL__
395
396 ccl_device_inline float2 interp(float2 a, float2 b, float t)
397 {
398         return a + t*(b - a);
399 }
400
401 #endif
402
403 /* Float3 Vector */
404
405 #ifndef __KERNEL_OPENCL__
406
407 ccl_device_inline float3 operator-(const float3 a)
408 {
409         return make_float3(-a.x, -a.y, -a.z);
410 }
411
412 ccl_device_inline float3 operator*(const float3 a, const float3 b)
413 {
414         return make_float3(a.x*b.x, a.y*b.y, a.z*b.z);
415 }
416
417 ccl_device_inline float3 operator*(const float3 a, float f)
418 {
419         return make_float3(a.x*f, a.y*f, a.z*f);
420 }
421
422 ccl_device_inline float3 operator*(float f, const float3 a)
423 {
424         return make_float3(a.x*f, a.y*f, a.z*f);
425 }
426
427 ccl_device_inline float3 operator/(float f, const float3 a)
428 {
429         return make_float3(f/a.x, f/a.y, f/a.z);
430 }
431
432 ccl_device_inline float3 operator/(const float3 a, float f)
433 {
434         float invf = 1.0f/f;
435         return make_float3(a.x*invf, a.y*invf, a.z*invf);
436 }
437
438 ccl_device_inline float3 operator/(const float3 a, const float3 b)
439 {
440         return make_float3(a.x/b.x, a.y/b.y, a.z/b.z);
441 }
442
443 ccl_device_inline float3 operator+(const float3 a, const float3 b)
444 {
445         return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
446 }
447
448 ccl_device_inline float3 operator-(const float3 a, const float3 b)
449 {
450         return make_float3(a.x-b.x, a.y-b.y, a.z-b.z);
451 }
452
453 ccl_device_inline float3 operator+=(float3& a, const float3 b)
454 {
455         return a = a + b;
456 }
457
458 ccl_device_inline float3 operator*=(float3& a, const float3 b)
459 {
460         return a = a * b;
461 }
462
463 ccl_device_inline float3 operator*=(float3& a, float f)
464 {
465         return a = a * f;
466 }
467
468 ccl_device_inline float3 operator/=(float3& a, const float3 b)
469 {
470         return a = a / b;
471 }
472
473 ccl_device_inline float3 operator/=(float3& a, float f)
474 {
475         float invf = 1.0f/f;
476         return a = a * invf;
477 }
478
479 ccl_device_inline float dot(const float3 a, const float3 b)
480 {
481 #if defined(__KERNEL_SSE41__) && defined(__KERNEL_SSE__)
482         return _mm_cvtss_f32(_mm_dp_ps(a, b, 0x7F));
483 #else   
484         return a.x*b.x + a.y*b.y + a.z*b.z;
485 #endif
486 }
487
488 ccl_device_inline float dot(const float4 a, const float4 b)
489 {
490 #if defined(__KERNEL_SSE41__) && defined(__KERNEL_SSE__)
491         return _mm_cvtss_f32(_mm_dp_ps(a, b, 0xFF));
492 #else   
493         return (a.x*b.x + a.y*b.y) + (a.z*b.z + a.w*b.w);
494 #endif
495 }
496
497 ccl_device_inline float3 cross(const float3 a, const float3 b)
498 {
499         float3 r = make_float3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x);
500         return r;
501 }
502
503 #endif
504
505 ccl_device_inline float len(const float3 a)
506 {
507 #if defined(__KERNEL_SSE41__) && defined(__KERNEL_SSE__)
508         return _mm_cvtss_f32(_mm_sqrt_ss(_mm_dp_ps(a.m128, a.m128, 0x7F)));
509 #else
510         return sqrtf(dot(a, a));
511 #endif
512 }
513
514 ccl_device_inline float len_squared(const float3 a)
515 {
516         return dot(a, a);
517 }
518
519 #ifndef __KERNEL_OPENCL__
520
521 ccl_device_inline float len_squared(const float4 a)
522 {
523         return dot(a, a);
524 }
525
526 ccl_device_inline float3 normalize(const float3 a)
527 {
528 #if defined(__KERNEL_SSE41__) && defined(__KERNEL_SSE__)
529         __m128 norm = _mm_sqrt_ps(_mm_dp_ps(a.m128, a.m128, 0x7F));
530         return _mm_div_ps(a.m128, norm);
531 #else
532         return a/len(a);
533 #endif
534 }
535
536 #endif
537
538 ccl_device_inline float3 normalize_len(const float3 a, float *t)
539 {
540         *t = len(a);
541         return a/(*t);
542 }
543
544 ccl_device_inline float3 safe_normalize(const float3 a)
545 {
546         float t = len(a);
547         return (t)? a/t: a;
548 }
549
550 #ifndef __KERNEL_OPENCL__
551
552 ccl_device_inline bool operator==(const float3 a, const float3 b)
553 {
554 #ifdef __KERNEL_SSE__
555         return (_mm_movemask_ps(_mm_cmpeq_ps(a.m128, b.m128)) & 7) == 7;
556 #else
557         return (a.x == b.x && a.y == b.y && a.z == b.z);
558 #endif
559 }
560
561 ccl_device_inline bool operator!=(const float3 a, const float3 b)
562 {
563         return !(a == b);
564 }
565
566 ccl_device_inline float3 min(float3 a, float3 b)
567 {
568 #ifdef __KERNEL_SSE__
569         return _mm_min_ps(a.m128, b.m128);
570 #else
571         return make_float3(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z));
572 #endif
573 }
574
575 ccl_device_inline float3 max(float3 a, float3 b)
576 {
577 #ifdef __KERNEL_SSE__
578         return _mm_max_ps(a.m128, b.m128);
579 #else
580         return make_float3(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z));
581 #endif
582 }
583
584 ccl_device_inline float3 clamp(float3 a, float3 mn, float3 mx)
585 {
586         return min(max(a, mn), mx);
587 }
588
589 ccl_device_inline float3 fabs(float3 a)
590 {
591 #ifdef __KERNEL_SSE__
592         __m128 mask = _mm_castsi128_ps(_mm_set1_epi32(0x7fffffff));
593         return _mm_and_ps(a.m128, mask);
594 #else
595         return make_float3(fabsf(a.x), fabsf(a.y), fabsf(a.z));
596 #endif
597 }
598
599 #endif
600
601 ccl_device_inline float3 float2_to_float3(const float2 a)
602 {
603         return make_float3(a.x, a.y, 0.0f);
604 }
605
606 ccl_device_inline float3 float4_to_float3(const float4 a)
607 {
608         return make_float3(a.x, a.y, a.z);
609 }
610
611 ccl_device_inline float4 float3_to_float4(const float3 a)
612 {
613         return make_float4(a.x, a.y, a.z, 1.0f);
614 }
615
616 #ifndef __KERNEL_GPU__
617
618 ccl_device_inline void print_float3(const char *label, const float3& a)
619 {
620         printf("%s: %.8f %.8f %.8f\n", label, (double)a.x, (double)a.y, (double)a.z);
621 }
622
623 ccl_device_inline float3 rcp(const float3& a)
624 {
625 #ifdef __KERNEL_SSE__
626         float4 r = _mm_rcp_ps(a.m128);
627         return _mm_sub_ps(_mm_add_ps(r, r), _mm_mul_ps(_mm_mul_ps(r, r), a));
628 #else
629         return make_float3(1.0f/a.x, 1.0f/a.y, 1.0f/a.z);
630 #endif
631 }
632
633 #endif
634
635 ccl_device_inline float3 interp(float3 a, float3 b, float t)
636 {
637         return a + t*(b - a);
638 }
639
640 ccl_device_inline bool is_zero(const float3 a)
641 {
642 #ifdef __KERNEL_SSE__
643         return a == make_float3(0.0f);
644 #else
645         return (a.x == 0.0f && a.y == 0.0f && a.z == 0.0f);
646 #endif
647 }
648
649 ccl_device_inline float reduce_add(const float3 a)
650 {
651         return (a.x + a.y + a.z);
652 }
653
654 ccl_device_inline float average(const float3 a)
655 {
656         return reduce_add(a)*(1.0f/3.0f);
657 }
658
659 /* Float4 Vector */
660
661 #ifdef __KERNEL_SSE__
662
663 template<size_t index_0, size_t index_1, size_t index_2, size_t index_3> __forceinline const float4 shuffle(const float4& b)
664 {
665         return _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(b), _MM_SHUFFLE(index_3, index_2, index_1, index_0)));
666 }
667
668 template<> __forceinline const float4 shuffle<0, 0, 2, 2>(const float4& b)
669 {
670         return _mm_moveldup_ps(b);
671 }
672
673 template<> __forceinline const float4 shuffle<1, 1, 3, 3>(const float4& b)
674 {
675         return _mm_movehdup_ps(b);
676 }
677
678 template<> __forceinline const float4 shuffle<0, 1, 0, 1>(const float4& b)
679 {
680         return _mm_castpd_ps(_mm_movedup_pd(_mm_castps_pd(b)));
681 }
682
683 #endif
684
685 #ifndef __KERNEL_OPENCL__
686
687 ccl_device_inline float4 operator-(const float4& a)
688 {
689 #ifdef __KERNEL_SSE__
690         __m128 mask = _mm_castsi128_ps(_mm_set1_epi32(0x80000000));
691         return _mm_xor_ps(a.m128, mask);
692 #else
693         return make_float4(-a.x, -a.y, -a.z, -a.w);
694 #endif
695 }
696
697 ccl_device_inline float4 operator*(const float4& a, const float4& b)
698 {
699 #ifdef __KERNEL_SSE__
700         return _mm_mul_ps(a.m128, b.m128);
701 #else
702         return make_float4(a.x*b.x, a.y*b.y, a.z*b.z, a.w*b.w);
703 #endif
704 }
705
706 ccl_device_inline float4 operator*(const float4& a, float f)
707 {
708 #ifdef __KERNEL_SSE__
709         return a * make_float4(f);
710 #else
711         return make_float4(a.x*f, a.y*f, a.z*f, a.w*f);
712 #endif
713 }
714
715 ccl_device_inline float4 operator*(float f, const float4& a)
716 {
717         return a * f;
718 }
719
720 ccl_device_inline float4 rcp(const float4& a)
721 {
722 #ifdef __KERNEL_SSE__
723         float4 r = _mm_rcp_ps(a.m128);
724         return _mm_sub_ps(_mm_add_ps(r, r), _mm_mul_ps(_mm_mul_ps(r, r), a));
725 #else
726         return make_float4(1.0f/a.x, 1.0f/a.y, 1.0f/a.z, 1.0f/a.w);
727 #endif
728 }
729
730 ccl_device_inline float4 operator/(const float4& a, float f)
731 {
732         return a * (1.0f/f);
733 }
734
735 ccl_device_inline float4 operator/(const float4& a, const float4& b)
736 {
737 #ifdef __KERNEL_SSE__
738         return a * rcp(b);
739 #else
740         return make_float4(a.x/b.x, a.y/b.y, a.z/b.z, a.w/b.w);
741 #endif
742
743 }
744
745 ccl_device_inline float4 operator+(const float4& a, const float4& b)
746 {
747 #ifdef __KERNEL_SSE__
748         return _mm_add_ps(a.m128, b.m128);
749 #else
750         return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
751 #endif
752 }
753
754 ccl_device_inline float4 operator-(const float4& a, const float4& b)
755 {
756 #ifdef __KERNEL_SSE__
757         return _mm_sub_ps(a.m128, b.m128);
758 #else
759         return make_float4(a.x-b.x, a.y-b.y, a.z-b.z, a.w-b.w);
760 #endif
761 }
762
763 ccl_device_inline float4 operator+=(float4& a, const float4& b)
764 {
765         return a = a + b;
766 }
767
768 ccl_device_inline float4 operator*=(float4& a, const float4& b)
769 {
770         return a = a * b;
771 }
772
773 ccl_device_inline float4 operator/=(float4& a, float f)
774 {
775         return a = a / f;
776 }
777
778 ccl_device_inline int4 operator<(const float4& a, const float4& b)
779 {
780 #ifdef __KERNEL_SSE__
781         return _mm_cvtps_epi32(_mm_cmplt_ps(a.m128, b.m128)); /* todo: avoid cvt */
782 #else
783         return make_int4(a.x < b.x, a.y < b.y, a.z < b.z, a.w < b.w);
784 #endif
785 }
786
787 ccl_device_inline int4 operator>=(float4 a, float4 b)
788 {
789 #ifdef __KERNEL_SSE__
790         return _mm_cvtps_epi32(_mm_cmpge_ps(a.m128, b.m128)); /* todo: avoid cvt */
791 #else
792         return make_int4(a.x >= b.x, a.y >= b.y, a.z >= b.z, a.w >= b.w);
793 #endif
794 }
795
796 ccl_device_inline int4 operator<=(const float4& a, const float4& b)
797 {
798 #ifdef __KERNEL_SSE__
799         return _mm_cvtps_epi32(_mm_cmple_ps(a.m128, b.m128)); /* todo: avoid cvt */
800 #else
801         return make_int4(a.x <= b.x, a.y <= b.y, a.z <= b.z, a.w <= b.w);
802 #endif
803 }
804
805 ccl_device_inline bool operator==(const float4 a, const float4 b)
806 {
807 #ifdef __KERNEL_SSE__
808         return (_mm_movemask_ps(_mm_cmpeq_ps(a.m128, b.m128)) & 15) == 15;
809 #else
810         return (a.x == b.x && a.y == b.y && a.z == b.z && a.w == b.w);
811 #endif
812 }
813
814 ccl_device_inline float4 cross(const float4& a, const float4& b)
815 {
816 #ifdef __KERNEL_SSE__
817         return (shuffle<1,2,0,0>(a)*shuffle<2,0,1,0>(b)) - (shuffle<2,0,1,0>(a)*shuffle<1,2,0,0>(b));
818 #else
819         return make_float4(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x, 0.0f);
820 #endif
821 }
822
823 ccl_device_inline bool is_zero(const float4& a)
824 {
825 #ifdef __KERNEL_SSE__
826         return a == make_float4(0.0f);
827 #else
828         return (a.x == 0.0f && a.y == 0.0f && a.z == 0.0f && a.w == 0.0f);
829 #endif
830 }
831
832 ccl_device_inline float reduce_add(const float4& a)
833 {
834 #ifdef __KERNEL_SSE__
835         float4 h = shuffle<1,0,3,2>(a) + a;
836         return _mm_cvtss_f32(shuffle<2,3,0,1>(h) + h); /* todo: efficiency? */
837 #else
838         return ((a.x + a.y) + (a.z + a.w));
839 #endif
840 }
841
842 ccl_device_inline float average(const float4& a)
843 {
844         return reduce_add(a) * 0.25f;
845 }
846
847 ccl_device_inline float len(const float4 a)
848 {
849         return sqrtf(dot(a, a));
850 }
851
852 ccl_device_inline float4 normalize(const float4 a)
853 {
854         return a/len(a);
855 }
856
857 ccl_device_inline float4 safe_normalize(const float4 a)
858 {
859         float t = len(a);
860         return (t)? a/t: a;
861 }
862
863 ccl_device_inline float4 min(float4 a, float4 b)
864 {
865 #ifdef __KERNEL_SSE__
866         return _mm_min_ps(a.m128, b.m128);
867 #else
868         return make_float4(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z), min(a.w, b.w));
869 #endif
870 }
871
872 ccl_device_inline float4 max(float4 a, float4 b)
873 {
874 #ifdef __KERNEL_SSE__
875         return _mm_max_ps(a.m128, b.m128);
876 #else
877         return make_float4(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z), max(a.w, b.w));
878 #endif
879 }
880
881 #endif
882
883 #ifndef __KERNEL_GPU__
884
885 ccl_device_inline float4 select(const int4& mask, const float4& a, const float4& b)
886 {
887 #ifdef __KERNEL_SSE__
888         return _mm_or_ps(_mm_and_ps(_mm_cvtepi32_ps(mask), a), _mm_andnot_ps(_mm_cvtepi32_ps(mask), b)); /* todo: avoid cvt */
889 #else
890         return make_float4((mask.x)? a.x: b.x, (mask.y)? a.y: b.y, (mask.z)? a.z: b.z, (mask.w)? a.w: b.w);
891 #endif
892 }
893
894 ccl_device_inline float4 reduce_min(const float4& a)
895 {
896 #ifdef __KERNEL_SSE__
897         float4 h = min(shuffle<1,0,3,2>(a), a);
898         return min(shuffle<2,3,0,1>(h), h);
899 #else
900         return make_float4(min(min(a.x, a.y), min(a.z, a.w)));
901 #endif
902 }
903
904 ccl_device_inline float4 reduce_max(const float4& a)
905 {
906 #ifdef __KERNEL_SSE__
907         float4 h = max(shuffle<1,0,3,2>(a), a);
908         return max(shuffle<2,3,0,1>(h), h);
909 #else
910         return make_float4(max(max(a.x, a.y), max(a.z, a.w)));
911 #endif
912 }
913
914 #if 0
915 ccl_device_inline float4 reduce_add(const float4& a)
916 {
917 #ifdef __KERNEL_SSE__
918         float4 h = shuffle<1,0,3,2>(a) + a;
919         return shuffle<2,3,0,1>(h) + h;
920 #else
921         return make_float4((a.x + a.y) + (a.z + a.w));
922 #endif
923 }
924 #endif
925
926 ccl_device_inline void print_float4(const char *label, const float4& a)
927 {
928         printf("%s: %.8f %.8f %.8f %.8f\n", label, (double)a.x, (double)a.y, (double)a.z, (double)a.w);
929 }
930
931 #endif
932
933 /* Int3 */
934
935 #ifndef __KERNEL_OPENCL__
936
937 ccl_device_inline int3 min(int3 a, int3 b)
938 {
939 #if defined(__KERNEL_SSE__) && defined(__KERNEL_SSE41__)
940         return _mm_min_epi32(a.m128, b.m128);
941 #else
942         return make_int3(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z));
943 #endif
944 }
945
946 ccl_device_inline int3 max(int3 a, int3 b)
947 {
948 #if defined(__KERNEL_SSE__) && defined(__KERNEL_SSE41__)
949         return _mm_max_epi32(a.m128, b.m128);
950 #else
951         return make_int3(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z));
952 #endif
953 }
954
955 ccl_device_inline int3 clamp(const int3& a, int mn, int mx)
956 {
957 #ifdef __KERNEL_SSE__
958         return min(max(a, make_int3(mn)), make_int3(mx));
959 #else
960         return make_int3(clamp(a.x, mn, mx), clamp(a.y, mn, mx), clamp(a.z, mn, mx));
961 #endif
962 }
963
964 ccl_device_inline int3 clamp(const int3& a, int3& mn, int mx)
965 {
966 #ifdef __KERNEL_SSE__
967         return min(max(a, mn), make_int3(mx));
968 #else
969         return make_int3(clamp(a.x, mn.x, mx), clamp(a.y, mn.y, mx), clamp(a.z, mn.z, mx));
970 #endif
971 }
972
973 #endif
974
975 #ifndef __KERNEL_GPU__
976
977 ccl_device_inline void print_int3(const char *label, const int3& a)
978 {
979         printf("%s: %d %d %d\n", label, a.x, a.y, a.z);
980 }
981
982 #endif
983
984 /* Int4 */
985
986 #ifndef __KERNEL_GPU__
987
988 ccl_device_inline int4 operator+(const int4& a, const int4& b)
989 {
990 #ifdef __KERNEL_SSE__
991         return _mm_add_epi32(a.m128, b.m128);
992 #else
993         return make_int4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
994 #endif
995 }
996
997 ccl_device_inline int4 operator+=(int4& a, const int4& b)
998 {
999         return a = a + b;
1000 }
1001
1002 ccl_device_inline int4 operator>>(const int4& a, int i)
1003 {
1004 #ifdef __KERNEL_SSE__
1005         return _mm_srai_epi32(a.m128, i);
1006 #else
1007         return make_int4(a.x >> i, a.y >> i, a.z >> i, a.w >> i);
1008 #endif
1009 }
1010
1011 ccl_device_inline int4 min(int4 a, int4 b)
1012 {
1013 #if defined(__KERNEL_SSE__) && defined(__KERNEL_SSE41__)
1014         return _mm_min_epi32(a.m128, b.m128);
1015 #else
1016         return make_int4(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z), min(a.w, b.w));
1017 #endif
1018 }
1019
1020 ccl_device_inline int4 max(int4 a, int4 b)
1021 {
1022 #if defined(__KERNEL_SSE__) && defined(__KERNEL_SSE41__)
1023         return _mm_max_epi32(a.m128, b.m128);
1024 #else
1025         return make_int4(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z), max(a.w, b.w));
1026 #endif
1027 }
1028
1029 ccl_device_inline int4 clamp(const int4& a, const int4& mn, const int4& mx)
1030 {
1031         return min(max(a, mn), mx);
1032 }
1033
1034 ccl_device_inline int4 select(const int4& mask, const int4& a, const int4& b)
1035 {
1036 #ifdef __KERNEL_SSE__
1037         __m128 m = _mm_cvtepi32_ps(mask);
1038         return _mm_castps_si128(_mm_or_ps(_mm_and_ps(m, _mm_castsi128_ps(a)), _mm_andnot_ps(m, _mm_castsi128_ps(b)))); /* todo: avoid cvt */
1039 #else
1040         return make_int4((mask.x)? a.x: b.x, (mask.y)? a.y: b.y, (mask.z)? a.z: b.z, (mask.w)? a.w: b.w);
1041 #endif
1042 }
1043
1044 ccl_device_inline void print_int4(const char *label, const int4& a)
1045 {
1046         printf("%s: %d %d %d %d\n", label, a.x, a.y, a.z, a.w);
1047 }
1048
1049 #endif
1050
1051 /* Int/Float conversion */
1052
1053 #ifndef __KERNEL_OPENCL__
1054
1055 ccl_device_inline int as_int(uint i)
1056 {
1057         union { uint ui; int i; } u;
1058         u.ui = i;
1059         return u.i;
1060 }
1061
1062 ccl_device_inline uint as_uint(int i)
1063 {
1064         union { uint ui; int i; } u;
1065         u.i = i;
1066         return u.ui;
1067 }
1068
1069 ccl_device_inline uint as_uint(float f)
1070 {
1071         union { uint i; float f; } u;
1072         u.f = f;
1073         return u.i;
1074 }
1075
1076 ccl_device_inline int __float_as_int(float f)
1077 {
1078         union { int i; float f; } u;
1079         u.f = f;
1080         return u.i;
1081 }
1082
1083 ccl_device_inline float __int_as_float(int i)
1084 {
1085         union { int i; float f; } u;
1086         u.i = i;
1087         return u.f;
1088 }
1089
1090 ccl_device_inline uint __float_as_uint(float f)
1091 {
1092         union { uint i; float f; } u;
1093         u.f = f;
1094         return u.i;
1095 }
1096
1097 ccl_device_inline float __uint_as_float(uint i)
1098 {
1099         union { uint i; float f; } u;
1100         u.i = i;
1101         return u.f;
1102 }
1103
1104 /* Interpolation */
1105
1106 template<class A, class B> A lerp(const A& a, const A& b, const B& t)
1107 {
1108         return (A)(a * ((B)1 - t) + b * t);
1109 }
1110
1111 /* Triangle */
1112
1113 ccl_device_inline float triangle_area(const float3 v1, const float3 v2, const float3 v3)
1114 {
1115         return len(cross(v3 - v2, v1 - v2))*0.5f;
1116 }
1117
1118 #endif
1119
1120 /* Orthonormal vectors */
1121
1122 ccl_device_inline void make_orthonormals(const float3 N, float3 *a, float3 *b)
1123 {
1124 #if 0
1125         if(fabsf(N.y) >= 0.999f) {
1126                 *a = make_float3(1, 0, 0);
1127                 *b = make_float3(0, 0, 1);
1128                 return;
1129         }
1130         if(fabsf(N.z) >= 0.999f) {
1131                 *a = make_float3(1, 0, 0);
1132                 *b = make_float3(0, 1, 0);
1133                 return;
1134         }
1135 #endif
1136
1137         if(N.x != N.y || N.x != N.z)
1138                 *a = make_float3(N.z-N.y, N.x-N.z, N.y-N.x);  //(1,1,1)x N
1139         else
1140                 *a = make_float3(N.z-N.y, N.x+N.z, -N.y-N.x);  //(-1,1,1)x N
1141
1142         *a = normalize(*a);
1143         *b = cross(N, *a);
1144 }
1145
1146 /* Color division */
1147
1148 ccl_device_inline float3 safe_invert_color(float3 a)
1149 {
1150         float x, y, z;
1151
1152         x = (a.x != 0.0f)? 1.0f/a.x: 0.0f;
1153         y = (a.y != 0.0f)? 1.0f/a.y: 0.0f;
1154         z = (a.z != 0.0f)? 1.0f/a.z: 0.0f;
1155
1156         return make_float3(x, y, z);
1157 }
1158
1159 ccl_device_inline float3 safe_divide_color(float3 a, float3 b)
1160 {
1161         float x, y, z;
1162
1163         x = (b.x != 0.0f)? a.x/b.x: 0.0f;
1164         y = (b.y != 0.0f)? a.y/b.y: 0.0f;
1165         z = (b.z != 0.0f)? a.z/b.z: 0.0f;
1166
1167         return make_float3(x, y, z);
1168 }
1169
1170 ccl_device_inline float3 safe_divide_even_color(float3 a, float3 b)
1171 {
1172         float x, y, z;
1173
1174         x = (b.x != 0.0f)? a.x/b.x: 0.0f;
1175         y = (b.y != 0.0f)? a.y/b.y: 0.0f;
1176         z = (b.z != 0.0f)? a.z/b.z: 0.0f;
1177
1178         /* try to get grey even if b is zero */
1179         if(b.x == 0.0f) {
1180                 if(b.y == 0.0f) {
1181                         x = z;
1182                         y = z;
1183                 }
1184                 else if(b.z == 0.0f) {
1185                         x = y;
1186                         z = y;
1187                 }
1188                 else
1189                         x = 0.5f*(y + z);
1190         }
1191         else if(b.y == 0.0f) {
1192                 if(b.z == 0.0f) {
1193                         y = x;
1194                         z = x;
1195                 }
1196                 else
1197                         y = 0.5f*(x + z);
1198         }
1199         else if(b.z == 0.0f) {
1200                 z = 0.5f*(x + y);
1201         }
1202
1203         return make_float3(x, y, z);
1204 }
1205
1206 /* Rotation of point around axis and angle */
1207
1208 ccl_device_inline float3 rotate_around_axis(float3 p, float3 axis, float angle)
1209 {
1210         float costheta = cosf(angle);
1211         float sintheta = sinf(angle);
1212         float3 r;
1213
1214         r.x = ((costheta + (1 - costheta) * axis.x * axis.x) * p.x) +
1215                 (((1 - costheta) * axis.x * axis.y - axis.z * sintheta) * p.y) +
1216                 (((1 - costheta) * axis.x * axis.z + axis.y * sintheta) * p.z);
1217
1218         r.y = (((1 - costheta) * axis.x * axis.y + axis.z * sintheta) * p.x) +
1219                 ((costheta + (1 - costheta) * axis.y * axis.y) * p.y) +
1220                 (((1 - costheta) * axis.y * axis.z - axis.x * sintheta) * p.z);
1221
1222         r.z = (((1 - costheta) * axis.x * axis.z - axis.y * sintheta) * p.x) +
1223                 (((1 - costheta) * axis.y * axis.z + axis.x * sintheta) * p.y) +
1224                 ((costheta + (1 - costheta) * axis.z * axis.z) * p.z);
1225
1226         return r;
1227 }
1228
1229 /* NaN-safe math ops */
1230
1231 ccl_device_inline float safe_sqrtf(float f)
1232 {
1233         return sqrtf(max(f, 0.0f));
1234 }
1235
1236 ccl_device float safe_asinf(float a)
1237 {
1238         return asinf(clamp(a, -1.0f, 1.0f));
1239 }
1240
1241 ccl_device float safe_acosf(float a)
1242 {
1243         return acosf(clamp(a, -1.0f, 1.0f));
1244 }
1245
1246 ccl_device float compatible_powf(float x, float y)
1247 {
1248 #ifdef __KERNEL_GPU__
1249         if(y == 0.0f) /* x^0 -> 1, including 0^0 */
1250                 return 1.0f;
1251
1252         /* GPU pow doesn't accept negative x, do manual checks here */
1253         if(x < 0.0f) {
1254                 if(fmodf(-y, 2.0f) == 0.0f)
1255                         return powf(-x, y);
1256                 else
1257                         return -powf(-x, y);
1258         }
1259         else if(x == 0.0f)
1260                 return 0.0f;
1261 #endif
1262         return powf(x, y);
1263 }
1264
1265 ccl_device float safe_powf(float a, float b)
1266 {
1267         if(UNLIKELY(a < 0.0f && b != float_to_int(b)))
1268                 return 0.0f;
1269
1270         return compatible_powf(a, b);
1271 }
1272
1273 ccl_device float safe_logf(float a, float b)
1274 {
1275         if(UNLIKELY(a < 0.0f || b < 0.0f))
1276                 return 0.0f;
1277
1278         return logf(a)/logf(b);
1279 }
1280
1281 ccl_device float safe_divide(float a, float b)
1282 {
1283         return (b != 0.0f)? a/b: 0.0f;
1284 }
1285
1286 ccl_device float safe_modulo(float a, float b)
1287 {
1288         return (b != 0.0f)? fmodf(a, b): 0.0f;
1289 }
1290
1291 /* Ray Intersection */
1292
1293 ccl_device bool ray_sphere_intersect(
1294         float3 ray_P, float3 ray_D, float ray_t,
1295         float3 sphere_P, float sphere_radius,
1296         float3 *isect_P, float *isect_t)
1297 {
1298         float3 d = sphere_P - ray_P;
1299         float radiussq = sphere_radius*sphere_radius;
1300         float tsq = dot(d, d);
1301
1302         if(tsq > radiussq) { /* ray origin outside sphere */
1303                 float tp = dot(d, ray_D);
1304
1305                 if(tp < 0.0f) /* dir points away from sphere */
1306                         return false;
1307
1308                 float dsq = tsq - tp*tp; /* pythagoras */
1309
1310                 if(dsq > radiussq) /* closest point on ray outside sphere */
1311                         return false;
1312
1313                 float t = tp - sqrtf(radiussq - dsq); /* pythagoras */
1314
1315                 if(t < ray_t) {
1316                         *isect_t = t;
1317                         *isect_P = ray_P + ray_D*t;
1318                         return true;
1319                 }
1320         }
1321
1322         return false;
1323 }
1324
1325 ccl_device bool ray_aligned_disk_intersect(
1326         float3 ray_P, float3 ray_D, float ray_t,
1327         float3 disk_P, float disk_radius,
1328         float3 *isect_P, float *isect_t)
1329 {
1330         /* aligned disk normal */
1331         float disk_t;
1332         float3 disk_N = normalize_len(ray_P - disk_P, &disk_t);
1333         float div = dot(ray_D, disk_N);
1334
1335         if(UNLIKELY(div == 0.0f))
1336                 return false;
1337
1338         /* compute t to intersection point */
1339         float t = -disk_t/div;
1340         if(t < 0.0f || t > ray_t)
1341                 return false;
1342         
1343         /* test if within radius */
1344         float3 P = ray_P + ray_D*t;
1345         if(len_squared(P - disk_P) > disk_radius*disk_radius)
1346                 return false;
1347
1348         *isect_P = P;
1349         *isect_t = t;
1350
1351         return true;
1352 }
1353
1354 ccl_device bool ray_triangle_intersect(
1355         float3 ray_P, float3 ray_D, float ray_t,
1356         float3 v0, float3 v1, float3 v2,
1357         float3 *isect_P, float *isect_t)
1358 {
1359         /* Calculate intersection */
1360         float3 e1 = v1 - v0;
1361         float3 e2 = v2 - v0;
1362         float3 s1 = cross(ray_D, e2);
1363
1364         const float divisor = dot(s1, e1);
1365         if(UNLIKELY(divisor == 0.0f))
1366                 return false;
1367
1368         const float invdivisor = 1.0f/divisor;
1369
1370         /* compute first barycentric coordinate */
1371         const float3 d = ray_P - v0;
1372         const float u = dot(d, s1)*invdivisor;
1373         if(u < 0.0f)
1374                 return false;
1375
1376         /* Compute second barycentric coordinate */
1377         const float3 s2 = cross(d, e1);
1378         const float v = dot(ray_D, s2)*invdivisor;
1379         if(v < 0.0f)
1380                 return false;
1381
1382         const float b0 = 1.0f - u - v;
1383         if(b0 < 0.0f)
1384                 return false;
1385
1386         /* compute t to intersection point */
1387         const float t = dot(e2, s2)*invdivisor;
1388         if(t < 0.0f || t > ray_t)
1389                 return false;
1390
1391         *isect_t = t;
1392         *isect_P = ray_P + ray_D*t;
1393
1394         return true;
1395 }
1396
1397 ccl_device bool ray_triangle_intersect_uv(
1398         float3 ray_P, float3 ray_D, float ray_t,
1399         float3 v0, float3 v1, float3 v2,
1400         float *isect_u, float *isect_v, float *isect_t)
1401 {
1402         /* Calculate intersection */
1403         float3 e1 = v1 - v0;
1404         float3 e2 = v2 - v0;
1405         float3 s1 = cross(ray_D, e2);
1406
1407         const float divisor = dot(s1, e1);
1408         if(UNLIKELY(divisor == 0.0f))
1409                 return false;
1410
1411         const float invdivisor = 1.0f/divisor;
1412
1413         /* compute first barycentric coordinate */
1414         const float3 d = ray_P - v0;
1415         const float u = dot(d, s1)*invdivisor;
1416         if(u < 0.0f)
1417                 return false;
1418
1419         /* Compute second barycentric coordinate */
1420         const float3 s2 = cross(d, e1);
1421         const float v = dot(ray_D, s2)*invdivisor;
1422         if(v < 0.0f)
1423                 return false;
1424
1425         const float b0 = 1.0f - u - v;
1426         if(b0 < 0.0f)
1427                 return false;
1428
1429         /* compute t to intersection point */
1430         const float t = dot(e2, s2)*invdivisor;
1431         if(t < 0.0f || t > ray_t)
1432                 return false;
1433
1434         *isect_u = u;
1435         *isect_v = v;
1436         *isect_t = t;
1437
1438         return true;
1439 }
1440
1441 ccl_device bool ray_quad_intersect(
1442         float3 ray_P, float3 ray_D, float ray_t,
1443         float3 quad_P, float3 quad_u, float3 quad_v,
1444         float3 *isect_P, float *isect_t)
1445 {
1446         float3 v0 = quad_P - quad_u*0.5f - quad_v*0.5f;
1447         float3 v1 = quad_P + quad_u*0.5f - quad_v*0.5f;
1448         float3 v2 = quad_P + quad_u*0.5f + quad_v*0.5f;
1449         float3 v3 = quad_P - quad_u*0.5f + quad_v*0.5f;
1450
1451         if(ray_triangle_intersect(ray_P, ray_D, ray_t, v0, v1, v2, isect_P, isect_t))
1452                 return true;
1453         else if(ray_triangle_intersect(ray_P, ray_D, ray_t, v0, v2, v3, isect_P, isect_t))
1454                 return true;
1455         
1456         return false;
1457 }
1458
1459 /* projections */
1460 ccl_device_inline float2 map_to_tube(const float3 co)
1461 {
1462         float len, u, v;
1463         len = sqrtf(co.x * co.x + co.y * co.y);
1464         if(len > 0.0f) {
1465                 u = (1.0f - (atan2f(co.x / len, co.y / len) / M_PI_F)) * 0.5f;
1466                 v = (co.x + 1.0f) * 0.5f;
1467         }
1468         else {
1469                 u = v = 0.0f;
1470         }
1471         return make_float2(u, v);
1472 }
1473
1474 ccl_device_inline float2 map_to_sphere(const float3 co)
1475 {
1476         float l = len(co);
1477         float u, v;
1478         if(l > 0.0f) {
1479                 if(UNLIKELY(co.x == 0.0f && co.y == 0.0f)) {
1480                         u = 0.0f;  /* othwise domain error */
1481                 }
1482                 else {
1483                         u = (1.0f - atan2f(co.x, co.y) / M_PI_F) / 2.0f;
1484                 }
1485                 v = 1.0f - safe_acosf(co.z / l) / M_PI_F;
1486         }
1487         else {
1488                 u = v = 0.0f;
1489         }
1490         return make_float2(u, v);
1491 }
1492
1493 ccl_device_inline int util_max_axis(float3 vec)
1494 {
1495         if(vec.x > vec.y) {
1496                 if(vec.x > vec.z)
1497                         return 0;
1498                 else
1499                         return 2;
1500         }
1501         else {
1502                 if(vec.y > vec.z)
1503                         return 1;
1504                 else
1505                         return 2;
1506         }
1507 }
1508
1509 CCL_NAMESPACE_END
1510
1511 #endif /* __UTIL_MATH_H__ */
1512