Cycles: merge of changes from tomato branch.
[blender-staging.git] / intern / cycles / util / util_math.h
1 /*
2  * Copyright 2011, Blender Foundation.
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
19 #ifndef __UTIL_MATH_H__
20 #define __UTIL_MATH_H__
21
22 /* Math
23  *
24  * Basic math functions on scalar and vector types. This header is used by
25  * both the kernel code when compiled as C++, and other C++ non-kernel code. */
26
27 #ifndef __KERNEL_OPENCL__
28
29 #define _USE_MATH_DEFINES
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 #ifndef M_PI_F
44 #define M_PI_F          ((float)3.14159265358979323846264338327950288)
45 #endif
46 #ifndef M_PI_2_F
47 #define M_PI_2_F        ((float)1.57079632679489661923132169163975144)
48 #endif
49 #ifndef M_PI_4_F
50 #define M_PI_4_F        ((float)0.785398163397448309615660845819875721)
51 #endif
52 #ifndef M_1_PI_F
53 #define M_1_PI_F        ((float)0.318309886183790671537767526745028724)
54 #endif
55 #ifndef M_2_PI_F
56 #define M_2_PI_F        ((float)0.636619772367581343075535053490057448)
57 #endif
58 #ifndef M_SQRT2_F
59 #define M_SQRT2_F       ((float)1.41421356237309504880)
60 #endif
61
62
63 /* Scalar */
64
65 #ifdef _WIN32
66
67 #ifndef __KERNEL_GPU__
68
69 #if(!defined(FREE_WINDOWS))
70 #define copysignf(x, y) ((float)_copysign(x, y))
71 #define hypotf(x, y) _hypotf(x, y)
72 #define isnan(x) _isnan(x)
73 #define isfinite(x) _finite(x)
74 #endif
75
76 #endif
77
78 #ifndef __KERNEL_OPENCL__
79
80 __device_inline float fmaxf(float a, float b)
81 {
82         return (a > b)? a: b;
83 }
84
85 __device_inline float fminf(float a, float b)
86 {
87         return (a < b)? a: b;
88 }
89
90 #endif
91
92 #endif
93
94 #ifndef __KERNEL_GPU__
95
96 __device_inline int max(int a, int b)
97 {
98         return (a > b)? a: b;
99 }
100
101 __device_inline int min(int a, int b)
102 {
103         return (a < b)? a: b;
104 }
105
106 __device_inline float max(float a, float b)
107 {
108         return (a > b)? a: b;
109 }
110
111 __device_inline float min(float a, float b)
112 {
113         return (a < b)? a: b;
114 }
115
116 __device_inline double max(double a, double b)
117 {
118         return (a > b)? a: b;
119 }
120
121 __device_inline double min(double a, double b)
122 {
123         return (a < b)? a: b;
124 }
125
126 #endif
127
128 __device_inline float min4(float a, float b, float c, float d)
129 {
130         return min(min(a, b), min(c, d));
131 }
132
133 __device_inline float max4(float a, float b, float c, float d)
134 {
135         return max(max(a, b), max(c, d));
136 }
137
138 #ifndef __KERNEL_OPENCL__
139
140 __device_inline int clamp(int a, int mn, int mx)
141 {
142         return min(max(a, mn), mx);
143 }
144
145 __device_inline float clamp(float a, float mn, float mx)
146 {
147         return min(max(a, mn), mx);
148 }
149
150 #endif
151
152 __device_inline float signf(float f)
153 {
154         return (f < 0.0f)? -1.0f: 1.0f;
155 }
156
157 __device_inline float nonzerof(float f, float eps)
158 {
159         if(fabsf(f) < eps)
160                 return signf(f)*eps;
161         else
162                 return f;
163 }
164
165 __device_inline float smoothstepf(float f)
166 {
167         float ff = f*f;
168         return (3.0f*ff - 2.0f*ff*f);
169 }
170
171 /* Float2 Vector */
172
173 #ifndef __KERNEL_OPENCL__
174
175 __device_inline bool is_zero(const float2 a)
176 {
177         return (a.x == 0.0f && a.y == 0.0f);
178 }
179
180 #endif
181
182 #ifndef __KERNEL_OPENCL__
183
184 __device_inline float average(const float2 a)
185 {
186         return (a.x + a.y)*(1.0f/2.0f);
187 }
188
189 #endif
190
191 #ifndef __KERNEL_OPENCL__
192
193 __device_inline float2 operator-(const float2 a)
194 {
195         return make_float2(-a.x, -a.y);
196 }
197
198 __device_inline float2 operator*(const float2 a, const float2 b)
199 {
200         return make_float2(a.x*b.x, a.y*b.y);
201 }
202
203 __device_inline float2 operator*(const float2 a, float f)
204 {
205         return make_float2(a.x*f, a.y*f);
206 }
207
208 __device_inline float2 operator*(float f, const float2 a)
209 {
210         return make_float2(a.x*f, a.y*f);
211 }
212
213 __device_inline float2 operator/(float f, const float2 a)
214 {
215         return make_float2(f/a.x, f/a.y);
216 }
217
218 __device_inline float2 operator/(const float2 a, float f)
219 {
220         float invf = 1.0f/f;
221         return make_float2(a.x*invf, a.y*invf);
222 }
223
224 __device_inline float2 operator/(const float2 a, const float2 b)
225 {
226         return make_float2(a.x/b.x, a.y/b.y);
227 }
228
229 __device_inline float2 operator+(const float2 a, const float2 b)
230 {
231         return make_float2(a.x+b.x, a.y+b.y);
232 }
233
234 __device_inline float2 operator-(const float2 a, const float2 b)
235 {
236         return make_float2(a.x-b.x, a.y-b.y);
237 }
238
239 __device_inline float2 operator+=(float2& a, const float2 b)
240 {
241         return a = a + b;
242 }
243
244 __device_inline float2 operator*=(float2& a, const float2 b)
245 {
246         return a = a * b;
247 }
248
249 __device_inline float2 operator*=(float2& a, float f)
250 {
251         return a = a * f;
252 }
253
254 __device_inline float2 operator/=(float2& a, const float2 b)
255 {
256         return a = a / b;
257 }
258
259 __device_inline float2 operator/=(float2& a, float f)
260 {
261         float invf = 1.0f/f;
262         return a = a * invf;
263 }
264
265
266 __device_inline float dot(const float2 a, const float2 b)
267 {
268         return a.x*b.x + a.y*b.y;
269 }
270
271 __device_inline float cross(const float2 a, const float2 b)
272 {
273         return (a.x*b.y - a.y*b.x);
274 }
275
276 #endif
277
278 #ifndef __KERNEL_OPENCL__
279
280 __device_inline bool operator==(const int2 a, const int2 b)
281 {
282         return (a.x == b.x && a.y == b.y);
283 }
284
285 __device_inline float len(const float2 a)
286 {
287         return sqrtf(dot(a, a));
288 }
289
290 __device_inline float2 normalize(const float2 a)
291 {
292         return a/len(a);
293 }
294
295 __device_inline float2 normalize_len(const float2 a, float *t)
296 {
297         *t = len(a);
298         return a/(*t);
299 }
300
301 __device_inline bool operator==(const float2 a, const float2 b)
302 {
303         return (a.x == b.x && a.y == b.y);
304 }
305
306 __device_inline bool operator!=(const float2 a, const float2 b)
307 {
308         return !(a == b);
309 }
310
311 __device_inline float2 min(float2 a, float2 b)
312 {
313         return make_float2(min(a.x, b.x), min(a.y, b.y));
314 }
315
316 __device_inline float2 max(float2 a, float2 b)
317 {
318         return make_float2(max(a.x, b.x), max(a.y, b.y));
319 }
320
321 __device_inline float2 clamp(float2 a, float2 mn, float2 mx)
322 {
323         return min(max(a, mn), mx);
324 }
325
326 __device_inline float2 fabs(float2 a)
327 {
328         return make_float2(fabsf(a.x), fabsf(a.y));
329 }
330
331 __device_inline float2 as_float2(const float4 a)
332 {
333         return make_float2(a.x, a.y);
334 }
335
336 #endif
337
338 #ifndef __KERNEL_GPU__
339
340 __device_inline void print_float2(const char *label, const float2& a)
341 {
342         printf("%s: %.8f %.8f\n", label, (double)a.x, (double)a.y);
343 }
344
345 #endif
346
347 #ifndef __KERNEL_OPENCL__
348
349 __device_inline float2 interp(float2 a, float2 b, float t)
350 {
351         return a + t*(b - a);
352 }
353
354 #endif
355
356 /* Float3 Vector */
357
358 #ifndef __KERNEL_OPENCL__
359
360 __device_inline float3 operator-(const float3 a)
361 {
362         return make_float3(-a.x, -a.y, -a.z);
363 }
364
365 __device_inline float3 operator*(const float3 a, const float3 b)
366 {
367         return make_float3(a.x*b.x, a.y*b.y, a.z*b.z);
368 }
369
370 __device_inline float3 operator*(const float3 a, float f)
371 {
372         return make_float3(a.x*f, a.y*f, a.z*f);
373 }
374
375 __device_inline float3 operator*(float f, const float3 a)
376 {
377         return make_float3(a.x*f, a.y*f, a.z*f);
378 }
379
380 __device_inline float3 operator/(float f, const float3 a)
381 {
382         return make_float3(f/a.x, f/a.y, f/a.z);
383 }
384
385 __device_inline float3 operator/(const float3 a, float f)
386 {
387         float invf = 1.0f/f;
388         return make_float3(a.x*invf, a.y*invf, a.z*invf);
389 }
390
391 __device_inline float3 operator/(const float3 a, const float3 b)
392 {
393         return make_float3(a.x/b.x, a.y/b.y, a.z/b.z);
394 }
395
396 __device_inline float3 operator+(const float3 a, const float3 b)
397 {
398         return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
399 }
400
401 __device_inline float3 operator-(const float3 a, const float3 b)
402 {
403         return make_float3(a.x-b.x, a.y-b.y, a.z-b.z);
404 }
405
406 __device_inline float3 operator+=(float3& a, const float3 b)
407 {
408         return a = a + b;
409 }
410
411 __device_inline float3 operator*=(float3& a, const float3 b)
412 {
413         return a = a * b;
414 }
415
416 __device_inline float3 operator*=(float3& a, float f)
417 {
418         return a = a * f;
419 }
420
421 __device_inline float3 operator/=(float3& a, const float3 b)
422 {
423         return a = a / b;
424 }
425
426 __device_inline float3 operator/=(float3& a, float f)
427 {
428         float invf = 1.0f/f;
429         return a = a * invf;
430 }
431
432 __device_inline float dot(const float3 a, const float3 b)
433 {
434         return a.x*b.x + a.y*b.y + a.z*b.z;
435 }
436
437 __device_inline float3 cross(const float3 a, const float3 b)
438 {
439         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);
440         return r;
441 }
442
443 #endif
444
445 __device_inline float len(const float3 a)
446 {
447         return sqrtf(dot(a, a));
448 }
449
450 __device_inline float len_squared(const float3 a)
451 {
452         return dot(a, a);
453 }
454
455 #ifndef __KERNEL_OPENCL__
456
457 __device_inline float3 normalize(const float3 a)
458 {
459         return a/len(a);
460 }
461
462 #endif
463
464 __device_inline float3 normalize_len(const float3 a, float *t)
465 {
466         *t = len(a);
467         return a/(*t);
468 }
469
470 #ifndef __KERNEL_OPENCL__
471
472 __device_inline bool operator==(const float3 a, const float3 b)
473 {
474 #ifdef __KERNEL_SSE__
475         return (_mm_movemask_ps(_mm_cmpeq_ps(a.m128, b.m128)) & 7) == 7;
476 #else
477         return (a.x == b.x && a.y == b.y && a.z == b.z);
478 #endif
479 }
480
481 __device_inline bool operator!=(const float3 a, const float3 b)
482 {
483         return !(a == b);
484 }
485
486 __device_inline float3 min(float3 a, float3 b)
487 {
488 #ifdef __KERNEL_SSE__
489         return _mm_min_ps(a.m128, b.m128);
490 #else
491         return make_float3(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z));
492 #endif
493 }
494
495 __device_inline float3 max(float3 a, float3 b)
496 {
497 #ifdef __KERNEL_SSE__
498         return _mm_max_ps(a.m128, b.m128);
499 #else
500         return make_float3(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z));
501 #endif
502 }
503
504 __device_inline float3 clamp(float3 a, float3 mn, float3 mx)
505 {
506         return min(max(a, mn), mx);
507 }
508
509 __device_inline float3 fabs(float3 a)
510 {
511 #ifdef __KERNEL_SSE__
512         __m128 mask = _mm_castsi128_ps(_mm_set1_epi32(0x7fffffff));
513         return _mm_and_ps(a.m128, mask);
514 #else
515         return make_float3(fabsf(a.x), fabsf(a.y), fabsf(a.z));
516 #endif
517 }
518
519 #endif
520
521 __device_inline float3 float2_to_float3(const float2 a)
522 {
523         return make_float3(a.x, a.y, 0.0f);
524 }
525
526 __device_inline float3 float4_to_float3(const float4 a)
527 {
528         return make_float3(a.x, a.y, a.z);
529 }
530
531 __device_inline float4 float3_to_float4(const float3 a)
532 {
533         return make_float4(a.x, a.y, a.z, 1.0f);
534 }
535
536 #ifndef __KERNEL_GPU__
537
538 __device_inline void print_float3(const char *label, const float3& a)
539 {
540         printf("%s: %.8f %.8f %.8f\n", label, (double)a.x, (double)a.y, (double)a.z);
541 }
542
543 __device_inline float3 rcp(const float3& a)
544 {
545 #ifdef __KERNEL_SSE__
546         float4 r = _mm_rcp_ps(a.m128);
547         return _mm_sub_ps(_mm_add_ps(r, r), _mm_mul_ps(_mm_mul_ps(r, r), a));
548 #else
549         return make_float3(1.0f/a.x, 1.0f/a.y, 1.0f/a.z);
550 #endif
551 }
552
553 #endif
554
555 __device_inline float3 interp(float3 a, float3 b, float t)
556 {
557         return a + t*(b - a);
558 }
559
560 __device_inline bool is_zero(const float3 a)
561 {
562 #ifdef __KERNEL_SSE__
563         return a == make_float3(0.0f);
564 #else
565         return (a.x == 0.0f && a.y == 0.0f && a.z == 0.0f);
566 #endif
567 }
568
569 __device_inline float reduce_add(const float3 a)
570 {
571 #ifdef __KERNEL_SSE__
572         return (a.x + a.y + a.z);
573 #else
574         return (a.x + a.y + a.z);
575 #endif
576 }
577
578 __device_inline float average(const float3 a)
579 {
580         return reduce_add(a)*(1.0f/3.0f);
581 }
582
583 /* Float4 Vector */
584
585 #ifdef __KERNEL_SSE__
586
587 template<size_t index_0, size_t index_1, size_t index_2, size_t index_3> __forceinline const float4 shuffle(const float4& b)
588 {
589         return _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(b), _MM_SHUFFLE(index_3, index_2, index_1, index_0)));
590 }
591
592 template<> __forceinline const float4 shuffle<0, 0, 2, 2>(const float4& b)
593 {
594         return _mm_moveldup_ps(b);
595 }
596
597 template<> __forceinline const float4 shuffle<1, 1, 3, 3>(const float4& b)
598 {
599         return _mm_movehdup_ps(b);
600 }
601
602 template<> __forceinline const float4 shuffle<0, 1, 0, 1>(const float4& b)
603 {
604         return _mm_castpd_ps(_mm_movedup_pd(_mm_castps_pd(b)));
605 }
606
607 #endif
608
609 #ifndef __KERNEL_OPENCL__
610
611 __device_inline float4 operator-(const float4& a)
612 {
613 #ifdef __KERNEL_SSE__
614         __m128 mask = _mm_castsi128_ps(_mm_set1_epi32(0x80000000));
615         return _mm_xor_ps(a.m128, mask);
616 #else
617         return make_float4(-a.x, -a.y, -a.z, -a.w);
618 #endif
619 }
620
621 __device_inline float4 operator*(const float4& a, const float4& b)
622 {
623 #ifdef __KERNEL_SSE__
624         return _mm_mul_ps(a.m128, b.m128);
625 #else
626         return make_float4(a.x*b.x, a.y*b.y, a.z*b.z, a.w*b.w);
627 #endif
628 }
629
630 __device_inline float4 operator*(const float4& a, float f)
631 {
632 #ifdef __KERNEL_SSE__
633         return a * make_float4(f);
634 #else
635         return make_float4(a.x*f, a.y*f, a.z*f, a.w*f);
636 #endif
637 }
638
639 __device_inline float4 operator*(float f, const float4& a)
640 {
641         return a * f;
642 }
643
644 __device_inline float4 rcp(const float4& a)
645 {
646 #ifdef __KERNEL_SSE__
647         float4 r = _mm_rcp_ps(a.m128);
648         return _mm_sub_ps(_mm_add_ps(r, r), _mm_mul_ps(_mm_mul_ps(r, r), a));
649 #else
650         return make_float4(1.0f/a.x, 1.0f/a.y, 1.0f/a.z, 1.0f/a.w);
651 #endif
652 }
653
654 __device_inline float4 operator/(const float4& a, float f)
655 {
656         return a * (1.0f/f);
657 }
658
659 __device_inline float4 operator/(const float4& a, const float4& b)
660 {
661 #ifdef __KERNEL_SSE__
662         return a * rcp(b);
663 #else
664         return make_float4(a.x/b.x, a.y/b.y, a.z/b.z, a.w/b.w);
665 #endif
666
667 }
668
669 __device_inline float4 operator+(const float4& a, const float4& b)
670 {
671 #ifdef __KERNEL_SSE__
672         return _mm_add_ps(a.m128, b.m128);
673 #else
674         return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
675 #endif
676 }
677
678 __device_inline float4 operator-(const float4& a, const float4& b)
679 {
680 #ifdef __KERNEL_SSE__
681         return _mm_sub_ps(a.m128, b.m128);
682 #else
683         return make_float4(a.x-b.x, a.y-b.y, a.z-b.z, a.w-b.w);
684 #endif
685 }
686
687 __device_inline float4 operator+=(float4& a, const float4& b)
688 {
689         return a = a + b;
690 }
691
692 __device_inline float4 operator*=(float4& a, const float4& b)
693 {
694         return a = a * b;
695 }
696
697 __device_inline float4 operator/=(float4& a, float f)
698 {
699         return a = a / f;
700 }
701
702 __device_inline int4 operator<(const float4& a, const float4& b)
703 {
704 #ifdef __KERNEL_SSE__
705         return _mm_cvtps_epi32(_mm_cmplt_ps(a.m128, b.m128)); /* todo: avoid cvt */
706 #else
707         return make_int4(a.x < b.x, a.y < b.y, a.z < b.z, a.w < b.w);
708 #endif
709 }
710
711 __device_inline int4 operator>=(float4 a, float4 b)
712 {
713 #ifdef __KERNEL_SSE__
714         return _mm_cvtps_epi32(_mm_cmpge_ps(a.m128, b.m128)); /* todo: avoid cvt */
715 #else
716         return make_int4(a.x >= b.x, a.y >= b.y, a.z >= b.z, a.w >= b.w);
717 #endif
718 }
719
720 __device_inline int4 operator<=(const float4& a, const float4& b)
721 {
722 #ifdef __KERNEL_SSE__
723         return _mm_cvtps_epi32(_mm_cmple_ps(a.m128, b.m128)); /* todo: avoid cvt */
724 #else
725         return make_int4(a.x <= b.x, a.y <= b.y, a.z <= b.z, a.w <= b.w);
726 #endif
727 }
728
729 __device_inline bool operator==(const float4 a, const float4 b)
730 {
731 #ifdef __KERNEL_SSE__
732         return (_mm_movemask_ps(_mm_cmpeq_ps(a.m128, b.m128)) & 15) == 15;
733 #else
734         return (a.x == b.x && a.y == b.y && a.z == b.z && a.w == b.w);
735 #endif
736 }
737
738 __device_inline float4 cross(const float4& a, const float4& b)
739 {
740 #ifdef __KERNEL_SSE__
741         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));
742 #else
743         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);
744 #endif
745 }
746
747 __device_inline bool is_zero(const float4& a)
748 {
749 #ifdef __KERNEL_SSE__
750         return a == make_float4(0.0f);
751 #else
752         return (a.x == 0.0f && a.y == 0.0f && a.z == 0.0f && a.w == 0.0f);
753 #endif
754 }
755
756 __device_inline float reduce_add(const float4& a)
757 {
758 #ifdef __KERNEL_SSE__
759         float4 h = shuffle<1,0,3,2>(a) + a;
760         return _mm_cvtss_f32(shuffle<2,3,0,1>(h) + h); /* todo: efficiency? */
761 #else
762         return ((a.x + a.y) + (a.z + a.w));
763 #endif
764 }
765
766 __device_inline float average(const float4& a)
767 {
768         return reduce_add(a) * 0.25f;
769 }
770
771 __device_inline float dot(const float4& a, const float4& b)
772 {
773         return reduce_add(a * b);
774 }
775
776 __device_inline float len(const float4 a)
777 {
778         return sqrtf(dot(a, a));
779 }
780
781 __device_inline float4 normalize(const float4 a)
782 {
783         return a/len(a);
784 }
785
786 __device_inline float4 min(float4 a, float4 b)
787 {
788 #ifdef __KERNEL_SSE__
789         return _mm_min_ps(a.m128, b.m128);
790 #else
791         return make_float4(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z), min(a.w, b.w));
792 #endif
793 }
794
795 __device_inline float4 max(float4 a, float4 b)
796 {
797 #ifdef __KERNEL_SSE__
798         return _mm_max_ps(a.m128, b.m128);
799 #else
800         return make_float4(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z), max(a.w, b.w));
801 #endif
802 }
803
804 #endif
805
806 #ifndef __KERNEL_GPU__
807
808 __device_inline float4 select(const int4& mask, const float4& a, const float4& b)
809 {
810 #ifdef __KERNEL_SSE__
811         /* blendv is sse4, and apparently broken on vs2008 */
812         return _mm_or_ps(_mm_and_ps(_mm_cvtepi32_ps(mask), a), _mm_andnot_ps(_mm_cvtepi32_ps(mask), b)); /* todo: avoid cvt */
813 #else
814         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);
815 #endif
816 }
817
818 __device_inline float4 reduce_min(const float4& a)
819 {
820 #ifdef __KERNEL_SSE__
821         float4 h = min(shuffle<1,0,3,2>(a), a);
822         return min(shuffle<2,3,0,1>(h), h);
823 #else
824         return make_float4(min(min(a.x, a.y), min(a.z, a.w)));
825 #endif
826 }
827
828 __device_inline float4 reduce_max(const float4& a)
829 {
830 #ifdef __KERNEL_SSE__
831         float4 h = max(shuffle<1,0,3,2>(a), a);
832         return max(shuffle<2,3,0,1>(h), h);
833 #else
834         return make_float4(max(max(a.x, a.y), max(a.z, a.w)));
835 #endif
836 }
837
838 #if 0
839 __device_inline float4 reduce_add(const float4& a)
840 {
841 #ifdef __KERNEL_SSE__
842         float4 h = shuffle<1,0,3,2>(a) + a;
843         return shuffle<2,3,0,1>(h) + h;
844 #else
845         return make_float4((a.x + a.y) + (a.z + a.w));
846 #endif
847 }
848 #endif
849
850 __device_inline void print_float4(const char *label, const float4& a)
851 {
852         printf("%s: %.8f %.8f %.8f %.8f\n", label, (double)a.x, (double)a.y, (double)a.z, (double)a.w);
853 }
854
855 #endif
856
857 /* Int3 */
858
859 #ifndef __KERNEL_OPENCL__
860
861 __device_inline int3 min(int3 a, int3 b)
862 {
863 #ifdef __KERNEL_SSE__
864         return _mm_min_epi32(a.m128, b.m128);
865 #else
866         return make_int3(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z));
867 #endif
868 }
869
870 __device_inline int3 max(int3 a, int3 b)
871 {
872 #ifdef __KERNEL_SSE__
873         return _mm_max_epi32(a.m128, b.m128);
874 #else
875         return make_int3(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z));
876 #endif
877 }
878
879 __device_inline int3 clamp(const int3& a, int mn, int mx)
880 {
881 #ifdef __KERNEL_SSE__
882         return min(max(a, make_int3(mn)), make_int3(mx));
883 #else
884         return make_int3(clamp(a.x, mn, mx), clamp(a.y, mn, mx), clamp(a.z, mn, mx));
885 #endif
886 }
887
888 __device_inline int3 clamp(const int3& a, int3& mn, int mx)
889 {
890 #ifdef __KERNEL_SSE__
891         return min(max(a, mn), make_int3(mx));
892 #else
893         return make_int3(clamp(a.x, mn.x, mx), clamp(a.y, mn.y, mx), clamp(a.z, mn.z, mx));
894 #endif
895 }
896
897 #endif
898
899 #ifndef __KERNEL_GPU__
900
901 __device_inline void print_int3(const char *label, const int3& a)
902 {
903         printf("%s: %d %d %d\n", label, a.x, a.y, a.z);
904 }
905
906 #endif
907
908 /* Int4 */
909
910 #ifndef __KERNEL_GPU__
911
912 __device_inline int4 operator+(const int4& a, const int4& b)
913 {
914 #ifdef __KERNEL_SSE__
915         return _mm_add_epi32(a.m128, b.m128);
916 #else
917         return make_int4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
918 #endif
919 }
920
921 __device_inline int4 operator+=(int4& a, const int4& b)
922 {
923         return a = a + b;
924 }
925
926 __device_inline int4 operator>>(const int4& a, int i)
927 {
928 #ifdef __KERNEL_SSE__
929         return _mm_srai_epi32(a.m128, i);
930 #else
931         return make_int4(a.x >> i, a.y >> i, a.z >> i, a.w >> i);
932 #endif
933 }
934
935 __device_inline int4 min(int4 a, int4 b)
936 {
937 #ifdef __KERNEL_SSE__
938         return _mm_min_epi32(a.m128, b.m128);
939 #else
940         return make_int4(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z), min(a.w, b.w));
941 #endif
942 }
943
944 __device_inline int4 max(int4 a, int4 b)
945 {
946 #ifdef __KERNEL_SSE__
947         return _mm_max_epi32(a.m128, b.m128);
948 #else
949         return make_int4(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z), max(a.w, b.w));
950 #endif
951 }
952
953 __device_inline int4 clamp(const int4& a, const int4& mn, const int4& mx)
954 {
955         return min(max(a, mn), mx);
956 }
957
958 __device_inline int4 select(const int4& mask, const int4& a, const int4& b)
959 {
960 #ifdef __KERNEL_SSE__
961         __m128 m = _mm_cvtepi32_ps(mask);
962         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 */
963 #else
964         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);
965 #endif
966 }
967
968 __device_inline void print_int4(const char *label, const int4& a)
969 {
970         printf("%s: %d %d %d %d\n", label, a.x, a.y, a.z, a.w);
971 }
972
973 #endif
974
975 /* Int/Float conversion */
976
977 #ifndef __KERNEL_OPENCL__
978
979 __device_inline unsigned int as_int(uint i)
980 {
981         union { unsigned int ui; int i; } u;
982         u.ui = i;
983         return u.i;
984 }
985
986 __device_inline unsigned int as_uint(int i)
987 {
988         union { unsigned int ui; int i; } u;
989         u.i = i;
990         return u.ui;
991 }
992
993 __device_inline unsigned int as_uint(float f)
994 {
995         union { unsigned int i; float f; } u;
996         u.f = f;
997         return u.i;
998 }
999
1000 __device_inline int __float_as_int(float f)
1001 {
1002         union { int i; float f; } u;
1003         u.f = f;
1004         return u.i;
1005 }
1006
1007 __device_inline float __int_as_float(int i)
1008 {
1009         union { int i; float f; } u;
1010         u.i = i;
1011         return u.f;
1012 }
1013
1014 __device_inline uint __float_as_uint(float f)
1015 {
1016         union { uint i; float f; } u;
1017         u.f = f;
1018         return u.i;
1019 }
1020
1021 __device_inline float __uint_as_float(uint i)
1022 {
1023         union { uint i; float f; } u;
1024         u.i = i;
1025         return u.f;
1026 }
1027
1028 /* Interpolation */
1029
1030 template<class A, class B> A lerp(const A& a, const A& b, const B& t)
1031 {
1032         return (A)(a * ((B)1 - t) + b * t);
1033 }
1034
1035 /* Triangle */
1036
1037 __device_inline float triangle_area(const float3 v1, const float3 v2, const float3 v3)
1038 {
1039         return len(cross(v3 - v2, v1 - v2))*0.5f;
1040 }
1041
1042 #endif
1043
1044 /* Orthonormal vectors */
1045
1046 __device_inline void make_orthonormals(const float3 N, float3 *a, float3 *b)
1047 {
1048         if(N.x != N.y || N.x != N.z)
1049                 *a = make_float3(N.z-N.y, N.x-N.z, N.y-N.x);  //(1,1,1)x N
1050         else
1051                 *a = make_float3(N.z-N.y, N.x+N.z, -N.y-N.x);  //(-1,1,1)x N
1052
1053         *a = normalize(*a);
1054         *b = cross(N, *a);
1055 }
1056
1057 /* Color division */
1058
1059 __device_inline float3 safe_divide_color(float3 a, float3 b)
1060 {
1061         float x, y, z;
1062
1063         x = (b.x != 0.0f)? a.x/b.x: 0.0f;
1064         y = (b.y != 0.0f)? a.y/b.y: 0.0f;
1065         z = (b.z != 0.0f)? a.z/b.z: 0.0f;
1066
1067         return make_float3(x, y, z);
1068 }
1069
1070 CCL_NAMESPACE_END
1071
1072 #endif /* __UTIL_MATH_H__ */
1073