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