src/compiler/dmd/arrayfloat.d

   1 /***************************
   2  * D programming language http://www.digitalmars.com/d/
   3  * Runtime support for float array operations.
   4  * Based on code originally written by Burton Radons.
   5  * Placed in public domain.
   6  */
   7
   8 module rt.arrayfloat;
   9
  10 private import util.cpuid;
  11
  12 version (unittest)
  13 {
  14     private import core.stdc.stdio : printf;
  15     /* This is so unit tests will test every CPU variant
  16      */
  17     int cpuid;
  18     const int CPUID_MAX = 5;
  19     bool mmx()      { return cpuid == 1 && util.cpuid.mmx(); }
  20     bool sse()      { return cpuid == 2 && util.cpuid.sse(); }
  21     bool sse2()     { return cpuid == 3 && util.cpuid.sse2(); }
  22     bool amd3dnow() { return cpuid == 4 && util.cpuid.amd3dnow(); }
  23 }
  24 else
  25 {
  26     alias util.cpuid.mmx mmx;
  27     alias util.cpuid.sse sse;
  28     alias util.cpuid.sse2 sse2;
  29     alias util.cpuid.amd3dnow amd3dnow;
  30 }
  31
  32 //version = log;
  33
  34 bool disjoint(T)(T[] a, T[] b)
  35 {
  36     return (a.ptr + a.length <= b.ptr || b.ptr + b.length <= a.ptr);
  37 }
  38
  39 alias float T;
  40
  41 extern (C):
  42
  43 /* ======================================================================== */
  44
  45 /***********************
  46  * Computes:
  47  *      a[] = b[] + c[]
  48  */
  49
  50 T[] _arraySliceSliceAddSliceAssign_f(T[] a, T[] c, T[] b)
  51 in
  52 {
  53         assert(a.length == b.length && b.length == c.length);
  54         assert(disjoint(a, b));
  55         assert(disjoint(a, c));
  56         assert(disjoint(b, c));
  57 }
  58 body
  59 {
  60     //printf("_arraySliceSliceAddSliceAssign_f()\n");
  61     auto aptr = a.ptr;
  62     auto aend = aptr + a.length;
  63     auto bptr = b.ptr;
  64     auto cptr = c.ptr;
  65
  66     version (D_InlineAsm_X86)
  67     {
  68         // SSE version is 834% faster
  69         if (sse() && b.length >= 16)
  70         {
  71             version (log) printf("\tsse unaligned\n");
  72             auto n = aptr + (b.length & ~15);
  73
  74             // Unaligned case
  75             asm
  76             {
  77                 mov EAX, bptr; // left operand
  78                 mov ECX, cptr; // right operand
  79                 mov ESI, aptr; // destination operand
  80                 mov EDI, n;    // end comparison
  81
  82                 align 8;
  83             startsseloopb:
  84                 movups XMM0, [EAX];
  85                 movups XMM1, [EAX+16];
  86                 movups XMM2, [EAX+32];
  87                 movups XMM3, [EAX+48];
  88                 add EAX, 64;
  89                 movups XMM4, [ECX];
  90                 movups XMM5, [ECX+16];
  91                 movups XMM6, [ECX+32];
  92                 movups XMM7, [ECX+48];
  93                 add ESI, 64;
  94                 addps XMM0, XMM4;
  95                 addps XMM1, XMM5;
  96                 addps XMM2, XMM6;
  97                 addps XMM3, XMM7;
  98                 add ECX, 64;
  99                 movups [ESI+ 0-64], XMM0;
 100                 movups [ESI+16-64], XMM1;
 101                 movups [ESI+32-64], XMM2;
 102                 movups [ESI+48-64], XMM3;
 103                 cmp ESI, EDI;
 104                 jb startsseloopb;
 105
 106                 mov aptr, ESI;
 107                 mov bptr, EAX;
 108                 mov cptr, ECX;
 109             }
 110         }
 111         else
 112         // 3DNow! version is only 13% faster
 113         if (amd3dnow() && b.length >= 8)
 114         {
 115             version (log) printf("\tamd3dnow\n");
 116             auto n = aptr + (b.length & ~7);
 117
 118             asm
 119             {
 120                 mov ESI, aptr; // destination operand
 121                 mov EDI, n;    // end comparison
 122                 mov EAX, bptr; // left operand
 123                 mov ECX, cptr; // right operand
 124
 125                 align 4;
 126             start3dnow:
 127                 movq MM0, [EAX];
 128                 movq MM1, [EAX+8];
 129                 movq MM2, [EAX+16];
 130                 movq MM3, [EAX+24];
 131                 pfadd MM0, [ECX];
 132                 pfadd MM1, [ECX+8];
 133                 pfadd MM2, [ECX+16];
 134                 pfadd MM3, [ECX+24];
 135                 movq [ESI], MM0;
 136                 movq [ESI+8], MM1;
 137                 movq [ESI+16], MM2;
 138                 movq [ESI+24], MM3;
 139                 add ECX, 32;
 140                 add ESI, 32;
 141                 add EAX, 32;
 142                 cmp ESI, EDI;
 143                 jb start3dnow;
 144
 145                 emms;
 146                 mov aptr, ESI;
 147                 mov bptr, EAX;
 148                 mov cptr, ECX;
 149             }
 150         }
 151     }
 152
 153     // Handle remainder
 154     version (log) if (aptr < aend) printf("\tbase\n");
 155     while (aptr < aend)
 156         *aptr++ = *bptr++ + *cptr++;
 157
 158     return a;
 159 }
 160
 161
 162 unittest
 163 {
 164     printf("_arraySliceSliceAddSliceAssign_f unittest\n");
 165     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
 166     {
 167         version (log) printf("    cpuid %d\n", cpuid);
 168
 169         for (int j = 0; j < 2; j++)
 170         {
 171             const int dim = 67;
 172             T[] a = new T[dim + j];     // aligned on 16 byte boundary
 173             a = a[j .. dim + j];        // misalign for second iteration
 174             T[] b = new T[dim + j];
 175             b = b[j .. dim + j];
 176             T[] c = new T[dim + j];
 177             c = c[j .. dim + j];
 178
 179             for (int i = 0; i < dim; i++)
 180             {   a[i] = cast(T)i;
 181                 b[i] = cast(T)(i + 7);
 182                 c[i] = cast(T)(i * 2);
 183             }
 184
 185             c[] = a[] + b[];
 186
 187             for (int i = 0; i < dim; i++)
 188             {
 189                 if (c[i] != cast(T)(a[i] + b[i]))
 190                 {
 191                     printf("[%d]: %g != %g + %g\n", i, c[i], a[i], b[i]);
 192                     assert(0);
 193                 }
 194             }
 195         }
 196     }
 197 }
 198
 199 /* ======================================================================== */
 200
 201 /***********************
 202  * Computes:
 203  *      a[] = b[] - c[]
 204  */
 205
 206 T[] _arraySliceSliceMinSliceAssign_f(T[] a, T[] c, T[] b)
 207 in
 208 {
 209         assert(a.length == b.length && b.length == c.length);
 210         assert(disjoint(a, b));
 211         assert(disjoint(a, c));
 212         assert(disjoint(b, c));
 213 }
 214 body
 215 {
 216     auto aptr = a.ptr;
 217     auto aend = aptr + a.length;
 218     auto bptr = b.ptr;
 219     auto cptr = c.ptr;
 220
 221     version (D_InlineAsm_X86)
 222     {
 223         // SSE version is 834% faster
 224         if (sse() && b.length >= 16)
 225         {
 226             auto n = aptr + (b.length & ~15);
 227
 228             // Unaligned case
 229             asm
 230             {
 231                 mov EAX, bptr; // left operand
 232                 mov ECX, cptr; // right operand
 233                 mov ESI, aptr; // destination operand
 234                 mov EDI, n;    // end comparison
 235
 236                 align 8;
 237             startsseloopb:
 238                 movups XMM0, [EAX];
 239                 movups XMM1, [EAX+16];
 240                 movups XMM2, [EAX+32];
 241                 movups XMM3, [EAX+48];
 242                 add EAX, 64;
 243                 movups XMM4, [ECX];
 244                 movups XMM5, [ECX+16];
 245                 movups XMM6, [ECX+32];
 246                 movups XMM7, [ECX+48];
 247                 add ESI, 64;
 248                 subps XMM0, XMM4;
 249                 subps XMM1, XMM5;
 250                 subps XMM2, XMM6;
 251                 subps XMM3, XMM7;
 252                 add ECX, 64;
 253                 movups [ESI+ 0-64], XMM0;
 254                 movups [ESI+16-64], XMM1;
 255                 movups [ESI+32-64], XMM2;
 256                 movups [ESI+48-64], XMM3;
 257                 cmp ESI, EDI;
 258                 jb startsseloopb;
 259
 260                 mov aptr, ESI;
 261                 mov bptr, EAX;
 262                 mov cptr, ECX;
 263             }
 264         }
 265         else
 266         // 3DNow! version is only 13% faster
 267         if (amd3dnow() && b.length >= 8)
 268         {
 269             auto n = aptr + (b.length & ~7);
 270
 271             asm
 272             {
 273                 mov ESI, aptr; // destination operand
 274                 mov EDI, n;    // end comparison
 275                 mov EAX, bptr; // left operand
 276                 mov ECX, cptr; // right operand
 277
 278                 align 4;
 279             start3dnow:
 280                 movq MM0, [EAX];
 281                 movq MM1, [EAX+8];
 282                 movq MM2, [EAX+16];
 283                 movq MM3, [EAX+24];
 284                 pfsub MM0, [ECX];
 285                 pfsub MM1, [ECX+8];
 286                 pfsub MM2, [ECX+16];
 287                 pfsub MM3, [ECX+24];
 288                 movq [ESI], MM0;
 289                 movq [ESI+8], MM1;
 290                 movq [ESI+16], MM2;
 291                 movq [ESI+24], MM3;
 292                 add ECX, 32;
 293                 add ESI, 32;
 294                 add EAX, 32;
 295                 cmp ESI, EDI;
 296                 jb start3dnow;
 297
 298                 emms;
 299                 mov aptr, ESI;
 300                 mov bptr, EAX;
 301                 mov cptr, ECX;
 302             }
 303         }
 304     }
 305
 306     // Handle remainder
 307     while (aptr < aend)
 308         *aptr++ = *bptr++ - *cptr++;
 309
 310     return a;
 311 }
 312
 313
 314 unittest
 315 {
 316     printf("_arraySliceSliceMinSliceAssign_f unittest\n");
 317     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
 318     {
 319         version (log) printf("    cpuid %d\n", cpuid);
 320
 321         for (int j = 0; j < 2; j++)
 322         {
 323             const int dim = 67;
 324             T[] a = new T[dim + j];     // aligned on 16 byte boundary
 325             a = a[j .. dim + j];        // misalign for second iteration
 326             T[] b = new T[dim + j];
 327             b = b[j .. dim + j];
 328             T[] c = new T[dim + j];
 329             c = c[j .. dim + j];
 330
 331             for (int i = 0; i < dim; i++)
 332             {   a[i] = cast(T)i;
 333                 b[i] = cast(T)(i + 7);
 334                 c[i] = cast(T)(i * 2);
 335             }
 336
 337             c[] = a[] - b[];
 338
 339             for (int i = 0; i < dim; i++)
 340             {
 341                 if (c[i] != cast(T)(a[i] - b[i]))
 342                 {
 343                     printf("[%d]: %g != %gd - %g\n", i, c[i], a[i], b[i]);
 344                     assert(0);
 345                 }
 346             }
 347         }
 348     }
 349 }
 350
 351 /* ======================================================================== */
 352
 353 /***********************
 354  * Computes:
 355  *      a[] = b[] + value
 356  */
 357
 358 T[] _arraySliceExpAddSliceAssign_f(T[] a, T value, T[] b)
 359 in
 360 {
 361     assert(a.length == b.length);
 362     assert(disjoint(a, b));
 363 }
 364 body
 365 {
 366     //printf("_arraySliceExpAddSliceAssign_f()\n");
 367     auto aptr = a.ptr;
 368     auto aend = aptr + a.length;
 369     auto bptr = b.ptr;
 370
 371     version (D_InlineAsm_X86)
 372     {
 373         // SSE version is 665% faster
 374         if (sse() && a.length >= 16)
 375         {
 376             auto n = aptr + (a.length & ~15);
 377
 378             // Unaligned case
 379             asm
 380             {
 381                 mov EAX, bptr;
 382                 mov ESI, aptr;
 383                 mov EDI, n;
 384                 movss XMM4, value;
 385                 shufps XMM4, XMM4, 0;
 386
 387                 align 8;
 388             startsseloop:
 389                 add ESI, 64;
 390                 movups XMM0, [EAX];
 391                 movups XMM1, [EAX+16];
 392                 movups XMM2, [EAX+32];
 393                 movups XMM3, [EAX+48];
 394                 add EAX, 64;
 395                 addps XMM0, XMM4;
 396                 addps XMM1, XMM4;
 397                 addps XMM2, XMM4;
 398                 addps XMM3, XMM4;
 399                 movups [ESI+ 0-64], XMM0;
 400                 movups [ESI+16-64], XMM1;
 401                 movups [ESI+32-64], XMM2;
 402                 movups [ESI+48-64], XMM3;
 403                 cmp ESI, EDI;
 404                 jb startsseloop;
 405
 406                 mov aptr, ESI;
 407                 mov bptr, EAX;
 408             }
 409         }
 410         else
 411         // 3DNow! version is 69% faster
 412         if (amd3dnow() && a.length >= 8)
 413         {
 414             auto n = aptr + (a.length & ~7);
 415
 416             ulong w = *cast(uint *) &value;
 417             ulong v = w | (w << 32L);
 418
 419             asm
 420             {
 421                 mov ESI, aptr;
 422                 mov EDI, n;
 423                 mov EAX, bptr;
 424                 movq MM4, qword ptr [v];
 425
 426                 align 8;
 427             start3dnow:
 428                 movq MM0, [EAX];
 429                 movq MM1, [EAX+8];
 430                 movq MM2, [EAX+16];
 431                 movq MM3, [EAX+24];
 432                 pfadd MM0, MM4;
 433                 pfadd MM1, MM4;
 434                 pfadd MM2, MM4;
 435                 pfadd MM3, MM4;
 436                 movq [ESI],    MM0;
 437                 movq [ESI+8],  MM1;
 438                 movq [ESI+16], MM2;
 439                 movq [ESI+24], MM3;
 440                 add ESI, 32;
 441                 add EAX, 32;
 442                 cmp ESI, EDI;
 443                 jb start3dnow;
 444
 445                 emms;
 446                 mov aptr, ESI;
 447                 mov bptr, EAX;
 448             }
 449         }
 450     }
 451
 452     while (aptr < aend)
 453         *aptr++ = *bptr++ + value;
 454
 455     return a;
 456 }
 457
 458 unittest
 459 {
 460     printf("_arraySliceExpAddSliceAssign_f unittest\n");
 461     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
 462     {
 463         version (log) printf("    cpuid %d\n", cpuid);
 464
 465         for (int j = 0; j < 2; j++)
 466         {
 467             const int dim = 67;
 468             T[] a = new T[dim + j];     // aligned on 16 byte boundary
 469             a = a[j .. dim + j];        // misalign for second iteration
 470             T[] b = new T[dim + j];
 471             b = b[j .. dim + j];
 472             T[] c = new T[dim + j];
 473             c = c[j .. dim + j];
 474
 475             for (int i = 0; i < dim; i++)
 476             {   a[i] = cast(T)i;
 477                 b[i] = cast(T)(i + 7);
 478                 c[i] = cast(T)(i * 2);
 479             }
 480
 481             c[] = a[] + 6;
 482
 483             for (int i = 0; i < dim; i++)
 484             {
 485                 if (c[i] != cast(T)(a[i] + 6))
 486                 {
 487                     printf("[%d]: %g != %g + 6\n", i, c[i], a[i]);
 488                     assert(0);
 489                 }
 490             }
 491         }
 492     }
 493 }
 494
 495 /* ======================================================================== */
 496
 497 /***********************
 498  * Computes:
 499  *      a[] += value
 500  */
 501
 502 T[] _arrayExpSliceAddass_f(T[] a, T value)
 503 {
 504     //printf("_arrayExpSliceAddass_f(a.length = %d, value = %Lg)\n", a.length, cast(real)value);
 505     auto aptr = a.ptr;
 506     auto aend = aptr + a.length;
 507
 508     version (D_InlineAsm_X86)
 509     {
 510         // SSE version is 302% faster
 511         if (sse() && a.length >= 16)
 512         {
 513             // align pointer
 514             auto n = cast(T*)((cast(uint)aptr + 15) & ~15);
 515             while (aptr < n)
 516                 *aptr++ += value;
 517             n = cast(T*)((cast(uint)aend) & ~15);
 518             if (aptr < n)
 519
 520             // Aligned case
 521             asm
 522             {
 523                 mov ESI, aptr;
 524                 mov EDI, n;
 525                 movss XMM4, value;
 526                 shufps XMM4, XMM4, 0;
 527
 528                 align 8;
 529             startsseloopa:
 530                 movaps XMM0, [ESI];
 531                 movaps XMM1, [ESI+16];
 532                 movaps XMM2, [ESI+32];
 533                 movaps XMM3, [ESI+48];
 534                 add ESI, 64;
 535                 addps XMM0, XMM4;
 536                 addps XMM1, XMM4;
 537                 addps XMM2, XMM4;
 538                 addps XMM3, XMM4;
 539                 movaps [ESI+ 0-64], XMM0;
 540                 movaps [ESI+16-64], XMM1;
 541                 movaps [ESI+32-64], XMM2;
 542                 movaps [ESI+48-64], XMM3;
 543                 cmp ESI, EDI;
 544                 jb startsseloopa;
 545
 546                 mov aptr, ESI;
 547             }
 548         }
 549         else
 550         // 3DNow! version is 63% faster
 551         if (amd3dnow() && a.length >= 8)
 552         {
 553             auto n = aptr + (a.length & ~7);
 554
 555             ulong w = *cast(uint *) &value;
 556             ulong v = w | (w << 32L);
 557
 558             asm
 559             {
 560                 mov ESI, dword ptr [aptr];
 561                 mov EDI, dword ptr [n];
 562                 movq MM4, qword ptr [v];
 563
 564                 align 8;
 565             start3dnow:
 566                 movq MM0, [ESI];
 567                 movq MM1, [ESI+8];
 568                 movq MM2, [ESI+16];
 569                 movq MM3, [ESI+24];
 570                 pfadd MM0, MM4;
 571                 pfadd MM1, MM4;
 572                 pfadd MM2, MM4;
 573                 pfadd MM3, MM4;
 574                 movq [ESI], MM0;
 575                 movq [ESI+8], MM1;
 576                 movq [ESI+16], MM2;
 577                 movq [ESI+24], MM3;
 578                 add ESI, 32;
 579                 cmp ESI, EDI;
 580                 jb start3dnow;
 581
 582                 emms;
 583                 mov dword ptr [aptr], ESI;
 584             }
 585         }
 586     }
 587
 588     while (aptr < aend)
 589         *aptr++ += value;
 590
 591     return a;
 592 }
 593
 594 unittest
 595 {
 596     printf("_arrayExpSliceAddass_f unittest\n");
 597     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
 598     {
 599         version (log) printf("    cpuid %d\n", cpuid);
 600
 601         for (int j = 0; j < 2; j++)
 602         {
 603             const int dim = 67;
 604             T[] a = new T[dim + j];     // aligned on 16 byte boundary
 605             a = a[j .. dim + j];        // misalign for second iteration
 606             T[] b = new T[dim + j];
 607             b = b[j .. dim + j];
 608             T[] c = new T[dim + j];
 609             c = c[j .. dim + j];
 610
 611             for (int i = 0; i < dim; i++)
 612             {   a[i] = cast(T)i;
 613                 b[i] = cast(T)(i + 7);
 614                 c[i] = cast(T)(i * 2);
 615             }
 616
 617             a[] = c[];
 618             c[] += 6;
 619
 620             for (int i = 0; i < dim; i++)
 621             {
 622                 if (c[i] != cast(T)(a[i] + 6))
 623                 {
 624                     printf("[%d]: %g != %g + 6\n", i, c[i], a[i]);
 625                     assert(0);
 626                 }
 627             }
 628         }
 629     }
 630 }
 631
 632 /* ======================================================================== */
 633
 634 /***********************
 635  * Computes:
 636  *      a[] += b[]
 637  */
 638
 639 T[] _arraySliceSliceAddass_f(T[] a, T[] b)
 640 in
 641 {
 642     assert (a.length == b.length);
 643     assert (disjoint(a, b));
 644 }
 645 body
 646 {
 647     //printf("_arraySliceSliceAddass_f()\n");
 648     auto aptr = a.ptr;
 649     auto aend = aptr + a.length;
 650     auto bptr = b.ptr;
 651
 652     version (D_InlineAsm_X86)
 653     {
 654         // SSE version is 468% faster
 655         if (sse() && a.length >= 16)
 656         {
 657             auto n = aptr + (a.length & ~15);
 658
 659             // Unaligned case
 660             asm
 661             {
 662                 mov ECX, bptr; // right operand
 663                 mov ESI, aptr; // destination operand
 664                 mov EDI, n; // end comparison
 665
 666                 align 8;
 667             startsseloopb:
 668                 movups XMM0, [ESI];
 669                 movups XMM1, [ESI+16];
 670                 movups XMM2, [ESI+32];
 671                 movups XMM3, [ESI+48];
 672                 add ESI, 64;
 673                 movups XMM4, [ECX];
 674                 movups XMM5, [ECX+16];
 675                 movups XMM6, [ECX+32];
 676                 movups XMM7, [ECX+48];
 677                 add ECX, 64;
 678                 addps XMM0, XMM4;
 679                 addps XMM1, XMM5;
 680                 addps XMM2, XMM6;
 681                 addps XMM3, XMM7;
 682                 movups [ESI+ 0-64], XMM0;
 683                 movups [ESI+16-64], XMM1;
 684                 movups [ESI+32-64], XMM2;
 685                 movups [ESI+48-64], XMM3;
 686                 cmp ESI, EDI;
 687                 jb startsseloopb;
 688
 689                 mov aptr, ESI;
 690                 mov bptr, ECX;
 691             }
 692         }
 693         else
 694         // 3DNow! version is 57% faster
 695         if (amd3dnow() && a.length >= 8)
 696         {
 697             auto n = aptr + (a.length & ~7);
 698
 699             asm
 700             {
 701                 mov ESI, dword ptr [aptr]; // destination operand
 702                 mov EDI, dword ptr [n];    // end comparison
 703                 mov ECX, dword ptr [bptr]; // right operand
 704
 705                 align 4;
 706             start3dnow:
 707                 movq MM0, [ESI];
 708                 movq MM1, [ESI+8];
 709                 movq MM2, [ESI+16];
 710                 movq MM3, [ESI+24];
 711                 pfadd MM0, [ECX];
 712                 pfadd MM1, [ECX+8];
 713                 pfadd MM2, [ECX+16];
 714                 pfadd MM3, [ECX+24];
 715                 movq [ESI], MM0;
 716                 movq [ESI+8], MM1;
 717                 movq [ESI+16], MM2;
 718                 movq [ESI+24], MM3;
 719                 add ESI, 32;
 720                 add ECX, 32;
 721                 cmp ESI, EDI;
 722                 jb start3dnow;
 723
 724                 emms;
 725                 mov dword ptr [aptr], ESI;
 726                 mov dword ptr [bptr], ECX;
 727             }
 728         }
 729     }
 730
 731     while (aptr < aend)
 732         *aptr++ += *bptr++;
 733
 734     return a;
 735 }
 736
 737 unittest
 738 {
 739     printf("_arraySliceSliceAddass_f unittest\n");
 740     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
 741     {
 742         version (log) printf("    cpuid %d\n", cpuid);
 743
 744         for (int j = 0; j < 2; j++)
 745         {
 746             const int dim = 67;
 747             T[] a = new T[dim + j];     // aligned on 16 byte boundary
 748             a = a[j .. dim + j];        // misalign for second iteration
 749             T[] b = new T[dim + j];
 750             b = b[j .. dim + j];
 751             T[] c = new T[dim + j];
 752             c = c[j .. dim + j];
 753
 754             for (int i = 0; i < dim; i++)
 755             {   a[i] = cast(T)i;
 756                 b[i] = cast(T)(i + 7);
 757                 c[i] = cast(T)(i * 2);
 758             }
 759
 760             a[] = c[];
 761             c[] += b[];
 762
 763             for (int i = 0; i < dim; i++)
 764             {
 765                 if (c[i] != cast(T)(a[i] + b[i]))
 766                 {
 767                     printf("[%d]: %g != %g + %g\n", i, c[i], a[i], b[i]);
 768                     assert(0);
 769                 }
 770             }
 771         }
 772     }
 773 }
 774
 775 /* ======================================================================== */
 776
 777 /***********************
 778  * Computes:
 779  *      a[] = b[] - value
 780  */
 781
 782 T[] _arraySliceExpMinSliceAssign_f(T[] a, T value, T[] b)
 783 in
 784 {
 785     assert (a.length == b.length);
 786     assert (disjoint(a, b));
 787 }
 788 body
 789 {
 790     //printf("_arraySliceExpMinSliceAssign_f()\n");
 791     auto aptr = a.ptr;
 792     auto aend = aptr + a.length;
 793     auto bptr = b.ptr;
 794
 795     version (D_InlineAsm_X86)
 796     {
 797         // SSE version is 622% faster
 798         if (sse() && a.length >= 16)
 799         {
 800             auto n = aptr + (a.length & ~15);
 801
 802             // Unaligned case
 803             asm
 804             {
 805                 mov EAX, bptr;
 806                 mov ESI, aptr;
 807                 mov EDI, n;
 808                 movss XMM4, value;
 809                 shufps XMM4, XMM4, 0;
 810
 811                 align 8;
 812             startsseloop:
 813                 add ESI, 64;
 814                 movups XMM0, [EAX];
 815                 movups XMM1, [EAX+16];
 816                 movups XMM2, [EAX+32];
 817                 movups XMM3, [EAX+48];
 818                 add EAX, 64;
 819                 subps XMM0, XMM4;
 820                 subps XMM1, XMM4;
 821                 subps XMM2, XMM4;
 822                 subps XMM3, XMM4;
 823                 movups [ESI+ 0-64], XMM0;
 824                 movups [ESI+16-64], XMM1;
 825                 movups [ESI+32-64], XMM2;
 826                 movups [ESI+48-64], XMM3;
 827                 cmp ESI, EDI;
 828                 jb startsseloop;
 829
 830                 mov aptr, ESI;
 831                 mov bptr, EAX;
 832             }
 833         }
 834         else
 835         // 3DNow! version is 67% faster
 836         if (amd3dnow() && a.length >= 8)
 837         {
 838             auto n = aptr + (a.length & ~7);
 839
 840             T[2] w;
 841
 842             w[0] = w[1] = value;
 843
 844             asm
 845             {
 846                 mov ESI, dword ptr [aptr];
 847                 mov EDI, dword ptr [n];
 848                 mov EAX, dword ptr [bptr];
 849                 movq MM4, qword ptr [w];
 850
 851                 align 8;
 852             start3dnow:
 853                 movq MM0, [EAX];
 854                 movq MM1, [EAX+8];
 855                 movq MM2, [EAX+16];
 856                 movq MM3, [EAX+24];
 857                 pfsub MM0, MM4;
 858                 pfsub MM1, MM4;
 859                 pfsub MM2, MM4;
 860                 pfsub MM3, MM4;
 861                 movq [ESI], MM0;
 862                 movq [ESI+8], MM1;
 863                 movq [ESI+16], MM2;
 864                 movq [ESI+24], MM3;
 865                 add ESI, 32;
 866                 add EAX, 32;
 867                 cmp ESI, EDI;
 868                 jb start3dnow;
 869
 870                 emms;
 871                 mov dword ptr [aptr], ESI;
 872                 mov dword ptr [bptr], EAX;
 873             }
 874         }
 875     }
 876
 877     while (aptr < aend)
 878         *aptr++ = *bptr++ - value;
 879
 880     return a;
 881 }
 882
 883 unittest
 884 {
 885     printf("_arraySliceExpMinSliceAssign_f unittest\n");
 886     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
 887     {
 888         version (log) printf("    cpuid %d\n", cpuid);
 889
 890         for (int j = 0; j < 2; j++)
 891         {
 892             const int dim = 67;
 893             T[] a = new T[dim + j];     // aligned on 16 byte boundary
 894             a = a[j .. dim + j];        // misalign for second iteration
 895             T[] b = new T[dim + j];
 896             b = b[j .. dim + j];
 897             T[] c = new T[dim + j];
 898             c = c[j .. dim + j];
 899
 900             for (int i = 0; i < dim; i++)
 901             {   a[i] = cast(T)i;
 902                 b[i] = cast(T)(i + 7);
 903                 c[i] = cast(T)(i * 2);
 904             }
 905
 906             c[] = a[] - 6;
 907
 908             for (int i = 0; i < dim; i++)
 909             {
 910                 if (c[i] != cast(T)(a[i] - 6))
 911                 {
 912                     printf("[%d]: %g != %g - 6\n", i, c[i], a[i]);
 913                     assert(0);
 914                 }
 915             }
 916         }
 917     }
 918 }
 919
 920 /* ======================================================================== */
 921
 922 /***********************
 923  * Computes:
 924  *      a[] = value - b[]
 925  */
 926
 927 T[] _arrayExpSliceMinSliceAssign_f(T[] a, T[] b, T value)
 928 in
 929 {
 930     assert (a.length == b.length);
 931     assert (disjoint(a, b));
 932 }
 933 body
 934 {
 935     //printf("_arrayExpSliceMinSliceAssign_f()\n");
 936     auto aptr = a.ptr;
 937     auto aend = aptr + a.length;
 938     auto bptr = b.ptr;
 939
 940     version (D_InlineAsm_X86)
 941     {
 942         // SSE version is 690% faster
 943         if (sse() && a.length >= 16)
 944         {
 945             auto n = aptr + (a.length & ~15);
 946
 947             // Unaligned case
 948             asm
 949             {
 950                 mov EAX, bptr;
 951                 mov ESI, aptr;
 952                 mov EDI, n;
 953                 movss XMM4, value;
 954                 shufps XMM4, XMM4, 0;
 955
 956                 align 8;
 957             startsseloop:
 958                 add ESI, 64;
 959                 movaps XMM5, XMM4;
 960                 movaps XMM6, XMM4;
 961                 movups XMM0, [EAX];
 962                 movups XMM1, [EAX+16];
 963                 movups XMM2, [EAX+32];
 964                 movups XMM3, [EAX+48];
 965                 add EAX, 64;
 966                 subps XMM5, XMM0;
 967                 subps XMM6, XMM1;
 968                 movups [ESI+ 0-64], XMM5;
 969                 movups [ESI+16-64], XMM6;
 970                 movaps XMM5, XMM4;
 971                 movaps XMM6, XMM4;
 972                 subps XMM5, XMM2;
 973                 subps XMM6, XMM3;
 974                 movups [ESI+32-64], XMM5;
 975                 movups [ESI+48-64], XMM6;
 976                 cmp ESI, EDI;
 977                 jb startsseloop;
 978
 979                 mov aptr, ESI;
 980                 mov bptr, EAX;
 981             }
 982         }
 983         else
 984         // 3DNow! version is 67% faster
 985         if (amd3dnow() && a.length >= 8)
 986         {
 987             auto n = aptr + (a.length & ~7);
 988
 989             ulong w = *cast(uint *) &value;
 990             ulong v = w | (w << 32L);
 991
 992             asm
 993             {
 994                 mov ESI, aptr;
 995                 mov EDI, n;
 996                 mov EAX, bptr;
 997                 movq MM4, qword ptr [v];
 998
 999                 align 8;
1000             start3dnow:
1001                 movq MM0, [EAX];
1002                 movq MM1, [EAX+8];
1003                 movq MM2, [EAX+16];
1004                 movq MM3, [EAX+24];
1005                 pfsubr MM0, MM4;
1006                 pfsubr MM1, MM4;
1007                 pfsubr MM2, MM4;
1008                 pfsubr MM3, MM4;
1009                 movq [ESI], MM0;
1010                 movq [ESI+8], MM1;
1011                 movq [ESI+16], MM2;
1012                 movq [ESI+24], MM3;
1013                 add ESI, 32;
1014                 add EAX, 32;
1015                 cmp ESI, EDI;
1016                 jb start3dnow;
1017
1018                 emms;
1019                 mov aptr, ESI;
1020                 mov bptr, EAX;
1021             }
1022         }
1023     }
1024
1025     while (aptr < aend)
1026         *aptr++ = value - *bptr++;
1027
1028     return a;
1029 }
1030
1031 unittest
1032 {
1033     printf("_arrayExpSliceMinSliceAssign_f unittest\n");
1034     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1035     {
1036         version (log) printf("    cpuid %d\n", cpuid);
1037
1038         for (int j = 0; j < 2; j++)
1039         {
1040             const int dim = 67;
1041             T[] a = new T[dim + j];     // aligned on 16 byte boundary
1042             a = a[j .. dim + j];        // misalign for second iteration
1043             T[] b = new T[dim + j];
1044             b = b[j .. dim + j];
1045             T[] c = new T[dim + j];
1046             c = c[j .. dim + j];
1047
1048             for (int i = 0; i < dim; i++)
1049             {   a[i] = cast(T)i;
1050                 b[i] = cast(T)(i + 7);
1051                 c[i] = cast(T)(i * 2);
1052             }
1053
1054             c[] = 6 - a[];
1055
1056             for (int i = 0; i < dim; i++)
1057             {
1058                 if (c[i] != cast(T)(6 - a[i]))
1059                 {
1060                     printf("[%d]: %g != 6 - %g\n", i, c[i], a[i]);
1061                     assert(0);
1062                 }
1063             }
1064         }
1065     }
1066 }
1067
1068 /* ======================================================================== */
1069
1070 /***********************
1071  * Computes:
1072  *      a[] -= value
1073  */
1074
1075 T[] _arrayExpSliceMinass_f(T[] a, T value)
1076 {
1077     //printf("_arrayExpSliceMinass_f(a.length = %d, value = %Lg)\n", a.length, cast(real)value);
1078     auto aptr = a.ptr;
1079     auto aend = aptr + a.length;
1080
1081     version (D_InlineAsm_X86)
1082     {
1083         // SSE version is 304% faster
1084         if (sse() && a.length >= 16)
1085         {
1086             // align pointer
1087             auto n = cast(T*)((cast(uint)aptr + 15) & ~15);
1088             while (aptr < n)
1089                 *aptr++ -= value;
1090             n = cast(T*)((cast(uint)aend) & ~15);
1091             if (aptr < n)
1092
1093             // Aligned case
1094             asm
1095             {
1096                 mov ESI, aptr;
1097                 mov EDI, n;
1098                 movss XMM4, value;
1099                 shufps XMM4, XMM4, 0;
1100
1101                 align 8;
1102             startsseloopa:
1103                 movaps XMM0, [ESI];
1104                 movaps XMM1, [ESI+16];
1105                 movaps XMM2, [ESI+32];
1106                 movaps XMM3, [ESI+48];
1107                 add ESI, 64;
1108                 subps XMM0, XMM4;
1109                 subps XMM1, XMM4;
1110                 subps XMM2, XMM4;
1111                 subps XMM3, XMM4;
1112                 movaps [ESI+ 0-64], XMM0;
1113                 movaps [ESI+16-64], XMM1;
1114                 movaps [ESI+32-64], XMM2;
1115                 movaps [ESI+48-64], XMM3;
1116                 cmp ESI, EDI;
1117                 jb startsseloopa;
1118
1119                 mov aptr, ESI;
1120             }
1121         }
1122         else
1123         // 3DNow! version is 63% faster
1124         if (amd3dnow() && a.length >= 8)
1125         {
1126             auto n = aptr + (a.length & ~7);
1127
1128             ulong w = *cast(uint *) &value;
1129             ulong v = w | (w << 32L);
1130
1131             asm
1132             {
1133                 mov ESI, dword ptr [aptr];
1134                 mov EDI, dword ptr [n];
1135                 movq MM4, qword ptr [v];
1136
1137                 align 8;
1138             start:
1139                 movq MM0, [ESI];
1140                 movq MM1, [ESI+8];
1141                 movq MM2, [ESI+16];
1142                 movq MM3, [ESI+24];
1143                 pfsub MM0, MM4;
1144                 pfsub MM1, MM4;
1145                 pfsub MM2, MM4;
1146                 pfsub MM3, MM4;
1147                 movq [ESI], MM0;
1148                 movq [ESI+8], MM1;
1149                 movq [ESI+16], MM2;
1150                 movq [ESI+24], MM3;
1151                 add ESI, 32;
1152                 cmp ESI, EDI;
1153                 jb start;
1154
1155                 emms;
1156                 mov dword ptr [aptr], ESI;
1157             }
1158         }
1159     }
1160
1161     while (aptr < aend)
1162         *aptr++ -= value;
1163
1164     return a;
1165 }
1166
1167 unittest
1168 {
1169     printf("_arrayExpSliceminass_f unittest\n");
1170     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1171     {
1172         version (log) printf("    cpuid %d\n", cpuid);
1173
1174         for (int j = 0; j < 2; j++)
1175         {
1176             const int dim = 67;
1177             T[] a = new T[dim + j];     // aligned on 16 byte boundary
1178             a = a[j .. dim + j];        // misalign for second iteration
1179             T[] b = new T[dim + j];
1180             b = b[j .. dim + j];
1181             T[] c = new T[dim + j];
1182             c = c[j .. dim + j];
1183
1184             for (int i = 0; i < dim; i++)
1185             {   a[i] = cast(T)i;
1186                 b[i] = cast(T)(i + 7);
1187                 c[i] = cast(T)(i * 2);
1188             }
1189
1190             a[] = c[];
1191             c[] -= 6;
1192
1193             for (int i = 0; i < dim; i++)
1194             {
1195                 if (c[i] != cast(T)(a[i] - 6))
1196                 {
1197                     printf("[%d]: %g != %g - 6\n", i, c[i], a[i]);
1198                     assert(0);
1199                 }
1200             }
1201         }
1202     }
1203 }
1204
1205 /* ======================================================================== */
1206
1207 /***********************
1208  * Computes:
1209  *      a[] -= b[]
1210  */
1211
1212 T[] _arraySliceSliceMinass_f(T[] a, T[] b)
1213 in
1214 {
1215     assert (a.length == b.length);
1216     assert (disjoint(a, b));
1217 }
1218 body
1219 {
1220     //printf("_arraySliceSliceMinass_f()\n");
1221     auto aptr = a.ptr;
1222     auto aend = aptr + a.length;
1223     auto bptr = b.ptr;
1224
1225     version (D_InlineAsm_X86)
1226     {
1227         // SSE version is 468% faster
1228         if (sse() && a.length >= 16)
1229         {
1230             auto n = aptr + (a.length & ~15);
1231
1232             // Unaligned case
1233             asm
1234             {
1235                 mov ECX, bptr; // right operand
1236                 mov ESI, aptr; // destination operand
1237                 mov EDI, n; // end comparison
1238
1239                 align 8;
1240             startsseloopb:
1241                 movups XMM0, [ESI];
1242                 movups XMM1, [ESI+16];
1243                 movups XMM2, [ESI+32];
1244                 movups XMM3, [ESI+48];
1245                 add ESI, 64;
1246                 movups XMM4, [ECX];
1247                 movups XMM5, [ECX+16];
1248                 movups XMM6, [ECX+32];
1249                 movups XMM7, [ECX+48];
1250                 add ECX, 64;
1251                 subps XMM0, XMM4;
1252                 subps XMM1, XMM5;
1253                 subps XMM2, XMM6;
1254                 subps XMM3, XMM7;
1255                 movups [ESI+ 0-64], XMM0;
1256                 movups [ESI+16-64], XMM1;
1257                 movups [ESI+32-64], XMM2;
1258                 movups [ESI+48-64], XMM3;
1259                 cmp ESI, EDI;
1260                 jb startsseloopb;
1261
1262                 mov aptr, ESI;
1263                 mov bptr, ECX;
1264             }
1265         }
1266         else
1267         // 3DNow! version is 57% faster
1268         if (amd3dnow() && a.length >= 8)
1269         {
1270             auto n = aptr + (a.length & ~7);
1271
1272             asm
1273             {
1274                 mov ESI, dword ptr [aptr]; // destination operand
1275                 mov EDI, dword ptr [n]; // end comparison
1276                 mov ECX, dword ptr [bptr]; // right operand
1277
1278                 align 4;
1279             start:
1280                 movq MM0, [ESI];
1281                 movq MM1, [ESI+8];
1282                 movq MM2, [ESI+16];
1283                 movq MM3, [ESI+24];
1284                 pfsub MM0, [ECX];
1285                 pfsub MM1, [ECX+8];
1286                 pfsub MM2, [ECX+16];
1287                 pfsub MM3, [ECX+24];
1288                 movq [ESI], MM0;
1289                 movq [ESI+8], MM1;
1290                 movq [ESI+16], MM2;
1291                 movq [ESI+24], MM3;
1292                 add ESI, 32;
1293                 add ECX, 32;
1294                 cmp ESI, EDI;
1295                 jb start;
1296
1297                 emms;
1298                 mov aptr, ESI;
1299                 mov bptr, ECX;
1300             }
1301         }
1302     }
1303
1304     while (aptr < aend)
1305         *aptr++ -= *bptr++;
1306
1307     return a;
1308 }
1309
1310 unittest
1311 {
1312     printf("_arrayExpSliceMinass_f unittest\n");
1313     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1314     {
1315         version (log) printf("    cpuid %d\n", cpuid);
1316
1317         for (int j = 0; j < 2; j++)
1318         {
1319             const int dim = 67;
1320             T[] a = new T[dim + j];     // aligned on 16 byte boundary
1321             a = a[j .. dim + j];        // misalign for second iteration
1322             T[] b = new T[dim + j];
1323             b = b[j .. dim + j];
1324             T[] c = new T[dim + j];
1325             c = c[j .. dim + j];
1326
1327             for (int i = 0; i < dim; i++)
1328             {   a[i] = cast(T)i;
1329                 b[i] = cast(T)(i + 7);
1330                 c[i] = cast(T)(i * 2);
1331             }
1332
1333             a[] = c[];
1334             c[] -= 6;
1335
1336             for (int i = 0; i < dim; i++)
1337             {
1338                 if (c[i] != cast(T)(a[i] - 6))
1339                 {
1340                     printf("[%d]: %g != %g - 6\n", i, c[i], a[i]);
1341                     assert(0);
1342                 }
1343             }
1344         }
1345     }
1346 }
1347
1348 /* ======================================================================== */
1349
1350 /***********************
1351  * Computes:
1352  *      a[] = b[] * value
1353  */
1354
1355 T[] _arraySliceExpMulSliceAssign_f(T[] a, T value, T[] b)
1356 in
1357 {
1358     assert(a.length == b.length);
1359     assert(disjoint(a, b));
1360 }
1361 body
1362 {
1363     //printf("_arraySliceExpMulSliceAssign_f()\n");
1364     auto aptr = a.ptr;
1365     auto aend = aptr + a.length;
1366     auto bptr = b.ptr;
1367
1368     version (D_InlineAsm_X86)
1369     {
1370         // SSE version is 607% faster
1371         if (sse() && a.length >= 16)
1372         {
1373             auto n = aptr + (a.length & ~15);
1374
1375             // Unaligned case
1376             asm
1377             {
1378                 mov EAX, bptr;
1379                 mov ESI, aptr;
1380                 mov EDI, n;
1381                 movss XMM4, value;
1382                 shufps XMM4, XMM4, 0;
1383
1384                 align 8;
1385             startsseloop:
1386                 add ESI, 64;
1387                 movups XMM0, [EAX];
1388                 movups XMM1, [EAX+16];
1389                 movups XMM2, [EAX+32];
1390                 movups XMM3, [EAX+48];
1391                 add EAX, 64;
1392                 mulps XMM0, XMM4;
1393                 mulps XMM1, XMM4;
1394                 mulps XMM2, XMM4;
1395                 mulps XMM3, XMM4;
1396                 movups [ESI+ 0-64], XMM0;
1397                 movups [ESI+16-64], XMM1;
1398                 movups [ESI+32-64], XMM2;
1399                 movups [ESI+48-64], XMM3;
1400                 cmp ESI, EDI;
1401                 jb startsseloop;
1402
1403                 mov aptr, ESI;
1404                 mov bptr, EAX;
1405             }
1406         }
1407         else
1408         // 3DNow! version is 69% faster
1409         if (amd3dnow() && a.length >= 8)
1410         {
1411             auto n = aptr + (a.length & ~7);
1412
1413             ulong w = *cast(uint *) &value;
1414             ulong v = w | (w << 32L);
1415
1416             asm
1417             {
1418                 mov ESI, dword ptr [aptr];
1419                 mov EDI, dword ptr [n];
1420                 mov EAX, dword ptr [bptr];
1421                 movq MM4, qword ptr [v];
1422
1423                 align 8;
1424             start:
1425                 movq MM0, [EAX];
1426                 movq MM1, [EAX+8];
1427                 movq MM2, [EAX+16];
1428                 movq MM3, [EAX+24];
1429                 pfmul MM0, MM4;
1430                 pfmul MM1, MM4;
1431                 pfmul MM2, MM4;
1432                 pfmul MM3, MM4;
1433                 movq [ESI], MM0;
1434                 movq [ESI+8], MM1;
1435                 movq [ESI+16], MM2;
1436                 movq [ESI+24], MM3;
1437                 add ESI, 32;
1438                 add EAX, 32;
1439                 cmp ESI, EDI;
1440                 jb start;
1441
1442                 emms;
1443                 mov aptr, ESI;
1444                 mov bptr, EAX;
1445             }
1446         }
1447     }
1448
1449     while (aptr < aend)
1450         *aptr++ = *bptr++ * value;
1451
1452     return a;
1453 }
1454
1455 unittest
1456 {
1457     printf("_arraySliceExpMulSliceAssign_f unittest\n");
1458     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1459     {
1460         version (log) printf("    cpuid %d\n", cpuid);
1461
1462         for (int j = 0; j < 2; j++)
1463         {
1464             const int dim = 67;
1465             T[] a = new T[dim + j];     // aligned on 16 byte boundary
1466             a = a[j .. dim + j];        // misalign for second iteration
1467             T[] b = new T[dim + j];
1468             b = b[j .. dim + j];
1469             T[] c = new T[dim + j];
1470             c = c[j .. dim + j];
1471
1472             for (int i = 0; i < dim; i++)
1473             {   a[i] = cast(T)i;
1474                 b[i] = cast(T)(i + 7);
1475                 c[i] = cast(T)(i * 2);
1476             }
1477
1478             c[] = a[] * 6;
1479
1480             for (int i = 0; i < dim; i++)
1481             {
1482                 if (c[i] != cast(T)(a[i] * 6))
1483                 {
1484                     printf("[%d]: %g != %g * 6\n", i, c[i], a[i]);
1485                     assert(0);
1486                 }
1487             }
1488         }
1489     }
1490 }
1491
1492 /* ======================================================================== */
1493
1494 /***********************
1495  * Computes:
1496  *      a[] = b[] * c[]
1497  */
1498
1499 T[] _arraySliceSliceMulSliceAssign_f(T[] a, T[] c, T[] b)
1500 in
1501 {
1502         assert(a.length == b.length && b.length == c.length);
1503         assert(disjoint(a, b));
1504         assert(disjoint(a, c));
1505         assert(disjoint(b, c));
1506 }
1507 body
1508 {
1509     //printf("_arraySliceSliceMulSliceAssign_f()\n");
1510     auto aptr = a.ptr;
1511     auto aend = aptr + a.length;
1512     auto bptr = b.ptr;
1513     auto cptr = c.ptr;
1514
1515     version (D_InlineAsm_X86)
1516     {
1517         // SSE version is 833% faster
1518         if (sse() && a.length >= 16)
1519         {
1520             auto n = aptr + (a.length & ~15);
1521
1522             // Unaligned case
1523             asm
1524             {
1525                 mov EAX, bptr; // left operand
1526                 mov ECX, cptr; // right operand
1527                 mov ESI, aptr; // destination operand
1528                 mov EDI, n; // end comparison
1529
1530                 align 8;
1531             startsseloopb:
1532                 movups XMM0, [EAX];
1533                 movups XMM1, [EAX+16];
1534                 movups XMM2, [EAX+32];
1535                 movups XMM3, [EAX+48];
1536                 add ESI, 64;
1537                 movups XMM4, [ECX];
1538                 movups XMM5, [ECX+16];
1539                 movups XMM6, [ECX+32];
1540                 movups XMM7, [ECX+48];
1541                 add EAX, 64;
1542                 mulps XMM0, XMM4;
1543                 mulps XMM1, XMM5;
1544                 mulps XMM2, XMM6;
1545                 mulps XMM3, XMM7;
1546                 add ECX, 64;
1547                 movups [ESI+ 0-64], XMM0;
1548                 movups [ESI+16-64], XMM1;
1549                 movups [ESI+32-64], XMM2;
1550                 movups [ESI+48-64], XMM3;
1551                 cmp ESI, EDI;
1552                 jb startsseloopb;
1553
1554                 mov aptr, ESI;
1555                 mov bptr, EAX;
1556                 mov cptr, ECX;
1557             }
1558         }
1559         else
1560         // 3DNow! version is only 13% faster
1561         if (amd3dnow() && a.length >= 8)
1562         {
1563             auto n = aptr + (a.length & ~7);
1564
1565             asm
1566             {
1567                 mov ESI, dword ptr [aptr]; // destination operand
1568                 mov EDI, dword ptr [n]; // end comparison
1569                 mov EAX, dword ptr [bptr]; // left operand
1570                 mov ECX, dword ptr [cptr]; // right operand
1571
1572                 align 4;
1573             start:
1574                 movq MM0, [EAX];
1575                 movq MM1, [EAX+8];
1576                 movq MM2, [EAX+16];
1577                 movq MM3, [EAX+24];
1578                 pfmul MM0, [ECX];
1579                 pfmul MM1, [ECX+8];
1580                 pfmul MM2, [ECX+16];
1581                 pfmul MM3, [ECX+24];
1582                 movq [ESI], MM0;
1583                 movq [ESI+8], MM1;
1584                 movq [ESI+16], MM2;
1585                 movq [ESI+24], MM3;
1586                 add ECX, 32;
1587                 add ESI, 32;
1588                 add EAX, 32;
1589                 cmp ESI, EDI;
1590                 jb start;
1591
1592                 emms;
1593                 mov aptr, ESI;
1594                 mov bptr, EAX;
1595                 mov cptr, ECX;
1596             }
1597         }
1598     }
1599
1600     while (aptr < aend)
1601         *aptr++ = *bptr++ * *cptr++;
1602
1603     return a;
1604 }
1605
1606 unittest
1607 {
1608     printf("_arraySliceSliceMulSliceAssign_f unittest\n");
1609     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1610     {
1611         version (log) printf("    cpuid %d\n", cpuid);
1612
1613         for (int j = 0; j < 2; j++)
1614         {
1615             const int dim = 67;
1616             T[] a = new T[dim + j];     // aligned on 16 byte boundary
1617             a = a[j .. dim + j];        // misalign for second iteration
1618             T[] b = new T[dim + j];
1619             b = b[j .. dim + j];
1620             T[] c = new T[dim + j];
1621             c = c[j .. dim + j];
1622
1623             for (int i = 0; i < dim; i++)
1624             {   a[i] = cast(T)i;
1625                 b[i] = cast(T)(i + 7);
1626                 c[i] = cast(T)(i * 2);
1627             }
1628
1629             c[] = a[] * b[];
1630
1631             for (int i = 0; i < dim; i++)
1632             {
1633                 if (c[i] != cast(T)(a[i] * b[i]))
1634                 {
1635                     printf("[%d]: %g != %g * %g\n", i, c[i], a[i], b[i]);
1636                     assert(0);
1637                 }
1638             }
1639         }
1640     }
1641 }
1642
1643 /* ======================================================================== */
1644
1645 /***********************
1646  * Computes:
1647  *      a[] *= value
1648  */
1649
1650 T[] _arrayExpSliceMulass_f(T[] a, T value)
1651 {
1652     //printf("_arrayExpSliceMulass_f(a.length = %d, value = %Lg)\n", a.length, cast(real)value);
1653     auto aptr = a.ptr;
1654     auto aend = aptr + a.length;
1655
1656     version (D_InlineAsm_X86)
1657     {
1658         // SSE version is 303% faster
1659         if (sse() && a.length >= 16)
1660         {
1661             // align pointer
1662             auto n = cast(T*)((cast(uint)aptr + 15) & ~15);
1663             while (aptr < n)
1664                 *aptr++ *= value;
1665             n = cast(T*)((cast(uint)aend) & ~15);
1666             if (aptr < n)
1667
1668             // Aligned case
1669             asm
1670             {
1671                 mov ESI, aptr;
1672                 mov EDI, n;
1673                 movss XMM4, value;
1674                 shufps XMM4, XMM4, 0;
1675
1676                 align 8;
1677             startsseloopa:
1678                 movaps XMM0, [ESI];
1679                 movaps XMM1, [ESI+16];
1680                 movaps XMM2, [ESI+32];
1681                 movaps XMM3, [ESI+48];
1682                 add ESI, 64;
1683                 mulps XMM0, XMM4;
1684                 mulps XMM1, XMM4;
1685                 mulps XMM2, XMM4;
1686                 mulps XMM3, XMM4;
1687                 movaps [ESI+ 0-64], XMM0;
1688                 movaps [ESI+16-64], XMM1;
1689                 movaps [ESI+32-64], XMM2;
1690                 movaps [ESI+48-64], XMM3;
1691                 cmp ESI, EDI;
1692                 jb startsseloopa;
1693
1694                 mov aptr, ESI;
1695             }
1696         }
1697         else
1698         // 3DNow! version is 63% faster
1699         if (amd3dnow() && a.length >= 8)
1700         {
1701             auto n = aptr + (a.length & ~7);
1702
1703             ulong w = *cast(uint *) &value;
1704             ulong v = w | (w << 32L);
1705
1706             asm
1707             {
1708                 mov ESI, dword ptr [aptr];
1709                 mov EDI, dword ptr [n];
1710                 movq MM4, qword ptr [v];
1711
1712                 align 8;
1713             start:
1714                 movq MM0, [ESI];
1715                 movq MM1, [ESI+8];
1716                 movq MM2, [ESI+16];
1717                 movq MM3, [ESI+24];
1718                 pfmul MM0, MM4;
1719                 pfmul MM1, MM4;
1720                 pfmul MM2, MM4;
1721                 pfmul MM3, MM4;
1722                 movq [ESI], MM0;
1723                 movq [ESI+8], MM1;
1724                 movq [ESI+16], MM2;
1725                 movq [ESI+24], MM3;
1726                 add ESI, 32;
1727                 cmp ESI, EDI;
1728                 jb start;
1729
1730                 emms;
1731                 mov dword ptr [aptr], ESI;
1732             }
1733         }
1734     }
1735
1736     while (aptr < aend)
1737         *aptr++ *= value;
1738
1739     return a;
1740 }
1741
1742 unittest
1743 {
1744     printf("_arrayExpSliceMulass_f unittest\n");
1745     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1746     {
1747         version (log) printf("    cpuid %d\n", cpuid);
1748
1749         for (int j = 0; j < 2; j++)
1750         {
1751             const int dim = 67;
1752             T[] a = new T[dim + j];     // aligned on 16 byte boundary
1753             a = a[j .. dim + j];        // misalign for second iteration
1754             T[] b = new T[dim + j];
1755             b = b[j .. dim + j];
1756             T[] c = new T[dim + j];
1757             c = c[j .. dim + j];
1758
1759             for (int i = 0; i < dim; i++)
1760             {   a[i] = cast(T)i;
1761                 b[i] = cast(T)(i + 7);
1762                 c[i] = cast(T)(i * 2);
1763             }
1764
1765             a[] = c[];
1766             c[] *= 6;
1767
1768             for (int i = 0; i < dim; i++)
1769             {
1770                 if (c[i] != cast(T)(a[i] * 6))
1771                 {
1772                     printf("[%d]: %g != %g * 6\n", i, c[i], a[i]);
1773                     assert(0);
1774                 }
1775             }
1776         }
1777     }
1778 }
1779
1780 /* ======================================================================== */
1781
1782 /***********************
1783  * Computes:
1784  *      a[] *= b[]
1785  */
1786
1787 T[] _arraySliceSliceMulass_f(T[] a, T[] b)
1788 in
1789 {
1790     assert (a.length == b.length);
1791     assert (disjoint(a, b));
1792 }
1793 body
1794 {
1795     //printf("_arraySliceSliceMulass_f()\n");
1796     auto aptr = a.ptr;
1797     auto aend = aptr + a.length;
1798     auto bptr = b.ptr;
1799
1800     version (D_InlineAsm_X86)
1801     {
1802         // SSE version is 525% faster
1803         if (sse() && a.length >= 16)
1804         {
1805             auto n = aptr + (a.length & ~15);
1806
1807             // Unaligned case
1808             asm
1809             {
1810                 mov ECX, bptr; // right operand
1811                 mov ESI, aptr; // destination operand
1812                 mov EDI, n; // end comparison
1813
1814                 align 8;
1815             startsseloopb:
1816                 movups XMM0, [ESI];
1817                 movups XMM1, [ESI+16];
1818                 movups XMM2, [ESI+32];
1819                 movups XMM3, [ESI+48];
1820                 add ESI, 64;
1821                 movups XMM4, [ECX];
1822                 movups XMM5, [ECX+16];
1823                 movups XMM6, [ECX+32];
1824                 movups XMM7, [ECX+48];
1825                 add ECX, 64;
1826                 mulps XMM0, XMM4;
1827                 mulps XMM1, XMM5;
1828                 mulps XMM2, XMM6;
1829                 mulps XMM3, XMM7;
1830                 movups [ESI+ 0-64], XMM0;
1831                 movups [ESI+16-64], XMM1;
1832                 movups [ESI+32-64], XMM2;
1833                 movups [ESI+48-64], XMM3;
1834                 cmp ESI, EDI;
1835                 jb startsseloopb;
1836
1837                 mov aptr, ESI;
1838                 mov bptr, ECX;
1839             }
1840         }
1841         else
1842         // 3DNow! version is 57% faster
1843         if (amd3dnow() && a.length >= 8)
1844         {
1845             auto n = aptr + (a.length & ~7);
1846
1847             asm
1848             {
1849                 mov ESI, dword ptr [aptr]; // destination operand
1850                 mov EDI, dword ptr [n]; // end comparison
1851                 mov ECX, dword ptr [bptr]; // right operand
1852
1853                 align 4;
1854             start:
1855                 movq MM0, [ESI];
1856                 movq MM1, [ESI+8];
1857                 movq MM2, [ESI+16];
1858                 movq MM3, [ESI+24];
1859                 pfmul MM0, [ECX];
1860                 pfmul MM1, [ECX+8];
1861                 pfmul MM2, [ECX+16];
1862                 pfmul MM3, [ECX+24];
1863                 movq [ESI], MM0;
1864                 movq [ESI+8], MM1;
1865                 movq [ESI+16], MM2;
1866                 movq [ESI+24], MM3;
1867                 add ESI, 32;
1868                 add ECX, 32;
1869                 cmp ESI, EDI;
1870                 jb start;
1871
1872                 emms;
1873                 mov dword ptr [aptr], ESI;
1874                 mov dword ptr [bptr], ECX;
1875             }
1876         }
1877     }
1878
1879     while (aptr < aend)
1880         *aptr++ *= *bptr++;
1881
1882     return a;
1883 }
1884
1885 unittest
1886 {
1887     printf("_arrayExpSliceMulass_f unittest\n");
1888     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
1889     {
1890         version (log) printf("    cpuid %d\n", cpuid);
1891
1892         for (int j = 0; j < 2; j++)
1893         {
1894             const int dim = 67;
1895             T[] a = new T[dim + j];     // aligned on 16 byte boundary
1896             a = a[j .. dim + j];        // misalign for second iteration
1897             T[] b = new T[dim + j];
1898             b = b[j .. dim + j];
1899             T[] c = new T[dim + j];
1900             c = c[j .. dim + j];
1901
1902             for (int i = 0; i < dim; i++)
1903             {   a[i] = cast(T)i;
1904                 b[i] = cast(T)(i + 7);
1905                 c[i] = cast(T)(i * 2);
1906             }
1907
1908             a[] = c[];
1909             c[] *= 6;
1910
1911             for (int i = 0; i < dim; i++)
1912             {
1913                 if (c[i] != cast(T)(a[i] * 6))
1914                 {
1915                     printf("[%d]: %g != %g * 6\n", i, c[i], a[i]);
1916                     assert(0);
1917                 }
1918             }
1919         }
1920     }
1921 }
1922
1923 /* ======================================================================== */
1924
1925 /***********************
1926  * Computes:
1927  *      a[] = b[] / value
1928  */
1929
1930 T[] _arraySliceExpDivSliceAssign_f(T[] a, T value, T[] b)
1931 in
1932 {
1933     assert(a.length == b.length);
1934     assert(disjoint(a, b));
1935 }
1936 body
1937 {
1938     //printf("_arraySliceExpDivSliceAssign_f()\n");
1939     auto aptr = a.ptr;
1940     auto aend = aptr + a.length;
1941     auto bptr = b.ptr;
1942
1943     /* Multiplying by the reciprocal is faster, but does
1944      * not produce as accurate an answer.
1945      */
1946     T recip = cast(T)1 / value;
1947
1948     version (D_InlineAsm_X86)
1949     {
1950         // SSE version is 587% faster
1951         if (sse() && a.length >= 16)
1952         {
1953             auto n = aptr + (a.length & ~15);
1954
1955             // Unaligned case
1956             asm
1957             {
1958                 mov EAX, bptr;
1959                 mov ESI, aptr;
1960                 mov EDI, n;
1961                 movss XMM4, recip;
1962                 //movss XMM4, value
1963                 //rcpss XMM4, XMM4
1964                 shufps XMM4, XMM4, 0;
1965
1966                 align 8;
1967             startsseloop:
1968                 add ESI, 64;
1969                 movups XMM0, [EAX];
1970                 movups XMM1, [EAX+16];
1971                 movups XMM2, [EAX+32];
1972                 movups XMM3, [EAX+48];
1973                 add EAX, 64;
1974                 mulps XMM0, XMM4;
1975                 mulps XMM1, XMM4;
1976                 mulps XMM2, XMM4;
1977                 mulps XMM3, XMM4;
1978                 //divps XMM0, XMM4;
1979                 //divps XMM1, XMM4;
1980                 //divps XMM2, XMM4;
1981                 //divps XMM3, XMM4;
1982                 movups [ESI+ 0-64], XMM0;
1983                 movups [ESI+16-64], XMM1;
1984                 movups [ESI+32-64], XMM2;
1985                 movups [ESI+48-64], XMM3;
1986                 cmp ESI, EDI;
1987                 jb startsseloop;
1988
1989                 mov aptr, ESI;
1990                 mov bptr, EAX;
1991             }
1992         }
1993         else
1994         // 3DNow! version is 72% faster
1995         if (amd3dnow() && a.length >= 8)
1996         {
1997             auto n = aptr + (a.length & ~7);
1998
1999             T[2] w = void;
2000
2001             w[0] = recip;
2002             w[1] = recip;
2003
2004             asm
2005             {
2006                 mov ESI, dword ptr [aptr];
2007                 mov EDI, dword ptr [n];
2008                 mov EAX, dword ptr [bptr];
2009                 movq MM4, qword ptr [w];
2010
2011                 align 8;
2012             start:
2013                 movq MM0, [EAX];
2014                 movq MM1, [EAX+8];
2015                 movq MM2, [EAX+16];
2016                 movq MM3, [EAX+24];
2017                 pfmul MM0, MM4;
2018                 pfmul MM1, MM4;
2019                 pfmul MM2, MM4;
2020                 pfmul MM3, MM4;
2021                 movq [ESI], MM0;
2022                 movq [ESI+8], MM1;
2023                 movq [ESI+16], MM2;
2024                 movq [ESI+24], MM3;
2025                 add ESI, 32;
2026                 add EAX, 32;
2027                 cmp ESI, EDI;
2028                 jb start;
2029
2030                 emms;
2031                 mov dword ptr [aptr], ESI;
2032                 mov dword ptr [bptr], EAX;
2033             }
2034         }
2035     }
2036
2037     while (aptr < aend)
2038         *aptr++ = *bptr++ * recip;
2039
2040     return a;
2041 }
2042
2043 unittest
2044 {
2045     printf("_arraySliceExpDivSliceAssign_f unittest\n");
2046     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
2047     {
2048         version (log) printf("    cpuid %d\n", cpuid);
2049
2050         for (int j = 0; j < 2; j++)
2051         {
2052             const int dim = 67;
2053             T[] a = new T[dim + j];     // aligned on 16 byte boundary
2054             a = a[j .. dim + j];        // misalign for second iteration
2055             T[] b = new T[dim + j];
2056             b = b[j .. dim + j];
2057             T[] c = new T[dim + j];
2058             c = c[j .. dim + j];
2059
2060             for (int i = 0; i < dim; i++)
2061             {   a[i] = cast(T)i;
2062                 b[i] = cast(T)(i + 7);
2063                 c[i] = cast(T)(i * 2);
2064             }
2065
2066             c[] = a[] / 8;
2067
2068             for (int i = 0; i < dim; i++)
2069             {
2070                 if (c[i] != cast(T)(a[i] / 8))
2071                 {
2072                     printf("[%d]: %g != %g / 8\n", i, c[i], a[i]);
2073                     assert(0);
2074                 }
2075             }
2076         }
2077     }
2078 }
2079
2080 /* ======================================================================== */
2081
2082 /***********************
2083  * Computes:
2084  *      a[] /= value
2085  */
2086
2087 T[] _arrayExpSliceDivass_f(T[] a, T value)
2088 {
2089     //printf("_arrayExpSliceDivass_f(a.length = %d, value = %Lg)\n", a.length, cast(real)value);
2090     auto aptr = a.ptr;
2091     auto aend = aptr + a.length;
2092
2093     /* Multiplying by the reciprocal is faster, but does
2094      * not produce as accurate an answer.
2095      */
2096     T recip = cast(T)1 / value;
2097
2098     version (D_InlineAsm_X86)
2099     {
2100         // SSE version is 245% faster
2101         if (sse() && a.length >= 16)
2102         {
2103             // align pointer
2104             auto n = cast(T*)((cast(uint)aptr + 15) & ~15);
2105             while (aptr < n)
2106                 *aptr++ *= recip;
2107             n = cast(T*)((cast(uint)aend) & ~15);
2108             if (aptr < n)
2109
2110             // Aligned case
2111             asm
2112             {
2113                 mov ESI, aptr;
2114                 mov EDI, n;
2115                 movss XMM4, recip;
2116                 //movss XMM4, value
2117                 //rcpss XMM4, XMM4
2118                 shufps XMM4, XMM4, 0;
2119
2120                 align 8;
2121             startsseloopa:
2122                 movaps XMM0, [ESI];
2123                 movaps XMM1, [ESI+16];
2124                 movaps XMM2, [ESI+32];
2125                 movaps XMM3, [ESI+48];
2126                 add ESI, 64;
2127                 mulps XMM0, XMM4;
2128                 mulps XMM1, XMM4;
2129                 mulps XMM2, XMM4;
2130                 mulps XMM3, XMM4;
2131                 //divps XMM0, XMM4;
2132                 //divps XMM1, XMM4;
2133                 //divps XMM2, XMM4;
2134                 //divps XMM3, XMM4;
2135                 movaps [ESI+ 0-64], XMM0;
2136                 movaps [ESI+16-64], XMM1;
2137                 movaps [ESI+32-64], XMM2;
2138                 movaps [ESI+48-64], XMM3;
2139                 cmp ESI, EDI;
2140                 jb startsseloopa;
2141
2142                 mov aptr, ESI;
2143             }
2144         }
2145         else
2146         // 3DNow! version is 57% faster
2147         if (amd3dnow() && a.length >= 8)
2148         {
2149             auto n = aptr + (a.length & ~7);
2150
2151             T[2] w = void;
2152
2153             w[0] = w[1] = recip;
2154
2155             asm
2156             {
2157                 mov ESI, dword ptr [aptr];
2158                 mov EDI, dword ptr [n];
2159                 movq MM4, qword ptr [w];
2160
2161                 align 8;
2162             start:
2163                 movq MM0, [ESI];
2164                 movq MM1, [ESI+8];
2165                 movq MM2, [ESI+16];
2166                 movq MM3, [ESI+24];
2167                 pfmul MM0, MM4;
2168                 pfmul MM1, MM4;
2169                 pfmul MM2, MM4;
2170                 pfmul MM3, MM4;
2171                 movq [ESI], MM0;
2172                 movq [ESI+8], MM1;
2173                 movq [ESI+16], MM2;
2174                 movq [ESI+24], MM3;
2175                 add ESI, 32;
2176                 cmp ESI, EDI;
2177                 jb start;
2178
2179                 emms;
2180                 mov dword ptr [aptr], ESI;
2181             }
2182         }
2183     }
2184
2185     while (aptr < aend)
2186         *aptr++ *= recip;
2187
2188     return a;
2189 }
2190
2191 unittest
2192 {
2193     printf("_arrayExpSliceDivass_f unittest\n");
2194     for (cpuid = 0; cpuid < CPUID_MAX; cpuid++)
2195     {
2196         version (log) printf("    cpuid %d\n", cpuid);
2197
2198         for (int j = 0; j < 2; j++)
2199         {
2200             const int dim = 67;
2201             T[] a = new T[dim + j];     // aligned on 16 byte boundary
2202             a = a[j .. dim + j];        // misalign for second iteration
2203             T[] b = new T[dim + j];
2204             b = b[j .. dim + j];
2205             T[] c = new T[dim + j];
2206             c = c[j .. dim + j];
2207
2208             for (int i = 0; i < dim; i++)
2209             {   a[i] = cast(T)i;
2210                 b[i] = cast(T)(i + 7);
2211                 c[i] = cast(T)(i * 2);
2212             }
2213
2214             a[] = c[];
2215             c[] /= 8;
2216
2217             for (int i = 0; i < dim; i++)
2218             {
2219                 if (c[i] != cast(T)(a[i] / 8))
2220                 {
2221                     printf("[%d]: %g != %g / 8\n", i, c[i], a[i]);
2222                     assert(0);
2223                 }
2224             }
2225         }
2226     }
2227 }
2228
2229
2230 /* ======================================================================== */
2231
2232 /***********************
2233  * Computes:
2234  *      a[] -= b[] * value
2235  */
2236
2237 T[] _arraySliceExpMulSliceMinass_f(T[] a, T value, T[] b)
2238 {
2239     return _arraySliceExpMulSliceAddass_f(a, -value, b);
2240 }
2241
2242 /***********************
2243  * Computes:
2244  *      a[] += b[] * value
2245  */
2246
2247 T[] _arraySliceExpMulSliceAddass_f(T[] a, T value, T[] b)
2248 in
2249 {
2250         assert(a.length == b.length);
2251         assert(disjoint(a, b));
2252 }
2253 body
2254 {
2255     auto aptr = a.ptr;
2256     auto aend = aptr + a.length;
2257     auto bptr = b.ptr;
2258
2259     // Handle remainder
2260     while (aptr < aend)
2261         *aptr++ += *bptr++ * value;
2262
2263     return a;
2264 }
2265
2266 unittest
2267 {
2268     printf("_arraySliceExpMulSliceAddass_f unittest\n");
2269
2270     cpuid = 1;
2271     {
2272         version (log) printf("    cpuid %d\n", cpuid);
2273
2274         for (int j = 0; j < 1; j++)
2275         {
2276             const int dim = 67;
2277             T[] a = new T[dim + j];     // aligned on 16 byte boundary
2278             a = a[j .. dim + j];        // misalign for second iteration
2279             T[] b = new T[dim + j];
2280             b = b[j .. dim + j];
2281             T[] c = new T[dim + j];
2282             c = c[j .. dim + j];
2283
2284             for (int i = 0; i < dim; i++)
2285             {   a[i] = cast(T)i;
2286                 b[i] = cast(T)(i + 7);
2287                 c[i] = cast(T)(i * 2);
2288             }
2289
2290             b[] = c[];
2291             c[] += a[] * 6;
2292
2293             for (int i = 0; i < dim; i++)
2294             {
2295                 //printf("[%d]: %g ?= %g + %g * 6\n", i, c[i], b[i], a[i]);
2296                 if (c[i] != cast(T)(b[i] + a[i] * 6))
2297                 {
2298                     printf("[%d]: %g ?= %g + %g * 6\n", i, c[i], b[i], a[i]);
2299                     assert(0);
2300                 }
2301             }
2302         }
2303     }
2304 }