src/compiler/dmd/arraydouble.d

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