sysdeps/ia64/fpu/s_log1p.S

   1 .file "log1p.s"
   2
   3 // Copyright (c) 2000, 2001, Intel Corporation
   4 // All rights reserved.
   5 //
   6 // Contributed 2/2/2000 by John Harrison, Ted Kubaska, Bob Norin, Shane Story,
   7 // and Ping Tak Peter Tang of the Computational Software Lab, Intel Corporation.
   8 //
   9 // WARRANTY DISCLAIMER
  10 //
  11 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  12 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  13 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  14 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
  15 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  16 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  17 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  18 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
  19 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
  20 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  21 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  22 //
  23 // Intel Corporation is the author of this code, and requests that all
  24 // problem reports or change requests be submitted to it directly at
  25 // http://developer.intel.com/opensource.
  26 //
  27 // History
  28 //==============================================================
  29 // 2/02/00  Initial version
  30 // 4/04/00  Unwind support added
  31 // 8/15/00  Bundle added after call to __libm_error_support to properly
  32 //          set [the previously overwritten] GR_Parameter_RESULT.
  33 //
  34 // *********************************************************************
  35 //
  36 // Function:   log1p(x) = ln(x+1), for double precision x values
  37 //
  38 // *********************************************************************
  39 //
  40 // Accuracy:   Very accurate for double precision values
  41 //
  42 // *********************************************************************
  43 //
  44 // Resources Used:
  45 //
  46 //    Floating-Point Registers: f8 (Input and Return Value)
  47 //                              f9,f33-f55,f99
  48 //
  49 //    General Purpose Registers:
  50 //      r32-r53
  51 //      r54-r57 (Used to pass arguments to error handling routine)
  52 //
  53 //    Predicate Registers:      p6-p15
  54 //
  55 // *********************************************************************
  56 //
  57 // IEEE Special Conditions:
  58 //
  59 //    Denormal  fault raised on denormal inputs
  60 //    Overflow exceptions cannot occur
  61 //    Underflow exceptions raised when appropriate for log1p
  62 //    (Error Handling Routine called for underflow)
  63 //    Inexact raised when appropriate by algorithm
  64 //
  65 //    log1p(inf) = inf
  66 //    log1p(-inf) = QNaN
  67 //    log1p(+/-0) = +/-0
  68 //    log1p(-1) =  -inf
  69 //    log1p(SNaN) = QNaN
  70 //    log1p(QNaN) = QNaN
  71 //    log1p(EM_special Values) = QNaN
  72 //
  73 // *********************************************************************
  74 //
  75 // Computation is based on the following kernel.
  76 //
  77 // ker_log_64( in_FR    :  X,
  78 //          in_FR    :  E,
  79 //          in_FR    :  Em1,
  80 //          in_GR    :  Expo_Range,
  81 //          out_FR   :  Y_hi,
  82 //          out_FR   :  Y_lo,
  83 //          out_FR   :  Scale,
  84 //          out_PR   :  Safe  )
  85 //
  86 // Overview
  87 //
  88 // The method consists of three cases.
  89 //
  90 // If   |X+Em1| < 2^(-80)       use case log1p_small;
  91 // elseif       |X+Em1| < 2^(-7)        use case log_near1;
  92 // else                         use case log_regular;
  93 //
  94 // Case log1p_small:
  95 //
  96 // log( 1 + (X+Em1) ) can be approximated by (X+Em1).
  97 //
  98 // Case log_near1:
  99 //
 100 //   log( 1 + (X+Em1) ) can be approximated by a simple polynomial
 101 //   in W = X+Em1. This polynomial resembles the truncated Taylor
 102 //   series W - W^/2 + W^3/3 - ...
 103 //
 104 // Case log_regular:
 105 //
 106 //   Here we use a table lookup method. The basic idea is that in
 107 //   order to compute log(Arg) for an argument Arg in [1,2), we
 108 //   construct a value G such that G*Arg is close to 1 and that
 109 //   log(1/G) is obtainable easily from a table of values calculated
 110 //   beforehand. Thus
 111 //
 112 //      log(Arg) = log(1/G) + log(G*Arg)
 113 //               = log(1/G) + log(1 + (G*Arg - 1))
 114 //
 115 //   Because |G*Arg - 1| is small, the second term on the right hand
 116 //   side can be approximated by a short polynomial. We elaborate
 117 //   this method in four steps.
 118 //
 119 //   Step 0: Initialization
 120 //
 121 //   We need to calculate log( E + X ). Obtain N, S_hi, S_lo such that
 122 //
 123 //      E + X = 2^N * ( S_hi + S_lo )   exactly
 124 //
 125 //   where S_hi in [1,2) and S_lo is a correction to S_hi in the sense
 126 //   that |S_lo| <= ulp(S_hi).
 127 //
 128 //   Step 1: Argument Reduction
 129 //
 130 //   Based on S_hi, obtain G_1, G_2, G_3 from a table and calculate
 131 //
 132 //      G := G_1 * G_2 * G_3
 133 //      r := (G * S_hi - 1)  + G * S_lo
 134 //
 135 //   These G_j's have the property that the product is exactly
 136 //   representable and that |r| < 2^(-12) as a result.
 137 //
 138 //   Step 2: Approximation
 139 //
 140 //
 141 //   log(1 + r) is approximated by a short polynomial poly(r).
 142 //
 143 //   Step 3: Reconstruction
 144 //
 145 //
 146 //   Finally, log( E + X ) is given by
 147 //
 148 //   log( E + X )   =   log( 2^N * (S_hi + S_lo) )
 149 //                 ~=~  N*log(2) + log(1/G) + log(1 + r)
 150 //                 ~=~  N*log(2) + log(1/G) + poly(r).
 151 //
 152 // **** Algorithm ****
 153 //
 154 // Case log1p_small:
 155 //
 156 // Although log(1 + (X+Em1)) is basically X+Em1, we would like to
 157 // preserve the inexactness nature as well as consistent behavior
 158 // under different rounding modes. Note that this case can only be
 159 // taken if E is set to be 1.0. In this case, Em1 is zero, and that
 160 // X can be very tiny and thus the final result can possibly underflow.
 161 // Thus, we compare X against a threshold that is dependent on the
 162 // input Expo_Range. If |X| is smaller than this threshold, we set
 163 // SAFE to be FALSE.
 164 //
 165 // The result is returned as Y_hi, Y_lo, and in the case of SAFE
 166 // is FALSE, an additional value Scale is also returned.
 167 //
 168 //      W    := X + Em1
 169 //      Threshold := Threshold_Table( Expo_Range )
 170 //      Tiny      := Tiny_Table( Expo_Range )
 171 //
 172 //      If ( |W| > Threshold ) then
 173 //         Y_hi  := W
 174 //         Y_lo  := -W*W
 175 //      Else
 176 //         Y_hi  := W
 177 //         Y_lo  := -Tiny
 178 //         Scale := 2^(-100)
 179 //         Safe  := FALSE
 180 //      EndIf
 181 //
 182 //
 183 // One may think that Y_lo should be -W*W/2; however, it does not matter
 184 // as Y_lo will be rounded off completely except for the correct effect in
 185 // directed rounding. Clearly -W*W is simplier to compute. Moreover,
 186 // because of the difference in exponent value, Y_hi + Y_lo or
 187 // Y_hi + Scale*Y_lo is always inexact.
 188 //
 189 // Case log_near1:
 190 //
 191 // Here we compute a simple polynomial. To exploit parallelism, we split
 192 // the polynomial into two portions.
 193 //
 194 //      W := X + Em1
 195 //      Wsq := W * W
 196 //      W4  := Wsq*Wsq
 197 //      W6  := W4*Wsq
 198 //      Y_hi := W + Wsq*(P_1 + W*(P_2 + W*(P_3 + W*P_4))
 199 //      Y_lo := W6*(P_5 + W*(P_6 + W*(P_7 + W*P_8)))
 200 //      set lsb(Y_lo) to be 1
 201 //
 202 // Case log_regular:
 203 //
 204 // We present the algorithm in four steps.
 205 //
 206 //   Step 0. Initialization
 207 //   ----------------------
 208 //
 209 //   Z := X + E
 210 //   N := unbaised exponent of Z
 211 //   S_hi := 2^(-N) * Z
 212 //   S_lo := 2^(-N) * { (max(X,E)-Z) + min(X,E) }
 213 //
 214 //   Note that S_lo is always 0 for the case E = 0.
 215 //
 216 //   Step 1. Argument Reduction
 217 //   --------------------------
 218 //
 219 //   Let
 220 //
 221 //      Z = 2^N * S_hi = 2^N * 1.d_1 d_2 d_3 ... d_63
 222 //
 223 //   We obtain G_1, G_2, G_3 by the following steps.
 224 //
 225 //
 226 //      Define          X_0 := 1.d_1 d_2 ... d_14. This is extracted
 227 //                      from S_hi.
 228 //
 229 //      Define          A_1 := 1.d_1 d_2 d_3 d_4. This is X_0 truncated
 230 //                      to lsb = 2^(-4).
 231 //
 232 //      Define          index_1 := [ d_1 d_2 d_3 d_4 ].
 233 //
 234 //      Fetch           Z_1 := (1/A_1) rounded UP in fixed point with
 235 //      fixed point     lsb = 2^(-15).
 236 //                      Z_1 looks like z_0.z_1 z_2 ... z_15
 237 //                      Note that the fetching is done using index_1.
 238 //                      A_1 is actually not needed in the implementation
 239 //                      and is used here only to explain how is the value
 240 //                      Z_1 defined.
 241 //
 242 //      Fetch           G_1 := (1/A_1) truncated to 21 sig. bits.
 243 //      floating pt.    Again, fetching is done using index_1. A_1
 244 //                      explains how G_1 is defined.
 245 //
 246 //      Calculate       X_1 := X_0 * Z_1 truncated to lsb = 2^(-14)
 247 //                           = 1.0 0 0 0 d_5 ... d_14
 248 //                      This is accomplised by integer multiplication.
 249 //                      It is proved that X_1 indeed always begin
 250 //                      with 1.0000 in fixed point.
 251 //
 252 //
 253 //      Define          A_2 := 1.0 0 0 0 d_5 d_6 d_7 d_8. This is X_1
 254 //                      truncated to lsb = 2^(-8). Similar to A_1,
 255 //                      A_2 is not needed in actual implementation. It
 256 //                      helps explain how some of the values are defined.
 257 //
 258 //      Define          index_2 := [ d_5 d_6 d_7 d_8 ].
 259 //
 260 //      Fetch           Z_2 := (1/A_2) rounded UP in fixed point with
 261 //      fixed point     lsb = 2^(-15). Fetch done using index_2.
 262 //                      Z_2 looks like z_0.z_1 z_2 ... z_15
 263 //
 264 //      Fetch           G_2 := (1/A_2) truncated to 21 sig. bits.
 265 //      floating pt.
 266 //
 267 //      Calculate       X_2 := X_1 * Z_2 truncated to lsb = 2^(-14)
 268 //                           = 1.0 0 0 0 0 0 0 0 d_9 d_10 ... d_14
 269 //                      This is accomplised by integer multiplication.
 270 //                      It is proved that X_2 indeed always begin
 271 //                      with 1.00000000 in fixed point.
 272 //
 273 //
 274 //      Define          A_3 := 1.0 0 0 0 0 0 0 0 d_9 d_10 d_11 d_12 d_13 1.
 275 //                      This is 2^(-14) + X_2 truncated to lsb = 2^(-13).
 276 //
 277 //      Define          index_3 := [ d_9 d_10 d_11 d_12 d_13 ].
 278 //
 279 //      Fetch           G_3 := (1/A_3) truncated to 21 sig. bits.
 280 //      floating pt.    Fetch is done using index_3.
 281 //
 282 //      Compute         G := G_1 * G_2 * G_3.
 283 //
 284 //      This is done exactly since each of G_j only has 21 sig. bits.
 285 //
 286 //      Compute
 287 //
 288 //              r := (G*S_hi - 1) + G*S_lo   using 2 FMA operations.
 289 //
 290 //      thus, r approximates G*(S_hi+S_lo) - 1 to within a couple of
 291 //      rounding errors.
 292 //
 293 //
 294 //  Step 2. Approximation
 295 //  ---------------------
 296 //
 297 //   This step computes an approximation to log( 1 + r ) where r is the
 298 //   reduced argument just obtained. It is proved that |r| <= 1.9*2^(-13);
 299 //   thus log(1+r) can be approximated by a short polynomial:
 300 //
 301 //      log(1+r) ~=~ poly = r + Q1 r^2 + ... + Q4 r^5
 302 //
 303 //
 304 //  Step 3. Reconstruction
 305 //  ----------------------
 306 //
 307 //   This step computes the desired result of log(X+E):
 308 //
 309 //      log(X+E)  =   log( 2^N * (S_hi + S_lo) )
 310 //                =   N*log(2) + log( S_hi + S_lo )
 311 //                =   N*log(2) + log(1/G) +
 312 //                    log(1 + C*(S_hi+S_lo) - 1 )
 313 //
 314 //   log(2), log(1/G_j) are stored as pairs of (single,double) numbers:
 315 //   log2_hi, log2_lo, log1byGj_hi, log1byGj_lo. The high parts are
 316 //   single-precision numbers and the low parts are double precision
 317 //   numbers. These have the property that
 318 //
 319 //      N*log2_hi + SUM ( log1byGj_hi )
 320 //
 321 //   is computable exactly in double-extended precision (64 sig. bits).
 322 //   Finally
 323 //
 324 //      Y_hi := N*log2_hi + SUM ( log1byGj_hi )
 325 //      Y_lo := poly_hi + [ poly_lo +
 326 //              ( SUM ( log1byGj_lo ) + N*log2_lo ) ]
 327 //      set lsb(Y_lo) to be 1
 328 //
 329
 330 #include "libm_support.h"
 331
 332 #ifdef _LIBC
 333 .rodata
 334 #else
 335 .data
 336 #endif
 337
 338 // P_7, P_6, P_5, P_4, P_3, P_2, and P_1
 339
 340 .align 64
 341 Constants_P:
 342 ASM_TYPE_DIRECTIVE(Constants_P,@object)
 343 data4  0xEFD62B15,0xE3936754,0x00003FFB,0x00000000
 344 data4  0xA5E56381,0x8003B271,0x0000BFFC,0x00000000
 345 data4  0x73282DB0,0x9249248C,0x00003FFC,0x00000000
 346 data4  0x47305052,0xAAAAAA9F,0x0000BFFC,0x00000000
 347 data4  0xCCD17FC9,0xCCCCCCCC,0x00003FFC,0x00000000
 348 data4  0x00067ED5,0x80000000,0x0000BFFD,0x00000000
 349 data4  0xAAAAAAAA,0xAAAAAAAA,0x00003FFD,0x00000000
 350 data4  0xFFFFFFFE,0xFFFFFFFF,0x0000BFFD,0x00000000
 351 ASM_SIZE_DIRECTIVE(Constants_P)
 352
 353 // log2_hi, log2_lo, Q_4, Q_3, Q_2, and Q_1
 354
 355 .align 64
 356 Constants_Q:
 357 ASM_TYPE_DIRECTIVE(Constants_Q,@object)
 358 data4  0x00000000,0xB1721800,0x00003FFE,0x00000000
 359 data4  0x4361C4C6,0x82E30865,0x0000BFE2,0x00000000
 360 data4  0x328833CB,0xCCCCCAF2,0x00003FFC,0x00000000
 361 data4  0xA9D4BAFB,0x80000077,0x0000BFFD,0x00000000
 362 data4  0xAAABE3D2,0xAAAAAAAA,0x00003FFD,0x00000000
 363 data4  0xFFFFDAB7,0xFFFFFFFF,0x0000BFFD,0x00000000
 364 ASM_SIZE_DIRECTIVE(Constants_Q)
 365
 366 // Z1 - 16 bit fixed, G1 and H1 - IEEE single
 367
 368 .align 64
 369 Constants_Z_G_H_h1:
 370 ASM_TYPE_DIRECTIVE(Constants_Z_G_H_h1,@object)
 371 data4  0x00008000,0x3F800000,0x00000000,0x00000000,0x00000000,0x00000000
 372 data4  0x00007879,0x3F70F0F0,0x3D785196,0x00000000,0x617D741C,0x3DA163A6
 373 data4  0x000071C8,0x3F638E38,0x3DF13843,0x00000000,0xCBD3D5BB,0x3E2C55E6
 374 data4  0x00006BCB,0x3F579430,0x3E2FF9A0,0x00000000,0xD86EA5E7,0xBE3EB0BF
 375 data4  0x00006667,0x3F4CCCC8,0x3E647FD6,0x00000000,0x86B12760,0x3E2E6A8C
 376 data4  0x00006187,0x3F430C30,0x3E8B3AE7,0x00000000,0x5C0739BA,0x3E47574C
 377 data4  0x00005D18,0x3F3A2E88,0x3EA30C68,0x00000000,0x13E8AF2F,0x3E20E30F
 378 data4  0x0000590C,0x3F321640,0x3EB9CEC8,0x00000000,0xF2C630BD,0xBE42885B
 379 data4  0x00005556,0x3F2AAAA8,0x3ECF9927,0x00000000,0x97E577C6,0x3E497F34
 380 data4  0x000051EC,0x3F23D708,0x3EE47FC5,0x00000000,0xA6B0A5AB,0x3E3E6A6E
 381 data4  0x00004EC5,0x3F1D89D8,0x3EF8947D,0x00000000,0xD328D9BE,0xBDF43E3C
 382 data4  0x00004BDB,0x3F17B420,0x3F05F3A1,0x00000000,0x0ADB090A,0x3E4094C3
 383 data4  0x00004925,0x3F124920,0x3F0F4303,0x00000000,0xFC1FE510,0xBE28FBB2
 384 data4  0x0000469F,0x3F0D3DC8,0x3F183EBF,0x00000000,0x10FDE3FA,0x3E3A7895
 385 data4  0x00004445,0x3F088888,0x3F20EC80,0x00000000,0x7CC8C98F,0x3E508CE5
 386 data4  0x00004211,0x3F042108,0x3F29516A,0x00000000,0xA223106C,0xBE534874
 387 ASM_SIZE_DIRECTIVE(Constants_Z_G_H_h1)
 388
 389 // Z2 - 16 bit fixed, G2 and H2 - IEEE single
 390
 391 .align 64
 392 Constants_Z_G_H_h2:
 393 ASM_TYPE_DIRECTIVE(Constants_Z_G_H_h2,@object)
 394 data4  0x00008000,0x3F800000,0x00000000,0x00000000,0x00000000,0x00000000
 395 data4  0x00007F81,0x3F7F00F8,0x3B7F875D,0x00000000,0x22C42273,0x3DB5A116
 396 data4  0x00007F02,0x3F7E03F8,0x3BFF015B,0x00000000,0x21F86ED3,0x3DE620CF
 397 data4  0x00007E85,0x3F7D08E0,0x3C3EE393,0x00000000,0x484F34ED,0xBDAFA07E
 398 data4  0x00007E08,0x3F7C0FC0,0x3C7E0586,0x00000000,0x3860BCF6,0xBDFE07F0
 399 data4  0x00007D8D,0x3F7B1880,0x3C9E75D2,0x00000000,0xA78093D6,0x3DEA370F
 400 data4  0x00007D12,0x3F7A2328,0x3CBDC97A,0x00000000,0x72A753D0,0x3DFF5791
 401 data4  0x00007C98,0x3F792FB0,0x3CDCFE47,0x00000000,0xA7EF896B,0x3DFEBE6C
 402 data4  0x00007C20,0x3F783E08,0x3CFC15D0,0x00000000,0x409ECB43,0x3E0CF156
 403 data4  0x00007BA8,0x3F774E38,0x3D0D874D,0x00000000,0xFFEF71DF,0xBE0B6F97
 404 data4  0x00007B31,0x3F766038,0x3D1CF49B,0x00000000,0x5D59EEE8,0xBE080483
 405 data4  0x00007ABB,0x3F757400,0x3D2C531D,0x00000000,0xA9192A74,0x3E1F91E9
 406 data4  0x00007A45,0x3F748988,0x3D3BA322,0x00000000,0xBF72A8CD,0xBE139A06
 407 data4  0x000079D1,0x3F73A0D0,0x3D4AE46F,0x00000000,0xF8FBA6CF,0x3E1D9202
 408 data4  0x0000795D,0x3F72B9D0,0x3D5A1756,0x00000000,0xBA796223,0xBE1DCCC4
 409 data4  0x000078EB,0x3F71D488,0x3D693B9D,0x00000000,0xB6B7C239,0xBE049391
 410 ASM_SIZE_DIRECTIVE(Constants_Z_G_H_h2)
 411
 412 // G3 and H3 - IEEE single and h3 -IEEE double
 413
 414 .align 64
 415 Constants_Z_G_H_h3:
 416 ASM_TYPE_DIRECTIVE(Constants_Z_G_H_h3,@object)
 417 data4  0x3F7FFC00,0x38800100,0x562224CD,0x3D355595
 418 data4  0x3F7FF400,0x39400480,0x06136FF6,0x3D8200A2
 419 data4  0x3F7FEC00,0x39A00640,0xE8DE9AF0,0x3DA4D68D
 420 data4  0x3F7FE400,0x39E00C41,0xB10238DC,0xBD8B4291
 421 data4  0x3F7FDC00,0x3A100A21,0x3B1952CA,0xBD89CCB8
 422 data4  0x3F7FD400,0x3A300F22,0x1DC46826,0xBDB10707
 423 data4  0x3F7FCC08,0x3A4FF51C,0xF43307DB,0x3DB6FCB9
 424 data4  0x3F7FC408,0x3A6FFC1D,0x62DC7872,0xBD9B7C47
 425 data4  0x3F7FBC10,0x3A87F20B,0x3F89154A,0xBDC3725E
 426 data4  0x3F7FB410,0x3A97F68B,0x62B9D392,0xBD93519D
 427 data4  0x3F7FAC18,0x3AA7EB86,0x0F21BD9D,0x3DC18441
 428 data4  0x3F7FA420,0x3AB7E101,0x2245E0A6,0xBDA64B95
 429 data4  0x3F7F9C20,0x3AC7E701,0xAABB34B8,0x3DB4B0EC
 430 data4  0x3F7F9428,0x3AD7DD7B,0x6DC40A7E,0x3D992337
 431 data4  0x3F7F8C30,0x3AE7D474,0x4F2083D3,0x3DC6E17B
 432 data4  0x3F7F8438,0x3AF7CBED,0x811D4394,0x3DAE314B
 433 data4  0x3F7F7C40,0x3B03E1F3,0xB08F2DB1,0xBDD46F21
 434 data4  0x3F7F7448,0x3B0BDE2F,0x6D34522B,0xBDDC30A4
 435 data4  0x3F7F6C50,0x3B13DAAA,0xB1F473DB,0x3DCB0070
 436 data4  0x3F7F6458,0x3B1BD766,0x6AD282FD,0xBDD65DDC
 437 data4  0x3F7F5C68,0x3B23CC5C,0xF153761A,0xBDCDAB83
 438 data4  0x3F7F5470,0x3B2BC997,0x341D0F8F,0xBDDADA40
 439 data4  0x3F7F4C78,0x3B33C711,0xEBC394E8,0x3DCD1BD7
 440 data4  0x3F7F4488,0x3B3BBCC6,0x52E3E695,0xBDC3532B
 441 data4  0x3F7F3C90,0x3B43BAC0,0xE846B3DE,0xBDA3961E
 442 data4  0x3F7F34A0,0x3B4BB0F4,0x785778D4,0xBDDADF06
 443 data4  0x3F7F2CA8,0x3B53AF6D,0xE55CE212,0x3DCC3ED1
 444 data4  0x3F7F24B8,0x3B5BA620,0x9E382C15,0xBDBA3103
 445 data4  0x3F7F1CC8,0x3B639D12,0x5C5AF197,0x3D635A0B
 446 data4  0x3F7F14D8,0x3B6B9444,0x71D34EFC,0xBDDCCB19
 447 data4  0x3F7F0CE0,0x3B7393BC,0x52CD7ADA,0x3DC74502
 448 data4  0x3F7F04F0,0x3B7B8B6D,0x7D7F2A42,0xBDB68F17
 449 ASM_SIZE_DIRECTIVE(Constants_Z_G_H_h3)
 450
 451 //
 452 //  Exponent Thresholds and Tiny Thresholds
 453 //  for 8, 11, 15, and 17 bit exponents
 454 //
 455 //  Expo_Range             Value
 456 //
 457 //  0 (8  bits)            2^(-126)
 458 //  1 (11 bits)            2^(-1022)
 459 //  2 (15 bits)            2^(-16382)
 460 //  3 (17 bits)            2^(-16382)
 461 //
 462 //  Tiny_Table
 463 //  ----------
 464 //  Expo_Range             Value
 465 //
 466 //  0 (8  bits)            2^(-16382)
 467 //  1 (11 bits)            2^(-16382)
 468 //  2 (15 bits)            2^(-16382)
 469 //  3 (17 bits)            2^(-16382)
 470 //
 471
 472 .align 64
 473 Constants_Threshold:
 474 ASM_TYPE_DIRECTIVE(Constants_Threshold,@object)
 475 data4  0x00000000,0x80000000,0x00003F81,0x00000000
 476 data4  0x00000000,0x80000000,0x00000001,0x00000000
 477 data4  0x00000000,0x80000000,0x00003C01,0x00000000
 478 data4  0x00000000,0x80000000,0x00000001,0x00000000
 479 data4  0x00000000,0x80000000,0x00000001,0x00000000
 480 data4  0x00000000,0x80000000,0x00000001,0x00000000
 481 data4  0x00000000,0x80000000,0x00000001,0x00000000
 482 data4  0x00000000,0x80000000,0x00000001,0x00000000
 483 ASM_SIZE_DIRECTIVE(Constants_Threshold)
 484
 485 .align 64
 486 Constants_1_by_LN10:
 487 ASM_TYPE_DIRECTIVE(Constants_1_by_LN10,@object)
 488 data4  0x37287195,0xDE5BD8A9,0x00003FFD,0x00000000
 489 data4  0xACCF70C8,0xD56EAABE,0x00003FBD,0x00000000
 490 ASM_SIZE_DIRECTIVE(Constants_1_by_LN10)
 491
 492 FR_Input_X = f8
 493 FR_Neg_One = f9
 494 FR_E       = f33
 495 FR_Em1     = f34
 496 FR_Y_hi    = f34
 497 // Shared with Em1
 498 FR_Y_lo    = f35
 499 FR_Scale   = f36
 500 FR_X_Prime = f37
 501 FR_Z       = f38
 502 FR_S_hi    = f38
 503 // Shared with Z
 504 FR_W       = f39
 505 FR_G       = f40
 506 FR_wsq     = f40
 507 // Shared with G
 508 FR_H       = f41
 509 FR_w4      = f41
 510 // Shared with H
 511 FR_h       = f42
 512 FR_w6      = f42
 513 // Shared with h
 514 FR_G_tmp   = f43
 515 FR_poly_lo = f43
 516 // Shared with G_tmp
 517 FR_P8      = f43
 518 // Shared with G_tmp
 519 FR_H_tmp   = f44
 520 FR_poly_hi = f44
 521   // Shared with H_tmp
 522 FR_P7      = f44
 523 // Shared with H_tmp
 524 FR_h_tmp   = f45
 525 FR_rsq     = f45
 526 // Shared with h_tmp
 527 FR_P6      = f45
 528 // Shared with h_tmp
 529 FR_abs_W   = f46
 530 FR_r       = f46
 531 // Shared with abs_W
 532 FR_AA      = f47
 533 FR_log2_hi = f47
 534 // Shared with AA
 535 FR_BB          = f48
 536 FR_log2_lo     = f48
 537 // Shared with BB
 538 FR_S_lo        = f49
 539 FR_two_negN    = f50
 540 FR_float_N     = f51
 541 FR_Q4          = f52
 542 FR_dummy       = f52
 543 // Shared with Q4
 544 FR_P4          = f52
 545 // Shared with Q4
 546 FR_Threshold    = f52
 547 // Shared with Q4
 548 FR_Q3          = f53
 549 FR_P3          = f53
 550 // Shared with Q3
 551 FR_Tiny        = f53
 552 // Shared with Q3
 553 FR_Q2          = f54
 554 FR_P2          = f54
 555 // Shared with Q2
 556 FR_1LN10_hi     = f54
 557 // Shared with Q2
 558 FR_Q1           = f55
 559 FR_P1           = f55
 560 // Shared with Q1
 561 FR_1LN10_lo     = f55
 562 // Shared with Q1
 563 FR_P5           = f98
 564 FR_SCALE        = f98
 565 FR_Output_X_tmp = f99
 566
 567 GR_Expo_Range   = r32
 568 GR_Table_Base   = r34
 569 GR_Table_Base1  = r35
 570 GR_Table_ptr    = r36
 571 GR_Index2       = r37
 572 GR_signif       = r38
 573 GR_X_0          = r39
 574 GR_X_1          = r40
 575 GR_X_2          = r41
 576 GR_Z_1          = r42
 577 GR_Z_2          = r43
 578 GR_N            = r44
 579 GR_Bias         = r45
 580 GR_M            = r46
 581 GR_ScaleN       = r47
 582 GR_Index3       = r48
 583 GR_Perturb      = r49
 584 GR_Table_Scale  = r50
 585
 586
 587 GR_SAVE_PFS     = r51
 588 GR_SAVE_B0      = r52
 589 GR_SAVE_GP      = r53
 590
 591 GR_Parameter_X       = r54
 592 GR_Parameter_Y       = r55
 593 GR_Parameter_RESULT  = r56
 594
 595 GR_Parameter_TAG = r57
 596
 597
 598 .section .text
 599 .proc log1p#
 600 .global log1p#
 601 .align 64
 602 log1p:
 603 #ifdef _LIBC
 604 .global __log1p
 605 __log1p:
 606 #endif
 607
 608 { .mfi
 609 alloc r32 = ar.pfs,0,22,4,0
 610 (p0)  fsub.s1 FR_Neg_One = f0,f1
 611 (p0)  cmp.eq.unc  p7, p0 = r0, r0
 612 }
 613
 614 { .mfi
 615 (p0)  cmp.ne.unc  p14, p0 = r0, r0
 616 (p0)  fnorm.s1 FR_X_Prime = FR_Input_X
 617 (p0)  cmp.eq.unc  p15, p0 = r0, r0 ;;
 618 }
 619
 620 { .mfi
 621       nop.m 999
 622 (p0)  fclass.m.unc p6, p0 =  FR_Input_X, 0x1E3
 623       nop.i 999
 624 }
 625 ;;
 626
 627 { .mfi
 628         nop.m 999
 629 (p0)  fclass.nm.unc p10, p0 =  FR_Input_X, 0x1FF
 630       nop.i 999
 631 }
 632 ;;
 633
 634 { .mfi
 635         nop.m 999
 636 (p0)  fcmp.eq.unc.s1 p9, p0 =  FR_Input_X, f0
 637       nop.i 999
 638 }
 639
 640 { .mfi
 641         nop.m 999
 642 (p0)  fadd FR_Em1 = f0,f0
 643         nop.i 999 ;;
 644 }
 645
 646 { .mfi
 647         nop.m 999
 648 (p0)  fadd FR_E = f0,f1
 649         nop.i 999 ;;
 650 }
 651
 652 { .mfi
 653         nop.m 999
 654 (p0)  fcmp.eq.unc.s1 p8, p0 =  FR_Input_X, FR_Neg_One
 655         nop.i 999
 656 }
 657
 658 { .mfi
 659         nop.m 999
 660 (p0)  fcmp.lt.unc.s1 p13, p0 =  FR_Input_X, FR_Neg_One
 661         nop.i 999
 662 }
 663
 664
 665 L(LOG_BEGIN):
 666
 667 { .mfi
 668         nop.m 999
 669 (p0)  fadd.s1 FR_Z = FR_X_Prime, FR_E
 670         nop.i 999
 671 }
 672
 673 { .mlx
 674         nop.m 999
 675 (p0)  movl GR_Table_Scale = 0x0000000000000018 ;;
 676 }
 677
 678 { .mmi
 679         nop.m 999
 680 //
 681 //    Create E = 1 and Em1 = 0
 682 //    Check for X == 0, meaning log(1+0)
 683 //    Check for X < -1, meaning log(negative)
 684 //    Check for X == -1, meaning log(0)
 685 //    Normalize x
 686 //    Identify NatVals, NaNs, Infs.
 687 //    Identify EM unsupporteds.
 688 //    Identify Negative values - us S1 so as
 689 //    not to raise denormal operand exception
 690 //    Set p15 to true for log1p
 691 //    Set p14 to false for log1p
 692 //    Set p7 true for log and log1p
 693 //
 694 (p0)  addl GR_Table_Base = @ltoff(Constants_Z_G_H_h1#),gp
 695       nop.i  999
 696 }
 697
 698 { .mfi
 699         nop.m 999
 700 (p0)  fmax.s1 FR_AA = FR_X_Prime, FR_E
 701         nop.i 999 ;;
 702 }
 703
 704 { .mfi
 705       ld8    GR_Table_Base = [GR_Table_Base]
 706 (p0)  fmin.s1 FR_BB = FR_X_Prime, FR_E
 707         nop.i 999
 708 }
 709
 710 { .mfb
 711         nop.m 999
 712 (p0)  fadd.s1 FR_W = FR_X_Prime, FR_Em1
 713 //
 714 //    Begin load of constants base
 715 //    FR_Z = Z = |x| + E
 716 //    FR_W = W = |x| + Em1
 717 //    AA = fmax(|x|,E)
 718 //    BB = fmin(|x|,E)
 719 //
 720 (p6)  br.cond.spnt L(LOG_64_special) ;;
 721 }
 722
 723 { .mib
 724         nop.m 999
 725         nop.i 999
 726 (p10) br.cond.spnt L(LOG_64_unsupported) ;;
 727 }
 728
 729 { .mib
 730         nop.m 999
 731         nop.i 999
 732 (p13) br.cond.spnt L(LOG_64_negative) ;;
 733 }
 734
 735 { .mib
 736 (p0)  getf.sig GR_signif = FR_Z
 737         nop.i 999
 738 (p9)  br.cond.spnt L(LOG_64_one) ;;
 739 }
 740
 741 { .mib
 742         nop.m 999
 743         nop.i 999
 744 (p8)  br.cond.spnt L(LOG_64_zero) ;;
 745 }
 746
 747 { .mfi
 748 (p0)  getf.exp GR_N =  FR_Z
 749 //
 750 //    Raise possible denormal operand exception
 751 //    Create Bias
 752 //
 753 //    This function computes ln( x + e )
 754 //    Input  FR 1: FR_X   = FR_Input_X
 755 //    Input  FR 2: FR_E   = FR_E
 756 //    Input  FR 3: FR_Em1 = FR_Em1
 757 //    Input  GR 1: GR_Expo_Range = GR_Expo_Range = 1
 758 //    Output FR 4: FR_Y_hi
 759 //    Output FR 5: FR_Y_lo
 760 //    Output FR 6: FR_Scale
 761 //    Output PR 7: PR_Safe
 762 //
 763 (p0)  fsub.s1 FR_S_lo = FR_AA, FR_Z
 764 //
 765 //    signif = getf.sig(Z)
 766 //    abs_W = fabs(w)
 767 //
 768 (p0)  extr.u GR_Table_ptr = GR_signif, 59, 4 ;;
 769 }
 770
 771 { .mfi
 772         nop.m 999
 773 (p0)  fmerge.se FR_S_hi =  f1,FR_Z
 774 (p0)  extr.u GR_X_0 = GR_signif, 49, 15
 775 }
 776
 777 { .mmi
 778       nop.m 999
 779 (p0)  addl GR_Table_Base1 = @ltoff(Constants_Z_G_H_h2#),gp
 780       nop.i 999
 781 }
 782 ;;
 783
 784 { .mlx
 785       ld8    GR_Table_Base1 = [GR_Table_Base1]
 786 (p0)  movl GR_Bias = 0x000000000000FFFF ;;
 787 }
 788
 789 { .mfi
 790         nop.m 999
 791 (p0)  fabs FR_abs_W =  FR_W
 792 (p0)  pmpyshr2.u GR_Table_ptr = GR_Table_ptr,GR_Table_Scale,0
 793 }
 794
 795 { .mfi
 796         nop.m 999
 797 //
 798 //    Branch out for special input values
 799 //
 800 (p0)  fcmp.lt.unc.s0 p8, p0 =  FR_Input_X, f0
 801         nop.i 999 ;;
 802 }
 803
 804 { .mfi
 805         nop.m 999
 806 //
 807 //    X_0 = extr.u(signif,49,15)
 808 //    Index1 = extr.u(signif,59,4)
 809 //
 810 (p0)  fadd.s1 FR_S_lo = FR_S_lo, FR_BB
 811         nop.i 999 ;;
 812 }
 813
 814 { .mii
 815         nop.m 999
 816         nop.i 999 ;;
 817 //
 818 //    Offset_to_Z1 = 24 * Index1
 819 //    For performance, don't use result
 820 //    for 3 or 4 cycles.
 821 //
 822 (p0)  add GR_Table_ptr = GR_Table_ptr, GR_Table_Base ;;
 823 }
 824 //
 825 //    Add Base to Offset for Z1
 826 //    Create Bias
 827
 828 { .mmi
 829 (p0)  ld4 GR_Z_1 = [GR_Table_ptr],4 ;;
 830 (p0)  ldfs  FR_G = [GR_Table_ptr],4
 831         nop.i 999 ;;
 832 }
 833
 834 { .mmi
 835 (p0)  ldfs  FR_H = [GR_Table_ptr],8 ;;
 836 (p0)  ldfd  FR_h = [GR_Table_ptr],0
 837 (p0)  pmpyshr2.u GR_X_1 = GR_X_0,GR_Z_1,15
 838 }
 839 //
 840 //    Load Z_1
 841 //    Get Base of Table2
 842 //
 843
 844 { .mfi
 845 (p0)  getf.exp GR_M = FR_abs_W
 846         nop.f 999
 847         nop.i 999 ;;
 848 }
 849
 850 { .mii
 851         nop.m 999
 852         nop.i 999 ;;
 853 //
 854 //    M = getf.exp(abs_W)
 855 //    S_lo = AA - Z
 856 //    X_1 = pmpyshr2(X_0,Z_1,15)
 857 //
 858 (p0)  sub GR_M = GR_M, GR_Bias ;;
 859 }
 860 //
 861 //    M = M - Bias
 862 //    Load G1
 863 //    N = getf.exp(Z)
 864 //
 865
 866 { .mii
 867 (p0)  cmp.gt.unc  p11, p0 =  -80, GR_M
 868 (p0)  cmp.gt.unc  p12, p0 =  -7, GR_M ;;
 869 (p0)  extr.u GR_Index2 = GR_X_1, 6, 4 ;;
 870 }
 871
 872 { .mib
 873         nop.m 999
 874 //
 875 //    if -80 > M, set p11
 876 //    Index2 = extr.u(X_1,6,4)
 877 //    if -7  > M, set p12
 878 //    Load H1
 879 //
 880 (p0)  pmpyshr2.u GR_Index2 = GR_Index2,GR_Table_Scale,0
 881 (p11) br.cond.spnt L(log1p_small) ;;
 882 }
 883
 884 { .mib
 885       nop.m 999
 886         nop.i 999
 887 (p12) br.cond.spnt L(log1p_near) ;;
 888 }
 889
 890 { .mii
 891 (p0)  sub GR_N = GR_N, GR_Bias
 892 //
 893 //    poly_lo = r * poly_lo
 894 //
 895 (p0)  add GR_Perturb = 0x1, r0 ;;
 896 (p0)  sub GR_ScaleN = GR_Bias, GR_N
 897 }
 898
 899 { .mii
 900 (p0)  setf.sig FR_float_N = GR_N
 901         nop.i 999 ;;
 902 //
 903 //    Prepare Index2 - pmpyshr2.u(X_1,Z_2,15)
 904 //    Load h1
 905 //    S_lo = S_lo + BB
 906 //    Branch for -80 > M
 907 //
 908 (p0)  add GR_Index2 = GR_Index2, GR_Table_Base1
 909 }
 910
 911 { .mmi
 912 (p0)  setf.exp FR_two_negN = GR_ScaleN
 913       nop.m 999
 914 (p0)  addl GR_Table_Base = @ltoff(Constants_Z_G_H_h3#),gp
 915 };;
 916
 917 //
 918 //    Index2 points to Z2
 919 //    Branch for -7 > M
 920 //
 921
 922 { .mmb
 923 (p0)  ld4 GR_Z_2 = [GR_Index2],4
 924       ld8 GR_Table_Base = [GR_Table_Base]
 925       nop.b 999 ;;
 926 }
 927 (p0)  nop.i 999
 928 //
 929 //    Load Z_2
 930 //    N = N - Bias
 931 //    Tablebase points to Table3
 932 //
 933
 934 { .mmi
 935 (p0)  ldfs  FR_G_tmp = [GR_Index2],4 ;;
 936 //
 937 //    Load G_2
 938 //    pmpyshr2  X_2= (X_1,Z_2,15)
 939 //    float_N = setf.sig(N)
 940 //    ScaleN = Bias - N
 941 //
 942 (p0)  ldfs  FR_H_tmp = [GR_Index2],8
 943         nop.i 999 ;;
 944 }
 945 //
 946 //    Load H_2
 947 //    two_negN = setf.exp(scaleN)
 948 //    G = G_1 * G_2
 949 //
 950
 951 { .mfi
 952 (p0)  ldfd  FR_h_tmp = [GR_Index2],0
 953         nop.f 999
 954 (p0)  pmpyshr2.u GR_X_2 = GR_X_1,GR_Z_2,15 ;;
 955 }
 956
 957 { .mii
 958         nop.m 999
 959 (p0)  extr.u GR_Index3 = GR_X_2, 1, 5 ;;
 960 //
 961 //    Load h_2
 962 //    H = H_1 + H_2
 963 //    h = h_1 + h_2
 964 //    Index3 = extr.u(X_2,1,5)
 965 //
 966 (p0)  shladd GR_Index3 = GR_Index3,4,GR_Table_Base
 967 }
 968
 969 { .mmi
 970         nop.m 999
 971         nop.m 999
 972 //
 973 //    float_N = fcvt.xf(float_N)
 974 //    load G3
 975 //
 976 (p0)  addl GR_Table_Base = @ltoff(Constants_Q#),gp ;;
 977 }
 978
 979 { .mfi
 980 ld8    GR_Table_Base = [GR_Table_Base]
 981 nop.f 999
 982 nop.i 999
 983 } ;;
 984
 985 { .mfi
 986 (p0)  ldfe FR_log2_hi = [GR_Table_Base],16
 987 (p0)  fmpy.s1 FR_S_lo = FR_S_lo, FR_two_negN
 988         nop.i 999 ;;
 989 }
 990
 991 { .mmf
 992         nop.m 999
 993 //
 994 //    G = G3 * G
 995 //    Load h3
 996 //    Load log2_hi
 997 //    H = H + H3
 998 //
 999 (p0)  ldfe FR_log2_lo = [GR_Table_Base],16
1000 (p0)  fmpy.s1 FR_G = FR_G, FR_G_tmp ;;
1001 }
1002
1003 { .mmf
1004 (p0)  ldfs  FR_G_tmp = [GR_Index3],4
1005 //
1006 //    h = h + h3
1007 //    r = G * S_hi + 1
1008 //    Load log2_lo
1009 //
1010 (p0)  ldfe FR_Q4 = [GR_Table_Base],16
1011 (p0)  fadd.s1 FR_h = FR_h, FR_h_tmp ;;
1012 }
1013
1014 { .mfi
1015 (p0)  ldfe FR_Q3 = [GR_Table_Base],16
1016 (p0)  fadd.s1 FR_H = FR_H, FR_H_tmp
1017         nop.i 999 ;;
1018 }
1019
1020 { .mmf
1021 (p0)  ldfs  FR_H_tmp = [GR_Index3],4
1022 (p0)  ldfe FR_Q2 = [GR_Table_Base],16
1023 //
1024 //    Comput Index for Table3
1025 //    S_lo = S_lo * two_negN
1026 //
1027 (p0)  fcvt.xf FR_float_N = FR_float_N ;;
1028 }
1029 //
1030 //    If S_lo == 0, set p8 false
1031 //    Load H3
1032 //    Load ptr to table of polynomial coeff.
1033 //
1034
1035 { .mmf
1036 (p0)  ldfd  FR_h_tmp = [GR_Index3],0
1037 (p0)  ldfe FR_Q1 = [GR_Table_Base],0
1038 (p0)  fcmp.eq.unc.s1 p0, p8 =  FR_S_lo, f0 ;;
1039 }
1040
1041 { .mfi
1042         nop.m 999
1043 (p0)  fmpy.s1 FR_G = FR_G, FR_G_tmp
1044         nop.i 999 ;;
1045 }
1046
1047 { .mfi
1048         nop.m 999
1049 (p0)  fadd.s1 FR_H = FR_H, FR_H_tmp
1050         nop.i 999 ;;
1051 }
1052
1053 { .mfi
1054         nop.m 999
1055 (p0)  fms.s1 FR_r = FR_G, FR_S_hi, f1
1056         nop.i 999
1057 }
1058
1059 { .mfi
1060         nop.m 999
1061 (p0)  fadd.s1 FR_h = FR_h, FR_h_tmp
1062         nop.i 999 ;;
1063 }
1064
1065 { .mfi
1066         nop.m 999
1067 (p0)  fma.s1 FR_Y_hi = FR_float_N, FR_log2_hi, FR_H
1068         nop.i 999 ;;
1069 }
1070
1071 { .mfi
1072         nop.m 999
1073 //
1074 //    Load Q4
1075 //    Load Q3
1076 //    Load Q2
1077 //    Load Q1
1078 //
1079 (p8) fma.s1 FR_r = FR_G, FR_S_lo, FR_r
1080         nop.i 999
1081 }
1082
1083 { .mfi
1084         nop.m 999
1085 //
1086 //    poly_lo = r * Q4 + Q3
1087 //    rsq = r* r
1088 //
1089 (p0)  fma.s1 FR_h = FR_float_N, FR_log2_lo, FR_h
1090         nop.i 999 ;;
1091 }
1092
1093 { .mfi
1094         nop.m 999
1095 //
1096 //    If (S_lo!=0) r = s_lo * G + r
1097 //
1098 (p0)  fma.s1 FR_poly_lo = FR_r, FR_Q4, FR_Q3
1099         nop.i 999
1100 }
1101 //
1102 //    Create a 0x00000....01
1103 //    poly_lo = poly_lo * rsq + h
1104 //
1105
1106 { .mfi
1107 (p0)  setf.sig FR_dummy = GR_Perturb
1108 (p0)  fmpy.s1 FR_rsq = FR_r, FR_r
1109         nop.i 999 ;;
1110 }
1111
1112 { .mfi
1113         nop.m 999
1114 //
1115 //    h = N * log2_lo + h
1116 //    Y_hi = n * log2_hi + H
1117 //
1118 (p0)  fma.s1 FR_poly_lo = FR_poly_lo, FR_r, FR_Q2
1119         nop.i 999
1120 }
1121
1122 { .mfi
1123         nop.m 999
1124 (p0)  fma.s1 FR_poly_hi = FR_Q1, FR_rsq, FR_r
1125         nop.i 999 ;;
1126 }
1127
1128 { .mfi
1129         nop.m 999
1130 //
1131 //    poly_lo = r * poly_o + Q2
1132 //    poly_hi = Q1 * rsq + r
1133 //
1134 (p0)  fmpy.s1 FR_poly_lo = FR_poly_lo, FR_r
1135         nop.i 999 ;;
1136 }
1137
1138 { .mfi
1139         nop.m 999
1140 (p0)  fma.s1 FR_poly_lo = FR_poly_lo, FR_rsq, FR_h
1141         nop.i 999 ;;
1142 }
1143
1144 { .mfb
1145         nop.m 999
1146 (p0)  fadd.s1 FR_Y_lo = FR_poly_hi, FR_poly_lo
1147 //
1148 //    Create the FR for a binary "or"
1149 //    Y_lo = poly_hi + poly_lo
1150 //
1151 // (p0)  for FR_dummy = FR_Y_lo,FR_dummy ;;
1152 //
1153 //    Turn the lsb of Y_lo ON
1154 //
1155 // (p0)  fmerge.se FR_Y_lo =  FR_Y_lo,FR_dummy ;;
1156 //
1157 //    Merge the new lsb into Y_lo, for alone doesn't
1158 //
1159 (p0)  br.cond.sptk L(LOG_main) ;;
1160 }
1161
1162
1163 L(log1p_near):
1164
1165 { .mmi
1166         nop.m 999
1167         nop.m 999
1168 //    /*******************************************************/
1169 //    /*********** Branch log1p_near  ************************/
1170 //    /*******************************************************/
1171 (p0)  addl GR_Table_Base = @ltoff(Constants_P#),gp ;;
1172 }
1173 //
1174 //    Load base address of poly. coeff.
1175 //
1176 {.mmi
1177       nop.m 999
1178       ld8    GR_Table_Base = [GR_Table_Base]
1179       nop.i 999
1180 };;
1181
1182 { .mmb
1183 (p0)  add GR_Table_ptr = 0x40,GR_Table_Base
1184 //
1185 //    Address tables with separate pointers
1186 //
1187 (p0)  ldfe FR_P8 = [GR_Table_Base],16
1188         nop.b 999 ;;
1189 }
1190
1191 { .mmb
1192 (p0)  ldfe FR_P4 = [GR_Table_ptr],16
1193 //
1194 //    Load P4
1195 //    Load P8
1196 //
1197 (p0)  ldfe FR_P7 = [GR_Table_Base],16
1198         nop.b 999 ;;
1199 }
1200
1201 { .mmf
1202 (p0)  ldfe FR_P3 = [GR_Table_ptr],16
1203 //
1204 //    Load P3
1205 //    Load P7
1206 //
1207 (p0)  ldfe FR_P6 = [GR_Table_Base],16
1208 (p0)  fmpy.s1 FR_wsq = FR_W, FR_W ;;
1209 }
1210
1211 { .mfi
1212 (p0)  ldfe FR_P2 = [GR_Table_ptr],16
1213         nop.f 999
1214         nop.i 999 ;;
1215 }
1216
1217 { .mfi
1218         nop.m 999
1219 (p0)  fma.s1 FR_Y_hi = FR_W, FR_P4, FR_P3
1220         nop.i 999
1221 }
1222 //
1223 //    Load P2
1224 //    Load P6
1225 //    Wsq = w * w
1226 //    Y_hi = p4 * w + p3
1227 //
1228
1229 { .mfi
1230 (p0)  ldfe FR_P5 = [GR_Table_Base],16
1231 (p0)  fma.s1 FR_Y_lo = FR_W, FR_P8, FR_P7
1232         nop.i 999 ;;
1233 }
1234
1235 { .mfi
1236 (p0)  ldfe FR_P1 = [GR_Table_ptr],16
1237 //
1238 //    Load P1
1239 //    Load P5
1240 //    Y_lo = p8 * w + P7
1241 //
1242 (p0)  fmpy.s1 FR_w4 = FR_wsq, FR_wsq
1243         nop.i 999 ;;
1244 }
1245
1246 { .mfi
1247         nop.m 999
1248 (p0)  fma.s1 FR_Y_hi = FR_W, FR_Y_hi, FR_P2
1249         nop.i 999
1250 }
1251
1252 { .mfi
1253         nop.m 999
1254 (p0)  fma.s1 FR_Y_lo = FR_W, FR_Y_lo, FR_P6
1255 (p0)  add GR_Perturb = 0x1, r0 ;;
1256 }
1257
1258 { .mfi
1259         nop.m 999
1260 //
1261 //    w4 = w2 * w2
1262 //    Y_hi = y_hi * w + p2
1263 //    Y_lo = y_lo * w + p6
1264 //    Create perturbation bit
1265 //
1266 (p0)  fmpy.s1 FR_w6 = FR_w4, FR_wsq
1267         nop.i 999 ;;
1268 }
1269
1270 { .mfi
1271         nop.m 999
1272 (p0)  fma.s1 FR_Y_hi = FR_W, FR_Y_hi, FR_P1
1273         nop.i 999
1274 }
1275 //
1276 //    Y_hi = y_hi * w + p1
1277 //    w6 = w4 * w2
1278 //
1279
1280 { .mfi
1281 (p0)  setf.sig FR_Q4 = GR_Perturb
1282 (p0)  fma.s1 FR_Y_lo = FR_W, FR_Y_lo, FR_P5
1283         nop.i 999 ;;
1284 }
1285
1286 { .mfi
1287         nop.m 999
1288 (p0)  fma.s1 FR_Y_hi = FR_wsq,FR_Y_hi, FR_W
1289         nop.i 999
1290 }
1291
1292 { .mfb
1293         nop.m 999
1294 //
1295 //    Y_hi = y_hi * wsq + w
1296 //    Y_lo = y_lo * w + p5
1297 //
1298 (p0)  fmpy.s1 FR_Y_lo = FR_w6, FR_Y_lo
1299 //
1300 //    Y_lo = y_lo * w6
1301 //
1302 // (p0)  for FR_dummy = FR_Y_lo,FR_dummy ;;
1303 //
1304 //    Set lsb on: Taken out to improve performance
1305 //
1306 // (p0)  fmerge.se FR_Y_lo =  FR_Y_lo,FR_dummy ;;
1307 //
1308 //    Make sure it's on in Y_lo also.  Taken out to improve
1309 //    performance
1310 //
1311 (p0)  br.cond.sptk L(LOG_main) ;;
1312 }
1313
1314
1315 L(log1p_small):
1316
1317 { .mmi
1318         nop.m 999
1319         nop.m 999
1320 //  /*******************************************************/
1321 //  /*********** Branch log1p_small  ***********************/
1322 //  /*******************************************************/
1323 (p0)  addl GR_Table_Base = @ltoff(Constants_Threshold#),gp
1324 }
1325
1326 { .mfi
1327         nop.m 999
1328 (p0)  mov FR_Em1 = FR_W
1329 (p0)  cmp.eq.unc  p7, p0 = r0, r0 ;;
1330 }
1331
1332 { .mlx
1333       ld8    GR_Table_Base = [GR_Table_Base]
1334 (p0)  movl GR_Expo_Range = 0x0000000000000002 ;;
1335 }
1336 //
1337 //    Set Safe to true
1338 //    Set Expo_Range = 0 for single
1339 //    Set Expo_Range = 2 for double
1340 //    Set Expo_Range = 4 for double-extended
1341 //
1342
1343 { .mmi
1344 (p0)  shladd GR_Table_Base = GR_Expo_Range,4,GR_Table_Base ;;
1345 (p0)  ldfe FR_Threshold = [GR_Table_Base],16
1346         nop.i 999
1347 }
1348
1349 { .mlx
1350         nop.m 999
1351 (p0)  movl GR_Bias = 0x000000000000FF9B ;;
1352 }
1353
1354 { .mfi
1355 (p0)  ldfe FR_Tiny = [GR_Table_Base],0
1356         nop.f 999
1357         nop.i 999 ;;
1358 }
1359
1360 { .mfi
1361         nop.m 999
1362 (p0)  fcmp.gt.unc.s1 p13, p12 =  FR_abs_W, FR_Threshold
1363         nop.i 999 ;;
1364 }
1365
1366 { .mfi
1367         nop.m 999
1368 (p13) fnmpy.s1 FR_Y_lo = FR_W, FR_W
1369         nop.i 999
1370 }
1371
1372 { .mfi
1373         nop.m 999
1374 (p13) fadd FR_SCALE = f0, f1
1375         nop.i 999 ;;
1376 }
1377
1378 { .mfi
1379         nop.m 999
1380 (p12) fsub.s1 FR_Y_lo = f0, FR_Tiny
1381 (p12) cmp.ne.unc  p7, p0 = r0, r0
1382 }
1383
1384 { .mfi
1385 (p12) setf.exp FR_SCALE = GR_Bias
1386         nop.f 999
1387         nop.i 999 ;;
1388 }
1389
1390 //
1391 //    Set p7 to SAFE = FALSE
1392 //    Set Scale = 2^-100
1393 //
1394 { .mfb
1395         nop.m 999
1396 (p0)  fma.d.s0 FR_Input_X = FR_Y_lo,FR_SCALE,FR_Y_hi
1397 (p0)  br.ret.sptk   b0
1398 }
1399 ;;
1400
1401 L(LOG_64_one):
1402
1403 { .mfb
1404         nop.m 999
1405 (p0)  fmpy.d.s0 FR_Input_X = FR_Input_X, f0
1406 (p0)  br.ret.sptk   b0
1407 }
1408 ;;
1409
1410 //
1411 //    Raise divide by zero for +/-0 input.
1412 //
1413 L(LOG_64_zero):
1414
1415 { .mfi
1416 (p0)  mov   GR_Parameter_TAG = 140
1417 //
1418 //    If we have log1p(0), return -Inf.
1419 //
1420 (p0)  fsub.s0 FR_Output_X_tmp = f0, f1
1421       nop.i 999 ;;
1422 }
1423 { .mfb
1424       nop.m 999
1425 (p0)  frcpa.s0 FR_Output_X_tmp, p8 =  FR_Output_X_tmp, f0
1426 (p0)  br.cond.sptk L(LOG_ERROR_Support) ;;
1427 }
1428
1429 L(LOG_64_special):
1430
1431 { .mfi
1432       nop.m 999
1433 //
1434 //    Return -Inf or value from handler.
1435 //
1436 (p0)  fclass.m.unc p7, p0 =  FR_Input_X, 0x1E1
1437       nop.i 999 ;;
1438 }
1439 { .mfb
1440       nop.m 999
1441 //
1442 //    Check for Natval, QNan, SNaN, +Inf
1443 //
1444 (p7)  fmpy.d.s0  f8 =  FR_Input_X, f1
1445 //
1446 //    For SNaN raise invalid and return QNaN.
1447 //    For QNaN raise invalid and return QNaN.
1448 //    For +Inf return +Inf.
1449 //
1450 (p7)  br.ret.sptk   b0
1451 }
1452 ;;
1453
1454 //
1455 //    For -Inf raise invalid and return QNaN.
1456 //
1457
1458 { .mfb
1459 (p0)  mov   GR_Parameter_TAG = 141
1460 (p0)  fmpy.d.s0  FR_Output_X_tmp =  FR_Input_X, f0
1461 (p0)  br.cond.sptk L(LOG_ERROR_Support) ;;
1462 }
1463
1464 //
1465 //    Report that log1p(-Inf) computed
1466 //
1467
1468 L(LOG_64_unsupported):
1469
1470 //
1471 //    Return generated NaN or other value .
1472 //
1473
1474 { .mfb
1475       nop.m 999
1476 (p0)  fmpy.d.s0 FR_Input_X = FR_Input_X, f0
1477 (p0)  br.ret.sptk   b0 ;;
1478 }
1479
1480 L(LOG_64_negative):
1481
1482 { .mfi
1483       nop.m 999
1484 //
1485 //    Deal with x < 0 in a special way
1486 //
1487 (p0)  frcpa.s0 FR_Output_X_tmp, p8 =  f0, f0
1488 //
1489 //    Deal with x < 0 in a special way - raise
1490 //    invalid and produce QNaN indefinite.
1491 //
1492 (p0)  mov   GR_Parameter_TAG = 141
1493 }
1494
1495 .endp log1p#
1496 ASM_SIZE_DIRECTIVE(log1p)
1497
1498 .proc __libm_error_region
1499 __libm_error_region:
1500 L(LOG_ERROR_Support):
1501 .prologue
1502
1503 // (1)
1504 { .mfi
1505         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
1506         nop.f 0
1507 .save   ar.pfs,GR_SAVE_PFS
1508         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
1509 }
1510 { .mfi
1511 .fframe 64
1512         add sp=-64,sp                          // Create new stack
1513         nop.f 0
1514         mov GR_SAVE_GP=gp                      // Save gp
1515 };;
1516
1517
1518 // (2)
1519 { .mmi
1520         stfd [GR_Parameter_Y] = f0,16         // STORE Parameter 2 on stack
1521         add GR_Parameter_X = 16,sp            // Parameter 1 address
1522 .save   b0, GR_SAVE_B0
1523         mov GR_SAVE_B0=b0                     // Save b0
1524 };;
1525
1526 .body
1527 // (3)
1528 { .mib
1529         stfd [GR_Parameter_X] =FR_Input_X               // STORE Parameter 1 on stack
1530         add   GR_Parameter_RESULT = 0,GR_Parameter_Y    // Parameter 3 address
1531         nop.b 0
1532 }
1533 { .mib
1534         stfd [GR_Parameter_Y] = FR_Output_X_tmp         // STORE Parameter 3 on stack
1535         add   GR_Parameter_Y = -16,GR_Parameter_Y
1536         br.call.sptk b0=__libm_error_support#           // Call error handling function
1537 };;
1538 { .mmi
1539         nop.m 0
1540         nop.m 0
1541         add   GR_Parameter_RESULT = 48,sp
1542 };;
1543
1544 // (4)
1545 { .mmi
1546         ldfd  FR_Input_X = [GR_Parameter_RESULT]       // Get return result off stack
1547 .restore sp
1548         add   sp = 64,sp                       // Restore stack pointer
1549         mov   b0 = GR_SAVE_B0                  // Restore return address
1550 };;
1551 { .mib
1552         mov   gp = GR_SAVE_GP                  // Restore gp
1553         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
1554         br.ret.sptk   b0
1555 };;
1556
1557 .endp __libm_error_region
1558 ASM_SIZE_DIRECTIVE(__libm_error_region)
1559
1560 .proc __libm_LOG_main
1561 __libm_LOG_main:
1562 L(LOG_main):
1563
1564 //
1565 //    kernel_log_64 computes ln(X + E)
1566 //
1567
1568 { .mfi
1569         nop.m 999
1570 (p7)  fadd.d.s0 FR_Input_X = FR_Y_lo,FR_Y_hi
1571         nop.i 999
1572 }
1573
1574 { .mmi
1575         nop.m 999
1576         nop.m 999
1577 (p14) addl GR_Table_Base = @ltoff(Constants_1_by_LN10#),gp ;;
1578 }
1579
1580 { .mmi
1581       nop.m 999
1582 (p14) ld8    GR_Table_Base = [GR_Table_Base]
1583       nop.i 999
1584 };;
1585
1586 { .mmi
1587 (p14) ldfe FR_1LN10_hi = [GR_Table_Base],16 ;;
1588 (p14) ldfe FR_1LN10_lo = [GR_Table_Base]
1589         nop.i 999 ;;
1590 }
1591
1592 { .mfi
1593         nop.m 999
1594 (p14) fmpy.s1 FR_Output_X_tmp = FR_Y_lo,FR_1LN10_hi
1595         nop.i 999 ;;
1596 }
1597
1598 { .mfi
1599         nop.m 999
1600 (p14) fma.s1  FR_Output_X_tmp = FR_Y_hi,FR_1LN10_lo,FR_Output_X_tmp
1601         nop.i 999 ;;
1602 }
1603
1604 { .mfb
1605         nop.m 999
1606 (p14) fma.d.s0 FR_Input_X = FR_Y_hi,FR_1LN10_hi,FR_Output_X_tmp
1607 (p0)  br.ret.sptk   b0 ;;
1608 }
1609 .endp __libm_LOG_main
1610 ASM_SIZE_DIRECTIVE(__libm_LOG_main)
1611
1612
1613 .type   __libm_error_support#,@function
1614 .global __libm_error_support#