1use itertools::izip;
2use pulp::{Arch, WithSimd};
3
4pub(crate) fn logaddexp(a: f64, b: f64) -> f64 {
5 if a == b {
6 return a + 2f64.ln();
7 }
8 let diff = a - b;
9 if diff > 0. {
10 a + (-diff).exp().ln_1p()
11 } else if diff < 0. {
12 b + diff.exp().ln_1p()
13 } else {
14 diff
16 }
17}
18
19struct Multiply<'a> {
20 x: &'a [f64],
21 y: &'a [f64],
22 out: &'a mut [f64],
23}
24
25impl<'a> WithSimd for Multiply<'a> {
26 type Output = ();
27
28 #[inline(always)]
29 fn with_simd<S: pulp::Simd>(self, simd: S) -> Self::Output {
30 let x = self.x;
31 let y = self.y;
32 let out = self.out;
33
34 let (x_out, x_tail) = S::as_simd_f64s(x);
35 let (y_out, y_tail) = S::as_simd_f64s(y);
36 let (out_head, out_tail) = S::as_mut_simd_f64s(out);
37
38 let (out_arrays, out_simd_tail) = pulp::as_arrays_mut::<4, _>(out_head);
39 let (x_arrays, x_simd_tail) = pulp::as_arrays::<4, _>(x_out);
40 let (y_arrays, y_simd_tail) = pulp::as_arrays::<4, _>(y_out);
41
42 izip!(out_arrays, x_arrays, y_arrays).for_each(
43 |([out0, out1, out2, out3], [x0, x1, x2, x3], [y0, y1, y2, y3])| {
44 *out0 = simd.mul_f64s(*x0, *y0);
45 *out1 = simd.mul_f64s(*x1, *y1);
46 *out2 = simd.mul_f64s(*x2, *y2);
47 *out3 = simd.mul_f64s(*x3, *y3);
48 },
49 );
50
51 izip!(
52 out_simd_tail.iter_mut(),
53 x_simd_tail.iter(),
54 y_simd_tail.iter()
55 )
56 .for_each(|(out, &x, &y)| {
57 *out = simd.mul_f64s(x, y);
58 });
59
60 izip!(out_tail.iter_mut(), x_tail.iter(), y_tail.iter()).for_each(|(out, &x, &y)| {
61 *out = x * y;
62 });
63 }
64}
65
66#[inline(never)]
67pub fn multiply(arch: Arch, x: &[f64], y: &[f64], out: &mut [f64]) {
68 arch.dispatch(Multiply { x, y, out })
69}
70
71struct ScalarProds2<'a> {
72 positive1: &'a [f64],
73 positive2: &'a [f64],
74 x: &'a [f64],
75 y: &'a [f64],
76}
77
78impl<'a> WithSimd for ScalarProds2<'a> {
79 type Output = (f64, f64);
80
81 #[inline(always)]
82 fn with_simd<S: pulp::Simd>(self, simd: S) -> Self::Output {
83 let positive1 = self.positive1;
84 let positive2 = self.positive2;
85 let x = self.x;
86 let y = self.y;
87
88 let (p1_out, p1_tail) = S::as_simd_f64s(positive1);
89 let (p2_out, p2_tail) = S::as_simd_f64s(positive2);
90 let (x_out, x_tail) = S::as_simd_f64s(x);
91 let (y_out, y_tail) = S::as_simd_f64s(y);
92
93 let mut s1_0 = simd.splat_f64s(0.0);
94 let mut s1_1 = simd.splat_f64s(0.0);
95 let mut s1_2 = simd.splat_f64s(0.0);
96 let mut s1_3 = simd.splat_f64s(0.0);
97 let mut s2_0 = simd.splat_f64s(0.0);
98 let mut s2_1 = simd.splat_f64s(0.0);
99 let mut s2_2 = simd.splat_f64s(0.0);
100 let mut s2_3 = simd.splat_f64s(0.0);
101
102 let (p1_out, p1_simd_tail) = pulp::as_arrays::<4, _>(p1_out);
103 let (p2_out, p2_simd_tail) = pulp::as_arrays::<4, _>(p2_out);
104 let (x_out, x_simd_tail) = pulp::as_arrays::<4, _>(x_out);
105 let (y_out, y_simd_tail) = pulp::as_arrays::<4, _>(y_out);
106
107 izip!(p1_out, p2_out, x_out, y_out).for_each(
108 |(
109 [p1_0, p1_1, p1_2, p1_3],
110 [p2_0, p2_1, p2_2, p2_3],
111 [x_0, x_1, x_2, x_3],
112 [y_0, y_1, y_2, y_3],
113 )| {
114 let sum0 = simd.add_f64s(*p1_0, *p2_0);
115 let sum1 = simd.add_f64s(*p1_1, *p2_1);
116 let sum2 = simd.add_f64s(*p1_2, *p2_2);
117 let sum3 = simd.add_f64s(*p1_3, *p2_3);
118 s1_0 = simd.mul_add_e_f64s(sum0, *x_0, s1_0);
119 s1_1 = simd.mul_add_e_f64s(sum1, *x_1, s1_1);
120 s1_2 = simd.mul_add_e_f64s(sum2, *x_2, s1_2);
121 s1_3 = simd.mul_add_e_f64s(sum3, *x_3, s1_3);
122 s2_0 = simd.mul_add_e_f64s(sum0, *y_0, s2_0);
123 s2_1 = simd.mul_add_e_f64s(sum1, *y_1, s2_1);
124 s2_2 = simd.mul_add_e_f64s(sum2, *y_2, s2_2);
125 s2_3 = simd.mul_add_e_f64s(sum3, *y_3, s2_3);
126 },
127 );
128
129 izip!(p1_simd_tail, p2_simd_tail, x_simd_tail, y_simd_tail).for_each(|(p1, p2, x, y)| {
130 let sum = simd.add_f64s(*p1, *p2);
131 s1_0 = simd.mul_add_e_f64s(sum, *x, s1_0);
132 s2_0 = simd.mul_add_e_f64s(sum, *y, s2_0);
133 });
134
135 let mut out = (
136 simd.reduce_sum_f64s(
137 simd.add_f64s(simd.add_f64s(s1_0, s1_1), simd.add_f64s(s1_2, s1_3)),
138 ),
139 simd.reduce_sum_f64s(
140 simd.add_f64s(simd.add_f64s(s2_0, s2_1), simd.add_f64s(s2_2, s2_3)),
141 ),
142 );
143
144 izip!(p1_tail.iter(), p2_tail.iter(), x_tail.iter(), y_tail.iter()).for_each(
145 |(p1, p2, x, y)| {
146 let sum = *p1 + *p2;
147 out.0 += sum * *x;
148 out.1 += sum * *y;
149 },
150 );
151 out
152 }
153}
154
155#[inline(never)]
156pub fn scalar_prods2(
157 arch: Arch,
158 positive1: &[f64],
159 positive2: &[f64],
160 x: &[f64],
161 y: &[f64],
162) -> (f64, f64) {
163 let n = positive1.len();
164
165 assert!(positive1.len() == n);
166 assert!(positive2.len() == n);
167 assert!(x.len() == n);
168 assert!(y.len() == n);
169
170 arch.dispatch(ScalarProds2 {
171 positive1,
172 positive2,
173 x,
174 y,
175 })
176}
177
178struct ScalarProds3<'a> {
179 positive1: &'a [f64],
180 negative1: &'a [f64],
181 positive2: &'a [f64],
182 x: &'a [f64],
183 y: &'a [f64],
184}
185
186impl<'a> WithSimd for ScalarProds3<'a> {
187 type Output = (f64, f64);
188
189 #[inline(always)]
190 fn with_simd<S: pulp::Simd>(self, simd: S) -> Self::Output {
191 let positive1 = self.positive1;
192 let negative1 = self.negative1;
193 let positive2 = self.positive2;
194 let x = self.x;
195 let y = self.y;
196
197 let (p1_out, p1_tail) = S::as_simd_f64s(positive1);
198 let (n1_out, n1_tail) = S::as_simd_f64s(negative1);
199 let (p2_out, p2_tail) = S::as_simd_f64s(positive2);
200 let (x_out, x_tail) = S::as_simd_f64s(x);
201 let (y_out, y_tail) = S::as_simd_f64s(y);
202
203 let mut s1_0 = simd.splat_f64s(0.0);
204 let mut s1_1 = simd.splat_f64s(0.0);
205 let mut s1_2 = simd.splat_f64s(0.0);
206 let mut s1_3 = simd.splat_f64s(0.0);
207 let mut s2_0 = simd.splat_f64s(0.0);
208 let mut s2_1 = simd.splat_f64s(0.0);
209 let mut s2_2 = simd.splat_f64s(0.0);
210 let mut s2_3 = simd.splat_f64s(0.0);
211
212 let (p1_out, p1_simd_tail) = pulp::as_arrays::<4, _>(p1_out);
213 let (n1_out, n1_simd_tail) = pulp::as_arrays::<4, _>(n1_out);
214 let (p2_out, p2_simd_tail) = pulp::as_arrays::<4, _>(p2_out);
215 let (x_out, x_simd_tail) = pulp::as_arrays::<4, _>(x_out);
216 let (y_out, y_simd_tail) = pulp::as_arrays::<4, _>(y_out);
217
218 izip!(p1_out, n1_out, p2_out, x_out, y_out).for_each(
219 |(
220 [p1_0, p1_1, p1_2, p1_3],
221 [n1_0, n1_1, n1_2, n1_3],
222 [p2_0, p2_1, p2_2, p2_3],
223 [x_0, x_1, x_2, x_3],
224 [y_0, y_1, y_2, y_3],
225 )| {
226 let sum0 = simd.sub_f64s(simd.add_f64s(*p1_0, *p2_0), *n1_0);
227 let sum1 = simd.sub_f64s(simd.add_f64s(*p1_1, *p2_1), *n1_1);
228 let sum2 = simd.sub_f64s(simd.add_f64s(*p1_2, *p2_2), *n1_2);
229 let sum3 = simd.sub_f64s(simd.add_f64s(*p1_3, *p2_3), *n1_3);
230 s1_0 = simd.mul_add_e_f64s(sum0, *x_0, s1_0);
231 s1_1 = simd.mul_add_e_f64s(sum1, *x_1, s1_1);
232 s1_2 = simd.mul_add_e_f64s(sum2, *x_2, s1_2);
233 s1_3 = simd.mul_add_e_f64s(sum3, *x_3, s1_3);
234 s2_0 = simd.mul_add_e_f64s(sum0, *y_0, s2_0);
235 s2_1 = simd.mul_add_e_f64s(sum1, *y_1, s2_1);
236 s2_2 = simd.mul_add_e_f64s(sum2, *y_2, s2_2);
237 s2_3 = simd.mul_add_e_f64s(sum3, *y_3, s2_3);
238 },
239 );
240
241 izip!(
242 p1_simd_tail,
243 n1_simd_tail,
244 p2_simd_tail,
245 x_simd_tail,
246 y_simd_tail
247 )
248 .for_each(|(p1, n1, p2, x, y)| {
249 let sum = simd.sub_f64s(simd.add_f64s(*p1, *p2), *n1);
250 s1_0 = simd.mul_add_e_f64s(sum, *x, s1_0);
251 s2_0 = simd.mul_add_e_f64s(sum, *y, s2_0);
252 });
253
254 let mut out = (
255 simd.reduce_sum_f64s(
256 simd.add_f64s(simd.add_f64s(s1_0, s1_1), simd.add_f64s(s1_2, s1_3)),
257 ),
258 simd.reduce_sum_f64s(
259 simd.add_f64s(simd.add_f64s(s2_0, s2_1), simd.add_f64s(s2_2, s2_3)),
260 ),
261 );
262
263 izip!(
264 p1_tail.iter(),
265 n1_tail.iter(),
266 p2_tail.iter(),
267 x_tail.iter(),
268 y_tail.iter()
269 )
270 .for_each(|(p1, n1, p2, x, y)| {
271 let sum = *p1 - *n1 + *p2;
272 out.0 += sum * *x;
273 out.1 += sum * *y;
274 });
275
276 out
277 }
278}
279
280#[inline(never)]
281pub fn scalar_prods3(
282 arch: Arch,
283 positive1: &[f64],
284 negative1: &[f64],
285 positive2: &[f64],
286 x: &[f64],
287 y: &[f64],
288) -> (f64, f64) {
289 let n = positive1.len();
290
291 assert!(positive1.len() == n);
292 assert!(positive2.len() == n);
293 assert!(negative1.len() == n);
294 assert!(x.len() == n);
295 assert!(y.len() == n);
296
297 arch.dispatch(ScalarProds3 {
298 positive1,
299 negative1,
300 positive2,
301 x,
302 y,
303 })
304}
305
306struct VectorDot<'a> {
307 x: &'a [f64],
308 y: &'a [f64],
309}
310
311impl<'a> WithSimd for VectorDot<'a> {
312 type Output = f64;
313
314 #[inline(always)]
315 fn with_simd<S: pulp::Simd>(self, simd: S) -> Self::Output {
316 let a = self.x;
317 let b = self.y;
318
319 assert!(a.len() == b.len());
320
321 let (x, x_tail) = S::as_simd_f64s(a);
322 let (y, y_tail) = S::as_simd_f64s(b);
323
324 let mut out0 = simd.splat_f64s(0f64);
325 let mut out1 = simd.splat_f64s(0f64);
326 let mut out2 = simd.splat_f64s(0f64);
327 let mut out3 = simd.splat_f64s(0f64);
328
329 let (x, x_simd_tail) = pulp::as_arrays::<4, _>(x);
330 let (y, y_simd_tail) = pulp::as_arrays::<4, _>(y);
331
332 izip!(x, y).for_each(|([x0, x1, x2, x3], [y0, y1, y2, y3])| {
333 out0 = simd.mul_add_e_f64s(*x0, *y0, out0);
334 out1 = simd.mul_add_e_f64s(*x1, *y1, out1);
335 out2 = simd.mul_add_e_f64s(*x2, *y2, out2);
336 out3 = simd.mul_add_e_f64s(*x3, *y3, out3);
337 });
338
339 izip!(x_simd_tail, y_simd_tail).for_each(|(&x, &y)| {
340 out0 = simd.mul_add_e_f64s(x, y, out0);
341 });
342
343 out0 = simd.add_f64s(out0, out1);
344 out1 = simd.add_f64s(out2, out3);
345 out0 = simd.add_f64s(out0, out1);
346 let mut result = simd.reduce_sum_f64s(out0);
347
348 x_tail.iter().zip(y_tail).for_each(|(&x, &y)| {
349 result += x * y;
350 });
351 result
352 }
353}
354
355pub fn vector_dot(arch: Arch, a: &[f64], b: &[f64]) -> f64 {
356 arch.dispatch(VectorDot { x: a, y: b })
357}
358
359struct Axpy<'a> {
360 x: &'a [f64],
361 y: &'a mut [f64],
362 a: f64,
363}
364
365impl<'a> WithSimd for Axpy<'a> {
366 type Output = ();
367
368 #[inline(always)]
369 fn with_simd<S: pulp::Simd>(self, simd: S) -> Self::Output {
370 let x = self.x;
371 let y = self.y;
372 let a = self.a;
373
374 let (x_out, x_tail) = S::as_simd_f64s(x);
375 let (y_out, y_tail) = S::as_mut_simd_f64s(y);
376
377 let a_splat = simd.splat_f64s(a);
378
379 let (x_arrays, x_simd_tail) = pulp::as_arrays::<4, _>(x_out);
380 let (y_arrays, y_simd_tail) = pulp::as_arrays_mut::<4, _>(y_out);
381
382 izip!(x_arrays, y_arrays).for_each(|([x0, x1, x2, x3], [y0, y1, y2, y3])| {
383 *y0 = simd.mul_add_e_f64s(a_splat, *x0, *y0);
384 *y1 = simd.mul_add_e_f64s(a_splat, *x1, *y1);
385 *y2 = simd.mul_add_e_f64s(a_splat, *x2, *y2);
386 *y3 = simd.mul_add_e_f64s(a_splat, *x3, *y3);
387 });
388
389 izip!(x_simd_tail.iter(), y_simd_tail.iter_mut()).for_each(|(&x, y)| {
390 *y = simd.mul_add_e_f64s(a_splat, x, *y);
391 });
392
393 izip!(x_tail.iter(), y_tail.iter_mut()).for_each(|(&x, y)| {
394 *y = a.mul_add(x, *y);
395 });
396 }
397}
398pub fn axpy(arch: Arch, x: &[f64], y: &mut [f64], a: f64) {
399 let n = x.len();
400 assert!(y.len() == n);
401
402 arch.dispatch(Axpy { x, y, a });
403}
404
405struct AxpyOut<'a> {
406 x: &'a [f64],
407 y: &'a [f64],
408 out: &'a mut [f64],
409 a: f64,
410}
411
412impl<'a> WithSimd for AxpyOut<'a> {
413 type Output = ();
414
415 #[inline(always)]
416 fn with_simd<S: pulp::Simd>(self, simd: S) -> Self::Output {
417 let x = self.x;
418 let y = self.y;
419 let out = self.out;
420 let a = self.a;
421
422 let (x_out, x_tail) = S::as_simd_f64s(x);
423 let (y_out, y_tail) = S::as_simd_f64s(y);
424 let (out_head, out_tail) = S::as_mut_simd_f64s(out);
425
426 let a_splat = simd.splat_f64s(a);
427
428 let (out_arrays, out_simd_tail) = pulp::as_arrays_mut::<4, _>(out_head);
429 let (x_arrays, x_simd_tail) = pulp::as_arrays::<4, _>(x_out);
430 let (y_arrays, y_simd_tail) = pulp::as_arrays::<4, _>(y_out);
431
432 izip!(out_arrays, x_arrays, y_arrays).for_each(
433 |([out0, out1, out2, out3], [x0, x1, x2, x3], [y0, y1, y2, y3])| {
434 *out0 = simd.mul_add_e_f64s(a_splat, *x0, *y0);
435 *out1 = simd.mul_add_e_f64s(a_splat, *x1, *y1);
436 *out2 = simd.mul_add_e_f64s(a_splat, *x2, *y2);
437 *out3 = simd.mul_add_e_f64s(a_splat, *x3, *y3);
438 },
439 );
440
441 izip!(
442 out_simd_tail.iter_mut(),
443 x_simd_tail.iter(),
444 y_simd_tail.iter()
445 )
446 .for_each(|(out, &x, &y)| {
447 *out = simd.mul_add_e_f64s(a_splat, x, y);
448 });
449
450 izip!(x_tail.iter(), y_tail.iter(), out_tail.iter_mut()).for_each(|(&x, &y, out)| {
451 *out = a.mul_add(x, y);
452 });
453 }
454}
455
456pub fn axpy_out(arch: Arch, x: &[f64], y: &[f64], a: f64, out: &mut [f64]) {
457 let n = x.len();
458 assert!(y.len() == n);
459 assert!(out.len() == n);
460
461 arch.dispatch(AxpyOut { x, y, out, a });
462}
463
464#[cfg(test)]
465mod tests {
466 use super::*;
467 use approx::assert_ulps_eq;
468 use pretty_assertions::assert_eq;
469 use proptest::prelude::*;
470
471 fn assert_approx_eq(a: f64, b: f64) {
472 if a.is_nan() && b.is_nan() | b.is_infinite() {
473 return;
474 }
475 if b.is_nan() && a.is_nan() | a.is_infinite() {
476 return;
477 }
478 assert_ulps_eq!(a, b, max_ulps = 8);
479 }
480
481 prop_compose! {
482 fn array2(maxsize: usize) (size in 0..maxsize) (
483 vec1 in prop::collection::vec(prop::num::f64::ANY, size),
484 vec2 in prop::collection::vec(prop::num::f64::ANY, size)
485 )
486 -> (Vec<f64>, Vec<f64>) {
487 (vec1, vec2)
488 }
489 }
490
491 prop_compose! {
492 fn array3(maxsize: usize) (size in 0..maxsize) (
493 vec1 in prop::collection::vec(prop::num::f64::ANY, size),
494 vec2 in prop::collection::vec(prop::num::f64::ANY, size),
495 vec3 in prop::collection::vec(prop::num::f64::ANY, size)
496 )
497 -> (Vec<f64>, Vec<f64>, Vec<f64>) {
498 (vec1, vec2, vec3)
499 }
500 }
501
502 prop_compose! {
503 fn array4(maxsize: usize) (size in 0..maxsize) (
504 vec1 in prop::collection::vec(prop::num::f64::ANY, size),
505 vec2 in prop::collection::vec(prop::num::f64::ANY, size),
506 vec3 in prop::collection::vec(prop::num::f64::ANY, size),
507 vec4 in prop::collection::vec(prop::num::f64::ANY, size)
508 )
509 -> (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>) {
510 (vec1, vec2, vec3, vec4)
511 }
512 }
513
514 prop_compose! {
515 fn array5(maxsize: usize) (size in 0..maxsize) (
516 vec1 in prop::collection::vec(prop::num::f64::ANY, size),
517 vec2 in prop::collection::vec(prop::num::f64::ANY, size),
518 vec3 in prop::collection::vec(prop::num::f64::ANY, size),
519 vec4 in prop::collection::vec(prop::num::f64::ANY, size),
520 vec5 in prop::collection::vec(prop::num::f64::ANY, size)
521 )
522 -> (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>) {
523 (vec1, vec2, vec3, vec4, vec5)
524 }
525 }
526
527 proptest! {
528 #[test]
529 fn check_logaddexp(x in -10f64..10f64, y in -10f64..10f64) {
530 let a = (x.exp() + y.exp()).ln();
531 let b = logaddexp(x, y);
532 let neginf = f64::NEG_INFINITY;
533 let nan = f64::NAN;
534 prop_assert!((a - b).abs() < 1e-10);
535 prop_assert_eq!(b, logaddexp(y, x));
536 prop_assert_eq!(x, logaddexp(x, neginf));
537 prop_assert_eq!(logaddexp(neginf, neginf), neginf);
538 prop_assert!(logaddexp(nan, x).is_nan());
539 }
540
541 #[test]
542 fn test_axpy((x, y) in array2(10), a in prop::num::f64::ANY) {
543 let arch = pulp::Arch::default();
544 let orig = y.clone();
545 let mut y = y.clone();
546 axpy(arch, &x[..], &mut y[..], a);
547 for ((&x, y), out) in x.iter().zip(orig).zip(y) {
548 assert_approx_eq(out, a.mul_add(x, y));
549 }
550 }
551
552 #[test]
553 fn test_scalar_prods2((x1, x2, y1, y2) in array4(10)) {
554 let arch = pulp::Arch::default();
555 let (p1, p2) = scalar_prods2(arch, &x1[..], &x2[..], &y1[..], &y2[..]);
556 let x1 = ndarray::Array1::from_vec(x1);
557 let x2 = ndarray::Array1::from_vec(x2);
558 let y1 = ndarray::Array1::from_vec(y1);
559 let y2 = ndarray::Array1::from_vec(y2);
560 assert_approx_eq(p1, (&x1 + &x2).dot(&y1));
561 assert_approx_eq(p2, (&x1 + &x2).dot(&y2));
562 }
563
564 #[test]
565 fn test_scalar_prods3((x1, x2, x3, y1, y2) in array5(10)) {
566 let arch = Arch::default();
567 let (p1, p2) = scalar_prods3(arch, &x1[..], &x2[..], &x3[..], &y1[..], &y2[..]);
568 let x1 = ndarray::Array1::from_vec(x1);
569 let x2 = ndarray::Array1::from_vec(x2);
570 let x3 = ndarray::Array1::from_vec(x3);
571 let y1 = ndarray::Array1::from_vec(y1);
572 let y2 = ndarray::Array1::from_vec(y2);
573 assert_approx_eq(p1, (&x1 - &x2 + &x3).dot(&y1));
574 assert_approx_eq(p2, (&x1 - &x2 + &x3).dot(&y2));
575 }
576
577 #[test]
578 fn test_axpy_out(a in prop::num::f64::ANY, (x, y, out) in array3(10)) {
579 let arch = Arch::default();
580 let mut out = out.clone();
581 axpy_out(arch, &x[..], &y[..], a, &mut out[..]);
582 let x = ndarray::Array1::from_vec(x);
583 let mut y = ndarray::Array1::from_vec(y);
584 y.scaled_add(a, &x);
585 for (&out1, out2) in out.iter().zip(y) {
586 assert_approx_eq(out1, out2);
587 }
588 }
589
590 #[test]
591 fn test_multiplty((x, y, out) in array3(10)) {
592 let arch = pulp::Arch::default();
593 let mut out = out.clone();
594 multiply(arch, &x[..], &y[..], &mut out[..]);
595 let x = ndarray::Array1::from_vec(x);
596 let y = ndarray::Array1::from_vec(y);
597 for (&out1, out2) in out.iter().zip(&x * &y) {
598 assert_approx_eq(out1, out2);
599 }
600 }
601
602 #[test]
603 fn test_vector_dot((x, y) in array2(10)) {
604 let arch = pulp::Arch::default();
605 let actual = vector_dot(arch, &x[..], &y[..]);
606 let x = ndarray::Array1::from_vec(x);
607 let y = ndarray::Array1::from_vec(y);
608 let expected = x.iter().zip(y.iter()).map(|(&x, &y)| x * y).sum();
609 assert_approx_eq(actual, expected);
610 }
611 }
612
613 #[test]
614 fn check_neginf() {
615 assert_eq!(logaddexp(f64::NEG_INFINITY, 2.), 2.);
616 assert_eq!(logaddexp(2., f64::NEG_INFINITY), 2.);
617 }
618}