src/gallium/auxiliary/gallivm/lp_bld_arit.c

   1 /**************************************************************************
   2  *
   3  * Copyright 2009-2010 VMware, Inc.
   4  * All Rights Reserved.
   5  *
   6  * Permission is hereby granted, free of charge, to any person obtaining a
   7  * copy of this software and associated documentation files (the
   8  * "Software"), to deal in the Software without restriction, including
   9  * without limitation the rights to use, copy, modify, merge, publish,
  10  * distribute, sub license, and/or sell copies of the Software, and to
  11  * permit persons to whom the Software is furnished to do so, subject to
  12  * the following conditions:
  13  *
  14  * The above copyright notice and this permission notice (including the
  15  * next paragraph) shall be included in all copies or substantial portions
  16  * of the Software.
  17  *
  18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  19  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  20  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NON-INFRINGEMENT.
  21  * IN NO EVENT SHALL VMWARE AND/OR ITS SUPPLIERS BE LIABLE FOR
  22  * ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  23  * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  24  * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  25  *
  26  **************************************************************************/
  27
  28
  29 /**
  30  * @file
  31  * Helper
  32  *
  33  * LLVM IR doesn't support all basic arithmetic operations we care about (most
  34  * notably min/max and saturated operations), and it is often necessary to
  35  * resort machine-specific intrinsics directly. The functions here hide all
  36  * these implementation details from the other modules.
  37  *
  38  * We also do simple expressions simplification here. Reasons are:
  39  * - it is very easy given we have all necessary information readily available
  40  * - LLVM optimization passes fail to simplify several vector expressions
  41  * - We often know value constraints which the optimization passes have no way
  42  *   of knowing, such as when source arguments are known to be in [0, 1] range.
  43  *
  44  * @author Jose Fonseca <jfonseca@vmware.com>
  45  */
  46
  47
  48 #include "util/u_memory.h"
  49 #include "util/u_debug.h"
  50 #include "util/u_math.h"
  51 #include "util/u_string.h"
  52 #include "util/u_cpu_detect.h"
  53
  54 #include "lp_bld_type.h"
  55 #include "lp_bld_const.h"
  56 #include "lp_bld_init.h"
  57 #include "lp_bld_intr.h"
  58 #include "lp_bld_logic.h"
  59 #include "lp_bld_pack.h"
  60 #include "lp_bld_debug.h"
  61 #include "lp_bld_arit.h"
  62
  63 #include "float.h"
  64
  65 #define EXP_POLY_DEGREE 5
  66
  67 #define LOG_POLY_DEGREE 4
  68
  69
  70 /**
  71  * Generate min(a, b)
  72  * No checks for special case values of a or b = 1 or 0 are done.
  73  */
  74 static LLVMValueRef
  75 lp_build_min_simple(struct lp_build_context *bld,
  76                     LLVMValueRef a,
  77                     LLVMValueRef b)
  78 {
  79    const struct lp_type type = bld->type;
  80    const char *intrinsic = NULL;
  81    unsigned intr_size = 0;
  82    LLVMValueRef cond;
  83
  84    assert(lp_check_value(type, a));
  85    assert(lp_check_value(type, b));
  86
  87    /* TODO: optimize the constant case */
  88
  89    if (type.floating && util_cpu_caps.has_sse) {
  90       if (type.width == 32) {
  91          if (type.length == 1) {
  92             intrinsic = "llvm.x86.sse.min.ss";
  93             intr_size = 128;
  94          }
  95          else if (type.length <= 4 || !util_cpu_caps.has_avx) {
  96             intrinsic = "llvm.x86.sse.min.ps";
  97             intr_size = 128;
  98          }
  99          else {
 100             intrinsic = "llvm.x86.avx.min.ps.256";
 101             intr_size = 256;
 102          }
 103       }
 104       if (type.width == 64 && util_cpu_caps.has_sse2) {
 105          if (type.length == 1) {
 106             intrinsic = "llvm.x86.sse2.min.sd";
 107             intr_size = 128;
 108          }
 109          else if (type.length == 2 || !util_cpu_caps.has_avx) {
 110             intrinsic = "llvm.x86.sse2.min.pd";
 111             intr_size = 128;
 112          }
 113          else {
 114             intrinsic = "llvm.x86.avx.min.pd.256";
 115             intr_size = 256;
 116          }
 117       }
 118    }
 119    else if (type.floating && util_cpu_caps.has_altivec) {
 120       if (type.width == 32 && type.length == 4) {
 121          intrinsic = "llvm.ppc.altivec.vminfp";
 122          intr_size = 128;
 123       }
 124    } else if (util_cpu_caps.has_sse2 && type.length >= 2) {
 125       intr_size = 128;
 126       if ((type.width == 8 || type.width == 16) &&
 127           (type.width * type.length <= 64) &&
 128           (gallivm_debug & GALLIVM_DEBUG_PERF)) {
 129          debug_printf("%s: inefficient code, bogus shuffle due to packing\n",
 130                       __FUNCTION__);
 131          }
 132       if (type.width == 8 && !type.sign) {
 133          intrinsic = "llvm.x86.sse2.pminu.b";
 134       }
 135       else if (type.width == 16 && type.sign) {
 136          intrinsic = "llvm.x86.sse2.pmins.w";
 137       }
 138       if (util_cpu_caps.has_sse4_1) {
 139          if (type.width == 8 && type.sign) {
 140             intrinsic = "llvm.x86.sse41.pminsb";
 141          }
 142          if (type.width == 16 && !type.sign) {
 143             intrinsic = "llvm.x86.sse41.pminuw";
 144          }
 145          if (type.width == 32 && !type.sign) {
 146             intrinsic = "llvm.x86.sse41.pminud";
 147         }
 148          if (type.width == 32 && type.sign) {
 149             intrinsic = "llvm.x86.sse41.pminsd";
 150          }
 151       }
 152    } else if (util_cpu_caps.has_altivec) {
 153      intr_size = 128;
 154      if (type.width == 8) {
 155        if (!type.sign) {
 156          intrinsic = "llvm.ppc.altivec.vminub";
 157        } else {
 158          intrinsic = "llvm.ppc.altivec.vminsb";
 159        }
 160      } else if (type.width == 16) {
 161        if (!type.sign) {
 162          intrinsic = "llvm.ppc.altivec.vminuh";
 163        } else {
 164          intrinsic = "llvm.ppc.altivec.vminsh";
 165        }
 166      } else if (type.width == 32) {
 167        if (!type.sign) {
 168          intrinsic = "llvm.ppc.altivec.vminuw";
 169        } else {
 170          intrinsic = "llvm.ppc.altivec.vminsw";
 171        }
 172      }
 173    }
 174
 175    if(intrinsic) {
 176       return lp_build_intrinsic_binary_anylength(bld->gallivm, intrinsic,
 177                                                  type,
 178                                                  intr_size, a, b);
 179    }
 180
 181    cond = lp_build_cmp(bld, PIPE_FUNC_LESS, a, b);
 182    return lp_build_select(bld, cond, a, b);
 183 }
 184
 185
 186 /**
 187  * Generate max(a, b)
 188  * No checks for special case values of a or b = 1 or 0 are done.
 189  */
 190 static LLVMValueRef
 191 lp_build_max_simple(struct lp_build_context *bld,
 192                     LLVMValueRef a,
 193                     LLVMValueRef b)
 194 {
 195    const struct lp_type type = bld->type;
 196    const char *intrinsic = NULL;
 197    unsigned intr_size = 0;
 198    LLVMValueRef cond;
 199
 200    assert(lp_check_value(type, a));
 201    assert(lp_check_value(type, b));
 202
 203    /* TODO: optimize the constant case */
 204
 205    if (type.floating && util_cpu_caps.has_sse) {
 206       if (type.width == 32) {
 207          if (type.length == 1) {
 208             intrinsic = "llvm.x86.sse.max.ss";
 209             intr_size = 128;
 210          }
 211          else if (type.length <= 4 || !util_cpu_caps.has_avx) {
 212             intrinsic = "llvm.x86.sse.max.ps";
 213             intr_size = 128;
 214          }
 215          else {
 216             intrinsic = "llvm.x86.avx.max.ps.256";
 217             intr_size = 256;
 218          }
 219       }
 220       if (type.width == 64 && util_cpu_caps.has_sse2) {
 221          if (type.length == 1) {
 222             intrinsic = "llvm.x86.sse2.max.sd";
 223             intr_size = 128;
 224          }
 225          else if (type.length == 2 || !util_cpu_caps.has_avx) {
 226             intrinsic = "llvm.x86.sse2.max.pd";
 227             intr_size = 128;
 228          }
 229          else {
 230             intrinsic = "llvm.x86.avx.max.pd.256";
 231             intr_size = 256;
 232          }
 233       }
 234    }
 235    else if (type.floating && util_cpu_caps.has_altivec) {
 236       if (type.width == 32 || type.length == 4) {
 237          intrinsic = "llvm.ppc.altivec.vmaxfp";
 238          intr_size = 128;
 239       }
 240    } else if (util_cpu_caps.has_sse2 && type.length >= 2) {
 241       intr_size = 128;
 242       if ((type.width == 8 || type.width == 16) &&
 243           (type.width * type.length <= 64) &&
 244           (gallivm_debug & GALLIVM_DEBUG_PERF)) {
 245          debug_printf("%s: inefficient code, bogus shuffle due to packing\n",
 246                       __FUNCTION__);
 247          }
 248       if (type.width == 8 && !type.sign) {
 249          intrinsic = "llvm.x86.sse2.pmaxu.b";
 250          intr_size = 128;
 251       }
 252       else if (type.width == 16 && type.sign) {
 253          intrinsic = "llvm.x86.sse2.pmaxs.w";
 254       }
 255       if (util_cpu_caps.has_sse4_1) {
 256          if (type.width == 8 && type.sign) {
 257             intrinsic = "llvm.x86.sse41.pmaxsb";
 258          }
 259          if (type.width == 16 && !type.sign) {
 260             intrinsic = "llvm.x86.sse41.pmaxuw";
 261          }
 262          if (type.width == 32 && !type.sign) {
 263             intrinsic = "llvm.x86.sse41.pmaxud";
 264         }
 265          if (type.width == 32 && type.sign) {
 266             intrinsic = "llvm.x86.sse41.pmaxsd";
 267          }
 268       }
 269    } else if (util_cpu_caps.has_altivec) {
 270      intr_size = 128;
 271      if (type.width == 8) {
 272        if (!type.sign) {
 273          intrinsic = "llvm.ppc.altivec.vmaxub";
 274        } else {
 275          intrinsic = "llvm.ppc.altivec.vmaxsb";
 276        }
 277      } else if (type.width == 16) {
 278        if (!type.sign) {
 279          intrinsic = "llvm.ppc.altivec.vmaxuh";
 280        } else {
 281          intrinsic = "llvm.ppc.altivec.vmaxsh";
 282        }
 283      } else if (type.width == 32) {
 284        if (!type.sign) {
 285          intrinsic = "llvm.ppc.altivec.vmaxuw";
 286        } else {
 287          intrinsic = "llvm.ppc.altivec.vmaxsw";
 288        }
 289      }
 290    }
 291
 292    if(intrinsic) {
 293       return lp_build_intrinsic_binary_anylength(bld->gallivm, intrinsic,
 294                                                  type,
 295                                                  intr_size, a, b);
 296    }
 297
 298    cond = lp_build_cmp(bld, PIPE_FUNC_GREATER, a, b);
 299    return lp_build_select(bld, cond, a, b);
 300 }
 301
 302
 303 /**
 304  * Generate 1 - a, or ~a depending on bld->type.
 305  */
 306 LLVMValueRef
 307 lp_build_comp(struct lp_build_context *bld,
 308               LLVMValueRef a)
 309 {
 310    LLVMBuilderRef builder = bld->gallivm->builder;
 311    const struct lp_type type = bld->type;
 312
 313    assert(lp_check_value(type, a));
 314
 315    if(a == bld->one)
 316       return bld->zero;
 317    if(a == bld->zero)
 318       return bld->one;
 319
 320    if(type.norm && !type.floating && !type.fixed && !type.sign) {
 321       if(LLVMIsConstant(a))
 322          return LLVMConstNot(a);
 323       else
 324          return LLVMBuildNot(builder, a, "");
 325    }
 326
 327    if(LLVMIsConstant(a))
 328       if (type.floating)
 329           return LLVMConstFSub(bld->one, a);
 330       else
 331           return LLVMConstSub(bld->one, a);
 332    else
 333       if (type.floating)
 334          return LLVMBuildFSub(builder, bld->one, a, "");
 335       else
 336          return LLVMBuildSub(builder, bld->one, a, "");
 337 }
 338
 339
 340 /**
 341  * Generate a + b
 342  */
 343 LLVMValueRef
 344 lp_build_add(struct lp_build_context *bld,
 345              LLVMValueRef a,
 346              LLVMValueRef b)
 347 {
 348    LLVMBuilderRef builder = bld->gallivm->builder;
 349    const struct lp_type type = bld->type;
 350    LLVMValueRef res;
 351
 352    assert(lp_check_value(type, a));
 353    assert(lp_check_value(type, b));
 354
 355    if(a == bld->zero)
 356       return b;
 357    if(b == bld->zero)
 358       return a;
 359    if(a == bld->undef || b == bld->undef)
 360       return bld->undef;
 361
 362    if(bld->type.norm) {
 363       const char *intrinsic = NULL;
 364
 365       if(a == bld->one || b == bld->one)
 366         return bld->one;
 367
 368       if (type.width * type.length == 128 &&
 369           !type.floating && !type.fixed) {
 370          if(util_cpu_caps.has_sse2) {
 371            if(type.width == 8)
 372              intrinsic = type.sign ? "llvm.x86.sse2.padds.b" : "llvm.x86.sse2.paddus.b";
 373            if(type.width == 16)
 374              intrinsic = type.sign ? "llvm.x86.sse2.padds.w" : "llvm.x86.sse2.paddus.w";
 375          } else if (util_cpu_caps.has_altivec) {
 376            if(type.width == 8)
 377               intrinsic = type.sign ? "llvm.ppc.altivec.vaddsbs" : "llvm.ppc.altivec.vaddubs";
 378            if(type.width == 16)
 379               intrinsic = type.sign ? "llvm.ppc.altivec.vaddsws" : "llvm.ppc.altivec.vadduws";
 380          }
 381       }
 382
 383       if(intrinsic)
 384          return lp_build_intrinsic_binary(builder, intrinsic, lp_build_vec_type(bld->gallivm, bld->type), a, b);
 385    }
 386
 387    if(LLVMIsConstant(a) && LLVMIsConstant(b))
 388       if (type.floating)
 389          res = LLVMConstFAdd(a, b);
 390       else
 391          res = LLVMConstAdd(a, b);
 392    else
 393       if (type.floating)
 394          res = LLVMBuildFAdd(builder, a, b, "");
 395       else
 396          res = LLVMBuildAdd(builder, a, b, "");
 397
 398    /* clamp to ceiling of 1.0 */
 399    if(bld->type.norm && (bld->type.floating || bld->type.fixed))
 400       res = lp_build_min_simple(bld, res, bld->one);
 401
 402    /* XXX clamp to floor of -1 or 0??? */
 403
 404    return res;
 405 }
 406
 407
 408 /** Return the scalar sum of the elements of a.
 409  * Should avoid this operation whenever possible.
 410  */
 411 LLVMValueRef
 412 lp_build_horizontal_add(struct lp_build_context *bld,
 413                         LLVMValueRef a)
 414 {
 415    LLVMBuilderRef builder = bld->gallivm->builder;
 416    const struct lp_type type = bld->type;
 417    LLVMValueRef index, res;
 418    unsigned i, length;
 419    LLVMValueRef shuffles1[LP_MAX_VECTOR_LENGTH / 2];
 420    LLVMValueRef shuffles2[LP_MAX_VECTOR_LENGTH / 2];
 421    LLVMValueRef vecres, elem2;
 422
 423    assert(lp_check_value(type, a));
 424
 425    if (type.length == 1) {
 426       return a;
 427    }
 428
 429    assert(!bld->type.norm);
 430
 431    /*
 432     * for byte vectors can do much better with psadbw.
 433     * Using repeated shuffle/adds here. Note with multiple vectors
 434     * this can be done more efficiently as outlined in the intel
 435     * optimization manual.
 436     * Note: could cause data rearrangement if used with smaller element
 437     * sizes.
 438     */
 439
 440    vecres = a;
 441    length = type.length / 2;
 442    while (length > 1) {
 443       LLVMValueRef vec1, vec2;
 444       for (i = 0; i < length; i++) {
 445          shuffles1[i] = lp_build_const_int32(bld->gallivm, i);
 446          shuffles2[i] = lp_build_const_int32(bld->gallivm, i + length);
 447       }
 448       vec1 = LLVMBuildShuffleVector(builder, vecres, vecres,
 449                                     LLVMConstVector(shuffles1, length), "");
 450       vec2 = LLVMBuildShuffleVector(builder, vecres, vecres,
 451                                     LLVMConstVector(shuffles2, length), "");
 452       if (type.floating) {
 453          vecres = LLVMBuildFAdd(builder, vec1, vec2, "");
 454       }
 455       else {
 456          vecres = LLVMBuildAdd(builder, vec1, vec2, "");
 457       }
 458       length = length >> 1;
 459    }
 460
 461    /* always have vector of size 2 here */
 462    assert(length == 1);
 463
 464    index = lp_build_const_int32(bld->gallivm, 0);
 465    res = LLVMBuildExtractElement(builder, vecres, index, "");
 466    index = lp_build_const_int32(bld->gallivm, 1);
 467    elem2 = LLVMBuildExtractElement(builder, vecres, index, "");
 468
 469    if (type.floating)
 470       res = LLVMBuildFAdd(builder, res, elem2, "");
 471     else
 472       res = LLVMBuildAdd(builder, res, elem2, "");
 473
 474    return res;
 475 }
 476
 477 /**
 478  * Return the horizontal sums of 4 float vectors as a float4 vector.
 479  * This uses the technique as outlined in Intel Optimization Manual.
 480  */
 481 static LLVMValueRef
 482 lp_build_horizontal_add4x4f(struct lp_build_context *bld,
 483                             LLVMValueRef src[4])
 484 {
 485    struct gallivm_state *gallivm = bld->gallivm;
 486    LLVMBuilderRef builder = gallivm->builder;
 487    LLVMValueRef shuffles[4];
 488    LLVMValueRef tmp[4];
 489    LLVMValueRef sumtmp[2], shuftmp[2];
 490
 491    /* lower half of regs */
 492    shuffles[0] = lp_build_const_int32(gallivm, 0);
 493    shuffles[1] = lp_build_const_int32(gallivm, 1);
 494    shuffles[2] = lp_build_const_int32(gallivm, 4);
 495    shuffles[3] = lp_build_const_int32(gallivm, 5);
 496    tmp[0] = LLVMBuildShuffleVector(builder, src[0], src[1],
 497                                    LLVMConstVector(shuffles, 4), "");
 498    tmp[2] = LLVMBuildShuffleVector(builder, src[2], src[3],
 499                                    LLVMConstVector(shuffles, 4), "");
 500
 501    /* upper half of regs */
 502    shuffles[0] = lp_build_const_int32(gallivm, 2);
 503    shuffles[1] = lp_build_const_int32(gallivm, 3);
 504    shuffles[2] = lp_build_const_int32(gallivm, 6);
 505    shuffles[3] = lp_build_const_int32(gallivm, 7);
 506    tmp[1] = LLVMBuildShuffleVector(builder, src[0], src[1],
 507                                    LLVMConstVector(shuffles, 4), "");
 508    tmp[3] = LLVMBuildShuffleVector(builder, src[2], src[3],
 509                                    LLVMConstVector(shuffles, 4), "");
 510
 511    sumtmp[0] = LLVMBuildFAdd(builder, tmp[0], tmp[1], "");
 512    sumtmp[1] = LLVMBuildFAdd(builder, tmp[2], tmp[3], "");
 513
 514    shuffles[0] = lp_build_const_int32(gallivm, 0);
 515    shuffles[1] = lp_build_const_int32(gallivm, 2);
 516    shuffles[2] = lp_build_const_int32(gallivm, 4);
 517    shuffles[3] = lp_build_const_int32(gallivm, 6);
 518    shuftmp[0] = LLVMBuildShuffleVector(builder, sumtmp[0], sumtmp[1],
 519                                        LLVMConstVector(shuffles, 4), "");
 520
 521    shuffles[0] = lp_build_const_int32(gallivm, 1);
 522    shuffles[1] = lp_build_const_int32(gallivm, 3);
 523    shuffles[2] = lp_build_const_int32(gallivm, 5);
 524    shuffles[3] = lp_build_const_int32(gallivm, 7);
 525    shuftmp[1] = LLVMBuildShuffleVector(builder, sumtmp[0], sumtmp[1],
 526                                        LLVMConstVector(shuffles, 4), "");
 527
 528    return LLVMBuildFAdd(builder, shuftmp[0], shuftmp[1], "");
 529 }
 530
 531
 532 /*
 533  * partially horizontally add 2-4 float vectors with length nx4,
 534  * i.e. only four adjacent values in each vector will be added,
 535  * assuming values are really grouped in 4 which also determines
 536  * output order.
 537  *
 538  * Return a vector of the same length as the initial vectors,
 539  * with the excess elements (if any) being undefined.
 540  * The element order is independent of number of input vectors.
 541  * For 3 vectors x0x1x2x3x4x5x6x7, y0y1y2y3y4y5y6y7, z0z1z2z3z4z5z6z7
 542  * the output order thus will be
 543  * sumx0-x3,sumy0-y3,sumz0-z3,undef,sumx4-x7,sumy4-y7,sumz4z7,undef
 544  */
 545 LLVMValueRef
 546 lp_build_hadd_partial4(struct lp_build_context *bld,
 547                        LLVMValueRef vectors[],
 548                        unsigned num_vecs)
 549 {
 550    struct gallivm_state *gallivm = bld->gallivm;
 551    LLVMBuilderRef builder = gallivm->builder;
 552    LLVMValueRef ret_vec;
 553    LLVMValueRef tmp[4];
 554    const char *intrinsic = NULL;
 555
 556    assert(num_vecs >= 2 && num_vecs <= 4);
 557    assert(bld->type.floating);
 558
 559    /* only use this with at least 2 vectors, as it is sort of expensive
 560     * (depending on cpu) and we always need two horizontal adds anyway,
 561     * so a shuffle/add approach might be better.
 562     */
 563
 564    tmp[0] = vectors[0];
 565    tmp[1] = vectors[1];
 566
 567    tmp[2] = num_vecs > 2 ? vectors[2] : vectors[0];
 568    tmp[3] = num_vecs > 3 ? vectors[3] : vectors[0];
 569
 570    if (util_cpu_caps.has_sse3 && bld->type.width == 32 &&
 571        bld->type.length == 4) {
 572       intrinsic = "llvm.x86.sse3.hadd.ps";
 573    }
 574    else if (util_cpu_caps.has_avx && bld->type.width == 32 &&
 575             bld->type.length == 8) {
 576       intrinsic = "llvm.x86.avx.hadd.ps.256";
 577    }
 578    if (intrinsic) {
 579       tmp[0] = lp_build_intrinsic_binary(builder, intrinsic,
 580                                        lp_build_vec_type(gallivm, bld->type),
 581                                        tmp[0], tmp[1]);
 582       if (num_vecs > 2) {
 583          tmp[1] = lp_build_intrinsic_binary(builder, intrinsic,
 584                                           lp_build_vec_type(gallivm, bld->type),
 585                                           tmp[2], tmp[3]);
 586       }
 587       else {
 588          tmp[1] = tmp[0];
 589       }
 590       return lp_build_intrinsic_binary(builder, intrinsic,
 591                                        lp_build_vec_type(gallivm, bld->type),
 592                                        tmp[0], tmp[1]);
 593    }
 594
 595    if (bld->type.length == 4) {
 596       ret_vec = lp_build_horizontal_add4x4f(bld, tmp);
 597    }
 598    else {
 599       LLVMValueRef partres[LP_MAX_VECTOR_LENGTH/4];
 600       unsigned j;
 601       unsigned num_iter = bld->type.length / 4;
 602       struct lp_type parttype = bld->type;
 603       parttype.length = 4;
 604       for (j = 0; j < num_iter; j++) {
 605          LLVMValueRef partsrc[4];
 606          unsigned i;
 607          for (i = 0; i < 4; i++) {
 608             partsrc[i] = lp_build_extract_range(gallivm, tmp[i], j*4, 4);
 609          }
 610          partres[j] = lp_build_horizontal_add4x4f(bld, partsrc);
 611       }
 612       ret_vec = lp_build_concat(gallivm, partres, parttype, num_iter);
 613    }
 614    return ret_vec;
 615 }
 616
 617 /**
 618  * Generate a - b
 619  */
 620 LLVMValueRef
 621 lp_build_sub(struct lp_build_context *bld,
 622              LLVMValueRef a,
 623              LLVMValueRef b)
 624 {
 625    LLVMBuilderRef builder = bld->gallivm->builder;
 626    const struct lp_type type = bld->type;
 627    LLVMValueRef res;
 628
 629    assert(lp_check_value(type, a));
 630    assert(lp_check_value(type, b));
 631
 632    if(b == bld->zero)
 633       return a;
 634    if(a == bld->undef || b == bld->undef)
 635       return bld->undef;
 636    if(a == b)
 637       return bld->zero;
 638
 639    if(bld->type.norm) {
 640       const char *intrinsic = NULL;
 641
 642       if(b == bld->one)
 643         return bld->zero;
 644
 645       if (type.width * type.length == 128 &&
 646           !type.floating && !type.fixed) {
 647          if (util_cpu_caps.has_sse2) {
 648            if(type.width == 8)
 649               intrinsic = type.sign ? "llvm.x86.sse2.psubs.b" : "llvm.x86.sse2.psubus.b";
 650            if(type.width == 16)
 651               intrinsic = type.sign ? "llvm.x86.sse2.psubs.w" : "llvm.x86.sse2.psubus.w";
 652          } else if (util_cpu_caps.has_altivec) {
 653            if(type.width == 8)
 654               intrinsic = type.sign ? "llvm.ppc.altivec.vsubsbs" : "llvm.ppc.altivec.vsububs";
 655            if(type.width == 16)
 656               intrinsic = type.sign ? "llvm.ppc.altivec.vsubsws" : "llvm.ppc.altivec.vsubuws";
 657          }
 658       }
 659
 660       if(intrinsic)
 661          return lp_build_intrinsic_binary(builder, intrinsic, lp_build_vec_type(bld->gallivm, bld->type), a, b);
 662    }
 663
 664    if(LLVMIsConstant(a) && LLVMIsConstant(b))
 665       if (type.floating)
 666          res = LLVMConstFSub(a, b);
 667       else
 668          res = LLVMConstSub(a, b);
 669    else
 670       if (type.floating)
 671          res = LLVMBuildFSub(builder, a, b, "");
 672       else
 673          res = LLVMBuildSub(builder, a, b, "");
 674
 675    if(bld->type.norm && (bld->type.floating || bld->type.fixed))
 676       res = lp_build_max_simple(bld, res, bld->zero);
 677
 678    return res;
 679 }
 680
 681
 682 /**
 683  * Normalized 8bit multiplication.
 684  *
 685  * - alpha plus one
 686  *
 687  *     makes the following approximation to the division (Sree)
 688  *
 689  *       a*b/255 ~= (a*(b + 1)) >> 256
 690  *
 691  *     which is the fastest method that satisfies the following OpenGL criteria
 692  *
 693  *       0*0 = 0 and 255*255 = 255
 694  *
 695  * - geometric series
 696  *
 697  *     takes the geometric series approximation to the division
 698  *
 699  *       t/255 = (t >> 8) + (t >> 16) + (t >> 24) ..
 700  *
 701  *     in this case just the first two terms to fit in 16bit arithmetic
 702  *
 703  *       t/255 ~= (t + (t >> 8)) >> 8
 704  *
 705  *     note that just by itself it doesn't satisfies the OpenGL criteria, as
 706  *     255*255 = 254, so the special case b = 255 must be accounted or roundoff
 707  *     must be used
 708  *
 709  * - geometric series plus rounding
 710  *
 711  *     when using a geometric series division instead of truncating the result
 712  *     use roundoff in the approximation (Jim Blinn)
 713  *
 714  *       t/255 ~= (t + (t >> 8) + 0x80) >> 8
 715  *
 716  *     achieving the exact results
 717  *
 718  * @sa Alvy Ray Smith, Image Compositing Fundamentals, Tech Memo 4, Aug 15, 1995,
 719  *     ftp://ftp.alvyray.com/Acrobat/4_Comp.pdf
 720  * @sa Michael Herf, The "double blend trick", May 2000,
 721  *     http://www.stereopsis.com/doubleblend.html
 722  */
 723 static LLVMValueRef
 724 lp_build_mul_u8n(struct gallivm_state *gallivm,
 725                  struct lp_type i16_type,
 726                  LLVMValueRef a, LLVMValueRef b)
 727 {
 728    LLVMBuilderRef builder = gallivm->builder;
 729    LLVMValueRef c8;
 730    LLVMValueRef ab;
 731
 732    assert(!i16_type.floating);
 733    assert(lp_check_value(i16_type, a));
 734    assert(lp_check_value(i16_type, b));
 735
 736    c8 = lp_build_const_int_vec(gallivm, i16_type, 8);
 737
 738 #if 0
 739
 740    /* a*b/255 ~= (a*(b + 1)) >> 256 */
 741    b = LLVMBuildAdd(builder, b, lp_build_const_int_vec(gallium, i16_type, 1), "");
 742    ab = LLVMBuildMul(builder, a, b, "");
 743
 744 #else
 745
 746    /* ab/255 ~= (ab + (ab >> 8) + 0x80) >> 8 */
 747    ab = LLVMBuildMul(builder, a, b, "");
 748    ab = LLVMBuildAdd(builder, ab, LLVMBuildLShr(builder, ab, c8, ""), "");
 749    ab = LLVMBuildAdd(builder, ab, lp_build_const_int_vec(gallivm, i16_type, 0x80), "");
 750
 751 #endif
 752
 753    ab = LLVMBuildLShr(builder, ab, c8, "");
 754
 755    return ab;
 756 }
 757
 758 /**
 759  * Normalized 16bit multiplication.
 760  *
 761  * Utilises same principle as above code.
 762  */
 763 static LLVMValueRef
 764 lp_build_mul_u16n(struct gallivm_state *gallivm,
 765                  struct lp_type i32_type,
 766                  LLVMValueRef a, LLVMValueRef b)
 767 {
 768    LLVMBuilderRef builder = gallivm->builder;
 769    LLVMValueRef c16;
 770    LLVMValueRef ab;
 771
 772    assert(!i32_type.floating);
 773    assert(lp_check_value(i32_type, a));
 774    assert(lp_check_value(i32_type, b));
 775
 776    c16 = lp_build_const_int_vec(gallivm, i32_type, 16);
 777
 778    /* ab/65535 ~= (ab + (ab >> 16) + 0x8000) >> 16 */
 779    ab = LLVMBuildMul(builder, a, b, "");
 780    ab = LLVMBuildAdd(builder, ab, LLVMBuildLShr(builder, ab, c16, ""), "");
 781    ab = LLVMBuildAdd(builder, ab, lp_build_const_int_vec(gallivm, i32_type, 0x8000), "");
 782
 783    ab = LLVMBuildLShr(builder, ab, c16, "");
 784
 785    return ab;
 786 }
 787
 788 /**
 789  * Generate a * b
 790  */
 791 LLVMValueRef
 792 lp_build_mul(struct lp_build_context *bld,
 793              LLVMValueRef a,
 794              LLVMValueRef b)
 795 {
 796    LLVMBuilderRef builder = bld->gallivm->builder;
 797    const struct lp_type type = bld->type;
 798    LLVMValueRef shift;
 799    LLVMValueRef res;
 800
 801    assert(lp_check_value(type, a));
 802    assert(lp_check_value(type, b));
 803
 804    if(a == bld->zero)
 805       return bld->zero;
 806    if(a == bld->one)
 807       return b;
 808    if(b == bld->zero)
 809       return bld->zero;
 810    if(b == bld->one)
 811       return a;
 812    if(a == bld->undef || b == bld->undef)
 813       return bld->undef;
 814
 815    if(!type.floating && !type.fixed && type.norm) {
 816       if(type.width == 8) {
 817          struct lp_type i16_type = lp_wider_type(type);
 818          LLVMValueRef al, ah, bl, bh, abl, abh, ab;
 819
 820          lp_build_unpack2(bld->gallivm, type, i16_type, a, &al, &ah);
 821          lp_build_unpack2(bld->gallivm, type, i16_type, b, &bl, &bh);
 822
 823          /* PMULLW, PSRLW, PADDW */
 824          abl = lp_build_mul_u8n(bld->gallivm, i16_type, al, bl);
 825          abh = lp_build_mul_u8n(bld->gallivm, i16_type, ah, bh);
 826
 827          ab = lp_build_pack2(bld->gallivm, i16_type, type, abl, abh);
 828
 829          return ab;
 830       }
 831
 832       if(type.width == 16) {
 833          struct lp_type i32_type = lp_wider_type(type);
 834          LLVMValueRef al, ah, bl, bh, abl, abh, ab;
 835
 836          lp_build_unpack2(bld->gallivm, type, i32_type, a, &al, &ah);
 837          lp_build_unpack2(bld->gallivm, type, i32_type, b, &bl, &bh);
 838
 839          /* PMULLW, PSRLW, PADDW */
 840          abl = lp_build_mul_u16n(bld->gallivm, i32_type, al, bl);
 841          abh = lp_build_mul_u16n(bld->gallivm, i32_type, ah, bh);
 842
 843          ab = lp_build_pack2(bld->gallivm, i32_type, type, abl, abh);
 844
 845          return ab;
 846       }
 847
 848       /* FIXME */
 849       assert(0);
 850    }
 851
 852    if(type.fixed)
 853       shift = lp_build_const_int_vec(bld->gallivm, type, type.width/2);
 854    else
 855       shift = NULL;
 856
 857    if(LLVMIsConstant(a) && LLVMIsConstant(b)) {
 858       if (type.floating)
 859          res = LLVMConstFMul(a, b);
 860       else
 861          res = LLVMConstMul(a, b);
 862       if(shift) {
 863          if(type.sign)
 864             res = LLVMConstAShr(res, shift);
 865          else
 866             res = LLVMConstLShr(res, shift);
 867       }
 868    }
 869    else {
 870       if (type.floating)
 871          res = LLVMBuildFMul(builder, a, b, "");
 872       else
 873          res = LLVMBuildMul(builder, a, b, "");
 874       if(shift) {
 875          if(type.sign)
 876             res = LLVMBuildAShr(builder, res, shift, "");
 877          else
 878             res = LLVMBuildLShr(builder, res, shift, "");
 879       }
 880    }
 881
 882    return res;
 883 }
 884
 885
 886 /**
 887  * Small vector x scale multiplication optimization.
 888  */
 889 LLVMValueRef
 890 lp_build_mul_imm(struct lp_build_context *bld,
 891                  LLVMValueRef a,
 892                  int b)
 893 {
 894    LLVMBuilderRef builder = bld->gallivm->builder;
 895    LLVMValueRef factor;
 896
 897    assert(lp_check_value(bld->type, a));
 898
 899    if(b == 0)
 900       return bld->zero;
 901
 902    if(b == 1)
 903       return a;
 904
 905    if(b == -1)
 906       return lp_build_negate(bld, a);
 907
 908    if(b == 2 && bld->type.floating)
 909       return lp_build_add(bld, a, a);
 910
 911    if(util_is_power_of_two(b)) {
 912       unsigned shift = ffs(b) - 1;
 913
 914       if(bld->type.floating) {
 915 #if 0
 916          /*
 917           * Power of two multiplication by directly manipulating the exponent.
 918           *
 919           * XXX: This might not be always faster, it will introduce a small error
 920           * for multiplication by zero, and it will produce wrong results
 921           * for Inf and NaN.
 922           */
 923          unsigned mantissa = lp_mantissa(bld->type);
 924          factor = lp_build_const_int_vec(bld->gallivm, bld->type, (unsigned long long)shift << mantissa);
 925          a = LLVMBuildBitCast(builder, a, lp_build_int_vec_type(bld->type), "");
 926          a = LLVMBuildAdd(builder, a, factor, "");
 927          a = LLVMBuildBitCast(builder, a, lp_build_vec_type(bld->gallivm, bld->type), "");
 928          return a;
 929 #endif
 930       }
 931       else {
 932          factor = lp_build_const_vec(bld->gallivm, bld->type, shift);
 933          return LLVMBuildShl(builder, a, factor, "");
 934       }
 935    }
 936
 937    factor = lp_build_const_vec(bld->gallivm, bld->type, (double)b);
 938    return lp_build_mul(bld, a, factor);
 939 }
 940
 941
 942 /**
 943  * Generate a / b
 944  */
 945 LLVMValueRef
 946 lp_build_div(struct lp_build_context *bld,
 947              LLVMValueRef a,
 948              LLVMValueRef b)
 949 {
 950    LLVMBuilderRef builder = bld->gallivm->builder;
 951    const struct lp_type type = bld->type;
 952
 953    assert(lp_check_value(type, a));
 954    assert(lp_check_value(type, b));
 955
 956    if(a == bld->zero)
 957       return bld->zero;
 958    if(a == bld->one)
 959       return lp_build_rcp(bld, b);
 960    if(b == bld->zero)
 961       return bld->undef;
 962    if(b == bld->one)
 963       return a;
 964    if(a == bld->undef || b == bld->undef)
 965       return bld->undef;
 966
 967    if(LLVMIsConstant(a) && LLVMIsConstant(b)) {
 968       if (type.floating)
 969          return LLVMConstFDiv(a, b);
 970       else if (type.sign)
 971          return LLVMConstSDiv(a, b);
 972       else
 973          return LLVMConstUDiv(a, b);
 974    }
 975
 976    if(((util_cpu_caps.has_sse && type.width == 32 && type.length == 4) ||
 977        (util_cpu_caps.has_avx && type.width == 32 && type.length == 8)) &&
 978       type.floating)
 979       return lp_build_mul(bld, a, lp_build_rcp(bld, b));
 980
 981    if (type.floating)
 982       return LLVMBuildFDiv(builder, a, b, "");
 983    else if (type.sign)
 984       return LLVMBuildSDiv(builder, a, b, "");
 985    else
 986       return LLVMBuildUDiv(builder, a, b, "");
 987 }
 988
 989
 990 /**
 991  * Linear interpolation -- without any checks.
 992  *
 993  * @sa http://www.stereopsis.com/doubleblend.html
 994  */
 995 static INLINE LLVMValueRef
 996 lp_build_lerp_simple(struct lp_build_context *bld,
 997                      LLVMValueRef x,
 998                      LLVMValueRef v0,
 999                      LLVMValueRef v1)
1000 {
1001    LLVMBuilderRef builder = bld->gallivm->builder;
1002    LLVMValueRef delta;
1003    LLVMValueRef res;
1004
1005    assert(lp_check_value(bld->type, x));
1006    assert(lp_check_value(bld->type, v0));
1007    assert(lp_check_value(bld->type, v1));
1008
1009    delta = lp_build_sub(bld, v1, v0);
1010
1011    res = lp_build_mul(bld, x, delta);
1012
1013    res = lp_build_add(bld, v0, res);
1014
1015    if (bld->type.fixed) {
1016       /* XXX: This step is necessary for lerping 8bit colors stored on 16bits,
1017        * but it will be wrong for other uses. Basically we need a more
1018        * powerful lp_type, capable of further distinguishing the values
1019        * interpretation from the value storage. */
1020       res = LLVMBuildAnd(builder, res, lp_build_const_int_vec(bld->gallivm, bld->type, (1 << bld->type.width/2) - 1), "");
1021    }
1022
1023    return res;
1024 }
1025
1026
1027 /**
1028  * Linear interpolation.
1029  */
1030 LLVMValueRef
1031 lp_build_lerp(struct lp_build_context *bld,
1032               LLVMValueRef x,
1033               LLVMValueRef v0,
1034               LLVMValueRef v1)
1035 {
1036    LLVMBuilderRef builder = bld->gallivm->builder;
1037    const struct lp_type type = bld->type;
1038    LLVMValueRef res;
1039
1040    assert(lp_check_value(type, x));
1041    assert(lp_check_value(type, v0));
1042    assert(lp_check_value(type, v1));
1043
1044    if (type.norm) {
1045       struct lp_type wide_type;
1046       struct lp_build_context wide_bld;
1047       LLVMValueRef xl, xh, v0l, v0h, v1l, v1h, resl, resh;
1048       LLVMValueRef shift;
1049
1050       assert(type.length >= 2);
1051       assert(!type.sign);
1052
1053       /*
1054        * Create a wider type, enough to hold the intermediate result of the
1055        * multiplication.
1056        */
1057       memset(&wide_type, 0, sizeof wide_type);
1058       wide_type.fixed  = TRUE;
1059       wide_type.width  = type.width*2;
1060       wide_type.length = type.length/2;
1061
1062       lp_build_context_init(&wide_bld, bld->gallivm, wide_type);
1063
1064       lp_build_unpack2(bld->gallivm, type, wide_type, x,  &xl,  &xh);
1065       lp_build_unpack2(bld->gallivm, type, wide_type, v0, &v0l, &v0h);
1066       lp_build_unpack2(bld->gallivm, type, wide_type, v1, &v1l, &v1h);
1067
1068       /*
1069        * Scale x from [0, 255] to [0, 256]
1070        */
1071
1072       shift = lp_build_const_int_vec(bld->gallivm, wide_type, type.width - 1);
1073
1074       xl = lp_build_add(&wide_bld, xl,
1075                         LLVMBuildAShr(builder, xl, shift, ""));
1076       xh = lp_build_add(&wide_bld, xh,
1077                         LLVMBuildAShr(builder, xh, shift, ""));
1078
1079       /*
1080        * Lerp both halves.
1081        */
1082
1083       resl = lp_build_lerp_simple(&wide_bld, xl, v0l, v1l);
1084       resh = lp_build_lerp_simple(&wide_bld, xh, v0h, v1h);
1085
1086       res = lp_build_pack2(bld->gallivm, wide_type, type, resl, resh);
1087    } else {
1088       res = lp_build_lerp_simple(bld, x, v0, v1);
1089    }
1090
1091    return res;
1092 }
1093
1094
1095 LLVMValueRef
1096 lp_build_lerp_2d(struct lp_build_context *bld,
1097                  LLVMValueRef x,
1098                  LLVMValueRef y,
1099                  LLVMValueRef v00,
1100                  LLVMValueRef v01,
1101                  LLVMValueRef v10,
1102                  LLVMValueRef v11)
1103 {
1104    LLVMValueRef v0 = lp_build_lerp(bld, x, v00, v01);
1105    LLVMValueRef v1 = lp_build_lerp(bld, x, v10, v11);
1106    return lp_build_lerp(bld, y, v0, v1);
1107 }
1108
1109
1110 /**
1111  * Generate min(a, b)
1112  * Do checks for special cases.
1113  */
1114 LLVMValueRef
1115 lp_build_min(struct lp_build_context *bld,
1116              LLVMValueRef a,
1117              LLVMValueRef b)
1118 {
1119    assert(lp_check_value(bld->type, a));
1120    assert(lp_check_value(bld->type, b));
1121
1122    if(a == bld->undef || b == bld->undef)
1123       return bld->undef;
1124
1125    if(a == b)
1126       return a;
1127
1128    if (bld->type.norm) {
1129       if (!bld->type.sign) {
1130          if (a == bld->zero || b == bld->zero) {
1131             return bld->zero;
1132          }
1133       }
1134       if(a == bld->one)
1135          return b;
1136       if(b == bld->one)
1137          return a;
1138    }
1139
1140    return lp_build_min_simple(bld, a, b);
1141 }
1142
1143
1144 /**
1145  * Generate max(a, b)
1146  * Do checks for special cases.
1147  */
1148 LLVMValueRef
1149 lp_build_max(struct lp_build_context *bld,
1150              LLVMValueRef a,
1151              LLVMValueRef b)
1152 {
1153    assert(lp_check_value(bld->type, a));
1154    assert(lp_check_value(bld->type, b));
1155
1156    if(a == bld->undef || b == bld->undef)
1157       return bld->undef;
1158
1159    if(a == b)
1160       return a;
1161
1162    if(bld->type.norm) {
1163       if(a == bld->one || b == bld->one)
1164          return bld->one;
1165       if (!bld->type.sign) {
1166          if (a == bld->zero) {
1167             return b;
1168          }
1169          if (b == bld->zero) {
1170             return a;
1171          }
1172       }
1173    }
1174
1175    return lp_build_max_simple(bld, a, b);
1176 }
1177
1178
1179 /**
1180  * Generate clamp(a, min, max)
1181  * Do checks for special cases.
1182  */
1183 LLVMValueRef
1184 lp_build_clamp(struct lp_build_context *bld,
1185                LLVMValueRef a,
1186                LLVMValueRef min,
1187                LLVMValueRef max)
1188 {
1189    assert(lp_check_value(bld->type, a));
1190    assert(lp_check_value(bld->type, min));
1191    assert(lp_check_value(bld->type, max));
1192
1193    a = lp_build_min(bld, a, max);
1194    a = lp_build_max(bld, a, min);
1195    return a;
1196 }
1197
1198
1199 /**
1200  * Generate abs(a)
1201  */
1202 LLVMValueRef
1203 lp_build_abs(struct lp_build_context *bld,
1204              LLVMValueRef a)
1205 {
1206    LLVMBuilderRef builder = bld->gallivm->builder;
1207    const struct lp_type type = bld->type;
1208    LLVMTypeRef vec_type = lp_build_vec_type(bld->gallivm, type);
1209
1210    assert(lp_check_value(type, a));
1211
1212    if(!type.sign)
1213       return a;
1214
1215    if(type.floating) {
1216       /* Mask out the sign bit */
1217       LLVMTypeRef int_vec_type = lp_build_int_vec_type(bld->gallivm, type);
1218       unsigned long long absMask = ~(1ULL << (type.width - 1));
1219       LLVMValueRef mask = lp_build_const_int_vec(bld->gallivm, type, ((unsigned long long) absMask));
1220       a = LLVMBuildBitCast(builder, a, int_vec_type, "");
1221       a = LLVMBuildAnd(builder, a, mask, "");
1222       a = LLVMBuildBitCast(builder, a, vec_type, "");
1223       return a;
1224    }
1225
1226    if(type.width*type.length == 128 && util_cpu_caps.has_ssse3) {
1227       switch(type.width) {
1228       case 8:
1229          return lp_build_intrinsic_unary(builder, "llvm.x86.ssse3.pabs.b.128", vec_type, a);
1230       case 16:
1231          return lp_build_intrinsic_unary(builder, "llvm.x86.ssse3.pabs.w.128", vec_type, a);
1232       case 32:
1233          return lp_build_intrinsic_unary(builder, "llvm.x86.ssse3.pabs.d.128", vec_type, a);
1234       }
1235    }
1236    else if (type.width*type.length == 256 && util_cpu_caps.has_ssse3 &&
1237             (gallivm_debug & GALLIVM_DEBUG_PERF) &&
1238             (type.width == 8 || type.width == 16 || type.width == 32)) {
1239       debug_printf("%s: inefficient code, should split vectors manually\n",
1240                    __FUNCTION__);
1241    }
1242
1243    return lp_build_max(bld, a, LLVMBuildNeg(builder, a, ""));
1244 }
1245
1246
1247 LLVMValueRef
1248 lp_build_negate(struct lp_build_context *bld,
1249                 LLVMValueRef a)
1250 {
1251    LLVMBuilderRef builder = bld->gallivm->builder;
1252
1253    assert(lp_check_value(bld->type, a));
1254
1255 #if HAVE_LLVM >= 0x0207
1256    if (bld->type.floating)
1257       a = LLVMBuildFNeg(builder, a, "");
1258    else
1259 #endif
1260       a = LLVMBuildNeg(builder, a, "");
1261
1262    return a;
1263 }
1264
1265
1266 /** Return -1, 0 or +1 depending on the sign of a */
1267 LLVMValueRef
1268 lp_build_sgn(struct lp_build_context *bld,
1269              LLVMValueRef a)
1270 {
1271    LLVMBuilderRef builder = bld->gallivm->builder;
1272    const struct lp_type type = bld->type;
1273    LLVMValueRef cond;
1274    LLVMValueRef res;
1275
1276    assert(lp_check_value(type, a));
1277
1278    /* Handle non-zero case */
1279    if(!type.sign) {
1280       /* if not zero then sign must be positive */
1281       res = bld->one;
1282    }
1283    else if(type.floating) {
1284       LLVMTypeRef vec_type;
1285       LLVMTypeRef int_type;
1286       LLVMValueRef mask;
1287       LLVMValueRef sign;
1288       LLVMValueRef one;
1289       unsigned long long maskBit = (unsigned long long)1 << (type.width - 1);
1290
1291       int_type = lp_build_int_vec_type(bld->gallivm, type);
1292       vec_type = lp_build_vec_type(bld->gallivm, type);
1293       mask = lp_build_const_int_vec(bld->gallivm, type, maskBit);
1294
1295       /* Take the sign bit and add it to 1 constant */
1296       sign = LLVMBuildBitCast(builder, a, int_type, "");
1297       sign = LLVMBuildAnd(builder, sign, mask, "");
1298       one = LLVMConstBitCast(bld->one, int_type);
1299       res = LLVMBuildOr(builder, sign, one, "");
1300       res = LLVMBuildBitCast(builder, res, vec_type, "");
1301    }
1302    else
1303    {
1304       /* signed int/norm/fixed point */
1305       /* could use psign with sse3 and appropriate vectors here */
1306       LLVMValueRef minus_one = lp_build_const_vec(bld->gallivm, type, -1.0);
1307       cond = lp_build_cmp(bld, PIPE_FUNC_GREATER, a, bld->zero);
1308       res = lp_build_select(bld, cond, bld->one, minus_one);
1309    }
1310
1311    /* Handle zero */
1312    cond = lp_build_cmp(bld, PIPE_FUNC_EQUAL, a, bld->zero);
1313    res = lp_build_select(bld, cond, bld->zero, res);
1314
1315    return res;
1316 }
1317
1318
1319 /**
1320  * Set the sign of float vector 'a' according to 'sign'.
1321  * If sign==0, return abs(a).
1322  * If sign==1, return -abs(a);
1323  * Other values for sign produce undefined results.
1324  */
1325 LLVMValueRef
1326 lp_build_set_sign(struct lp_build_context *bld,
1327                   LLVMValueRef a, LLVMValueRef sign)
1328 {
1329    LLVMBuilderRef builder = bld->gallivm->builder;
1330    const struct lp_type type = bld->type;
1331    LLVMTypeRef int_vec_type = lp_build_int_vec_type(bld->gallivm, type);
1332    LLVMTypeRef vec_type = lp_build_vec_type(bld->gallivm, type);
1333    LLVMValueRef shift = lp_build_const_int_vec(bld->gallivm, type, type.width - 1);
1334    LLVMValueRef mask = lp_build_const_int_vec(bld->gallivm, type,
1335                              ~((unsigned long long) 1 << (type.width - 1)));
1336    LLVMValueRef val, res;
1337
1338    assert(type.floating);
1339    assert(lp_check_value(type, a));
1340
1341    /* val = reinterpret_cast<int>(a) */
1342    val = LLVMBuildBitCast(builder, a, int_vec_type, "");
1343    /* val = val & mask */
1344    val = LLVMBuildAnd(builder, val, mask, "");
1345    /* sign = sign << shift */
1346    sign = LLVMBuildShl(builder, sign, shift, "");
1347    /* res = val | sign */
1348    res = LLVMBuildOr(builder, val, sign, "");
1349    /* res = reinterpret_cast<float>(res) */
1350    res = LLVMBuildBitCast(builder, res, vec_type, "");
1351
1352    return res;
1353 }
1354
1355
1356 /**
1357  * Convert vector of (or scalar) int to vector of (or scalar) float.
1358  */
1359 LLVMValueRef
1360 lp_build_int_to_float(struct lp_build_context *bld,
1361                       LLVMValueRef a)
1362 {
1363    LLVMBuilderRef builder = bld->gallivm->builder;
1364    const struct lp_type type = bld->type;
1365    LLVMTypeRef vec_type = lp_build_vec_type(bld->gallivm, type);
1366
1367    assert(type.floating);
1368
1369    return LLVMBuildSIToFP(builder, a, vec_type, "");
1370 }
1371
1372 static boolean
1373 sse41_rounding_available(const struct lp_type type)
1374 {
1375    if ((util_cpu_caps.has_sse4_1 &&
1376        (type.length == 1 || type.width*type.length == 128)) ||
1377        (util_cpu_caps.has_avx && type.width*type.length == 256))
1378       return TRUE;
1379
1380    return FALSE;
1381 }
1382
1383 enum lp_build_round_sse41_mode
1384 {
1385    LP_BUILD_ROUND_SSE41_NEAREST = 0,
1386    LP_BUILD_ROUND_SSE41_FLOOR = 1,
1387    LP_BUILD_ROUND_SSE41_CEIL = 2,
1388    LP_BUILD_ROUND_SSE41_TRUNCATE = 3
1389 };
1390
1391
1392 /**
1393  * Helper for SSE4.1's ROUNDxx instructions.
1394  *
1395  * NOTE: In the SSE4.1's nearest mode, if two values are equally close, the
1396  * result is the even value.  That is, rounding 2.5 will be 2.0, and not 3.0.
1397  */
1398 static INLINE LLVMValueRef
1399 lp_build_round_sse41(struct lp_build_context *bld,
1400                      LLVMValueRef a,
1401                      enum lp_build_round_sse41_mode mode)
1402 {
1403    LLVMBuilderRef builder = bld->gallivm->builder;
1404    const struct lp_type type = bld->type;
1405    LLVMTypeRef i32t = LLVMInt32TypeInContext(bld->gallivm->context);
1406    const char *intrinsic;
1407    LLVMValueRef res;
1408
1409    assert(type.floating);
1410
1411    assert(lp_check_value(type, a));
1412    assert(util_cpu_caps.has_sse4_1);
1413
1414    if (type.length == 1) {
1415       LLVMTypeRef vec_type;
1416       LLVMValueRef undef;
1417       LLVMValueRef args[3];
1418       LLVMValueRef index0 = LLVMConstInt(i32t, 0, 0);
1419
1420       switch(type.width) {
1421       case 32:
1422          intrinsic = "llvm.x86.sse41.round.ss";
1423          break;
1424       case 64:
1425          intrinsic = "llvm.x86.sse41.round.sd";
1426          break;
1427       default:
1428          assert(0);
1429          return bld->undef;
1430       }
1431
1432       vec_type = LLVMVectorType(bld->elem_type, 4);
1433
1434       undef = LLVMGetUndef(vec_type);
1435
1436       args[0] = undef;
1437       args[1] = LLVMBuildInsertElement(builder, undef, a, index0, "");
1438       args[2] = LLVMConstInt(i32t, mode, 0);
1439
1440       res = lp_build_intrinsic(builder, intrinsic,
1441                                vec_type, args, Elements(args));
1442
1443       res = LLVMBuildExtractElement(builder, res, index0, "");
1444    }
1445    else {
1446       if (type.width * type.length == 128) {
1447          switch(type.width) {
1448          case 32:
1449             intrinsic = "llvm.x86.sse41.round.ps";
1450             break;
1451          case 64:
1452             intrinsic = "llvm.x86.sse41.round.pd";
1453             break;
1454          default:
1455             assert(0);
1456             return bld->undef;
1457          }
1458       }
1459       else {
1460          assert(type.width * type.length == 256);
1461          assert(util_cpu_caps.has_avx);
1462
1463          switch(type.width) {
1464          case 32:
1465             intrinsic = "llvm.x86.avx.round.ps.256";
1466             break;
1467          case 64:
1468             intrinsic = "llvm.x86.avx.round.pd.256";
1469             break;
1470          default:
1471             assert(0);
1472             return bld->undef;
1473          }
1474       }
1475
1476       res = lp_build_intrinsic_binary(builder, intrinsic,
1477                                       bld->vec_type, a,
1478                                       LLVMConstInt(i32t, mode, 0));
1479    }
1480
1481    return res;
1482 }
1483
1484
1485 static INLINE LLVMValueRef
1486 lp_build_iround_nearest_sse2(struct lp_build_context *bld,
1487                              LLVMValueRef a)
1488 {
1489    LLVMBuilderRef builder = bld->gallivm->builder;
1490    const struct lp_type type = bld->type;
1491    LLVMTypeRef i32t = LLVMInt32TypeInContext(bld->gallivm->context);
1492    LLVMTypeRef ret_type = lp_build_int_vec_type(bld->gallivm, type);
1493    const char *intrinsic;
1494    LLVMValueRef res;
1495
1496    assert(type.floating);
1497    /* using the double precision conversions is a bit more complicated */
1498    assert(type.width == 32);
1499
1500    assert(lp_check_value(type, a));
1501    assert(util_cpu_caps.has_sse2);
1502
1503    /* This is relying on MXCSR rounding mode, which should always be nearest. */
1504    if (type.length == 1) {
1505       LLVMTypeRef vec_type;
1506       LLVMValueRef undef;
1507       LLVMValueRef arg;
1508       LLVMValueRef index0 = LLVMConstInt(i32t, 0, 0);
1509
1510       vec_type = LLVMVectorType(bld->elem_type, 4);
1511
1512       intrinsic = "llvm.x86.sse.cvtss2si";
1513
1514       undef = LLVMGetUndef(vec_type);
1515
1516       arg = LLVMBuildInsertElement(builder, undef, a, index0, "");
1517
1518       res = lp_build_intrinsic_unary(builder, intrinsic,
1519                                      ret_type, arg);
1520    }
1521    else {
1522       if (type.width* type.length == 128) {
1523          intrinsic = "llvm.x86.sse2.cvtps2dq";
1524       }
1525       else {
1526          assert(type.width*type.length == 256);
1527          assert(util_cpu_caps.has_avx);
1528
1529          intrinsic = "llvm.x86.avx.cvt.ps2dq.256";
1530       }
1531       res = lp_build_intrinsic_unary(builder, intrinsic,
1532                                      ret_type, a);
1533    }
1534
1535    return res;
1536 }
1537
1538
1539 /**
1540  * Return the integer part of a float (vector) value (== round toward zero).
1541  * The returned value is a float (vector).
1542  * Ex: trunc(-1.5) = -1.0
1543  */
1544 LLVMValueRef
1545 lp_build_trunc(struct lp_build_context *bld,
1546                LLVMValueRef a)
1547 {
1548    LLVMBuilderRef builder = bld->gallivm->builder;
1549    const struct lp_type type = bld->type;
1550
1551    assert(type.floating);
1552    assert(lp_check_value(type, a));
1553
1554    if (sse41_rounding_available(type)) {
1555       return lp_build_round_sse41(bld, a, LP_BUILD_ROUND_SSE41_TRUNCATE);
1556    }
1557    else {
1558       LLVMTypeRef vec_type = lp_build_vec_type(bld->gallivm, type);
1559       LLVMTypeRef int_vec_type = lp_build_int_vec_type(bld->gallivm, type);
1560       LLVMValueRef res;
1561       res = LLVMBuildFPToSI(builder, a, int_vec_type, "");
1562       res = LLVMBuildSIToFP(builder, res, vec_type, "");
1563       return res;
1564    }
1565 }
1566
1567
1568 /**
1569  * Return float (vector) rounded to nearest integer (vector).  The returned
1570  * value is a float (vector).
1571  * Ex: round(0.9) = 1.0
1572  * Ex: round(-1.5) = -2.0
1573  */
1574 LLVMValueRef
1575 lp_build_round(struct lp_build_context *bld,
1576                LLVMValueRef a)
1577 {
1578    LLVMBuilderRef builder = bld->gallivm->builder;
1579    const struct lp_type type = bld->type;
1580
1581    assert(type.floating);
1582    assert(lp_check_value(type, a));
1583
1584    if (sse41_rounding_available(type)) {
1585       return lp_build_round_sse41(bld, a, LP_BUILD_ROUND_SSE41_NEAREST);
1586    }
1587    else {
1588       LLVMTypeRef vec_type = lp_build_vec_type(bld->gallivm, type);
1589       LLVMValueRef res;
1590       res = lp_build_iround(bld, a);
1591       res = LLVMBuildSIToFP(builder, res, vec_type, "");
1592       return res;
1593    }
1594 }
1595
1596
1597 /**
1598  * Return floor of float (vector), result is a float (vector)
1599  * Ex: floor(1.1) = 1.0
1600  * Ex: floor(-1.1) = -2.0
1601  */
1602 LLVMValueRef
1603 lp_build_floor(struct lp_build_context *bld,
1604                LLVMValueRef a)
1605 {
1606    LLVMBuilderRef builder = bld->gallivm->builder;
1607    const struct lp_type type = bld->type;
1608
1609    assert(type.floating);
1610    assert(lp_check_value(type, a));
1611
1612    if (sse41_rounding_available(type)) {
1613       return lp_build_round_sse41(bld, a, LP_BUILD_ROUND_SSE41_FLOOR);
1614    }
1615    else {
1616       LLVMTypeRef vec_type = lp_build_vec_type(bld->gallivm, type);
1617       LLVMValueRef res;
1618       res = lp_build_ifloor(bld, a);
1619       res = LLVMBuildSIToFP(builder, res, vec_type, "");
1620       return res;
1621    }
1622 }
1623
1624
1625 /**
1626  * Return ceiling of float (vector), returning float (vector).
1627  * Ex: ceil( 1.1) = 2.0
1628  * Ex: ceil(-1.1) = -1.0
1629  */
1630 LLVMValueRef
1631 lp_build_ceil(struct lp_build_context *bld,
1632               LLVMValueRef a)
1633 {
1634    LLVMBuilderRef builder = bld->gallivm->builder;
1635    const struct lp_type type = bld->type;
1636
1637    assert(type.floating);
1638    assert(lp_check_value(type, a));
1639
1640    if (sse41_rounding_available(type)) {
1641       return lp_build_round_sse41(bld, a, LP_BUILD_ROUND_SSE41_CEIL);
1642    }
1643    else {
1644       LLVMTypeRef vec_type = lp_build_vec_type(bld->gallivm, type);
1645       LLVMValueRef res;
1646       res = lp_build_iceil(bld, a);
1647       res = LLVMBuildSIToFP(builder, res, vec_type, "");
1648       return res;
1649    }
1650 }
1651
1652
1653 /**
1654  * Return fractional part of 'a' computed as a - floor(a)
1655  * Typically used in texture coord arithmetic.
1656  */
1657 LLVMValueRef
1658 lp_build_fract(struct lp_build_context *bld,
1659                LLVMValueRef a)
1660 {
1661    assert(bld->type.floating);
1662    return lp_build_sub(bld, a, lp_build_floor(bld, a));
1663 }
1664
1665
1666 /**
1667  * Prevent returning a fractional part of 1.0 for very small negative values of
1668  * 'a' by clamping against 0.99999(9).
1669  */
1670 static inline LLVMValueRef
1671 clamp_fract(struct lp_build_context *bld, LLVMValueRef fract)
1672 {
1673    LLVMValueRef max;
1674
1675    /* this is the largest number smaller than 1.0 representable as float */
1676    max = lp_build_const_vec(bld->gallivm, bld->type,
1677                             1.0 - 1.0/(1LL << (lp_mantissa(bld->type) + 1)));
1678    return lp_build_min(bld, fract, max);
1679 }
1680
1681
1682 /**
1683  * Same as lp_build_fract, but guarantees that the result is always smaller
1684  * than one.
1685  */
1686 LLVMValueRef
1687 lp_build_fract_safe(struct lp_build_context *bld,
1688                     LLVMValueRef a)
1689 {
1690    return clamp_fract(bld, lp_build_fract(bld, a));
1691 }
1692
1693
1694 /**
1695  * Return the integer part of a float (vector) value (== round toward zero).
1696  * The returned value is an integer (vector).
1697  * Ex: itrunc(-1.5) = -1
1698  */
1699 LLVMValueRef
1700 lp_build_itrunc(struct lp_build_context *bld,
1701                 LLVMValueRef a)
1702 {
1703    LLVMBuilderRef builder = bld->gallivm->builder;
1704    const struct lp_type type = bld->type;
1705    LLVMTypeRef int_vec_type = lp_build_int_vec_type(bld->gallivm, type);
1706
1707    assert(type.floating);
1708    assert(lp_check_value(type, a));
1709
1710    return LLVMBuildFPToSI(builder, a, int_vec_type, "");
1711 }
1712
1713
1714 /**
1715  * Return float (vector) rounded to nearest integer (vector).  The returned
1716  * value is an integer (vector).
1717  * Ex: iround(0.9) = 1
1718  * Ex: iround(-1.5) = -2
1719  */
1720 LLVMValueRef
1721 lp_build_iround(struct lp_build_context *bld,
1722                 LLVMValueRef a)
1723 {
1724    LLVMBuilderRef builder = bld->gallivm->builder;
1725    const struct lp_type type = bld->type;
1726    LLVMTypeRef int_vec_type = bld->int_vec_type;
1727    LLVMValueRef res;
1728
1729    assert(type.floating);
1730
1731    assert(lp_check_value(type, a));
1732
1733    if ((util_cpu_caps.has_sse2 &&
1734        ((type.width == 32) && (type.length == 1 || type.length == 4))) ||
1735        (util_cpu_caps.has_avx && type.width == 32 && type.length == 8)) {
1736       return lp_build_iround_nearest_sse2(bld, a);
1737    }
1738    if (sse41_rounding_available(type)) {
1739       res = lp_build_round_sse41(bld, a, LP_BUILD_ROUND_SSE41_NEAREST);
1740    }
1741    else {
1742       LLVMValueRef half;
1743
1744       half = lp_build_const_vec(bld->gallivm, type, 0.5);
1745
1746       if (type.sign) {
1747          LLVMTypeRef vec_type = bld->vec_type;
1748          LLVMValueRef mask = lp_build_const_int_vec(bld->gallivm, type,
1749                                     (unsigned long long)1 << (type.width - 1));
1750          LLVMValueRef sign;
1751
1752          /* get sign bit */
1753          sign = LLVMBuildBitCast(builder, a, int_vec_type, "");
1754          sign = LLVMBuildAnd(builder, sign, mask, "");
1755
1756          /* sign * 0.5 */
1757          half = LLVMBuildBitCast(builder, half, int_vec_type, "");
1758          half = LLVMBuildOr(builder, sign, half, "");
1759          half = LLVMBuildBitCast(builder, half, vec_type, "");
1760       }
1761
1762       res = LLVMBuildFAdd(builder, a, half, "");
1763    }
1764
1765    res = LLVMBuildFPToSI(builder, res, int_vec_type, "");
1766
1767    return res;
1768 }
1769
1770
1771 /**
1772  * Return floor of float (vector), result is an int (vector)
1773  * Ex: ifloor(1.1) = 1.0
1774  * Ex: ifloor(-1.1) = -2.0
1775  */
1776 LLVMValueRef
1777 lp_build_ifloor(struct lp_build_context *bld,
1778                 LLVMValueRef a)
1779 {
1780    LLVMBuilderRef builder = bld->gallivm->builder;
1781    const struct lp_type type = bld->type;
1782    LLVMTypeRef int_vec_type = bld->int_vec_type;
1783    LLVMValueRef res;
1784
1785    assert(type.floating);
1786    assert(lp_check_value(type, a));
1787
1788    res = a;
1789    if (type.sign) {
1790       if (sse41_rounding_available(type)) {
1791          res = lp_build_round_sse41(bld, a, LP_BUILD_ROUND_SSE41_FLOOR);
1792       }
1793       else {
1794          /* Take the sign bit and add it to 1 constant */
1795          LLVMTypeRef vec_type = bld->vec_type;
1796          unsigned mantissa = lp_mantissa(type);
1797          LLVMValueRef mask = lp_build_const_int_vec(bld->gallivm, type,
1798                                   (unsigned long long)1 << (type.width - 1));
1799          LLVMValueRef sign;
1800          LLVMValueRef offset;
1801
1802          /* sign = a < 0 ? ~0 : 0 */
1803          sign = LLVMBuildBitCast(builder, a, int_vec_type, "");
1804          sign = LLVMBuildAnd(builder, sign, mask, "");
1805          sign = LLVMBuildAShr(builder, sign,
1806                               lp_build_const_int_vec(bld->gallivm, type,
1807                                                      type.width - 1),
1808                               "ifloor.sign");
1809
1810          /* offset = -0.99999(9)f */
1811          offset = lp_build_const_vec(bld->gallivm, type,
1812                                      -(double)(((unsigned long long)1 << mantissa) - 10)/((unsigned long long)1 << mantissa));
1813          offset = LLVMConstBitCast(offset, int_vec_type);
1814
1815          /* offset = a < 0 ? offset : 0.0f */
1816          offset = LLVMBuildAnd(builder, offset, sign, "");
1817          offset = LLVMBuildBitCast(builder, offset, vec_type, "ifloor.offset");
1818
1819          res = LLVMBuildFAdd(builder, res, offset, "ifloor.res");
1820       }
1821    }
1822
1823    /* round to nearest (toward zero) */
1824    res = LLVMBuildFPToSI(builder, res, int_vec_type, "ifloor.res");
1825
1826    return res;
1827 }
1828
1829
1830 /**
1831  * Return ceiling of float (vector), returning int (vector).
1832  * Ex: iceil( 1.1) = 2
1833  * Ex: iceil(-1.1) = -1
1834  */
1835 LLVMValueRef
1836 lp_build_iceil(struct lp_build_context *bld,
1837                LLVMValueRef a)
1838 {
1839    LLVMBuilderRef builder = bld->gallivm->builder;
1840    const struct lp_type type = bld->type;
1841    LLVMTypeRef int_vec_type = bld->int_vec_type;
1842    LLVMValueRef res;
1843
1844    assert(type.floating);
1845    assert(lp_check_value(type, a));
1846
1847    if (sse41_rounding_available(type)) {
1848       res = lp_build_round_sse41(bld, a, LP_BUILD_ROUND_SSE41_CEIL);
1849    }
1850    else {
1851       LLVMTypeRef vec_type = bld->vec_type;
1852       unsigned mantissa = lp_mantissa(type);
1853       LLVMValueRef offset;
1854
1855       /* offset = 0.99999(9)f */
1856       offset = lp_build_const_vec(bld->gallivm, type,
1857                                   (double)(((unsigned long long)1 << mantissa) - 10)/((unsigned long long)1 << mantissa));
1858
1859       if (type.sign) {
1860          LLVMValueRef mask = lp_build_const_int_vec(bld->gallivm, type,
1861                                 (unsigned long long)1 << (type.width - 1));
1862          LLVMValueRef sign;
1863
1864          /* sign = a < 0 ? 0 : ~0 */
1865          sign = LLVMBuildBitCast(builder, a, int_vec_type, "");
1866          sign = LLVMBuildAnd(builder, sign, mask, "");
1867          sign = LLVMBuildAShr(builder, sign,
1868                               lp_build_const_int_vec(bld->gallivm, type,
1869                                                      type.width - 1),
1870                               "iceil.sign");
1871          sign = LLVMBuildNot(builder, sign, "iceil.not");
1872
1873          /* offset = a < 0 ? 0.0 : offset */
1874          offset = LLVMConstBitCast(offset, int_vec_type);
1875          offset = LLVMBuildAnd(builder, offset, sign, "");
1876          offset = LLVMBuildBitCast(builder, offset, vec_type, "iceil.offset");
1877       }
1878
1879       res = LLVMBuildFAdd(builder, a, offset, "iceil.res");
1880    }
1881
1882    /* round to nearest (toward zero) */
1883    res = LLVMBuildFPToSI(builder, res, int_vec_type, "iceil.res");
1884
1885    return res;
1886 }
1887
1888
1889 /**
1890  * Combined ifloor() & fract().
1891  *
1892  * Preferred to calling the functions separately, as it will ensure that the
1893  * strategy (floor() vs ifloor()) that results in less redundant work is used.
1894  */
1895 void
1896 lp_build_ifloor_fract(struct lp_build_context *bld,
1897                       LLVMValueRef a,
1898                       LLVMValueRef *out_ipart,
1899                       LLVMValueRef *out_fpart)
1900 {
1901    LLVMBuilderRef builder = bld->gallivm->builder;
1902    const struct lp_type type = bld->type;
1903    LLVMValueRef ipart;
1904
1905    assert(type.floating);
1906    assert(lp_check_value(type, a));
1907
1908    if (sse41_rounding_available(type)) {
1909       /*
1910        * floor() is easier.
1911        */
1912
1913       ipart = lp_build_floor(bld, a);
1914       *out_fpart = LLVMBuildFSub(builder, a, ipart, "fpart");
1915       *out_ipart = LLVMBuildFPToSI(builder, ipart, bld->int_vec_type, "ipart");
1916    }
1917    else {
1918       /*
1919        * ifloor() is easier.
1920        */
1921
1922       *out_ipart = lp_build_ifloor(bld, a);
1923       ipart = LLVMBuildSIToFP(builder, *out_ipart, bld->vec_type, "ipart");
1924       *out_fpart = LLVMBuildFSub(builder, a, ipart, "fpart");
1925    }
1926 }
1927
1928
1929 /**
1930  * Same as lp_build_ifloor_fract, but guarantees that the fractional part is
1931  * always smaller than one.
1932  */
1933 void
1934 lp_build_ifloor_fract_safe(struct lp_build_context *bld,
1935                            LLVMValueRef a,
1936                            LLVMValueRef *out_ipart,
1937                            LLVMValueRef *out_fpart)
1938 {
1939    lp_build_ifloor_fract(bld, a, out_ipart, out_fpart);
1940    *out_fpart = clamp_fract(bld, *out_fpart);
1941 }
1942
1943
1944 LLVMValueRef
1945 lp_build_sqrt(struct lp_build_context *bld,
1946               LLVMValueRef a)
1947 {
1948    LLVMBuilderRef builder = bld->gallivm->builder;
1949    const struct lp_type type = bld->type;
1950    LLVMTypeRef vec_type = lp_build_vec_type(bld->gallivm, type);
1951    char intrinsic[32];
1952
1953    assert(lp_check_value(type, a));
1954
1955    /* TODO: optimize the constant case */
1956
1957    assert(type.floating);
1958    if (type.length == 1) {
1959       util_snprintf(intrinsic, sizeof intrinsic, "llvm.sqrt.f%u", type.width);
1960    }
1961    else {
1962       util_snprintf(intrinsic, sizeof intrinsic, "llvm.sqrt.v%uf%u", type.length, type.width);
1963    }
1964
1965    return lp_build_intrinsic_unary(builder, intrinsic, vec_type, a);
1966 }
1967
1968
1969 /**
1970  * Do one Newton-Raphson step to improve reciprocate precision:
1971  *
1972  *   x_{i+1} = x_i * (2 - a * x_i)
1973  *
1974  * XXX: Unfortunately this won't give IEEE-754 conformant results for 0 or
1975  * +/-Inf, giving NaN instead.  Certain applications rely on this behavior,
1976  * such as Google Earth, which does RCP(RSQRT(0.0) when drawing the Earth's
1977  * halo. It would be necessary to clamp the argument to prevent this.
1978  *
1979  * See also:
1980  * - http://en.wikipedia.org/wiki/Division_(digital)#Newton.E2.80.93Raphson_division
1981  * - http://softwarecommunity.intel.com/articles/eng/1818.htm
1982  */
1983 static INLINE LLVMValueRef
1984 lp_build_rcp_refine(struct lp_build_context *bld,
1985                     LLVMValueRef a,
1986                     LLVMValueRef rcp_a)
1987 {
1988    LLVMBuilderRef builder = bld->gallivm->builder;
1989    LLVMValueRef two = lp_build_const_vec(bld->gallivm, bld->type, 2.0);
1990    LLVMValueRef res;
1991
1992    res = LLVMBuildFMul(builder, a, rcp_a, "");
1993    res = LLVMBuildFSub(builder, two, res, "");
1994    res = LLVMBuildFMul(builder, rcp_a, res, "");
1995
1996    return res;
1997 }
1998
1999
2000 LLVMValueRef
2001 lp_build_rcp(struct lp_build_context *bld,
2002              LLVMValueRef a)
2003 {
2004    LLVMBuilderRef builder = bld->gallivm->builder;
2005    const struct lp_type type = bld->type;
2006
2007    assert(lp_check_value(type, a));
2008
2009    if(a == bld->zero)
2010       return bld->undef;
2011    if(a == bld->one)
2012       return bld->one;
2013    if(a == bld->undef)
2014       return bld->undef;
2015
2016    assert(type.floating);
2017
2018    if(LLVMIsConstant(a))
2019       return LLVMConstFDiv(bld->one, a);
2020
2021    /*
2022     * We don't use RCPPS because:
2023     * - it only has 10bits of precision
2024     * - it doesn't even get the reciprocate of 1.0 exactly
2025     * - doing Newton-Rapshon steps yields wrong (NaN) values for 0.0 or Inf
2026     * - for recent processors the benefit over DIVPS is marginal, a case
2027     *   dependent
2028     *
2029     * We could still use it on certain processors if benchmarks show that the
2030     * RCPPS plus necessary workarounds are still preferrable to DIVPS; or for
2031     * particular uses that require less workarounds.
2032     */
2033
2034    if (FALSE && ((util_cpu_caps.has_sse && type.width == 32 && type.length == 4) ||
2035          (util_cpu_caps.has_avx && type.width == 32 && type.length == 8))){
2036       const unsigned num_iterations = 0;
2037       LLVMValueRef res;
2038       unsigned i;
2039       const char *intrinsic = NULL;
2040
2041       if (type.length == 4) {
2042          intrinsic = "llvm.x86.sse.rcp.ps";
2043       }
2044       else {
2045          intrinsic = "llvm.x86.avx.rcp.ps.256";
2046       }
2047
2048       res = lp_build_intrinsic_unary(builder, intrinsic, bld->vec_type, a);
2049
2050       for (i = 0; i < num_iterations; ++i) {
2051          res = lp_build_rcp_refine(bld, a, res);
2052       }
2053
2054       return res;
2055    }
2056
2057    return LLVMBuildFDiv(builder, bld->one, a, "");
2058 }
2059
2060
2061 /**
2062  * Do one Newton-Raphson step to improve rsqrt precision:
2063  *
2064  *   x_{i+1} = 0.5 * x_i * (3.0 - a * x_i * x_i)
2065  *
2066  * See also Intel 64 and IA-32 Architectures Optimization Manual.
2067  */
2068 static INLINE LLVMValueRef
2069 lp_build_rsqrt_refine(struct lp_build_context *bld,
2070                       LLVMValueRef a,
2071                       LLVMValueRef rsqrt_a)
2072 {
2073    LLVMBuilderRef builder = bld->gallivm->builder;
2074    LLVMValueRef half = lp_build_const_vec(bld->gallivm, bld->type, 0.5);
2075    LLVMValueRef three = lp_build_const_vec(bld->gallivm, bld->type, 3.0);
2076    LLVMValueRef res;
2077
2078    res = LLVMBuildFMul(builder, rsqrt_a, rsqrt_a, "");
2079    res = LLVMBuildFMul(builder, a, res, "");
2080    res = LLVMBuildFSub(builder, three, res, "");
2081    res = LLVMBuildFMul(builder, rsqrt_a, res, "");
2082    res = LLVMBuildFMul(builder, half, res, "");
2083
2084    return res;
2085 }
2086
2087
2088 /**
2089  * Generate 1/sqrt(a).
2090  * Result is undefined for values < 0, infinity for +0.
2091  */
2092 LLVMValueRef
2093 lp_build_rsqrt(struct lp_build_context *bld,
2094                LLVMValueRef a)
2095 {
2096    LLVMBuilderRef builder = bld->gallivm->builder;
2097    const struct lp_type type = bld->type;
2098
2099    assert(lp_check_value(type, a));
2100
2101    assert(type.floating);
2102
2103    /*
2104     * This should be faster but all denormals will end up as infinity.
2105     */
2106    if (0 && ((util_cpu_caps.has_sse && type.width == 32 && type.length == 4) ||
2107         (util_cpu_caps.has_avx && type.width == 32 && type.length == 8))) {
2108       const unsigned num_iterations = 1;
2109       LLVMValueRef res;
2110       unsigned i;
2111       const char *intrinsic = NULL;
2112
2113       if (type.length == 4) {
2114          intrinsic = "llvm.x86.sse.rsqrt.ps";
2115       }
2116       else {
2117          intrinsic = "llvm.x86.avx.rsqrt.ps.256";
2118       }
2119       if (num_iterations) {
2120          /*
2121           * Newton-Raphson will result in NaN instead of infinity for zero,
2122           * and NaN instead of zero for infinity.
2123           * Also, need to ensure rsqrt(1.0) == 1.0.
2124           * All numbers smaller than FLT_MIN will result in +infinity
2125           * (rsqrtps treats all denormals as zero).
2126           */
2127          /*
2128           * Certain non-c99 compilers don't know INFINITY and might not support
2129           * hacks to evaluate it at compile time neither.
2130           */
2131          const unsigned posinf_int = 0x7F800000;
2132          LLVMValueRef cmp;
2133          LLVMValueRef flt_min = lp_build_const_vec(bld->gallivm, type, FLT_MIN);
2134          LLVMValueRef inf = lp_build_const_int_vec(bld->gallivm, type, posinf_int);
2135
2136          inf = LLVMBuildBitCast(builder, inf, lp_build_vec_type(bld->gallivm, type), "");
2137
2138          res = lp_build_intrinsic_unary(builder, intrinsic, bld->vec_type, a);
2139
2140          for (i = 0; i < num_iterations; ++i) {
2141             res = lp_build_rsqrt_refine(bld, a, res);
2142          }
2143          cmp = lp_build_compare(bld->gallivm, type, PIPE_FUNC_LESS, a, flt_min);
2144          res = lp_build_select(bld, cmp, inf, res);
2145          cmp = lp_build_compare(bld->gallivm, type, PIPE_FUNC_EQUAL, a, inf);
2146          res = lp_build_select(bld, cmp, bld->zero, res);
2147          cmp = lp_build_compare(bld->gallivm, type, PIPE_FUNC_EQUAL, a, bld->one);
2148          res = lp_build_select(bld, cmp, bld->one, res);
2149       }
2150       else {
2151          /* rsqrt(1.0) != 1.0 here */
2152          res = lp_build_intrinsic_unary(builder, intrinsic, bld->vec_type, a);
2153
2154       }
2155
2156       return res;
2157    }
2158
2159    return lp_build_rcp(bld, lp_build_sqrt(bld, a));
2160 }
2161
2162
2163 /**
2164  * Generate sin(a) using SSE2
2165  */
2166 LLVMValueRef
2167 lp_build_sin(struct lp_build_context *bld,
2168              LLVMValueRef a)
2169 {
2170    struct gallivm_state *gallivm = bld->gallivm;
2171    LLVMBuilderRef builder = gallivm->builder;
2172    struct lp_type int_type = lp_int_type(bld->type);
2173    LLVMBuilderRef b = builder;
2174
2175    /*
2176     *  take the absolute value,
2177     *  x = _mm_and_ps(x, *(v4sf*)_ps_inv_sign_mask);
2178     */
2179
2180    LLVMValueRef inv_sig_mask = lp_build_const_int_vec(gallivm, bld->type, ~0x80000000);
2181    LLVMValueRef a_v4si = LLVMBuildBitCast(b, a, bld->int_vec_type, "a_v4si");
2182
2183    LLVMValueRef absi = LLVMBuildAnd(b, a_v4si, inv_sig_mask, "absi");
2184    LLVMValueRef x_abs = LLVMBuildBitCast(b, absi, bld->vec_type, "x_abs");
2185
2186    /*
2187     * extract the sign bit (upper one)
2188     * sign_bit = _mm_and_ps(sign_bit, *(v4sf*)_ps_sign_mask);
2189     */
2190    LLVMValueRef sig_mask = lp_build_const_int_vec(gallivm, bld->type, 0x80000000);
2191    LLVMValueRef sign_bit_i = LLVMBuildAnd(b, a_v4si, sig_mask, "sign_bit_i");
2192
2193    /*
2194     * scale by 4/Pi
2195     * y = _mm_mul_ps(x, *(v4sf*)_ps_cephes_FOPI);
2196     */
2197
2198    LLVMValueRef FOPi = lp_build_const_vec(gallivm, bld->type, 1.27323954473516);
2199    LLVMValueRef scale_y = LLVMBuildFMul(b, x_abs, FOPi, "scale_y");
2200
2201    /*
2202     * store the integer part of y in mm0
2203     * emm2 = _mm_cvttps_epi32(y);
2204     */
2205
2206    LLVMValueRef emm2_i = LLVMBuildFPToSI(b, scale_y, bld->int_vec_type, "emm2_i");
2207
2208    /*
2209     * j=(j+1) & (~1) (see the cephes sources)
2210     * emm2 = _mm_add_epi32(emm2, *(v4si*)_pi32_1);
2211     */
2212
2213    LLVMValueRef all_one = lp_build_const_int_vec(gallivm, bld->type, 1);
2214    LLVMValueRef emm2_add =  LLVMBuildAdd(b, emm2_i, all_one, "emm2_add");
2215    /*
2216     * emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_inv1);
2217     */
2218    LLVMValueRef inv_one = lp_build_const_int_vec(gallivm, bld->type, ~1);
2219    LLVMValueRef emm2_and =  LLVMBuildAnd(b, emm2_add, inv_one, "emm2_and");
2220
2221    /*
2222     * y = _mm_cvtepi32_ps(emm2);
2223     */
2224    LLVMValueRef y_2 = LLVMBuildSIToFP(b, emm2_and, bld->vec_type, "y_2");
2225
2226    /* get the swap sign flag
2227     * emm0 = _mm_and_si128(emm2, *(v4si*)_pi32_4);
2228     */
2229    LLVMValueRef pi32_4 = lp_build_const_int_vec(gallivm, bld->type, 4);
2230    LLVMValueRef emm0_and =  LLVMBuildAnd(b, emm2_add, pi32_4, "emm0_and");
2231
2232    /*
2233     * emm2 = _mm_slli_epi32(emm0, 29);
2234     */
2235    LLVMValueRef const_29 = lp_build_const_int_vec(gallivm, bld->type, 29);
2236    LLVMValueRef swap_sign_bit = LLVMBuildShl(b, emm0_and, const_29, "swap_sign_bit");
2237
2238    /*
2239     * get the polynom selection mask
2240     * there is one polynom for 0 <= x <= Pi/4
2241     * and another one for Pi/4<x<=Pi/2
2242     * Both branches will be computed.
2243     *
2244     * emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_2);
2245     * emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
2246     */
2247
2248    LLVMValueRef pi32_2 = lp_build_const_int_vec(gallivm, bld->type, 2);
2249    LLVMValueRef emm2_3 =  LLVMBuildAnd(b, emm2_and, pi32_2, "emm2_3");
2250    LLVMValueRef poly_mask = lp_build_compare(gallivm,
2251                                              int_type, PIPE_FUNC_EQUAL,
2252                                              emm2_3, lp_build_const_int_vec(gallivm, bld->type, 0));
2253    /*
2254     *   sign_bit = _mm_xor_ps(sign_bit, swap_sign_bit);
2255     */
2256    LLVMValueRef sign_bit_1 =  LLVMBuildXor(b, sign_bit_i, swap_sign_bit, "sign_bit");
2257
2258    /*
2259     * _PS_CONST(minus_cephes_DP1, -0.78515625);
2260     * _PS_CONST(minus_cephes_DP2, -2.4187564849853515625e-4);
2261     * _PS_CONST(minus_cephes_DP3, -3.77489497744594108e-8);
2262     */
2263    LLVMValueRef DP1 = lp_build_const_vec(gallivm, bld->type, -0.78515625);
2264    LLVMValueRef DP2 = lp_build_const_vec(gallivm, bld->type, -2.4187564849853515625e-4);
2265    LLVMValueRef DP3 = lp_build_const_vec(gallivm, bld->type, -3.77489497744594108e-8);
2266
2267    /*
2268     * The magic pass: "Extended precision modular arithmetic"
2269     * x = ((x - y * DP1) - y * DP2) - y * DP3;
2270     * xmm1 = _mm_mul_ps(y, xmm1);
2271     * xmm2 = _mm_mul_ps(y, xmm2);
2272     * xmm3 = _mm_mul_ps(y, xmm3);
2273     */
2274    LLVMValueRef xmm1 = LLVMBuildFMul(b, y_2, DP1, "xmm1");
2275    LLVMValueRef xmm2 = LLVMBuildFMul(b, y_2, DP2, "xmm2");
2276    LLVMValueRef xmm3 = LLVMBuildFMul(b, y_2, DP3, "xmm3");
2277
2278    /*
2279     * x = _mm_add_ps(x, xmm1);
2280     * x = _mm_add_ps(x, xmm2);
2281     * x = _mm_add_ps(x, xmm3);
2282     */
2283
2284    LLVMValueRef x_1 = LLVMBuildFAdd(b, x_abs, xmm1, "x_1");
2285    LLVMValueRef x_2 = LLVMBuildFAdd(b, x_1, xmm2, "x_2");
2286    LLVMValueRef x_3 = LLVMBuildFAdd(b, x_2, xmm3, "x_3");
2287
2288    /*
2289     * Evaluate the first polynom  (0 <= x <= Pi/4)
2290     *
2291     * z = _mm_mul_ps(x,x);
2292     */
2293    LLVMValueRef z = LLVMBuildFMul(b, x_3, x_3, "z");
2294
2295    /*
2296     * _PS_CONST(coscof_p0,  2.443315711809948E-005);
2297     * _PS_CONST(coscof_p1, -1.388731625493765E-003);
2298     * _PS_CONST(coscof_p2,  4.166664568298827E-002);
2299     */
2300    LLVMValueRef coscof_p0 = lp_build_const_vec(gallivm, bld->type, 2.443315711809948E-005);
2301    LLVMValueRef coscof_p1 = lp_build_const_vec(gallivm, bld->type, -1.388731625493765E-003);
2302    LLVMValueRef coscof_p2 = lp_build_const_vec(gallivm, bld->type, 4.166664568298827E-002);
2303
2304    /*
2305     * y = *(v4sf*)_ps_coscof_p0;
2306     * y = _mm_mul_ps(y, z);
2307     */
2308    LLVMValueRef y_3 = LLVMBuildFMul(b, z, coscof_p0, "y_3");
2309    LLVMValueRef y_4 = LLVMBuildFAdd(b, y_3, coscof_p1, "y_4");
2310    LLVMValueRef y_5 = LLVMBuildFMul(b, y_4, z, "y_5");
2311    LLVMValueRef y_6 = LLVMBuildFAdd(b, y_5, coscof_p2, "y_6");
2312    LLVMValueRef y_7 = LLVMBuildFMul(b, y_6, z, "y_7");
2313    LLVMValueRef y_8 = LLVMBuildFMul(b, y_7, z, "y_8");
2314
2315
2316    /*
2317     * tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5);
2318     * y = _mm_sub_ps(y, tmp);
2319     * y = _mm_add_ps(y, *(v4sf*)_ps_1);
2320     */
2321    LLVMValueRef half = lp_build_const_vec(gallivm, bld->type, 0.5);
2322    LLVMValueRef tmp = LLVMBuildFMul(b, z, half, "tmp");
2323    LLVMValueRef y_9 = LLVMBuildFSub(b, y_8, tmp, "y_8");
2324    LLVMValueRef one = lp_build_const_vec(gallivm, bld->type, 1.0);
2325    LLVMValueRef y_10 = LLVMBuildFAdd(b, y_9, one, "y_9");
2326
2327    /*
2328     * _PS_CONST(sincof_p0, -1.9515295891E-4);
2329     * _PS_CONST(sincof_p1,  8.3321608736E-3);
2330     * _PS_CONST(sincof_p2, -1.6666654611E-1);
2331     */
2332    LLVMValueRef sincof_p0 = lp_build_const_vec(gallivm, bld->type, -1.9515295891E-4);
2333    LLVMValueRef sincof_p1 = lp_build_const_vec(gallivm, bld->type, 8.3321608736E-3);
2334    LLVMValueRef sincof_p2 = lp_build_const_vec(gallivm, bld->type, -1.6666654611E-1);
2335
2336    /*
2337     * Evaluate the second polynom  (Pi/4 <= x <= 0)
2338     *
2339     * y2 = *(v4sf*)_ps_sincof_p0;
2340     * y2 = _mm_mul_ps(y2, z);
2341     * y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p1);
2342     * y2 = _mm_mul_ps(y2, z);
2343     * y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p2);
2344     * y2 = _mm_mul_ps(y2, z);
2345     * y2 = _mm_mul_ps(y2, x);
2346     * y2 = _mm_add_ps(y2, x);
2347     */
2348
2349    LLVMValueRef y2_3 = LLVMBuildFMul(b, z, sincof_p0, "y2_3");
2350    LLVMValueRef y2_4 = LLVMBuildFAdd(b, y2_3, sincof_p1, "y2_4");
2351    LLVMValueRef y2_5 = LLVMBuildFMul(b, y2_4, z, "y2_5");
2352    LLVMValueRef y2_6 = LLVMBuildFAdd(b, y2_5, sincof_p2, "y2_6");
2353    LLVMValueRef y2_7 = LLVMBuildFMul(b, y2_6, z, "y2_7");
2354    LLVMValueRef y2_8 = LLVMBuildFMul(b, y2_7, x_3, "y2_8");
2355    LLVMValueRef y2_9 = LLVMBuildFAdd(b, y2_8, x_3, "y2_9");
2356
2357    /*
2358     * select the correct result from the two polynoms
2359     * xmm3 = poly_mask;
2360     * y2 = _mm_and_ps(xmm3, y2); //, xmm3);
2361     * y = _mm_andnot_ps(xmm3, y);
2362     * y = _mm_add_ps(y,y2);
2363     */
2364    LLVMValueRef y2_i = LLVMBuildBitCast(b, y2_9, bld->int_vec_type, "y2_i");
2365    LLVMValueRef y_i = LLVMBuildBitCast(b, y_10, bld->int_vec_type, "y_i");
2366    LLVMValueRef y2_and = LLVMBuildAnd(b, y2_i, poly_mask, "y2_and");
2367    LLVMValueRef inv = lp_build_const_int_vec(gallivm, bld->type, ~0);
2368    LLVMValueRef poly_mask_inv = LLVMBuildXor(b, poly_mask, inv, "poly_mask_inv");
2369    LLVMValueRef y_and = LLVMBuildAnd(b, y_i, poly_mask_inv, "y_and");
2370    LLVMValueRef y_combine = LLVMBuildAdd(b, y_and, y2_and, "y_combine");
2371
2372    /*
2373     * update the sign
2374     * y = _mm_xor_ps(y, sign_bit);
2375     */
2376    LLVMValueRef y_sign = LLVMBuildXor(b, y_combine, sign_bit_1, "y_sin");
2377    LLVMValueRef y_result = LLVMBuildBitCast(b, y_sign, bld->vec_type, "y_result");
2378    return y_result;
2379 }
2380
2381
2382 /**
2383  * Generate cos(a) using SSE2
2384  */
2385 LLVMValueRef
2386 lp_build_cos(struct lp_build_context *bld,
2387              LLVMValueRef a)
2388 {
2389    struct gallivm_state *gallivm = bld->gallivm;
2390    LLVMBuilderRef builder = gallivm->builder;
2391    struct lp_type int_type = lp_int_type(bld->type);
2392    LLVMBuilderRef b = builder;
2393
2394    /*
2395     *  take the absolute value,
2396     *  x = _mm_and_ps(x, *(v4sf*)_ps_inv_sign_mask);
2397     */
2398
2399    LLVMValueRef inv_sig_mask = lp_build_const_int_vec(gallivm, bld->type, ~0x80000000);
2400    LLVMValueRef a_v4si = LLVMBuildBitCast(b, a, bld->int_vec_type, "a_v4si");
2401
2402    LLVMValueRef absi = LLVMBuildAnd(b, a_v4si, inv_sig_mask, "absi");
2403    LLVMValueRef x_abs = LLVMBuildBitCast(b, absi, bld->vec_type, "x_abs");
2404
2405    /*
2406     * scale by 4/Pi
2407     * y = _mm_mul_ps(x, *(v4sf*)_ps_cephes_FOPI);
2408     */
2409
2410    LLVMValueRef FOPi = lp_build_const_vec(gallivm, bld->type, 1.27323954473516);
2411    LLVMValueRef scale_y = LLVMBuildFMul(b, x_abs, FOPi, "scale_y");
2412
2413    /*
2414     * store the integer part of y in mm0
2415     * emm2 = _mm_cvttps_epi32(y);
2416     */
2417
2418    LLVMValueRef emm2_i = LLVMBuildFPToSI(b, scale_y, bld->int_vec_type, "emm2_i");
2419
2420    /*
2421     * j=(j+1) & (~1) (see the cephes sources)
2422     * emm2 = _mm_add_epi32(emm2, *(v4si*)_pi32_1);
2423     */
2424
2425    LLVMValueRef all_one = lp_build_const_int_vec(gallivm, bld->type, 1);
2426    LLVMValueRef emm2_add =  LLVMBuildAdd(b, emm2_i, all_one, "emm2_add");
2427    /*
2428     * emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_inv1);
2429     */
2430    LLVMValueRef inv_one = lp_build_const_int_vec(gallivm, bld->type, ~1);
2431    LLVMValueRef emm2_and =  LLVMBuildAnd(b, emm2_add, inv_one, "emm2_and");
2432
2433    /*
2434     * y = _mm_cvtepi32_ps(emm2);
2435     */
2436    LLVMValueRef y_2 = LLVMBuildSIToFP(b, emm2_and, bld->vec_type, "y_2");
2437
2438
2439    /*
2440     * emm2 = _mm_sub_epi32(emm2, *(v4si*)_pi32_2);
2441     */
2442    LLVMValueRef const_2 = lp_build_const_int_vec(gallivm, bld->type, 2);
2443    LLVMValueRef emm2_2 = LLVMBuildSub(b, emm2_and, const_2, "emm2_2");
2444
2445
2446    /* get the swap sign flag
2447     * emm0 = _mm_andnot_si128(emm2, *(v4si*)_pi32_4);
2448     */
2449    LLVMValueRef inv = lp_build_const_int_vec(gallivm, bld->type, ~0);
2450    LLVMValueRef emm0_not = LLVMBuildXor(b, emm2_2, inv, "emm0_not");
2451    LLVMValueRef pi32_4 = lp_build_const_int_vec(gallivm, bld->type, 4);
2452    LLVMValueRef emm0_and =  LLVMBuildAnd(b, emm0_not, pi32_4, "emm0_and");
2453
2454    /*
2455     * emm2 = _mm_slli_epi32(emm0, 29);
2456     */
2457    LLVMValueRef const_29 = lp_build_const_int_vec(gallivm, bld->type, 29);
2458    LLVMValueRef sign_bit = LLVMBuildShl(b, emm0_and, const_29, "sign_bit");
2459
2460    /*
2461     * get the polynom selection mask
2462     * there is one polynom for 0 <= x <= Pi/4
2463     * and another one for Pi/4<x<=Pi/2
2464     * Both branches will be computed.
2465     *
2466     * emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_2);
2467     * emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
2468     */
2469
2470    LLVMValueRef pi32_2 = lp_build_const_int_vec(gallivm, bld->type, 2);
2471    LLVMValueRef emm2_3 =  LLVMBuildAnd(b, emm2_2, pi32_2, "emm2_3");
2472    LLVMValueRef poly_mask = lp_build_compare(gallivm,
2473                                              int_type, PIPE_FUNC_EQUAL,
2474                                              emm2_3, lp_build_const_int_vec(gallivm, bld->type, 0));
2475
2476    /*
2477     * _PS_CONST(minus_cephes_DP1, -0.78515625);
2478     * _PS_CONST(minus_cephes_DP2, -2.4187564849853515625e-4);
2479     * _PS_CONST(minus_cephes_DP3, -3.77489497744594108e-8);
2480     */
2481    LLVMValueRef DP1 = lp_build_const_vec(gallivm, bld->type, -0.78515625);
2482    LLVMValueRef DP2 = lp_build_const_vec(gallivm, bld->type, -2.4187564849853515625e-4);
2483    LLVMValueRef DP3 = lp_build_const_vec(gallivm, bld->type, -3.77489497744594108e-8);
2484
2485    /*
2486     * The magic pass: "Extended precision modular arithmetic"
2487     * x = ((x - y * DP1) - y * DP2) - y * DP3;
2488     * xmm1 = _mm_mul_ps(y, xmm1);
2489     * xmm2 = _mm_mul_ps(y, xmm2);
2490     * xmm3 = _mm_mul_ps(y, xmm3);
2491     */
2492    LLVMValueRef xmm1 = LLVMBuildFMul(b, y_2, DP1, "xmm1");
2493    LLVMValueRef xmm2 = LLVMBuildFMul(b, y_2, DP2, "xmm2");
2494    LLVMValueRef xmm3 = LLVMBuildFMul(b, y_2, DP3, "xmm3");
2495
2496    /*
2497     * x = _mm_add_ps(x, xmm1);
2498     * x = _mm_add_ps(x, xmm2);
2499     * x = _mm_add_ps(x, xmm3);
2500     */
2501
2502    LLVMValueRef x_1 = LLVMBuildFAdd(b, x_abs, xmm1, "x_1");
2503    LLVMValueRef x_2 = LLVMBuildFAdd(b, x_1, xmm2, "x_2");
2504    LLVMValueRef x_3 = LLVMBuildFAdd(b, x_2, xmm3, "x_3");
2505
2506    /*
2507     * Evaluate the first polynom  (0 <= x <= Pi/4)
2508     *
2509     * z = _mm_mul_ps(x,x);
2510     */
2511    LLVMValueRef z = LLVMBuildFMul(b, x_3, x_3, "z");
2512
2513    /*
2514     * _PS_CONST(coscof_p0,  2.443315711809948E-005);
2515     * _PS_CONST(coscof_p1, -1.388731625493765E-003);
2516     * _PS_CONST(coscof_p2,  4.166664568298827E-002);
2517     */
2518    LLVMValueRef coscof_p0 = lp_build_const_vec(gallivm, bld->type, 2.443315711809948E-005);
2519    LLVMValueRef coscof_p1 = lp_build_const_vec(gallivm, bld->type, -1.388731625493765E-003);
2520    LLVMValueRef coscof_p2 = lp_build_const_vec(gallivm, bld->type, 4.166664568298827E-002);
2521
2522    /*
2523     * y = *(v4sf*)_ps_coscof_p0;
2524     * y = _mm_mul_ps(y, z);
2525     */
2526    LLVMValueRef y_3 = LLVMBuildFMul(b, z, coscof_p0, "y_3");
2527    LLVMValueRef y_4 = LLVMBuildFAdd(b, y_3, coscof_p1, "y_4");
2528    LLVMValueRef y_5 = LLVMBuildFMul(b, y_4, z, "y_5");
2529    LLVMValueRef y_6 = LLVMBuildFAdd(b, y_5, coscof_p2, "y_6");
2530    LLVMValueRef y_7 = LLVMBuildFMul(b, y_6, z, "y_7");
2531    LLVMValueRef y_8 = LLVMBuildFMul(b, y_7, z, "y_8");
2532
2533
2534    /*
2535     * tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5);
2536     * y = _mm_sub_ps(y, tmp);
2537     * y = _mm_add_ps(y, *(v4sf*)_ps_1);
2538     */
2539    LLVMValueRef half = lp_build_const_vec(gallivm, bld->type, 0.5);
2540    LLVMValueRef tmp = LLVMBuildFMul(b, z, half, "tmp");
2541    LLVMValueRef y_9 = LLVMBuildFSub(b, y_8, tmp, "y_8");
2542    LLVMValueRef one = lp_build_const_vec(gallivm, bld->type, 1.0);
2543    LLVMValueRef y_10 = LLVMBuildFAdd(b, y_9, one, "y_9");
2544
2545    /*
2546     * _PS_CONST(sincof_p0, -1.9515295891E-4);
2547     * _PS_CONST(sincof_p1,  8.3321608736E-3);
2548     * _PS_CONST(sincof_p2, -1.6666654611E-1);
2549     */
2550    LLVMValueRef sincof_p0 = lp_build_const_vec(gallivm, bld->type, -1.9515295891E-4);
2551    LLVMValueRef sincof_p1 = lp_build_const_vec(gallivm, bld->type, 8.3321608736E-3);
2552    LLVMValueRef sincof_p2 = lp_build_const_vec(gallivm, bld->type, -1.6666654611E-1);
2553
2554    /*
2555     * Evaluate the second polynom  (Pi/4 <= x <= 0)
2556     *
2557     * y2 = *(v4sf*)_ps_sincof_p0;
2558     * y2 = _mm_mul_ps(y2, z);
2559     * y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p1);
2560     * y2 = _mm_mul_ps(y2, z);
2561     * y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p2);
2562     * y2 = _mm_mul_ps(y2, z);
2563     * y2 = _mm_mul_ps(y2, x);
2564     * y2 = _mm_add_ps(y2, x);
2565     */
2566
2567    LLVMValueRef y2_3 = LLVMBuildFMul(b, z, sincof_p0, "y2_3");
2568    LLVMValueRef y2_4 = LLVMBuildFAdd(b, y2_3, sincof_p1, "y2_4");
2569    LLVMValueRef y2_5 = LLVMBuildFMul(b, y2_4, z, "y2_5");
2570    LLVMValueRef y2_6 = LLVMBuildFAdd(b, y2_5, sincof_p2, "y2_6");
2571    LLVMValueRef y2_7 = LLVMBuildFMul(b, y2_6, z, "y2_7");
2572    LLVMValueRef y2_8 = LLVMBuildFMul(b, y2_7, x_3, "y2_8");
2573    LLVMValueRef y2_9 = LLVMBuildFAdd(b, y2_8, x_3, "y2_9");
2574
2575    /*
2576     * select the correct result from the two polynoms
2577     * xmm3 = poly_mask;
2578     * y2 = _mm_and_ps(xmm3, y2); //, xmm3);
2579     * y = _mm_andnot_ps(xmm3, y);
2580     * y = _mm_add_ps(y,y2);
2581     */
2582    LLVMValueRef y2_i = LLVMBuildBitCast(b, y2_9, bld->int_vec_type, "y2_i");
2583    LLVMValueRef y_i = LLVMBuildBitCast(b, y_10, bld->int_vec_type, "y_i");
2584    LLVMValueRef y2_and = LLVMBuildAnd(b, y2_i, poly_mask, "y2_and");
2585    LLVMValueRef poly_mask_inv = LLVMBuildXor(b, poly_mask, inv, "poly_mask_inv");
2586    LLVMValueRef y_and = LLVMBuildAnd(b, y_i, poly_mask_inv, "y_and");
2587    LLVMValueRef y_combine = LLVMBuildAdd(b, y_and, y2_and, "y_combine");
2588
2589    /*
2590     * update the sign
2591     * y = _mm_xor_ps(y, sign_bit);
2592     */
2593    LLVMValueRef y_sign = LLVMBuildXor(b, y_combine, sign_bit, "y_sin");
2594    LLVMValueRef y_result = LLVMBuildBitCast(b, y_sign, bld->vec_type, "y_result");
2595    return y_result;
2596 }
2597
2598
2599 /**
2600  * Generate pow(x, y)
2601  */
2602 LLVMValueRef
2603 lp_build_pow(struct lp_build_context *bld,
2604              LLVMValueRef x,
2605              LLVMValueRef y)
2606 {
2607    /* TODO: optimize the constant case */
2608    if (gallivm_debug & GALLIVM_DEBUG_PERF &&
2609        LLVMIsConstant(x) && LLVMIsConstant(y)) {
2610       debug_printf("%s: inefficient/imprecise constant arithmetic\n",
2611                    __FUNCTION__);
2612    }
2613
2614    return lp_build_exp2(bld, lp_build_mul(bld, lp_build_log2(bld, x), y));
2615 }
2616
2617
2618 /**
2619  * Generate exp(x)
2620  */
2621 LLVMValueRef
2622 lp_build_exp(struct lp_build_context *bld,
2623              LLVMValueRef x)
2624 {
2625    /* log2(e) = 1/log(2) */
2626    LLVMValueRef log2e = lp_build_const_vec(bld->gallivm, bld->type,
2627                                            1.4426950408889634);
2628
2629    assert(lp_check_value(bld->type, x));
2630
2631    return lp_build_exp2(bld, lp_build_mul(bld, log2e, x));
2632 }
2633
2634
2635 /**
2636  * Generate log(x)
2637  */
2638 LLVMValueRef
2639 lp_build_log(struct lp_build_context *bld,
2640              LLVMValueRef x)
2641 {
2642    /* log(2) */
2643    LLVMValueRef log2 = lp_build_const_vec(bld->gallivm, bld->type,
2644                                           0.69314718055994529);
2645
2646    assert(lp_check_value(bld->type, x));
2647
2648    return lp_build_mul(bld, log2, lp_build_log2(bld, x));
2649 }
2650
2651
2652 /**
2653  * Generate polynomial.
2654  * Ex:  coeffs[0] + x * coeffs[1] + x^2 * coeffs[2].
2655  */
2656 static LLVMValueRef
2657 lp_build_polynomial(struct lp_build_context *bld,
2658                     LLVMValueRef x,
2659                     const double *coeffs,
2660                     unsigned num_coeffs)
2661 {
2662    const struct lp_type type = bld->type;
2663    LLVMValueRef even = NULL, odd = NULL;
2664    LLVMValueRef x2;
2665    unsigned i;
2666
2667    assert(lp_check_value(bld->type, x));
2668
2669    /* TODO: optimize the constant case */
2670    if (gallivm_debug & GALLIVM_DEBUG_PERF &&
2671        LLVMIsConstant(x)) {
2672       debug_printf("%s: inefficient/imprecise constant arithmetic\n",
2673                    __FUNCTION__);
2674    }
2675
2676    /*
2677     * Calculate odd and even terms seperately to decrease data dependency
2678     * Ex:
2679     *     c[0] + x^2 * c[2] + x^4 * c[4] ...
2680     *     + x * (c[1] + x^2 * c[3] + x^4 * c[5]) ...
2681     */
2682    x2 = lp_build_mul(bld, x, x);
2683
2684    for (i = num_coeffs; i--; ) {
2685       LLVMValueRef coeff;
2686
2687       coeff = lp_build_const_vec(bld->gallivm, type, coeffs[i]);
2688
2689       if (i % 2 == 0) {
2690          if (even)
2691             even = lp_build_add(bld, coeff, lp_build_mul(bld, x2, even));
2692          else
2693             even = coeff;
2694       } else {
2695          if (odd)
2696             odd = lp_build_add(bld, coeff, lp_build_mul(bld, x2, odd));
2697          else
2698             odd = coeff;
2699       }
2700    }
2701
2702    if (odd)
2703       return lp_build_add(bld, lp_build_mul(bld, odd, x), even);
2704    else if (even)
2705       return even;
2706    else
2707       return bld->undef;
2708 }
2709
2710
2711 /**
2712  * Minimax polynomial fit of 2**x, in range [0, 1[
2713  */
2714 const double lp_build_exp2_polynomial[] = {
2715 #if EXP_POLY_DEGREE == 5
2716    0.999999925063526176901,
2717    0.693153073200168932794,
2718    0.240153617044375388211,
2719    0.0558263180532956664775,
2720    0.00898934009049466391101,
2721    0.00187757667519147912699
2722 #elif EXP_POLY_DEGREE == 4
2723    1.00000259337069434683,
2724    0.693003834469974940458,
2725    0.24144275689150793076,
2726    0.0520114606103070150235,
2727    0.0135341679161270268764
2728 #elif EXP_POLY_DEGREE == 3
2729    0.999925218562710312959,
2730    0.695833540494823811697,
2731    0.226067155427249155588,
2732    0.0780245226406372992967
2733 #elif EXP_POLY_DEGREE == 2
2734    1.00172476321474503578,
2735    0.657636275736077639316,
2736    0.33718943461968720704
2737 #else
2738 #error
2739 #endif
2740 };
2741
2742
2743 void
2744 lp_build_exp2_approx(struct lp_build_context *bld,
2745                      LLVMValueRef x,
2746                      LLVMValueRef *p_exp2_int_part,
2747                      LLVMValueRef *p_frac_part,
2748                      LLVMValueRef *p_exp2)
2749 {
2750    LLVMBuilderRef builder = bld->gallivm->builder;
2751    const struct lp_type type = bld->type;
2752    LLVMTypeRef vec_type = lp_build_vec_type(bld->gallivm, type);
2753    LLVMValueRef ipart = NULL;
2754    LLVMValueRef fpart = NULL;
2755    LLVMValueRef expipart = NULL;
2756    LLVMValueRef expfpart = NULL;
2757    LLVMValueRef res = NULL;
2758
2759    assert(lp_check_value(bld->type, x));
2760
2761    if(p_exp2_int_part || p_frac_part || p_exp2) {
2762       /* TODO: optimize the constant case */
2763       if (gallivm_debug & GALLIVM_DEBUG_PERF &&
2764           LLVMIsConstant(x)) {
2765          debug_printf("%s: inefficient/imprecise constant arithmetic\n",
2766                       __FUNCTION__);
2767       }
2768
2769       assert(type.floating && type.width == 32);
2770
2771       x = lp_build_min(bld, x, lp_build_const_vec(bld->gallivm, type,  129.0));
2772       x = lp_build_max(bld, x, lp_build_const_vec(bld->gallivm, type, -126.99999));
2773
2774       /* ipart = floor(x) */
2775       /* fpart = x - ipart */
2776       lp_build_ifloor_fract(bld, x, &ipart, &fpart);
2777    }
2778
2779    if(p_exp2_int_part || p_exp2) {
2780       /* expipart = (float) (1 << ipart) */
2781       expipart = LLVMBuildAdd(builder, ipart,
2782                               lp_build_const_int_vec(bld->gallivm, type, 127), "");
2783       expipart = LLVMBuildShl(builder, expipart,
2784                               lp_build_const_int_vec(bld->gallivm, type, 23), "");
2785       expipart = LLVMBuildBitCast(builder, expipart, vec_type, "");
2786    }
2787
2788    if(p_exp2) {
2789       expfpart = lp_build_polynomial(bld, fpart, lp_build_exp2_polynomial,
2790                                      Elements(lp_build_exp2_polynomial));
2791
2792       res = LLVMBuildFMul(builder, expipart, expfpart, "");
2793    }
2794
2795    if(p_exp2_int_part)
2796       *p_exp2_int_part = expipart;
2797
2798    if(p_frac_part)
2799       *p_frac_part = fpart;
2800
2801    if(p_exp2)
2802       *p_exp2 = res;
2803 }
2804
2805
2806 LLVMValueRef
2807 lp_build_exp2(struct lp_build_context *bld,
2808               LLVMValueRef x)
2809 {
2810    LLVMValueRef res;
2811    lp_build_exp2_approx(bld, x, NULL, NULL, &res);
2812    return res;
2813 }
2814
2815
2816 /**
2817  * Extract the exponent of a IEEE-754 floating point value.
2818  *
2819  * Optionally apply an integer bias.
2820  *
2821  * Result is an integer value with
2822  *
2823  *   ifloor(log2(x)) + bias
2824  */
2825 LLVMValueRef
2826 lp_build_extract_exponent(struct lp_build_context *bld,
2827                           LLVMValueRef x,
2828                           int bias)
2829 {
2830    LLVMBuilderRef builder = bld->gallivm->builder;
2831    const struct lp_type type = bld->type;
2832    unsigned mantissa = lp_mantissa(type);
2833    LLVMValueRef res;
2834
2835    assert(type.floating);
2836
2837    assert(lp_check_value(bld->type, x));
2838
2839    x = LLVMBuildBitCast(builder, x, bld->int_vec_type, "");
2840
2841    res = LLVMBuildLShr(builder, x,
2842                        lp_build_const_int_vec(bld->gallivm, type, mantissa), "");
2843    res = LLVMBuildAnd(builder, res,
2844                       lp_build_const_int_vec(bld->gallivm, type, 255), "");
2845    res = LLVMBuildSub(builder, res,
2846                       lp_build_const_int_vec(bld->gallivm, type, 127 - bias), "");
2847
2848    return res;
2849 }
2850
2851
2852 /**
2853  * Extract the mantissa of the a floating.
2854  *
2855  * Result is a floating point value with
2856  *
2857  *   x / floor(log2(x))
2858  */
2859 LLVMValueRef
2860 lp_build_extract_mantissa(struct lp_build_context *bld,
2861                           LLVMValueRef x)
2862 {
2863    LLVMBuilderRef builder = bld->gallivm->builder;
2864    const struct lp_type type = bld->type;
2865    unsigned mantissa = lp_mantissa(type);
2866    LLVMValueRef mantmask = lp_build_const_int_vec(bld->gallivm, type,
2867                                                   (1ULL << mantissa) - 1);
2868    LLVMValueRef one = LLVMConstBitCast(bld->one, bld->int_vec_type);
2869    LLVMValueRef res;
2870
2871    assert(lp_check_value(bld->type, x));
2872
2873    assert(type.floating);
2874
2875    x = LLVMBuildBitCast(builder, x, bld->int_vec_type, "");
2876
2877    /* res = x / 2**ipart */
2878    res = LLVMBuildAnd(builder, x, mantmask, "");
2879    res = LLVMBuildOr(builder, res, one, "");
2880    res = LLVMBuildBitCast(builder, res, bld->vec_type, "");
2881
2882    return res;
2883 }
2884
2885
2886
2887 /**
2888  * Minimax polynomial fit of log2((1.0 + sqrt(x))/(1.0 - sqrt(x)))/sqrt(x) ,for x in range of [0, 1/9[
2889  * These coefficients can be generate with
2890  * http://www.boost.org/doc/libs/1_36_0/libs/math/doc/sf_and_dist/html/math_toolkit/toolkit/internals2/minimax.html
2891  */
2892 const double lp_build_log2_polynomial[] = {
2893 #if LOG_POLY_DEGREE == 5
2894    2.88539008148777786488L,
2895    0.961796878841293367824L,
2896    0.577058946784739859012L,
2897    0.412914355135828735411L,
2898    0.308591899232910175289L,
2899    0.352376952300281371868L,
2900 #elif LOG_POLY_DEGREE == 4
2901    2.88539009343309178325L,
2902    0.961791550404184197881L,
2903    0.577440339438736392009L,
2904    0.403343858251329912514L,
2905    0.406718052498846252698L,
2906 #elif LOG_POLY_DEGREE == 3
2907    2.88538959748872753838L,
2908    0.961932915889597772928L,
2909    0.571118517972136195241L,
2910    0.493997535084709500285L,
2911 #else
2912 #error
2913 #endif
2914 };
2915
2916 /**
2917  * See http://www.devmaster.net/forums/showthread.php?p=43580
2918  * http://en.wikipedia.org/wiki/Logarithm#Calculation
2919  * http://www.nezumi.demon.co.uk/consult/logx.htm
2920  */
2921 void
2922 lp_build_log2_approx(struct lp_build_context *bld,
2923                      LLVMValueRef x,
2924                      LLVMValueRef *p_exp,
2925                      LLVMValueRef *p_floor_log2,
2926                      LLVMValueRef *p_log2)
2927 {
2928    LLVMBuilderRef builder = bld->gallivm->builder;
2929    const struct lp_type type = bld->type;
2930    LLVMTypeRef vec_type = lp_build_vec_type(bld->gallivm, type);
2931    LLVMTypeRef int_vec_type = lp_build_int_vec_type(bld->gallivm, type);
2932
2933    LLVMValueRef expmask = lp_build_const_int_vec(bld->gallivm, type, 0x7f800000);
2934    LLVMValueRef mantmask = lp_build_const_int_vec(bld->gallivm, type, 0x007fffff);
2935    LLVMValueRef one = LLVMConstBitCast(bld->one, int_vec_type);
2936
2937    LLVMValueRef i = NULL;
2938    LLVMValueRef y = NULL;
2939    LLVMValueRef z = NULL;
2940    LLVMValueRef exp = NULL;
2941    LLVMValueRef mant = NULL;
2942    LLVMValueRef logexp = NULL;
2943    LLVMValueRef logmant = NULL;
2944    LLVMValueRef res = NULL;
2945
2946    assert(lp_check_value(bld->type, x));
2947
2948    if(p_exp || p_floor_log2 || p_log2) {
2949       /* TODO: optimize the constant case */
2950       if (gallivm_debug & GALLIVM_DEBUG_PERF &&
2951           LLVMIsConstant(x)) {
2952          debug_printf("%s: inefficient/imprecise constant arithmetic\n",
2953                       __FUNCTION__);
2954       }
2955
2956       assert(type.floating && type.width == 32);
2957
2958       /*
2959        * We don't explicitly handle denormalized numbers. They will yield a
2960        * result in the neighbourhood of -127, which appears to be adequate
2961        * enough.
2962        */
2963
2964       i = LLVMBuildBitCast(builder, x, int_vec_type, "");
2965
2966       /* exp = (float) exponent(x) */
2967       exp = LLVMBuildAnd(builder, i, expmask, "");
2968    }
2969
2970    if(p_floor_log2 || p_log2) {
2971       logexp = LLVMBuildLShr(builder, exp, lp_build_const_int_vec(bld->gallivm, type, 23), "");
2972       logexp = LLVMBuildSub(builder, logexp, lp_build_const_int_vec(bld->gallivm, type, 127), "");
2973       logexp = LLVMBuildSIToFP(builder, logexp, vec_type, "");
2974    }
2975
2976    if(p_log2) {
2977       /* mant = 1 + (float) mantissa(x) */
2978       mant = LLVMBuildAnd(builder, i, mantmask, "");
2979       mant = LLVMBuildOr(builder, mant, one, "");
2980       mant = LLVMBuildBitCast(builder, mant, vec_type, "");
2981
2982       /* y = (mant - 1) / (mant + 1) */
2983       y = lp_build_div(bld,
2984          lp_build_sub(bld, mant, bld->one),
2985          lp_build_add(bld, mant, bld->one)
2986       );
2987
2988       /* z = y^2 */
2989       z = lp_build_mul(bld, y, y);
2990
2991       /* compute P(z) */
2992       logmant = lp_build_polynomial(bld, z, lp_build_log2_polynomial,
2993                                     Elements(lp_build_log2_polynomial));
2994
2995       /* logmant = y * P(z) */
2996       logmant = lp_build_mul(bld, y, logmant);
2997
2998       res = lp_build_add(bld, logmant, logexp);
2999    }
3000
3001    if(p_exp) {
3002       exp = LLVMBuildBitCast(builder, exp, vec_type, "");
3003       *p_exp = exp;
3004    }
3005
3006    if(p_floor_log2)
3007       *p_floor_log2 = logexp;
3008
3009    if(p_log2)
3010       *p_log2 = res;
3011 }
3012
3013
3014 LLVMValueRef
3015 lp_build_log2(struct lp_build_context *bld,
3016               LLVMValueRef x)
3017 {
3018    LLVMValueRef res;
3019    lp_build_log2_approx(bld, x, NULL, NULL, &res);
3020    return res;
3021 }
3022
3023
3024 /**
3025  * Faster (and less accurate) log2.
3026  *
3027  *    log2(x) = floor(log2(x)) - 1 + x / 2**floor(log2(x))
3028  *
3029  * Piece-wise linear approximation, with exact results when x is a
3030  * power of two.
3031  *
3032  * See http://www.flipcode.com/archives/Fast_log_Function.shtml
3033  */
3034 LLVMValueRef
3035 lp_build_fast_log2(struct lp_build_context *bld,
3036                    LLVMValueRef x)
3037 {
3038    LLVMBuilderRef builder = bld->gallivm->builder;
3039    LLVMValueRef ipart;
3040    LLVMValueRef fpart;
3041
3042    assert(lp_check_value(bld->type, x));
3043
3044    assert(bld->type.floating);
3045
3046    /* ipart = floor(log2(x)) - 1 */
3047    ipart = lp_build_extract_exponent(bld, x, -1);
3048    ipart = LLVMBuildSIToFP(builder, ipart, bld->vec_type, "");
3049
3050    /* fpart = x / 2**ipart */
3051    fpart = lp_build_extract_mantissa(bld, x);
3052
3053    /* ipart + fpart */
3054    return LLVMBuildFAdd(builder, ipart, fpart, "");
3055 }
3056
3057
3058 /**
3059  * Fast implementation of iround(log2(x)).
3060  *
3061  * Not an approximation -- it should give accurate results all the time.
3062  */
3063 LLVMValueRef
3064 lp_build_ilog2(struct lp_build_context *bld,
3065                LLVMValueRef x)
3066 {
3067    LLVMBuilderRef builder = bld->gallivm->builder;
3068    LLVMValueRef sqrt2 = lp_build_const_vec(bld->gallivm, bld->type, M_SQRT2);
3069    LLVMValueRef ipart;
3070
3071    assert(bld->type.floating);
3072
3073    assert(lp_check_value(bld->type, x));
3074
3075    /* x * 2^(0.5)   i.e., add 0.5 to the log2(x) */
3076    x = LLVMBuildFMul(builder, x, sqrt2, "");
3077
3078    /* ipart = floor(log2(x) + 0.5)  */
3079    ipart = lp_build_extract_exponent(bld, x, 0);
3080
3081    return ipart;
3082 }
3083
3084 LLVMValueRef
3085 lp_build_mod(struct lp_build_context *bld,
3086              LLVMValueRef x,
3087              LLVMValueRef y)
3088 {
3089    LLVMBuilderRef builder = bld->gallivm->builder;
3090    LLVMValueRef res;
3091    const struct lp_type type = bld->type;
3092
3093    assert(lp_check_value(type, x));
3094    assert(lp_check_value(type, y));
3095
3096    if (type.floating)
3097       res = LLVMBuildFRem(builder, x, y, "");
3098    else if (type.sign)
3099       res = LLVMBuildSRem(builder, x, y, "");
3100    else
3101       res = LLVMBuildURem(builder, x, y, "");
3102    return res;
3103 }