arraydouble.d

Revision 441, 43.1 kB (checked in by sean, 14 years ago)
Moved boilerplate copyright blurb out of DDoc block.
Property svn:eol-style set to `native`

Line
1	/**
2	* Contains SSE2 and MMX versions of certain operations for double.
3	*
4	* Copyright: Copyright Digital Mars 2008 - 2010.
5	* License: <a href="http://www.boost.org/LICENSE_1_0.txt">Boost License 1.0</a>.
6	* Authors: Walter Bright, based on code originally written by Burton Radons
7	*/
8
9	/* Copyright Digital Mars 2008 - 2010.
10	* Distributed under the Boost Software License, Version 1.0.
11	* (See accompanying file LICENSE_1_0.txt or copy at
12	* http://www.boost.org/LICENSE_1_0.txt)
13	*/
14	module rt.arraydouble;
15
16	private import core.cpuid;
17
18	version (unittest)
19	{
20	private import core.stdc.stdio : printf;
21	/* This is so unit tests will test every CPU variant
22	*/
23	int cpuid;
24	const int CPUID_MAX = 5;
25	bool mmx() { return cpuid == 1 && core.cpuid.mmx(); }
26	bool sse() { return cpuid == 2 && core.cpuid.sse(); }
27	bool sse2() { return cpuid == 3 && core.cpuid.sse2(); }
28	bool amd3dnow() { return cpuid == 4 && core.cpuid.amd3dnow(); }
29	}
30	else
31	{
32	alias core.cpuid.mmx mmx;
33	alias core.cpuid.sse sse;
34	alias core.cpuid.sse2 sse2;
35	alias core.cpuid.amd3dnow amd3dnow;
36	}
37
38	//version = log;
39
40	bool disjoint(T)(T[] a, T[] b)
41	{
42	return (a.ptr + a.length <= b.ptr \|\| b.ptr + b.length <= a.ptr);
43	}
44
45	/* Performance figures measured by Burton Radons
46	*/
47
48	alias double T;
49
50	extern (C):
51
52	/* ======================================================================== */
53
54	/***********************
55	* Computes:
56	* a[] = b[] + c[]
57	*/
58
59	T[] _arraySliceSliceAddSliceAssign_d(T[] a, T[] c, T[] b)
60	in
61	{
62	assert(a.length == b.length && b.length == c.length);
63	assert(disjoint(a, b));
64	assert(disjoint(a, c));
65	assert(disjoint(b, c));
66	}
67	body
68	{
69	auto aptr = a.ptr;
70	auto aend = aptr + a.length;
71	auto bptr = b.ptr;
72	auto cptr = c.ptr;
73
74	version (D_InlineAsm_X86)
75	{
76	// SSE2 version is 333% faster
77	if (sse2() && b.length >= 16)
78	{
79	auto n = aptr + (b.length & ~15);
80
81	// Unaligned case
82	asm
83	{
84	mov EAX, bptr; // left operand
85	mov ECX, cptr; // right operand
86	mov ESI, aptr; // destination operand
87	mov EDI, n; // end comparison
88
89	align 8;
90	startsseloopb:
91	movupd XMM0, [EAX];
92	movupd XMM1, [EAX+16];
93	movupd XMM2, [EAX+32];
94	movupd XMM3, [EAX+48];
95	add EAX, 64;
96	movupd XMM4, [ECX];
97	movupd XMM5, [ECX+16];
98	movupd XMM6, [ECX+32];
99	movupd XMM7, [ECX+48];
100	add ESI, 64;
101	addpd XMM0, XMM4;
102	addpd XMM1, XMM5;
103	addpd XMM2, XMM6;
104	addpd XMM3, XMM7;
105	add ECX, 64;
106	movupd [ESI+ 0-64], XMM0;
107	movupd [ESI+16-64], XMM1;
108	movupd [ESI+32-64], XMM2;
109	movupd [ESI+48-64], XMM3;
110	cmp ESI, EDI;
111	jb startsseloopb;
112
113	mov aptr, ESI;
114	mov bptr, EAX;
115	mov cptr, ECX;
116	}
117	}
118	}
119
120	// Handle remainder
121	while (aptr < aend)
122	aptr++ = bptr++ + *cptr++;
123
124	return a;
125	}
126
127
128	unittest
129	{
130	printf("_arraySliceSliceAddSliceAssign_d unittest\n");
131	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
132	{
133	version (log) printf(" cpuid %d\n", cpuid);
134
135	for (int j = 0; j < 2; j++)
136	{
137	const int dim = 67;
138	T[] a = new T[dim + j]; // aligned on 16 byte boundary
139	a = a[j .. dim + j]; // misalign for second iteration
140	T[] b = new T[dim + j];
141	b = b[j .. dim + j];
142	T[] c = new T[dim + j];
143	c = c[j .. dim + j];
144
145	for (int i = 0; i < dim; i++)
146	{ a[i] = cast(T)i;
147	b[i] = cast(T)(i + 7);
148	c[i] = cast(T)(i * 2);
149	}
150
151	c[] = a[] + b[];
152
153	for (int i = 0; i < dim; i++)
154	{
155	if (c[i] != cast(T)(a[i] + b[i]))
156	{
157	printf("[%d]: %g != %g + %g\n", i, c[i], a[i], b[i]);
158	assert(0);
159	}
160	}
161	}
162	}
163	}
164
165	/* ======================================================================== */
166
167	/***********************
168	* Computes:
169	* a[] = b[] - c[]
170	*/
171
172	T[] _arraySliceSliceMinSliceAssign_d(T[] a, T[] c, T[] b)
173	in
174	{
175	assert(a.length == b.length && b.length == c.length);
176	assert(disjoint(a, b));
177	assert(disjoint(a, c));
178	assert(disjoint(b, c));
179	}
180	body
181	{
182	auto aptr = a.ptr;
183	auto aend = aptr + a.length;
184	auto bptr = b.ptr;
185	auto cptr = c.ptr;
186
187	version (D_InlineAsm_X86)
188	{
189	// SSE2 version is 324% faster
190	if (sse2() && b.length >= 8)
191	{
192	auto n = aptr + (b.length & ~7);
193
194	// Unaligned case
195	asm
196	{
197	mov EAX, bptr; // left operand
198	mov ECX, cptr; // right operand
199	mov ESI, aptr; // destination operand
200	mov EDI, n; // end comparison
201
202	align 8;
203	startsseloopb:
204	movupd XMM0, [EAX];
205	movupd XMM1, [EAX+16];
206	movupd XMM2, [EAX+32];
207	movupd XMM3, [EAX+48];
208	add EAX, 64;
209	movupd XMM4, [ECX];
210	movupd XMM5, [ECX+16];
211	movupd XMM6, [ECX+32];
212	movupd XMM7, [ECX+48];
213	add ESI, 64;
214	subpd XMM0, XMM4;
215	subpd XMM1, XMM5;
216	subpd XMM2, XMM6;
217	subpd XMM3, XMM7;
218	add ECX, 64;
219	movupd [ESI+ 0-64], XMM0;
220	movupd [ESI+16-64], XMM1;
221	movupd [ESI+32-64], XMM2;
222	movupd [ESI+48-64], XMM3;
223	cmp ESI, EDI;
224	jb startsseloopb;
225
226	mov aptr, ESI;
227	mov bptr, EAX;
228	mov cptr, ECX;
229	}
230	}
231	}
232
233	// Handle remainder
234	while (aptr < aend)
235	aptr++ = bptr++ - *cptr++;
236
237	return a;
238	}
239
240
241	unittest
242	{
243	printf("_arraySliceSliceMinSliceAssign_d unittest\n");
244	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
245	{
246	version (log) printf(" cpuid %d\n", cpuid);
247
248	for (int j = 0; j < 2; j++)
249	{
250	const int dim = 67;
251	T[] a = new T[dim + j]; // aligned on 16 byte boundary
252	a = a[j .. dim + j]; // misalign for second iteration
253	T[] b = new T[dim + j];
254	b = b[j .. dim + j];
255	T[] c = new T[dim + j];
256	c = c[j .. dim + j];
257
258	for (int i = 0; i < dim; i++)
259	{ a[i] = cast(T)i;
260	b[i] = cast(T)(i + 7);
261	c[i] = cast(T)(i * 2);
262	}
263
264	c[] = a[] - b[];
265
266	for (int i = 0; i < dim; i++)
267	{
268	if (c[i] != cast(T)(a[i] - b[i]))
269	{
270	printf("[%d]: %g != %g - %g\n", i, c[i], a[i], b[i]);
271	assert(0);
272	}
273	}
274	}
275	}
276	}
277
278
279	/* ======================================================================== */
280
281	/***********************
282	* Computes:
283	* a[] = b[] + value
284	*/
285
286	T[] _arraySliceExpAddSliceAssign_d(T[] a, T value, T[] b)
287	in
288	{
289	assert(a.length == b.length);
290	assert(disjoint(a, b));
291	}
292	body
293	{
294	//printf("_arraySliceExpAddSliceAssign_d()\n");
295	auto aptr = a.ptr;
296	auto aend = aptr + a.length;
297	auto bptr = b.ptr;
298
299	version (D_InlineAsm_X86)
300	{
301	// SSE2 version is 305% faster
302	if (sse2() && a.length >= 8)
303	{
304	auto n = aptr + (a.length & ~7);
305
306	// Unaligned case
307	asm
308	{
309	mov EAX, bptr;
310	mov ESI, aptr;
311	mov EDI, n;
312	movsd XMM4, value;
313	shufpd XMM4, XMM4, 0;
314
315	align 8;
316	startsseloop:
317	add ESI, 64;
318	movupd XMM0, [EAX];
319	movupd XMM1, [EAX+16];
320	movupd XMM2, [EAX+32];
321	movupd XMM3, [EAX+48];
322	add EAX, 64;
323	addpd XMM0, XMM4;
324	addpd XMM1, XMM4;
325	addpd XMM2, XMM4;
326	addpd XMM3, XMM4;
327	movupd [ESI+ 0-64], XMM0;
328	movupd [ESI+16-64], XMM1;
329	movupd [ESI+32-64], XMM2;
330	movupd [ESI+48-64], XMM3;
331	cmp ESI, EDI;
332	jb startsseloop;
333
334	mov aptr, ESI;
335	mov bptr, EAX;
336	}
337	}
338	}
339
340	while (aptr < aend)
341	aptr++ = bptr++ + value;
342
343	return a;
344	}
345
346	unittest
347	{
348	printf("_arraySliceExpAddSliceAssign_d unittest\n");
349	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
350	{
351	version (log) printf(" cpuid %d\n", cpuid);
352
353	for (int j = 0; j < 2; j++)
354	{
355	const int dim = 67;
356	T[] a = new T[dim + j]; // aligned on 16 byte boundary
357	a = a[j .. dim + j]; // misalign for second iteration
358	T[] b = new T[dim + j];
359	b = b[j .. dim + j];
360	T[] c = new T[dim + j];
361	c = c[j .. dim + j];
362
363	for (int i = 0; i < dim; i++)
364	{ a[i] = cast(T)i;
365	b[i] = cast(T)(i + 7);
366	c[i] = cast(T)(i * 2);
367	}
368
369	c[] = a[] + 6;
370
371	for (int i = 0; i < dim; i++)
372	{
373	if (c[i] != cast(T)(a[i] + 6))
374	{
375	printf("[%d]: %g != %g + 6\n", i, c[i], a[i]);
376	assert(0);
377	}
378	}
379	}
380	}
381	}
382
383	/* ======================================================================== */
384
385	/***********************
386	* Computes:
387	* a[] += value
388	*/
389
390	T[] _arrayExpSliceAddass_d(T[] a, T value)
391	{
392	//printf("_arrayExpSliceAddass_d(a.length = %d, value = %Lg)\n", a.length, cast(real)value);
393	auto aptr = a.ptr;
394	auto aend = aptr + a.length;
395
396	version (D_InlineAsm_X86)
397	{
398	// SSE2 version is 114% faster
399	if (sse2() && a.length >= 8)
400	{
401	auto n = aptr + (a.length & ~7);
402	if (aptr < n)
403
404	// Unaligned case
405	asm
406	{
407	mov ESI, aptr;
408	mov EDI, n;
409	movsd XMM4, value;
410	shufpd XMM4, XMM4, 0;
411
412	align 8;
413	startsseloopa:
414	movupd XMM0, [ESI];
415	movupd XMM1, [ESI+16];
416	movupd XMM2, [ESI+32];
417	movupd XMM3, [ESI+48];
418	add ESI, 64;
419	addpd XMM0, XMM4;
420	addpd XMM1, XMM4;
421	addpd XMM2, XMM4;
422	addpd XMM3, XMM4;
423	movupd [ESI+ 0-64], XMM0;
424	movupd [ESI+16-64], XMM1;
425	movupd [ESI+32-64], XMM2;
426	movupd [ESI+48-64], XMM3;
427	cmp ESI, EDI;
428	jb startsseloopa;
429
430	mov aptr, ESI;
431	}
432	}
433	}
434
435	while (aptr < aend)
436	*aptr++ += value;
437
438	return a;
439	}
440
441	unittest
442	{
443	printf("_arrayExpSliceAddass_d unittest\n");
444	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
445	{
446	version (log) printf(" cpuid %d\n", cpuid);
447
448	for (int j = 0; j < 2; j++)
449	{
450	const int dim = 67;
451	T[] a = new T[dim + j]; // aligned on 16 byte boundary
452	a = a[j .. dim + j]; // misalign for second iteration
453	T[] b = new T[dim + j];
454	b = b[j .. dim + j];
455	T[] c = new T[dim + j];
456	c = c[j .. dim + j];
457
458	for (int i = 0; i < dim; i++)
459	{ a[i] = cast(T)i;
460	b[i] = cast(T)(i + 7);
461	c[i] = cast(T)(i * 2);
462	}
463
464	a[] = c[];
465	c[] += 6;
466
467	for (int i = 0; i < dim; i++)
468	{
469	if (c[i] != cast(T)(a[i] + 6))
470	{
471	printf("[%d]: %g != %g + 6\n", i, c[i], a[i]);
472	assert(0);
473	}
474	}
475	}
476	}
477	}
478
479	/* ======================================================================== */
480
481	/***********************
482	* Computes:
483	* a[] += b[]
484	*/
485
486	T[] _arraySliceSliceAddass_d(T[] a, T[] b)
487	in
488	{
489	assert (a.length == b.length);
490	assert (disjoint(a, b));
491	}
492	body
493	{
494	//printf("_arraySliceSliceAddass_d()\n");
495	auto aptr = a.ptr;
496	auto aend = aptr + a.length;
497	auto bptr = b.ptr;
498
499	version (D_InlineAsm_X86)
500	{
501	// SSE2 version is 183% faster
502	if (sse2() && a.length >= 8)
503	{
504	auto n = aptr + (a.length & ~7);
505
506	// Unaligned case
507	asm
508	{
509	mov ECX, bptr; // right operand
510	mov ESI, aptr; // destination operand
511	mov EDI, n; // end comparison
512
513	align 8;
514	startsseloopb:
515	movupd XMM0, [ESI];
516	movupd XMM1, [ESI+16];
517	movupd XMM2, [ESI+32];
518	movupd XMM3, [ESI+48];
519	add ESI, 64;
520	movupd XMM4, [ECX];
521	movupd XMM5, [ECX+16];
522	movupd XMM6, [ECX+32];
523	movupd XMM7, [ECX+48];
524	add ECX, 64;
525	addpd XMM0, XMM4;
526	addpd XMM1, XMM5;
527	addpd XMM2, XMM6;
528	addpd XMM3, XMM7;
529	movupd [ESI+ 0-64], XMM0;
530	movupd [ESI+16-64], XMM1;
531	movupd [ESI+32-64], XMM2;
532	movupd [ESI+48-64], XMM3;
533	cmp ESI, EDI;
534	jb startsseloopb;
535
536	mov aptr, ESI;
537	mov bptr, ECX;
538	}
539	}
540	}
541
542	while (aptr < aend)
543	aptr++ += bptr++;
544
545	return a;
546	}
547
548	unittest
549	{
550	printf("_arraySliceSliceAddass_d unittest\n");
551	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
552	{
553	version (log) printf(" cpuid %d\n", cpuid);
554
555	for (int j = 0; j < 2; j++)
556	{
557	const int dim = 67;
558	T[] a = new T[dim + j]; // aligned on 16 byte boundary
559	a = a[j .. dim + j]; // misalign for second iteration
560	T[] b = new T[dim + j];
561	b = b[j .. dim + j];
562	T[] c = new T[dim + j];
563	c = c[j .. dim + j];
564
565	for (int i = 0; i < dim; i++)
566	{ a[i] = cast(T)i;
567	b[i] = cast(T)(i + 7);
568	c[i] = cast(T)(i * 2);
569	}
570
571	a[] = c[];
572	c[] += b[];
573
574	for (int i = 0; i < dim; i++)
575	{
576	if (c[i] != cast(T)(a[i] + b[i]))
577	{
578	printf("[%d]: %g != %g + %g\n", i, c[i], a[i], b[i]);
579	assert(0);
580	}
581	}
582	}
583	}
584	}
585
586	/* ======================================================================== */
587
588	/***********************
589	* Computes:
590	* a[] = b[] - value
591	*/
592
593	T[] _arraySliceExpMinSliceAssign_d(T[] a, T value, T[] b)
594	in
595	{
596	assert (a.length == b.length);
597	assert (disjoint(a, b));
598	}
599	body
600	{
601	//printf("_arraySliceExpMinSliceAssign_d()\n");
602	auto aptr = a.ptr;
603	auto aend = aptr + a.length;
604	auto bptr = b.ptr;
605
606	version (D_InlineAsm_X86)
607	{
608	// SSE2 version is 305% faster
609	if (sse2() && a.length >= 8)
610	{
611	auto n = aptr + (a.length & ~7);
612
613	// Unaligned case
614	asm
615	{
616	mov EAX, bptr;
617	mov ESI, aptr;
618	mov EDI, n;
619	movsd XMM4, value;
620	shufpd XMM4, XMM4, 0;
621
622	align 8;
623	startsseloop:
624	add ESI, 64;
625	movupd XMM0, [EAX];
626	movupd XMM1, [EAX+16];
627	movupd XMM2, [EAX+32];
628	movupd XMM3, [EAX+48];
629	add EAX, 64;
630	subpd XMM0, XMM4;
631	subpd XMM1, XMM4;
632	subpd XMM2, XMM4;
633	subpd XMM3, XMM4;
634	movupd [ESI+ 0-64], XMM0;
635	movupd [ESI+16-64], XMM1;
636	movupd [ESI+32-64], XMM2;
637	movupd [ESI+48-64], XMM3;
638	cmp ESI, EDI;
639	jb startsseloop;
640
641	mov aptr, ESI;
642	mov bptr, EAX;
643	}
644	}
645	}
646
647	while (aptr < aend)
648	aptr++ = bptr++ - value;
649
650	return a;
651	}
652
653	unittest
654	{
655	printf("_arraySliceExpMinSliceAssign_d unittest\n");
656	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
657	{
658	version (log) printf(" cpuid %d\n", cpuid);
659
660	for (int j = 0; j < 2; j++)
661	{
662	const int dim = 67;
663	T[] a = new T[dim + j]; // aligned on 16 byte boundary
664	a = a[j .. dim + j]; // misalign for second iteration
665	T[] b = new T[dim + j];
666	b = b[j .. dim + j];
667	T[] c = new T[dim + j];
668	c = c[j .. dim + j];
669
670	for (int i = 0; i < dim; i++)
671	{ a[i] = cast(T)i;
672	b[i] = cast(T)(i + 7);
673	c[i] = cast(T)(i * 2);
674	}
675
676	c[] = a[] - 6;
677
678	for (int i = 0; i < dim; i++)
679	{
680	if (c[i] != cast(T)(a[i] - 6))
681	{
682	printf("[%d]: %g != %g - 6\n", i, c[i], a[i]);
683	assert(0);
684	}
685	}
686	}
687	}
688	}
689
690	/* ======================================================================== */
691
692	/***********************
693	* Computes:
694	* a[] = value - b[]
695	*/
696
697	T[] _arrayExpSliceMinSliceAssign_d(T[] a, T[] b, T value)
698	in
699	{
700	assert (a.length == b.length);
701	assert (disjoint(a, b));
702	}
703	body
704	{
705	//printf("_arrayExpSliceMinSliceAssign_d()\n");
706	auto aptr = a.ptr;
707	auto aend = aptr + a.length;
708	auto bptr = b.ptr;
709
710	version (D_InlineAsm_X86)
711	{
712	// SSE2 version is 66% faster
713	if (sse2() && a.length >= 8)
714	{
715	auto n = aptr + (a.length & ~7);
716
717	// Unaligned case
718	asm
719	{
720	mov EAX, bptr;
721	mov ESI, aptr;
722	mov EDI, n;
723	movsd XMM4, value;
724	shufpd XMM4, XMM4, 0;
725
726	align 8;
727	startsseloop:
728	add ESI, 64;
729	movapd XMM5, XMM4;
730	movapd XMM6, XMM4;
731	movupd XMM0, [EAX];
732	movupd XMM1, [EAX+16];
733	movupd XMM2, [EAX+32];
734	movupd XMM3, [EAX+48];
735	add EAX, 64;
736	subpd XMM5, XMM0;
737	subpd XMM6, XMM1;
738	movupd [ESI+ 0-64], XMM5;
739	movupd [ESI+16-64], XMM6;
740	movapd XMM5, XMM4;
741	movapd XMM6, XMM4;
742	subpd XMM5, XMM2;
743	subpd XMM6, XMM3;
744	movupd [ESI+32-64], XMM5;
745	movupd [ESI+48-64], XMM6;
746	cmp ESI, EDI;
747	jb startsseloop;
748
749	mov aptr, ESI;
750	mov bptr, EAX;
751	}
752	}
753	}
754
755	while (aptr < aend)
756	aptr++ = value - bptr++;
757
758	return a;
759	}
760
761	unittest
762	{
763	printf("_arrayExpSliceMinSliceAssign_d unittest\n");
764	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
765	{
766	version (log) printf(" cpuid %d\n", cpuid);
767
768	for (int j = 0; j < 2; j++)
769	{
770	const int dim = 67;
771	T[] a = new T[dim + j]; // aligned on 16 byte boundary
772	a = a[j .. dim + j]; // misalign for second iteration
773	T[] b = new T[dim + j];
774	b = b[j .. dim + j];
775	T[] c = new T[dim + j];
776	c = c[j .. dim + j];
777
778	for (int i = 0; i < dim; i++)
779	{ a[i] = cast(T)i;
780	b[i] = cast(T)(i + 7);
781	c[i] = cast(T)(i * 2);
782	}
783
784	c[] = 6 - a[];
785
786	for (int i = 0; i < dim; i++)
787	{
788	if (c[i] != cast(T)(6 - a[i]))
789	{
790	printf("[%d]: %g != 6 - %g\n", i, c[i], a[i]);
791	assert(0);
792	}
793	}
794	}
795	}
796	}
797
798	/* ======================================================================== */
799
800	/***********************
801	* Computes:
802	* a[] -= value
803	*/
804
805	T[] _arrayExpSliceMinass_d(T[] a, T value)
806	{
807	//printf("_arrayExpSliceMinass_d(a.length = %d, value = %Lg)\n", a.length, cast(real)value);
808	auto aptr = a.ptr;
809	auto aend = aptr + a.length;
810
811	version (D_InlineAsm_X86)
812	{
813	// SSE2 version is 115% faster
814	if (sse2() && a.length >= 8)
815	{
816	auto n = aptr + (a.length & ~7);
817	if (aptr < n)
818
819	// Unaligned case
820	asm
821	{
822	mov ESI, aptr;
823	mov EDI, n;
824	movsd XMM4, value;
825	shufpd XMM4, XMM4, 0;
826
827	align 8;
828	startsseloopa:
829	movupd XMM0, [ESI];
830	movupd XMM1, [ESI+16];
831	movupd XMM2, [ESI+32];
832	movupd XMM3, [ESI+48];
833	add ESI, 64;
834	subpd XMM0, XMM4;
835	subpd XMM1, XMM4;
836	subpd XMM2, XMM4;
837	subpd XMM3, XMM4;
838	movupd [ESI+ 0-64], XMM0;
839	movupd [ESI+16-64], XMM1;
840	movupd [ESI+32-64], XMM2;
841	movupd [ESI+48-64], XMM3;
842	cmp ESI, EDI;
843	jb startsseloopa;
844
845	mov aptr, ESI;
846	}
847	}
848	}
849
850	while (aptr < aend)
851	*aptr++ -= value;
852
853	return a;
854	}
855
856	unittest
857	{
858	printf("_arrayExpSliceMinass_d unittest\n");
859	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
860	{
861	version (log) printf(" cpuid %d\n", cpuid);
862
863	for (int j = 0; j < 2; j++)
864	{
865	const int dim = 67;
866	T[] a = new T[dim + j]; // aligned on 16 byte boundary
867	a = a[j .. dim + j]; // misalign for second iteration
868	T[] b = new T[dim + j];
869	b = b[j .. dim + j];
870	T[] c = new T[dim + j];
871	c = c[j .. dim + j];
872
873	for (int i = 0; i < dim; i++)
874	{ a[i] = cast(T)i;
875	b[i] = cast(T)(i + 7);
876	c[i] = cast(T)(i * 2);
877	}
878
879	a[] = c[];
880	c[] -= 6;
881
882	for (int i = 0; i < dim; i++)
883	{
884	if (c[i] != cast(T)(a[i] - 6))
885	{
886	printf("[%d]: %g != %g - 6\n", i, c[i], a[i]);
887	assert(0);
888	}
889	}
890	}
891	}
892	}
893
894	/* ======================================================================== */
895
896	/***********************
897	* Computes:
898	* a[] -= b[]
899	*/
900
901	T[] _arraySliceSliceMinass_d(T[] a, T[] b)
902	in
903	{
904	assert (a.length == b.length);
905	assert (disjoint(a, b));
906	}
907	body
908	{
909	//printf("_arraySliceSliceMinass_d()\n");
910	auto aptr = a.ptr;
911	auto aend = aptr + a.length;
912	auto bptr = b.ptr;
913
914	version (D_InlineAsm_X86)
915	{
916	// SSE2 version is 183% faster
917	if (sse2() && a.length >= 8)
918	{
919	auto n = aptr + (a.length & ~7);
920
921	// Unaligned case
922	asm
923	{
924	mov ECX, bptr; // right operand
925	mov ESI, aptr; // destination operand
926	mov EDI, n; // end comparison
927
928	align 8;
929	startsseloopb:
930	movupd XMM0, [ESI];
931	movupd XMM1, [ESI+16];
932	movupd XMM2, [ESI+32];
933	movupd XMM3, [ESI+48];
934	add ESI, 64;
935	movupd XMM4, [ECX];
936	movupd XMM5, [ECX+16];
937	movupd XMM6, [ECX+32];
938	movupd XMM7, [ECX+48];
939	add ECX, 64;
940	subpd XMM0, XMM4;
941	subpd XMM1, XMM5;
942	subpd XMM2, XMM6;
943	subpd XMM3, XMM7;
944	movupd [ESI+ 0-64], XMM0;
945	movupd [ESI+16-64], XMM1;
946	movupd [ESI+32-64], XMM2;
947	movupd [ESI+48-64], XMM3;
948	cmp ESI, EDI;
949	jb startsseloopb;
950
951	mov aptr, ESI;
952	mov bptr, ECX;
953	}
954	}
955	}
956
957	while (aptr < aend)
958	aptr++ -= bptr++;
959
960	return a;
961	}
962
963	unittest
964	{
965	printf("_arrayExpSliceMinass_d unittest\n");
966	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
967	{
968	version (log) printf(" cpuid %d\n", cpuid);
969
970	for (int j = 0; j < 2; j++)
971	{
972	const int dim = 67;
973	T[] a = new T[dim + j]; // aligned on 16 byte boundary
974	a = a[j .. dim + j]; // misalign for second iteration
975	T[] b = new T[dim + j];
976	b = b[j .. dim + j];
977	T[] c = new T[dim + j];
978	c = c[j .. dim + j];
979
980	for (int i = 0; i < dim; i++)
981	{ a[i] = cast(T)i;
982	b[i] = cast(T)(i + 7);
983	c[i] = cast(T)(i * 2);
984	}
985
986	a[] = c[];
987	c[] -= 6;
988
989	for (int i = 0; i < dim; i++)
990	{
991	if (c[i] != cast(T)(a[i] - 6))
992	{
993	printf("[%d]: %g != %g - 6\n", i, c[i], a[i]);
994	assert(0);
995	}
996	}
997	}
998	}
999	}
1000
1001	/* ======================================================================== */
1002
1003	/***********************
1004	* Computes:
1005	* a[] = b[] * value
1006	*/
1007
1008	T[] _arraySliceExpMulSliceAssign_d(T[] a, T value, T[] b)
1009	in
1010	{
1011	assert(a.length == b.length);
1012	assert(disjoint(a, b));
1013	}
1014	body
1015	{
1016	//printf("_arraySliceExpMulSliceAssign_d()\n");
1017	auto aptr = a.ptr;
1018	auto aend = aptr + a.length;
1019	auto bptr = b.ptr;
1020
1021	version (D_InlineAsm_X86)
1022	{
1023	// SSE2 version is 304% faster
1024	if (sse2() && a.length >= 8)
1025	{
1026	auto n = aptr + (a.length & ~7);
1027
1028	// Unaligned case
1029	asm
1030	{
1031	mov EAX, bptr;
1032	mov ESI, aptr;
1033	mov EDI, n;
1034	movsd XMM4, value;
1035	shufpd XMM4, XMM4, 0;
1036
1037	align 8;
1038	startsseloop:
1039	add ESI, 64;
1040	movupd XMM0, [EAX];
1041	movupd XMM1, [EAX+16];
1042	movupd XMM2, [EAX+32];
1043	movupd XMM3, [EAX+48];
1044	add EAX, 64;
1045	mulpd XMM0, XMM4;
1046	mulpd XMM1, XMM4;
1047	mulpd XMM2, XMM4;
1048	mulpd XMM3, XMM4;
1049	movupd [ESI+ 0-64], XMM0;
1050	movupd [ESI+16-64], XMM1;
1051	movupd [ESI+32-64], XMM2;
1052	movupd [ESI+48-64], XMM3;
1053	cmp ESI, EDI;
1054	jb startsseloop;
1055
1056	mov aptr, ESI;
1057	mov bptr, EAX;
1058	}
1059	}
1060	}
1061
1062	while (aptr < aend)
1063	aptr++ = bptr++ * value;
1064
1065	return a;
1066	}
1067
1068	unittest
1069	{
1070	printf("_arraySliceExpMulSliceAssign_d unittest\n");
1071	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1072	{
1073	version (log) printf(" cpuid %d\n", cpuid);
1074
1075	for (int j = 0; j < 2; j++)
1076	{
1077	const int dim = 67;
1078	T[] a = new T[dim + j]; // aligned on 16 byte boundary
1079	a = a[j .. dim + j]; // misalign for second iteration
1080	T[] b = new T[dim + j];
1081	b = b[j .. dim + j];
1082	T[] c = new T[dim + j];
1083	c = c[j .. dim + j];
1084
1085	for (int i = 0; i < dim; i++)
1086	{ a[i] = cast(T)i;
1087	b[i] = cast(T)(i + 7);
1088	c[i] = cast(T)(i * 2);
1089	}
1090
1091	c[] = a[] * 6;
1092
1093	for (int i = 0; i < dim; i++)
1094	{
1095	if (c[i] != cast(T)(a[i] * 6))
1096	{
1097	printf("[%d]: %g != %g * 6\n", i, c[i], a[i]);
1098	assert(0);
1099	}
1100	}
1101	}
1102	}
1103	}
1104
1105	/* ======================================================================== */
1106
1107	/***********************
1108	* Computes:
1109	* a[] = b[] * c[]
1110	*/
1111
1112	T[] _arraySliceSliceMulSliceAssign_d(T[] a, T[] c, T[] b)
1113	in
1114	{
1115	assert(a.length == b.length && b.length == c.length);
1116	assert(disjoint(a, b));
1117	assert(disjoint(a, c));
1118	assert(disjoint(b, c));
1119	}
1120	body
1121	{
1122	//printf("_arraySliceSliceMulSliceAssign_d()\n");
1123	auto aptr = a.ptr;
1124	auto aend = aptr + a.length;
1125	auto bptr = b.ptr;
1126	auto cptr = c.ptr;
1127
1128	version (D_InlineAsm_X86)
1129	{
1130	// SSE2 version is 329% faster
1131	if (sse2() && a.length >= 8)
1132	{
1133	auto n = aptr + (a.length & ~7);
1134
1135	// Unaligned case
1136	asm
1137	{
1138	mov EAX, bptr; // left operand
1139	mov ECX, cptr; // right operand
1140	mov ESI, aptr; // destination operand
1141	mov EDI, n; // end comparison
1142
1143	align 8;
1144	startsseloopb:
1145	movupd XMM0, [EAX];
1146	movupd XMM1, [EAX+16];
1147	movupd XMM2, [EAX+32];
1148	movupd XMM3, [EAX+48];
1149	add ESI, 64;
1150	movupd XMM4, [ECX];
1151	movupd XMM5, [ECX+16];
1152	movupd XMM6, [ECX+32];
1153	movupd XMM7, [ECX+48];
1154	add EAX, 64;
1155	mulpd XMM0, XMM4;
1156	mulpd XMM1, XMM5;
1157	mulpd XMM2, XMM6;
1158	mulpd XMM3, XMM7;
1159	add ECX, 64;
1160	movupd [ESI+ 0-64], XMM0;
1161	movupd [ESI+16-64], XMM1;
1162	movupd [ESI+32-64], XMM2;
1163	movupd [ESI+48-64], XMM3;
1164	cmp ESI, EDI;
1165	jb startsseloopb;
1166
1167	mov aptr, ESI;
1168	mov bptr, EAX;
1169	mov cptr, ECX;
1170	}
1171	}
1172	}
1173
1174	while (aptr < aend)
1175	aptr++ = bptr++ * *cptr++;
1176
1177	return a;
1178	}
1179
1180	unittest
1181	{
1182	printf("_arraySliceSliceMulSliceAssign_d unittest\n");
1183	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1184	{
1185	version (log) printf(" cpuid %d\n", cpuid);
1186
1187	for (int j = 0; j < 2; j++)
1188	{
1189	const int dim = 67;
1190	T[] a = new T[dim + j]; // aligned on 16 byte boundary
1191	a = a[j .. dim + j]; // misalign for second iteration
1192	T[] b = new T[dim + j];
1193	b = b[j .. dim + j];
1194	T[] c = new T[dim + j];
1195	c = c[j .. dim + j];
1196
1197	for (int i = 0; i < dim; i++)
1198	{ a[i] = cast(T)i;
1199	b[i] = cast(T)(i + 7);
1200	c[i] = cast(T)(i * 2);
1201	}
1202
1203	c[] = a[] * b[];
1204
1205	for (int i = 0; i < dim; i++)
1206	{
1207	if (c[i] != cast(T)(a[i] * b[i]))
1208	{
1209	printf("[%d]: %g != %g * %g\n", i, c[i], a[i], b[i]);
1210	assert(0);
1211	}
1212	}
1213	}
1214	}
1215	}
1216
1217	/* ======================================================================== */
1218
1219	/***********************
1220	* Computes:
1221	* a[] *= value
1222	*/
1223
1224	T[] _arrayExpSliceMulass_d(T[] a, T value)
1225	{
1226	//printf("_arrayExpSliceMulass_d(a.length = %d, value = %Lg)\n", a.length, cast(real)value);
1227	auto aptr = a.ptr;
1228	auto aend = aptr + a.length;
1229
1230	version (D_InlineAsm_X86)
1231	{
1232	// SSE2 version is 109% faster
1233	if (sse2() && a.length >= 8)
1234	{
1235	auto n = aptr + (a.length & ~7);
1236	if (aptr < n)
1237
1238	// Unaligned case
1239	asm
1240	{
1241	mov ESI, aptr;
1242	mov EDI, n;
1243	movsd XMM4, value;
1244	shufpd XMM4, XMM4, 0;
1245
1246	align 8;
1247	startsseloopa:
1248	movupd XMM0, [ESI];
1249	movupd XMM1, [ESI+16];
1250	movupd XMM2, [ESI+32];
1251	movupd XMM3, [ESI+48];
1252	add ESI, 64;
1253	mulpd XMM0, XMM4;
1254	mulpd XMM1, XMM4;
1255	mulpd XMM2, XMM4;
1256	mulpd XMM3, XMM4;
1257	movupd [ESI+ 0-64], XMM0;
1258	movupd [ESI+16-64], XMM1;
1259	movupd [ESI+32-64], XMM2;
1260	movupd [ESI+48-64], XMM3;
1261	cmp ESI, EDI;
1262	jb startsseloopa;
1263
1264	mov aptr, ESI;
1265	}
1266	}
1267	}
1268
1269	while (aptr < aend)
1270	aptr++ = value;
1271
1272	return a;
1273	}
1274
1275	unittest
1276	{
1277	printf("_arrayExpSliceMulass_d unittest\n");
1278	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1279	{
1280	version (log) printf(" cpuid %d\n", cpuid);
1281
1282	for (int j = 0; j < 2; j++)
1283	{
1284	const int dim = 67;
1285	T[] a = new T[dim + j]; // aligned on 16 byte boundary
1286	a = a[j .. dim + j]; // misalign for second iteration
1287	T[] b = new T[dim + j];
1288	b = b[j .. dim + j];
1289	T[] c = new T[dim + j];
1290	c = c[j .. dim + j];
1291
1292	for (int i = 0; i < dim; i++)
1293	{ a[i] = cast(T)i;
1294	b[i] = cast(T)(i + 7);
1295	c[i] = cast(T)(i * 2);
1296	}
1297
1298	a[] = c[];
1299	c[] *= 6;
1300
1301	for (int i = 0; i < dim; i++)
1302	{
1303	if (c[i] != cast(T)(a[i] * 6))
1304	{
1305	printf("[%d]: %g != %g * 6\n", i, c[i], a[i]);
1306	assert(0);
1307	}
1308	}
1309	}
1310	}
1311	}
1312
1313	/* ======================================================================== */
1314
1315	/***********************
1316	* Computes:
1317	* a[] *= b[]
1318	*/
1319
1320	T[] _arraySliceSliceMulass_d(T[] a, T[] b)
1321	in
1322	{
1323	assert (a.length == b.length);
1324	assert (disjoint(a, b));
1325	}
1326	body
1327	{
1328	//printf("_arraySliceSliceMulass_d()\n");
1329	auto aptr = a.ptr;
1330	auto aend = aptr + a.length;
1331	auto bptr = b.ptr;
1332
1333	version (D_InlineAsm_X86)
1334	{
1335	// SSE2 version is 205% faster
1336	if (sse2() && a.length >= 8)
1337	{
1338	auto n = aptr + (a.length & ~7);
1339
1340	// Unaligned case
1341	asm
1342	{
1343	mov ECX, bptr; // right operand
1344	mov ESI, aptr; // destination operand
1345	mov EDI, n; // end comparison
1346
1347	align 8;
1348	startsseloopb:
1349	movupd XMM0, [ESI];
1350	movupd XMM1, [ESI+16];
1351	movupd XMM2, [ESI+32];
1352	movupd XMM3, [ESI+48];
1353	add ESI, 64;
1354	movupd XMM4, [ECX];
1355	movupd XMM5, [ECX+16];
1356	movupd XMM6, [ECX+32];
1357	movupd XMM7, [ECX+48];
1358	add ECX, 64;
1359	mulpd XMM0, XMM4;
1360	mulpd XMM1, XMM5;
1361	mulpd XMM2, XMM6;
1362	mulpd XMM3, XMM7;
1363	movupd [ESI+ 0-64], XMM0;
1364	movupd [ESI+16-64], XMM1;
1365	movupd [ESI+32-64], XMM2;
1366	movupd [ESI+48-64], XMM3;
1367	cmp ESI, EDI;
1368	jb startsseloopb;
1369
1370	mov aptr, ESI;
1371	mov bptr, ECX;
1372	}
1373	}
1374	}
1375
1376	while (aptr < aend)
1377	aptr++ = *bptr++;
1378
1379	return a;
1380	}
1381
1382	unittest
1383	{
1384	printf("_arrayExpSliceMulass_d unittest\n");
1385	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1386	{
1387	version (log) printf(" cpuid %d\n", cpuid);
1388
1389	for (int j = 0; j < 2; j++)
1390	{
1391	const int dim = 67;
1392	T[] a = new T[dim + j]; // aligned on 16 byte boundary
1393	a = a[j .. dim + j]; // misalign for second iteration
1394	T[] b = new T[dim + j];
1395	b = b[j .. dim + j];
1396	T[] c = new T[dim + j];
1397	c = c[j .. dim + j];
1398
1399	for (int i = 0; i < dim; i++)
1400	{ a[i] = cast(T)i;
1401	b[i] = cast(T)(i + 7);
1402	c[i] = cast(T)(i * 2);
1403	}
1404
1405	a[] = c[];
1406	c[] *= 6;
1407
1408	for (int i = 0; i < dim; i++)
1409	{
1410	if (c[i] != cast(T)(a[i] * 6))
1411	{
1412	printf("[%d]: %g != %g * 6\n", i, c[i], a[i]);
1413	assert(0);
1414	}
1415	}
1416	}
1417	}
1418	}
1419
1420	/* ======================================================================== */
1421
1422	/***********************
1423	* Computes:
1424	* a[] = b[] / value
1425	*/
1426
1427	T[] _arraySliceExpDivSliceAssign_d(T[] a, T value, T[] b)
1428	in
1429	{
1430	assert(a.length == b.length);
1431	assert(disjoint(a, b));
1432	}
1433	body
1434	{
1435	//printf("_arraySliceExpDivSliceAssign_d()\n");
1436	auto aptr = a.ptr;
1437	auto aend = aptr + a.length;
1438	auto bptr = b.ptr;
1439
1440	/* Multiplying by the reciprocal is faster, but does
1441	* not produce as accurate an answer.
1442	*/
1443	T recip = cast(T)1 / value;
1444
1445	version (D_InlineAsm_X86)
1446	{
1447	// SSE2 version is 299% faster
1448	if (sse2() && a.length >= 8)
1449	{
1450	auto n = aptr + (a.length & ~7);
1451
1452	// Unaligned case
1453	asm
1454	{
1455	mov EAX, bptr;
1456	mov ESI, aptr;
1457	mov EDI, n;
1458	movsd XMM4, recip;
1459	//movsd XMM4, value
1460	//rcpsd XMM4, XMM4
1461	shufpd XMM4, XMM4, 0;
1462
1463	align 8;
1464	startsseloop:
1465	add ESI, 64;
1466	movupd XMM0, [EAX];
1467	movupd XMM1, [EAX+16];
1468	movupd XMM2, [EAX+32];
1469	movupd XMM3, [EAX+48];
1470	add EAX, 64;
1471	mulpd XMM0, XMM4;
1472	mulpd XMM1, XMM4;
1473	mulpd XMM2, XMM4;
1474	mulpd XMM3, XMM4;
1475	//divpd XMM0, XMM4;
1476	//divpd XMM1, XMM4;
1477	//divpd XMM2, XMM4;
1478	//divpd XMM3, XMM4;
1479	movupd [ESI+ 0-64], XMM0;
1480	movupd [ESI+16-64], XMM1;
1481	movupd [ESI+32-64], XMM2;
1482	movupd [ESI+48-64], XMM3;
1483	cmp ESI, EDI;
1484	jb startsseloop;
1485
1486	mov aptr, ESI;
1487	mov bptr, EAX;
1488	}
1489	}
1490	}
1491
1492	while (aptr < aend)
1493	{
1494	aptr++ = bptr++ / value;
1495	//aptr++ = bptr++ * recip;
1496	}
1497
1498	return a;
1499	}
1500
1501	unittest
1502	{
1503	printf("_arraySliceExpDivSliceAssign_d unittest\n");
1504	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1505	{
1506	version (log) printf(" cpuid %d\n", cpuid);
1507
1508	for (int j = 0; j < 2; j++)
1509	{
1510	const int dim = 67;
1511	T[] a = new T[dim + j]; // aligned on 16 byte boundary
1512	a = a[j .. dim + j]; // misalign for second iteration
1513	T[] b = new T[dim + j];
1514	b = b[j .. dim + j];
1515	T[] c = new T[dim + j];
1516	c = c[j .. dim + j];
1517
1518	for (int i = 0; i < dim; i++)
1519	{ a[i] = cast(T)i;
1520	b[i] = cast(T)(i + 7);
1521	c[i] = cast(T)(i * 2);
1522	}
1523
1524	c[] = a[] / 8;
1525
1526	for (int i = 0; i < dim; i++)
1527	{
1528	//printf("[%d]: %g ?= %g / 8\n", i, c[i], a[i]);
1529	if (c[i] != cast(T)(a[i] / 8))
1530	{
1531	printf("[%d]: %g != %g / 8\n", i, c[i], a[i]);
1532	assert(0);
1533	}
1534	}
1535	}
1536	}
1537	}
1538
1539	/* ======================================================================== */
1540
1541	/***********************
1542	* Computes:
1543	* a[] /= value
1544	*/
1545
1546	T[] _arrayExpSliceDivass_d(T[] a, T value)
1547	{
1548	//printf("_arrayExpSliceDivass_d(a.length = %d, value = %Lg)\n", a.length, cast(real)value);
1549	auto aptr = a.ptr;
1550	auto aend = aptr + a.length;
1551
1552	/* Multiplying by the reciprocal is faster, but does
1553	* not produce as accurate an answer.
1554	*/
1555	T recip = cast(T)1 / value;
1556
1557	version (D_InlineAsm_X86)
1558	{
1559	// SSE2 version is 65% faster
1560	if (sse2() && a.length >= 8)
1561	{
1562	auto n = aptr + (a.length & ~7);
1563
1564	// Unaligned case
1565	asm
1566	{
1567	mov ESI, aptr;
1568	mov EDI, n;
1569	movsd XMM4, recip;
1570	//movsd XMM4, value
1571	//rcpsd XMM4, XMM4
1572	shufpd XMM4, XMM4, 0;
1573
1574	align 8;
1575	startsseloopa:
1576	movupd XMM0, [ESI];
1577	movupd XMM1, [ESI+16];
1578	movupd XMM2, [ESI+32];
1579	movupd XMM3, [ESI+48];
1580	add ESI, 64;
1581	mulpd XMM0, XMM4;
1582	mulpd XMM1, XMM4;
1583	mulpd XMM2, XMM4;
1584	mulpd XMM3, XMM4;
1585	//divpd XMM0, XMM4;
1586	//divpd XMM1, XMM4;
1587	//divpd XMM2, XMM4;
1588	//divpd XMM3, XMM4;
1589	movupd [ESI+ 0-64], XMM0;
1590	movupd [ESI+16-64], XMM1;
1591	movupd [ESI+32-64], XMM2;
1592	movupd [ESI+48-64], XMM3;
1593	cmp ESI, EDI;
1594	jb startsseloopa;
1595
1596	mov aptr, ESI;
1597	}
1598	}
1599	}
1600
1601	while (aptr < aend)
1602	aptr++ = recip;
1603
1604	return a;
1605	}
1606
1607
1608	unittest
1609	{
1610	printf("_arrayExpSliceDivass_d unittest\n");
1611	for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1612	{
1613	version (log) printf(" cpuid %d\n", cpuid);
1614
1615	for (int j = 0; j < 2; j++)
1616	{
1617	const int dim = 67;
1618	T[] a = new T[dim + j]; // aligned on 16 byte boundary
1619	a = a[j .. dim + j]; // misalign for second iteration
1620	T[] b = new T[dim + j];
1621	b = b[j .. dim + j];
1622	T[] c = new T[dim + j];
1623	c = c[j .. dim + j];
1624
1625	for (int i = 0; i < dim; i++)
1626	{ a[i] = cast(T)i;
1627	b[i] = cast(T)(i + 7);
1628	c[i] = cast(T)(i * 2);
1629	}
1630
1631	a[] = c[];
1632	c[] /= 8;
1633
1634	for (int i = 0; i < dim; i++)
1635	{
1636	if (c[i] != cast(T)(a[i] / 8))
1637	{
1638	printf("[%d]: %g != %g / 8\n", i, c[i], a[i]);
1639	assert(0);
1640	}
1641	}
1642	}
1643	}
1644	}
1645
1646
1647	/* ======================================================================== */
1648
1649	/***********************
1650	* Computes:
1651	* a[] -= b[] * value
1652	*/
1653
1654	T[] _arraySliceExpMulSliceMinass_d(T[] a, T value, T[] b)
1655	{
1656	return _arraySliceExpMulSliceAddass_d(a, -value, b);
1657	}
1658
1659	/***********************
1660	* Computes:
1661	* a[] += b[] * value
1662	*/
1663
1664	T[] _arraySliceExpMulSliceAddass_d(T[] a, T value, T[] b)
1665	in
1666	{
1667	assert(a.length == b.length);
1668	assert(disjoint(a, b));
1669	}
1670	body
1671	{
1672	auto aptr = a.ptr;
1673	auto aend = aptr + a.length;
1674	auto bptr = b.ptr;
1675
1676	// Handle remainder
1677	while (aptr < aend)
1678	aptr++ += bptr++ * value;
1679
1680	return a;
1681	}
1682
1683	unittest
1684	{
1685	printf("_arraySliceExpMulSliceAddass_d unittest\n");
1686
1687	cpuid = 1;
1688	{
1689	version (log) printf(" cpuid %d\n", cpuid);
1690
1691	for (int j = 0; j < 1; j++)
1692	{
1693	const int dim = 67;
1694	T[] a = new T[dim + j]; // aligned on 16 byte boundary
1695	a = a[j .. dim + j]; // misalign for second iteration
1696	T[] b = new T[dim + j];
1697	b = b[j .. dim + j];
1698	T[] c = new T[dim + j];
1699	c = c[j .. dim + j];
1700
1701	for (int i = 0; i < dim; i++)
1702	{ a[i] = cast(T)i;
1703	b[i] = cast(T)(i + 7);
1704	c[i] = cast(T)(i * 2);
1705	}
1706
1707	b[] = c[];
1708	c[] += a[] * 6;
1709
1710	for (int i = 0; i < dim; i++)
1711	{
1712	//printf("[%d]: %g ?= %g + %g * 6\n", i, c[i], b[i], a[i]);
1713	if (c[i] != cast(T)(b[i] + a[i] * 6))
1714	{
1715	printf("[%d]: %g ?= %g + %g * 6\n", i, c[i], b[i], a[i]);
1716	assert(0);
1717	}
1718	}
1719	}
1720	}
1721	}

Note: See TracBrowser for help on using the browser.

druntime

root/trunk/src/rt/arraydouble.d

Download in other formats: