| blob | 4bd92c3b267cda0f53d626d94dd266406a73e3b1 |
1 .file "libm_lgammaf.s"
4 // Copyright (c) 2002 - 2005, Intel Corporation
5 // All rights reserved.
6 //
7 // Contributed 2002 by the Intel Numerics Group, 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.
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://www.intel.com/software/products/opensource/libraries/num.htm.
39 //
40 //*********************************************************************
41 //
42 // History:
43 // 01/10/02 Initial version
44 // 01/25/02 Corrected parameter store, load, and tag for __libm_error_support
45 // 02/01/02 Added support of SIGN(GAMMA(x)) calculation
46 // 05/20/02 Cleaned up namespace and sf0 syntax
47 // 09/16/02 Improved accuracy on intervals reduced to [1;1.25]
48 // 10/21/02 Now it returns SIGN(GAMMA(x))=-1 for negative zero
49 // 02/10/03 Reordered header: .section, .global, .proc, .align
50 // 07/22/03 Reformatted some data tables
51 // 03/31/05 Reformatted delimiters between data tables
52 //
53 //*********************************************************************
54 //
55 //*********************************************************************
56 //
57 // Function: __libm_lgammaf(float x, int* signgam, int szsigngam)
58 // computes the principle value of the logarithm of the GAMMA function
59 // of x. Signum of GAMMA(x) is stored to memory starting at the address
60 // specified by the signgam.
61 //
62 //*********************************************************************
63 //
64 // Resources Used:
65 //
66 // Floating-Point Registers: f6-f15
67 // f32-f97
68 //
69 // General Purpose Registers:
70 // r8-r11
71 // r14-r30
72 // r32-r36
73 // r37-r40 (Used to pass arguments to error handling routine)
74 //
75 // Predicate Registers: p6-p15
76 //
77 //*********************************************************************
78 //
79 // IEEE Special Conditions:
80 //
81 // lgamma(+inf) = +inf
82 // lgamma(-inf) = +inf
83 // lgamma(+/-0) = +inf
84 // lgamma(x<0, x - integer) = +inf
85 // lgamma(SNaN) = QNaN
86 // lgamma(QNaN) = QNaN
87 //
88 //*********************************************************************
89 //
90 // Overview
91 //
92 // The method consists of three cases.
93 //
94 // If 2^13 <= x < OVERFLOW_BOUNDARY use case lgammaf_pstirling;
95 // else if 1 < x < 2^13 use case lgammaf_regular;
96 // else if -9 < x < 1 use case lgammaf_negrecursion;
97 // else if -2^13 < x < -9 use case lgammaf_negpoly;
98 // else if x < -2^13 use case lgammaf_negstirling;
99 // else if x is close to negative
100 // roots of ln(GAMMA(x)) use case lgammaf_negroots;
101 //
102 //
103 // Case 2^13 <= x < OVERFLOW_BOUNDARY
104 // ----------------------------------
105 // Here we use algorithm based on the Stirling formula:
106 // ln(GAMMA(x)) = ln(sqrt(2*Pi)) + (x-0.5)*ln(x) - x
107 //
108 // Case 1 < x < 2^13
109 // -----------------
110 // To calculate ln(GAMMA(x)) for such arguments we use polynomial
111 // approximation on following intervals: [1.0; 1.25), [1.25; 1.5),
112 // [1.5, 1.75), [1.75; 2), [2; 4), [2^i; 2^(i+1)), i=1..8
113 //
114 // Following variants of approximation and argument reduction are used:
115 // 1. [1.0; 1.25)
116 // ln(GAMMA(x)) ~ (x-1.0)*P7(x)
117 //
118 // 2. [1.25; 1.5)
119 // ln(GAMMA(x)) ~ ln(GAMMA(x0))+(x-x0)*P8(x-x0),
120 // where x0 - point of local minimum on [1;2] rounded to nearest double
121 // precision number.
122 //
123 // 3. [1.5; 1.75)
124 // ln(GAMMA(x)) ~ P8(x)
125 //
126 // 4. [1.75; 2.0)
127 // ln(GAMMA(x)) ~ (x-2)*P7(x)
128 //
129 // 5. [2; 4)
130 // ln(GAMMA(x)) ~ (x-2)*P10(x)
131 //
132 // 6. [2^i; 2^(i+1)), i=2..8
133 // ln(GAMMA(x)) ~ P10((x-2^i)/2^i)
134 //
135 // Case -9 < x < 1
136 // ---------------
137 // Here we use the recursive formula:
138 // ln(GAMMA(x)) = ln(GAMMA(x+1)) - ln(x)
139 //
140 // Using this formula we reduce argument to base interval [1.0; 2.0]
141 //
142 // Case -2^13 < x < -9
143 // --------------------
144 // Here we use the formula:
145 // ln(GAMMA(x)) = ln(Pi/(|x|*GAMMA(|x|)*sin(Pi*|x|))) =
146 // = -ln(|x|) - ln((GAMMA(|x|)) - ln(sin(Pi*r)/(Pi*r)) - ln(|r|)
147 // where r = x - rounded_to_nearest(x), i.e |r| <= 0.5 and
148 // ln(sin(Pi*r)/(Pi*r)) is approximated by 8-degree polynomial of r^2
149 //
150 // Case x < -2^13
151 // --------------
152 // Here we use algorithm based on the Stirling formula:
153 // ln(GAMMA(x)) = -ln(sqrt(2*Pi)) + (|x|-0.5)ln(x) - |x| -
154 // - ln(sin(Pi*r)/(Pi*r)) - ln(|r|)
155 // where r = x - rounded_to_nearest(x).
156 //
157 // Neighbourhoods of negative roots
158 // --------------------------------
159 // Here we use polynomial approximation
160 // ln(GAMMA(x-x0)) = ln(GAMMA(x0)) + (x-x0)*P14(x-x0),
161 // where x0 is a root of ln(GAMMA(x)) rounded to nearest double
162 // precision number.
163 //
164 //
165 // Claculation of logarithm
166 // ------------------------
167 // Consider x = 2^N * xf so
168 // ln(x) = ln(frcpa(x)*x/frcpa(x))
169 // = ln(1/frcpa(x)) + ln(frcpa(x)*x)
170 //
171 // frcpa(x) = 2^(-N) * frcpa(xf)
172 //
173 // ln(1/frcpa(x)) = -ln(2^(-N)) - ln(frcpa(xf))
174 // = N*ln(2) - ln(frcpa(xf))
175 // = N*ln(2) + ln(1/frcpa(xf))
176 //
177 // ln(x) = ln(1/frcpa(x)) + ln(frcpa(x)*x) =
178 // = N*ln(2) + ln(1/frcpa(xf)) + ln(frcpa(x)*x)
179 // = N*ln(2) + T + ln(frcpa(x)*x)
180 //
181 // Let r = 1 - frcpa(x)*x, note that r is quite small by
182 // absolute value so
183 //
184 // ln(x) = N*ln(2) + T + ln(1+r) ~ N*ln(2) + T + Series(r),
185 // where T - is precomputed tabular value,
186 // Series(r) = (P3*r + P2)*r^2 + (P1*r + 1)
187 //
188 //*********************************************************************
190 GR_TAG = r8
191 GR_ad_Data = r8
192 GR_ad_Co = r9
193 GR_ad_SignGam = r10
194 GR_ad_Ce = r10
195 GR_SignExp = r11
197 GR_ad_C650 = r14
198 GR_ad_RootCo = r14
199 GR_ad_C0 = r15
200 GR_Dx = r15
201 GR_Ind = r16
202 GR_Offs = r17
203 GR_IntNum = r17
204 GR_ExpBias = r18
205 GR_ExpMask = r19
206 GR_Ind4T = r20
207 GR_RootInd = r20
208 GR_Sig = r21
209 GR_Exp = r22
210 GR_PureExp = r23
211 GR_ad_C43 = r24
212 GR_StirlBound = r25
213 GR_ad_T = r25
214 GR_IndX8 = r25
215 GR_Neg2 = r25
216 GR_2xDx = r25
217 GR_SingBound = r26
218 GR_IndX2 = r26
219 GR_Neg4 = r26
220 GR_ad_RootCe = r26
221 GR_Arg = r27
222 GR_ExpOf2 = r28
223 GR_fff7 = r28
224 GR_Root = r28
225 GR_ReqBound = r28
226 GR_N = r29
227 GR_ad_Root = r30
228 GR_ad_OvfBound = r30
229 GR_SignOfGamma = r31
231 GR_SAVE_B0 = r33
232 GR_SAVE_PFS = r34
233 GR_SAVE_GP = r35
234 GR_SAVE_SP = r36
236 GR_Parameter_X = r37
237 GR_Parameter_Y = r38
238 GR_Parameter_RESULT = r39
239 GR_Parameter_TAG = r40
241 //*********************************************************************
243 FR_X = f10
244 FR_Y = f1 // lgammaf is single argument function
245 FR_RESULT = f8
247 FR_x = f6
248 FR_x2 = f7
250 FR_x3 = f9
251 FR_x4 = f10
252 FR_xm2 = f11
253 FR_w = f11
254 FR_w2 = f12
255 FR_Q32 = f13
256 FR_Q10 = f14
257 FR_InvX = f15
259 FR_NormX = f32
261 FR_A0 = f33
262 FR_A1 = f34
263 FR_A2 = f35
264 FR_A3 = f36
265 FR_A4 = f37
266 FR_A5 = f38
267 FR_A6 = f39
268 FR_A7 = f40
269 FR_A8 = f41
270 FR_A9 = f42
271 FR_A10 = f43
273 FR_int_N = f44
274 FR_P3 = f45
275 FR_P2 = f46
276 FR_P1 = f47
277 FR_LocalMin = f48
278 FR_Ln2 = f49
279 FR_05 = f50
280 FR_LnSqrt2Pi = f51
281 FR_3 = f52
282 FR_r = f53
283 FR_r2 = f54
284 FR_T = f55
285 FR_N = f56
286 FR_xm05 = f57
287 FR_int_Ln = f58
288 FR_P32 = f59
289 FR_P10 = f60
291 FR_Xf = f61
292 FR_InvXf = f62
293 FR_rf = f63
294 FR_rf2 = f64
295 FR_Tf = f65
296 FR_Nf = f66
297 FR_xm05f = f67
298 FR_P32f = f68
299 FR_P10f = f69
300 FR_Lnf = f70
301 FR_Xf2 = f71
302 FR_Xf4 = f72
303 FR_Xf8 = f73
304 FR_Ln = f74
305 FR_xx = f75
306 FR_Root = f75
307 FR_Req = f76
308 FR_1pXf = f77
310 FR_S16 = f78
311 FR_R3 = f78
312 FR_S14 = f79
313 FR_R2 = f79
314 FR_S12 = f80
315 FR_R1 = f80
316 FR_S10 = f81
317 FR_R0 = f81
318 FR_S8 = f82
319 FR_rx = f82
320 FR_S6 = f83
321 FR_rx2 = f84
322 FR_S4 = f84
323 FR_S2 = f85
325 FR_Xp1 = f86
326 FR_Xp2 = f87
327 FR_Xp3 = f88
328 FR_Xp4 = f89
329 FR_Xp5 = f90
330 FR_Xp6 = f91
331 FR_Xp7 = f92
332 FR_Xp8 = f93
333 FR_OverflowBound = f93
335 FR_2 = f94
336 FR_tmp = f95
337 FR_int_Ntrunc = f96
338 FR_Ntrunc = f97
340 //*********************************************************************
342 RODATA
343 .align 32
344 LOCAL_OBJECT_START(lgammaf_data)
345 log_table_1:
346 data8 0xbfd0001008f39d59 // P3
347 data8 0x3fd5556073e0c45a // P2
348 data8 0x3fe62e42fefa39ef // ln(2)
349 data8 0x3fe0000000000000 // 0.5
350 //
351 data8 0x3F60040155D5889E //ln(1/frcpa(1+ 0/256)
352 data8 0x3F78121214586B54 //ln(1/frcpa(1+ 1/256)
353 data8 0x3F841929F96832F0 //ln(1/frcpa(1+ 2/256)
354 data8 0x3F8C317384C75F06 //ln(1/frcpa(1+ 3/256)
355 data8 0x3F91A6B91AC73386 //ln(1/frcpa(1+ 4/256)
356 data8 0x3F95BA9A5D9AC039 //ln(1/frcpa(1+ 5/256)
357 data8 0x3F99D2A8074325F4 //ln(1/frcpa(1+ 6/256)
358 data8 0x3F9D6B2725979802 //ln(1/frcpa(1+ 7/256)
359 data8 0x3FA0C58FA19DFAAA //ln(1/frcpa(1+ 8/256)
360 data8 0x3FA2954C78CBCE1B //ln(1/frcpa(1+ 9/256)
361 data8 0x3FA4A94D2DA96C56 //ln(1/frcpa(1+ 10/256)
362 data8 0x3FA67C94F2D4BB58 //ln(1/frcpa(1+ 11/256)
363 data8 0x3FA85188B630F068 //ln(1/frcpa(1+ 12/256)
364 data8 0x3FAA6B8ABE73AF4C //ln(1/frcpa(1+ 13/256)
365 data8 0x3FAC441E06F72A9E //ln(1/frcpa(1+ 14/256)
366 data8 0x3FAE1E6713606D07 //ln(1/frcpa(1+ 15/256)
367 data8 0x3FAFFA6911AB9301 //ln(1/frcpa(1+ 16/256)
368 data8 0x3FB0EC139C5DA601 //ln(1/frcpa(1+ 17/256)
369 data8 0x3FB1DBD2643D190B //ln(1/frcpa(1+ 18/256)
370 data8 0x3FB2CC7284FE5F1C //ln(1/frcpa(1+ 19/256)
371 data8 0x3FB3BDF5A7D1EE64 //ln(1/frcpa(1+ 20/256)
372 data8 0x3FB4B05D7AA012E0 //ln(1/frcpa(1+ 21/256)
373 data8 0x3FB580DB7CEB5702 //ln(1/frcpa(1+ 22/256)
374 data8 0x3FB674F089365A7A //ln(1/frcpa(1+ 23/256)
375 data8 0x3FB769EF2C6B568D //ln(1/frcpa(1+ 24/256)
376 data8 0x3FB85FD927506A48 //ln(1/frcpa(1+ 25/256)
377 data8 0x3FB9335E5D594989 //ln(1/frcpa(1+ 26/256)
378 data8 0x3FBA2B0220C8E5F5 //ln(1/frcpa(1+ 27/256)
379 data8 0x3FBB0004AC1A86AC //ln(1/frcpa(1+ 28/256)
380 data8 0x3FBBF968769FCA11 //ln(1/frcpa(1+ 29/256)
381 data8 0x3FBCCFEDBFEE13A8 //ln(1/frcpa(1+ 30/256)
382 data8 0x3FBDA727638446A2 //ln(1/frcpa(1+ 31/256)
383 data8 0x3FBEA3257FE10F7A //ln(1/frcpa(1+ 32/256)
384 data8 0x3FBF7BE9FEDBFDE6 //ln(1/frcpa(1+ 33/256)
385 data8 0x3FC02AB352FF25F4 //ln(1/frcpa(1+ 34/256)
386 data8 0x3FC097CE579D204D //ln(1/frcpa(1+ 35/256)
387 data8 0x3FC1178E8227E47C //ln(1/frcpa(1+ 36/256)
388 data8 0x3FC185747DBECF34 //ln(1/frcpa(1+ 37/256)
389 data8 0x3FC1F3B925F25D41 //ln(1/frcpa(1+ 38/256)
390 data8 0x3FC2625D1E6DDF57 //ln(1/frcpa(1+ 39/256)
391 data8 0x3FC2D1610C86813A //ln(1/frcpa(1+ 40/256)
392 data8 0x3FC340C59741142E //ln(1/frcpa(1+ 41/256)
393 data8 0x3FC3B08B6757F2A9 //ln(1/frcpa(1+ 42/256)
394 data8 0x3FC40DFB08378003 //ln(1/frcpa(1+ 43/256)
395 data8 0x3FC47E74E8CA5F7C //ln(1/frcpa(1+ 44/256)
396 data8 0x3FC4EF51F6466DE4 //ln(1/frcpa(1+ 45/256)
397 data8 0x3FC56092E02BA516 //ln(1/frcpa(1+ 46/256)
398 data8 0x3FC5D23857CD74D5 //ln(1/frcpa(1+ 47/256)
399 data8 0x3FC6313A37335D76 //ln(1/frcpa(1+ 48/256)
400 data8 0x3FC6A399DABBD383 //ln(1/frcpa(1+ 49/256)
401 data8 0x3FC70337DD3CE41B //ln(1/frcpa(1+ 50/256)
402 data8 0x3FC77654128F6127 //ln(1/frcpa(1+ 51/256)
403 data8 0x3FC7E9D82A0B022D //ln(1/frcpa(1+ 52/256)
404 data8 0x3FC84A6B759F512F //ln(1/frcpa(1+ 53/256)
405 data8 0x3FC8AB47D5F5A310 //ln(1/frcpa(1+ 54/256)
406 data8 0x3FC91FE49096581B //ln(1/frcpa(1+ 55/256)
407 data8 0x3FC981634011AA75 //ln(1/frcpa(1+ 56/256)
408 data8 0x3FC9F6C407089664 //ln(1/frcpa(1+ 57/256)
409 data8 0x3FCA58E729348F43 //ln(1/frcpa(1+ 58/256)
410 data8 0x3FCABB55C31693AD //ln(1/frcpa(1+ 59/256)
411 data8 0x3FCB1E104919EFD0 //ln(1/frcpa(1+ 60/256)
412 data8 0x3FCB94EE93E367CB //ln(1/frcpa(1+ 61/256)
413 data8 0x3FCBF851C067555F //ln(1/frcpa(1+ 62/256)
414 data8 0x3FCC5C0254BF23A6 //ln(1/frcpa(1+ 63/256)
415 data8 0x3FCCC000C9DB3C52 //ln(1/frcpa(1+ 64/256)
416 data8 0x3FCD244D99C85674 //ln(1/frcpa(1+ 65/256)
417 data8 0x3FCD88E93FB2F450 //ln(1/frcpa(1+ 66/256)
418 data8 0x3FCDEDD437EAEF01 //ln(1/frcpa(1+ 67/256)
419 data8 0x3FCE530EFFE71012 //ln(1/frcpa(1+ 68/256)
420 data8 0x3FCEB89A1648B971 //ln(1/frcpa(1+ 69/256)
421 data8 0x3FCF1E75FADF9BDE //ln(1/frcpa(1+ 70/256)
422 data8 0x3FCF84A32EAD7C35 //ln(1/frcpa(1+ 71/256)
423 data8 0x3FCFEB2233EA07CD //ln(1/frcpa(1+ 72/256)
424 data8 0x3FD028F9C7035C1C //ln(1/frcpa(1+ 73/256)
425 data8 0x3FD05C8BE0D9635A //ln(1/frcpa(1+ 74/256)
426 data8 0x3FD085EB8F8AE797 //ln(1/frcpa(1+ 75/256)
427 data8 0x3FD0B9C8E32D1911 //ln(1/frcpa(1+ 76/256)
428 data8 0x3FD0EDD060B78081 //ln(1/frcpa(1+ 77/256)
429 data8 0x3FD122024CF0063F //ln(1/frcpa(1+ 78/256)
430 data8 0x3FD14BE2927AECD4 //ln(1/frcpa(1+ 79/256)
431 data8 0x3FD180618EF18ADF //ln(1/frcpa(1+ 80/256)
432 data8 0x3FD1B50BBE2FC63B //ln(1/frcpa(1+ 81/256)
433 data8 0x3FD1DF4CC7CF242D //ln(1/frcpa(1+ 82/256)
434 data8 0x3FD214456D0EB8D4 //ln(1/frcpa(1+ 83/256)
435 data8 0x3FD23EC5991EBA49 //ln(1/frcpa(1+ 84/256)
436 data8 0x3FD2740D9F870AFB //ln(1/frcpa(1+ 85/256)
437 data8 0x3FD29ECDABCDFA04 //ln(1/frcpa(1+ 86/256)
438 data8 0x3FD2D46602ADCCEE //ln(1/frcpa(1+ 87/256)
439 data8 0x3FD2FF66B04EA9D4 //ln(1/frcpa(1+ 88/256)
440 data8 0x3FD335504B355A37 //ln(1/frcpa(1+ 89/256)
441 data8 0x3FD360925EC44F5D //ln(1/frcpa(1+ 90/256)
442 data8 0x3FD38BF1C3337E75 //ln(1/frcpa(1+ 91/256)
443 data8 0x3FD3C25277333184 //ln(1/frcpa(1+ 92/256)
444 data8 0x3FD3EDF463C1683E //ln(1/frcpa(1+ 93/256)
445 data8 0x3FD419B423D5E8C7 //ln(1/frcpa(1+ 94/256)
446 data8 0x3FD44591E0539F49 //ln(1/frcpa(1+ 95/256)
447 data8 0x3FD47C9175B6F0AD //ln(1/frcpa(1+ 96/256)
448 data8 0x3FD4A8B341552B09 //ln(1/frcpa(1+ 97/256)
449 data8 0x3FD4D4F3908901A0 //ln(1/frcpa(1+ 98/256)
450 data8 0x3FD501528DA1F968 //ln(1/frcpa(1+ 99/256)
451 data8 0x3FD52DD06347D4F6 //ln(1/frcpa(1+ 100/256)
452 data8 0x3FD55A6D3C7B8A8A //ln(1/frcpa(1+ 101/256)
453 data8 0x3FD5925D2B112A59 //ln(1/frcpa(1+ 102/256)
454 data8 0x3FD5BF406B543DB2 //ln(1/frcpa(1+ 103/256)
455 data8 0x3FD5EC433D5C35AE //ln(1/frcpa(1+ 104/256)
456 data8 0x3FD61965CDB02C1F //ln(1/frcpa(1+ 105/256)
457 data8 0x3FD646A84935B2A2 //ln(1/frcpa(1+ 106/256)
458 data8 0x3FD6740ADD31DE94 //ln(1/frcpa(1+ 107/256)
459 data8 0x3FD6A18DB74A58C5 //ln(1/frcpa(1+ 108/256)
460 data8 0x3FD6CF31058670EC //ln(1/frcpa(1+ 109/256)
461 data8 0x3FD6F180E852F0BA //ln(1/frcpa(1+ 110/256)
462 data8 0x3FD71F5D71B894F0 //ln(1/frcpa(1+ 111/256)
463 data8 0x3FD74D5AEFD66D5C //ln(1/frcpa(1+ 112/256)
464 data8 0x3FD77B79922BD37E //ln(1/frcpa(1+ 113/256)
465 data8 0x3FD7A9B9889F19E2 //ln(1/frcpa(1+ 114/256)
466 data8 0x3FD7D81B037EB6A6 //ln(1/frcpa(1+ 115/256)
467 data8 0x3FD8069E33827231 //ln(1/frcpa(1+ 116/256)
468 data8 0x3FD82996D3EF8BCB //ln(1/frcpa(1+ 117/256)
469 data8 0x3FD85855776DCBFB //ln(1/frcpa(1+ 118/256)
470 data8 0x3FD8873658327CCF //ln(1/frcpa(1+ 119/256)
471 data8 0x3FD8AA75973AB8CF //ln(1/frcpa(1+ 120/256)
472 data8 0x3FD8D992DC8824E5 //ln(1/frcpa(1+ 121/256)
473 data8 0x3FD908D2EA7D9512 //ln(1/frcpa(1+ 122/256)
474 data8 0x3FD92C59E79C0E56 //ln(1/frcpa(1+ 123/256)
475 data8 0x3FD95BD750EE3ED3 //ln(1/frcpa(1+ 124/256)
476 data8 0x3FD98B7811A3EE5B //ln(1/frcpa(1+ 125/256)
477 data8 0x3FD9AF47F33D406C //ln(1/frcpa(1+ 126/256)
478 data8 0x3FD9DF270C1914A8 //ln(1/frcpa(1+ 127/256)
479 data8 0x3FDA0325ED14FDA4 //ln(1/frcpa(1+ 128/256)
480 data8 0x3FDA33440224FA79 //ln(1/frcpa(1+ 129/256)
481 data8 0x3FDA57725E80C383 //ln(1/frcpa(1+ 130/256)
482 data8 0x3FDA87D0165DD199 //ln(1/frcpa(1+ 131/256)
483 data8 0x3FDAAC2E6C03F896 //ln(1/frcpa(1+ 132/256)
484 data8 0x3FDADCCC6FDF6A81 //ln(1/frcpa(1+ 133/256)
485 data8 0x3FDB015B3EB1E790 //ln(1/frcpa(1+ 134/256)
486 data8 0x3FDB323A3A635948 //ln(1/frcpa(1+ 135/256)
487 data8 0x3FDB56FA04462909 //ln(1/frcpa(1+ 136/256)
488 data8 0x3FDB881AA659BC93 //ln(1/frcpa(1+ 137/256)
489 data8 0x3FDBAD0BEF3DB165 //ln(1/frcpa(1+ 138/256)
490 data8 0x3FDBD21297781C2F //ln(1/frcpa(1+ 139/256)
491 data8 0x3FDC039236F08819 //ln(1/frcpa(1+ 140/256)
492 data8 0x3FDC28CB1E4D32FD //ln(1/frcpa(1+ 141/256)
493 data8 0x3FDC4E19B84723C2 //ln(1/frcpa(1+ 142/256)
494 data8 0x3FDC7FF9C74554C9 //ln(1/frcpa(1+ 143/256)
495 data8 0x3FDCA57B64E9DB05 //ln(1/frcpa(1+ 144/256)
496 data8 0x3FDCCB130A5CEBB0 //ln(1/frcpa(1+ 145/256)
497 data8 0x3FDCF0C0D18F326F //ln(1/frcpa(1+ 146/256)
498 data8 0x3FDD232075B5A201 //ln(1/frcpa(1+ 147/256)
499 data8 0x3FDD490246DEFA6B //ln(1/frcpa(1+ 148/256)
500 data8 0x3FDD6EFA918D25CD //ln(1/frcpa(1+ 149/256)
501 data8 0x3FDD9509707AE52F //ln(1/frcpa(1+ 150/256)
502 data8 0x3FDDBB2EFE92C554 //ln(1/frcpa(1+ 151/256)
503 data8 0x3FDDEE2F3445E4AF //ln(1/frcpa(1+ 152/256)
504 data8 0x3FDE148A1A2726CE //ln(1/frcpa(1+ 153/256)
505 data8 0x3FDE3AFC0A49FF40 //ln(1/frcpa(1+ 154/256)
506 data8 0x3FDE6185206D516E //ln(1/frcpa(1+ 155/256)
507 data8 0x3FDE882578823D52 //ln(1/frcpa(1+ 156/256)
508 data8 0x3FDEAEDD2EAC990C //ln(1/frcpa(1+ 157/256)
509 data8 0x3FDED5AC5F436BE3 //ln(1/frcpa(1+ 158/256)
510 data8 0x3FDEFC9326D16AB9 //ln(1/frcpa(1+ 159/256)
511 data8 0x3FDF2391A2157600 //ln(1/frcpa(1+ 160/256)
512 data8 0x3FDF4AA7EE03192D //ln(1/frcpa(1+ 161/256)
513 data8 0x3FDF71D627C30BB0 //ln(1/frcpa(1+ 162/256)
514 data8 0x3FDF991C6CB3B379 //ln(1/frcpa(1+ 163/256)
515 data8 0x3FDFC07ADA69A910 //ln(1/frcpa(1+ 164/256)
516 data8 0x3FDFE7F18EB03D3E //ln(1/frcpa(1+ 165/256)
517 data8 0x3FE007C053C5002E //ln(1/frcpa(1+ 166/256)
518 data8 0x3FE01B942198A5A1 //ln(1/frcpa(1+ 167/256)
519 data8 0x3FE02F74400C64EB //ln(1/frcpa(1+ 168/256)
520 data8 0x3FE04360BE7603AD //ln(1/frcpa(1+ 169/256)
521 data8 0x3FE05759AC47FE34 //ln(1/frcpa(1+ 170/256)
522 data8 0x3FE06B5F1911CF52 //ln(1/frcpa(1+ 171/256)
523 data8 0x3FE078BF0533C568 //ln(1/frcpa(1+ 172/256)
524 data8 0x3FE08CD9687E7B0E //ln(1/frcpa(1+ 173/256)
525 data8 0x3FE0A10074CF9019 //ln(1/frcpa(1+ 174/256)
526 data8 0x3FE0B5343A234477 //ln(1/frcpa(1+ 175/256)
527 data8 0x3FE0C974C89431CE //ln(1/frcpa(1+ 176/256)
528 data8 0x3FE0DDC2305B9886 //ln(1/frcpa(1+ 177/256)
529 data8 0x3FE0EB524BAFC918 //ln(1/frcpa(1+ 178/256)
530 data8 0x3FE0FFB54213A476 //ln(1/frcpa(1+ 179/256)
531 data8 0x3FE114253DA97D9F //ln(1/frcpa(1+ 180/256)
532 data8 0x3FE128A24F1D9AFF //ln(1/frcpa(1+ 181/256)
533 data8 0x3FE1365252BF0865 //ln(1/frcpa(1+ 182/256)
534 data8 0x3FE14AE558B4A92D //ln(1/frcpa(1+ 183/256)
535 data8 0x3FE15F85A19C765B //ln(1/frcpa(1+ 184/256)
536 data8 0x3FE16D4D38C119FA //ln(1/frcpa(1+ 185/256)
537 data8 0x3FE18203C20DD133 //ln(1/frcpa(1+ 186/256)
538 data8 0x3FE196C7BC4B1F3B //ln(1/frcpa(1+ 187/256)
539 data8 0x3FE1A4A738B7A33C //ln(1/frcpa(1+ 188/256)
540 data8 0x3FE1B981C0C9653D //ln(1/frcpa(1+ 189/256)
541 data8 0x3FE1CE69E8BB106B //ln(1/frcpa(1+ 190/256)
542 data8 0x3FE1DC619DE06944 //ln(1/frcpa(1+ 191/256)
543 data8 0x3FE1F160A2AD0DA4 //ln(1/frcpa(1+ 192/256)
544 data8 0x3FE2066D7740737E //ln(1/frcpa(1+ 193/256)
545 data8 0x3FE2147DBA47A394 //ln(1/frcpa(1+ 194/256)
546 data8 0x3FE229A1BC5EBAC3 //ln(1/frcpa(1+ 195/256)
547 data8 0x3FE237C1841A502E //ln(1/frcpa(1+ 196/256)
548 data8 0x3FE24CFCE6F80D9A //ln(1/frcpa(1+ 197/256)
549 data8 0x3FE25B2C55CD5762 //ln(1/frcpa(1+ 198/256)
550 data8 0x3FE2707F4D5F7C41 //ln(1/frcpa(1+ 199/256)
551 data8 0x3FE285E0842CA384 //ln(1/frcpa(1+ 200/256)
552 data8 0x3FE294294708B773 //ln(1/frcpa(1+ 201/256)
553 data8 0x3FE2A9A2670AFF0C //ln(1/frcpa(1+ 202/256)
554 data8 0x3FE2B7FB2C8D1CC1 //ln(1/frcpa(1+ 203/256)
555 data8 0x3FE2C65A6395F5F5 //ln(1/frcpa(1+ 204/256)
556 data8 0x3FE2DBF557B0DF43 //ln(1/frcpa(1+ 205/256)
557 data8 0x3FE2EA64C3F97655 //ln(1/frcpa(1+ 206/256)
558 data8 0x3FE3001823684D73 //ln(1/frcpa(1+ 207/256)
559 data8 0x3FE30E97E9A8B5CD //ln(1/frcpa(1+ 208/256)
560 data8 0x3FE32463EBDD34EA //ln(1/frcpa(1+ 209/256)
561 data8 0x3FE332F4314AD796 //ln(1/frcpa(1+ 210/256)
562 data8 0x3FE348D90E7464D0 //ln(1/frcpa(1+ 211/256)
563 data8 0x3FE35779F8C43D6E //ln(1/frcpa(1+ 212/256)
564 data8 0x3FE36621961A6A99 //ln(1/frcpa(1+ 213/256)
565 data8 0x3FE37C299F3C366A //ln(1/frcpa(1+ 214/256)
566 data8 0x3FE38AE2171976E7 //ln(1/frcpa(1+ 215/256)
567 data8 0x3FE399A157A603E7 //ln(1/frcpa(1+ 216/256)
568 data8 0x3FE3AFCCFE77B9D1 //ln(1/frcpa(1+ 217/256)
569 data8 0x3FE3BE9D503533B5 //ln(1/frcpa(1+ 218/256)
570 data8 0x3FE3CD7480B4A8A3 //ln(1/frcpa(1+ 219/256)
571 data8 0x3FE3E3C43918F76C //ln(1/frcpa(1+ 220/256)
572 data8 0x3FE3F2ACB27ED6C7 //ln(1/frcpa(1+ 221/256)
573 data8 0x3FE4019C2125CA93 //ln(1/frcpa(1+ 222/256)
574 data8 0x3FE4181061389722 //ln(1/frcpa(1+ 223/256)
575 data8 0x3FE42711518DF545 //ln(1/frcpa(1+ 224/256)
576 data8 0x3FE436194E12B6BF //ln(1/frcpa(1+ 225/256)
577 data8 0x3FE445285D68EA69 //ln(1/frcpa(1+ 226/256)
578 data8 0x3FE45BCC464C893A //ln(1/frcpa(1+ 227/256)
579 data8 0x3FE46AED21F117FC //ln(1/frcpa(1+ 228/256)
580 data8 0x3FE47A1527E8A2D3 //ln(1/frcpa(1+ 229/256)
581 data8 0x3FE489445EFFFCCC //ln(1/frcpa(1+ 230/256)
582 data8 0x3FE4A018BCB69835 //ln(1/frcpa(1+ 231/256)
583 data8 0x3FE4AF5A0C9D65D7 //ln(1/frcpa(1+ 232/256)
584 data8 0x3FE4BEA2A5BDBE87 //ln(1/frcpa(1+ 233/256)
585 data8 0x3FE4CDF28F10AC46 //ln(1/frcpa(1+ 234/256)
586 data8 0x3FE4DD49CF994058 //ln(1/frcpa(1+ 235/256)
587 data8 0x3FE4ECA86E64A684 //ln(1/frcpa(1+ 236/256)
588 data8 0x3FE503C43CD8EB68 //ln(1/frcpa(1+ 237/256)
589 data8 0x3FE513356667FC57 //ln(1/frcpa(1+ 238/256)
590 data8 0x3FE522AE0738A3D8 //ln(1/frcpa(1+ 239/256)
591 data8 0x3FE5322E26867857 //ln(1/frcpa(1+ 240/256)
592 data8 0x3FE541B5CB979809 //ln(1/frcpa(1+ 241/256)
593 data8 0x3FE55144FDBCBD62 //ln(1/frcpa(1+ 242/256)
594 data8 0x3FE560DBC45153C7 //ln(1/frcpa(1+ 243/256)
595 data8 0x3FE5707A26BB8C66 //ln(1/frcpa(1+ 244/256)
596 data8 0x3FE587F60ED5B900 //ln(1/frcpa(1+ 245/256)
597 data8 0x3FE597A7977C8F31 //ln(1/frcpa(1+ 246/256)
598 data8 0x3FE5A760D634BB8B //ln(1/frcpa(1+ 247/256)
599 data8 0x3FE5B721D295F10F //ln(1/frcpa(1+ 248/256)
600 data8 0x3FE5C6EA94431EF9 //ln(1/frcpa(1+ 249/256)
601 data8 0x3FE5D6BB22EA86F6 //ln(1/frcpa(1+ 250/256)
602 data8 0x3FE5E6938645D390 //ln(1/frcpa(1+ 251/256)
603 data8 0x3FE5F673C61A2ED2 //ln(1/frcpa(1+ 252/256)
604 data8 0x3FE6065BEA385926 //ln(1/frcpa(1+ 253/256)
605 data8 0x3FE6164BFA7CC06B //ln(1/frcpa(1+ 254/256)
606 data8 0x3FE62643FECF9743 //ln(1/frcpa(1+ 255/256)
607 //
608 // [2;4)
609 data8 0xBEB2CC7A38B9355F,0x3F035F2D1833BF4C // A10,A9
610 data8 0xBFF51BAA7FD27785,0x3FFC9D5D5B6CDEFF // A2,A1
611 data8 0xBF421676F9CB46C7,0x3F7437F2FA1436C6 // A8,A7
612 data8 0xBFD7A7041DE592FE,0x3FE9F107FEE8BD29 // A4,A3
613 // [4;8)
614 data8 0x3F6BBBD68451C0CD,0xBF966EC3272A16F7 // A10,A9
615 data8 0x40022A24A39AD769,0x4014190EDF49C8C5 // A2,A1
616 data8 0x3FB130FD016EE241,0xBFC151B46E635248 // A8,A7
617 data8 0x3FDE8F611965B5FE,0xBFEB5110EB265E3D // A4,A3
618 // [8;16)
619 data8 0x3F736EF93508626A,0xBF9FE5DBADF58AF1 // A10,A9
620 data8 0x40110A9FC5192058,0x40302008A6F96B29 // A2,A1
621 data8 0x3FB8E74E0CE1E4B5,0xBFC9B5DA78873656 // A8,A7
622 data8 0x3FE99D0DF10022DC,0xBFF829C0388F9484 // A4,A3
623 // [16;32)
624 data8 0x3F7FFF9D6D7E9269,0xBFAA780A249AEDB1 // A10,A9
625 data8 0x402082A807AEA080,0x4045ED9868408013 // A2,A1
626 data8 0x3FC4E1E54C2F99B7,0xBFD5DE2D6FFF1490 // A8,A7
627 data8 0x3FF75FC89584AE87,0xC006B4BADD886CAE // A4,A3
628 // [32;64)
629 data8 0x3F8CE54375841A5F,0xBFB801ABCFFA1BE2 // A10,A9
630 data8 0x403040A8B1815BDA,0x405B99A917D24B7A // A2,A1
631 data8 0x3FD30CAB81BFFA03,0xBFE41AEF61ECF48B // A8,A7
632 data8 0x400650CC136BEC43,0xC016022046E8292B // A4,A3
633 // [64;128)
634 data8 0x3F9B69BD22CAA8B8,0xBFC6D48875B7A213 // A10,A9
635 data8 0x40402028CCAA2F6D,0x40709AACEB3CBE0F // A2,A1
636 data8 0x3FE22C6A5924761E,0xBFF342F5F224523D // A8,A7
637 data8 0x4015CD405CCA331F,0xC025AAD10482C769 // A4,A3
638 // [128;256)
639 data8 0x3FAAAD9CD0E40D06,0xBFD63FC8505D80CB // A10,A9
640 data8 0x40501008D56C2648,0x408364794B0F4376 // A2,A1
641 data8 0x3FF1BE0126E00284,0xC002D8E3F6F7F7CA // A8,A7
642 data8 0x40258C757E95D860,0xC0357FA8FD398011 // A4,A3
643 // [256;512)
644 data8 0x3FBA4DAC59D49FEB,0xBFE5F476D1C43A77 // A10,A9
645 data8 0x40600800D890C7C6,0x40962C42AAEC8EF0 // A2,A1
646 data8 0x40018680ECF19B89,0xC012A3EB96FB7BA4 // A8,A7
647 data8 0x40356C4CDD3B60F9,0xC0456A34BF18F440 // A4,A3
648 // [512;1024)
649 data8 0x3FCA1B54F6225A5A,0xBFF5CD67BA10E048 // A10,A9
650 data8 0x407003FED94C58C2,0x40A8F30B4ACBCD22 // A2,A1
651 data8 0x40116A135EB66D8C,0xC022891B1CED527E // A8,A7
652 data8 0x40455C4617FDD8BC,0xC0555F82729E59C4 // A4,A3
653 // [1024;2048)
654 data8 0x3FD9FFF9095C6EC9,0xC005B88CB25D76C9 // A10,A9
655 data8 0x408001FE58FA734D,0x40BBB953BAABB0F3 // A2,A1
656 data8 0x40215B2F9FEB5D87,0xC0327B539DEA5058 // A8,A7
657 data8 0x40555444B3E8D64D,0xC0655A2B26F9FC8A // A4,A3
658 // [2048;4096)
659 data8 0x3FE9F065A1C3D6B1,0xC015ACF6FAE8D78D // A10,A9
660 data8 0x409000FE383DD2B7,0x40CE7F5C1E8BCB8B // A2,A1
661 data8 0x40315324E5DB2EBE,0xC04274194EF70D18 // A8,A7
662 data8 0x4065504353FF2207,0xC075577FE1BFE7B6 // A4,A3
663 // [4096;8192)
664 data8 0x3FF9E6FBC6B1C70D,0xC025A62DAF76F85D // A10,A9
665 data8 0x40A0007E2F61EBE8,0x40E0A2A23FB5F6C3 // A2,A1
666 data8 0x40414E9BC0A0141A,0xC0527030F2B69D43 // A8,A7
667 data8 0x40754E417717B45B,0xC085562A447258E5 // A4,A3
668 //
669 data8 0xbfdffffffffaea15 // P1
670 data8 0x3FDD8B618D5AF8FE // point of local minimum on [1;2]
671 data8 0x3FED67F1C864BEB5 // ln(sqrt(2*Pi))
672 data8 0x4008000000000000 // 3.0
673 //
674 data8 0xBF9E1C289FB224AB,0x3FBF7422445C9460 // A6,A5
675 data8 0xBFF01E76D66F8D8A // A0
676 data8 0xBFE2788CFC6F91DA // A1 [1.0;1.25)
677 data8 0x3FCB8CC69000EB5C,0xBFD41997A0C2C641 // A6,A5
678 data8 0x3FFCAB0BFA0EA462 // A0
679 data8 0xBFBF19B9BCC38A42 // A0 [1.25;1.5)
680 data8 0x3FD51EE4DE0A364C,0xBFE00D7F98A16E4B // A6,A5
681 data8 0x40210CE1F327E9E4 // A0
682 data8 0x4001DB08F9DFA0CC // A0 [1.5;1.75)
683 data8 0x3FE24F606742D252,0xBFEC81D7D12574EC // A6,A5
684 data8 0x403BE636A63A9C27 // A0
685 data8 0x4000A0CB38D6CF0A // A0 [1.75;2.0)
686 data8 0x3FF1029A9DD542B4,0xBFFAD37C209D3B25 // A6,A5
687 data8 0x405385E6FD9BE7EA // A0
688 data8 0x478895F1C0000000 // Overflow boundary
689 data8 0x400062D97D26B523,0xC00A03E1529FF023 // A6,A5
690 data8 0x4069204C51E566CE // A0
691 data8 0x0000000000000000 // pad
692 data8 0x40101476B38FD501,0xC0199DE7B387C0FC // A6,A5
693 data8 0x407EB8DAEC83D759 // A0
694 data8 0x0000000000000000 // pad
695 data8 0x401FDB008D65125A,0xC0296B506E665581 // A6,A5
696 data8 0x409226D93107EF66 // A0
697 data8 0x0000000000000000 // pad
698 data8 0x402FB3EAAF3E7B2D,0xC039521142AD8E0D // A6,A5
699 data8 0x40A4EFA4F072792E // A0
700 data8 0x0000000000000000 // pad
701 data8 0x403FA024C66B2563,0xC0494569F250E691 // A6,A5
702 data8 0x40B7B747C9235BB8 // A0
703 data8 0x0000000000000000 // pad
704 data8 0x404F9607D6DA512C,0xC0593F0B2EDDB4BC // A6,A5
705 data8 0x40CA7E29C5F16DE2 // A0
706 data8 0x0000000000000000 // pad
707 data8 0x405F90C5F613D98D,0xC0693BD130E50AAF // A6,A5
708 data8 0x40DD4495238B190C // A0
709 data8 0x0000000000000000 // pad
710 //
711 // polynomial approximation of ln(sin(Pi*x)/(Pi*x)), |x| <= 0.5
712 data8 0xBFD58731A486E820,0xBFA4452CC28E15A9 // S16,S14
713 data8 0xBFD013F6E1B86C4F,0xBFD5B3F19F7A341F // S8,S6
714 data8 0xBFC86A0D5252E778,0xBFC93E08C9EE284B // S12,S10
715 data8 0xBFE15132555C9EDD,0xBFFA51A662480E35 // S4,S2
716 //
717 // [1.0;1.25)
718 data8 0xBFA697D6775F48EA,0x3FB9894B682A98E7 // A9,A8
719 data8 0xBFCA8969253CFF55,0x3FD15124EFB35D9D // A5,A4
720 data8 0xBFC1B00158AB719D,0x3FC5997D04E7F1C1 // A7,A6
721 data8 0xBFD9A4D50BAFF989,0x3FEA51A661F5176A // A3,A2
722 // [1.25;1.5)
723 data8 0x3F838E0D35A6171A,0xBF831BBBD61313B7 // A8,A7
724 data8 0x3FB08B40196425D0,0xBFC2E427A53EB830 // A4,A3
725 data8 0x3F9285DDDC20D6C3,0xBFA0C90C9C223044 // A6,A5
726 data8 0x3FDEF72BC8F5287C,0x3D890B3DAEBC1DFC // A2,A1
727 // [1.5;1.75)
728 data8 0x3F65D5A7EB31047F,0xBFA44EAC9BFA7FDE // A8,A7
729 data8 0x40051FEFE7A663D8,0xC012A5CFE00A2522 // A4,A3
730 data8 0x3FD0E1583AB00E08,0xBFF084AF95883BA5 // A6,A5
731 data8 0x40185982877AE0A2,0xC015F83DB73B57B7 // A2,A1
732 // [1.75;2.0)
733 data8 0x3F4A9222032EB39A,0xBF8CBC9587EEA5A3 // A8,A7
734 data8 0x3FF795400783BE49,0xC00851BC418B8A25 // A4,A3
735 data8 0x3FBBC992783E8C5B,0xBFDFA67E65E89B29 // A6,A5
736 data8 0x4012B408F02FAF88,0xC013284CE7CB0C39 // A2,A1
737 //
738 // roots
739 data8 0xC003A7FC9600F86C // -2.4570247382208005860
740 data8 0xC009260DBC9E59AF // -3.1435808883499798405
741 data8 0xC005FB410A1BD901 // -2.7476826467274126919
742 data8 0xC00FA471547C2FE5 // -3.9552942848585979085
743 //
744 // polynomial approximation of ln(GAMMA(x)) near roots
745 // near -2.4570247382208005860
746 data8 0x3FF694A6058D9592,0x40136EEBB003A92B // R3,R2
747 data8 0x3FF83FE966AF5360,0x3C90323B6D1FE86D // R1,R0
748 // near -3.1435808883499798405
749 data8 0x405C11371268DA38,0x4039D4D2977D2C23 // R3,R2
750 data8 0x401F20A65F2FAC62,0x3CDE9605E3AE7A62 // R1,R0
751 // near -2.7476826467274126919
752 data8 0xC034185AC31314FF,0x4023267F3C28DFE3 // R3,R2
753 data8 0xBFFEA12DA904B194,0x3CA8FB8530BA7689 // R1,R0
754 // near -2.7476826467274126919
755 data8 0xC0AD25359E70C888,0x406F76DEAEA1B8C6 // R3,R2
756 data8 0xC034B99D966C5644,0xBCBDDC0336980B58 // R1,R0
757 LOCAL_OBJECT_END(lgammaf_data)
759 //*********************************************************************
761 .section .text
762 GLOBAL_LIBM_ENTRY(__libm_lgammaf)
763 { .mfi
764 getf.exp GR_SignExp = f8
765 frcpa.s1 FR_InvX,p0 = f1,f8
766 mov GR_ExpOf2 = 0x10000
767 }
768 { .mfi
769 addl GR_ad_Data = @ltoff(lgammaf_data),gp
770 fcvt.fx.s1 FR_int_N = f8
771 mov GR_ExpMask = 0x1ffff
772 };;
773 { .mfi
774 getf.sig GR_Sig = f8
775 fclass.m p13,p0 = f8,0x1EF // is x NaTVal, NaN,
776 // +/-0, +/-INF or +/-deno?
777 mov GR_ExpBias = 0xffff
778 }
779 { .mfi
780 ld8 GR_ad_Data = [GR_ad_Data]
781 fma.s1 FR_Xp1 = f8,f1,f1
782 mov GR_StirlBound = 0x1000C
783 };;
784 { .mfi
785 setf.exp FR_2 = GR_ExpOf2
786 fmerge.se FR_x = f1,f8
787 dep.z GR_Ind = GR_SignExp,3,4
788 }
789 { .mfi
790 cmp.eq p8,p0 = GR_SignExp,GR_ExpBias
791 fcvt.fx.trunc.s1 FR_int_Ntrunc = f8
792 and GR_Exp = GR_ExpMask,GR_SignExp
793 };;
794 { .mfi
795 add GR_ad_C650 = 0xB20,GR_ad_Data
796 fcmp.lt.s1 p14,p15 = f8,f0
797 extr.u GR_Ind4T = GR_Sig,55,8
798 }
799 { .mfb
800 sub GR_PureExp = GR_Exp,GR_ExpBias
801 fnorm.s1 FR_NormX = f8
802 // jump if x is NaTVal, NaN, +/-0, +/-INF or +/-deno
803 (p13) br.cond.spnt lgammaf_spec
804 };;
805 lgammaf_core:
806 { .mfi
807 ldfpd FR_P1,FR_LocalMin = [GR_ad_C650],16
808 fms.s1 FR_xm2 = f8,f1,f1
809 add GR_ad_Co = 0x820,GR_ad_Data
810 }
811 { .mib
812 ldfpd FR_P3,FR_P2 = [GR_ad_Data],16
813 cmp.ltu p9,p0 = GR_SignExp,GR_ExpBias
814 // jump if x is from the interval [1; 2)
815 (p8) br.cond.spnt lgammaf_1_2
816 };;
817 { .mfi
818 setf.sig FR_int_Ln = GR_PureExp
819 fms.s1 FR_r = FR_InvX,f8,f1
820 shladd GR_ad_Co = GR_Ind,3,GR_ad_Co
821 }
822 { .mib
823 ldfpd FR_LnSqrt2Pi,FR_3 = [GR_ad_C650],16
824 cmp.lt p13,p12 = GR_Exp,GR_StirlBound
825 // jump if x is from the interval (0; 1)
826 (p9) br.cond.spnt lgammaf_0_1
827 };;
828 { .mfi
829 ldfpd FR_Ln2,FR_05 = [GR_ad_Data],16
830 fma.s1 FR_Xp2 = f1,f1,FR_Xp1 // (x+2)
831 shladd GR_ad_C650 = GR_Ind,2,GR_ad_C650
832 }
833 { .mfi
834 add GR_ad_Ce = 0x20,GR_ad_Co
835 nop.f 0
836 add GR_ad_C43 = 0x30,GR_ad_Co
837 };;
838 { .mfi
839 // load coefficients of polynomial approximation
840 // of ln(GAMMA(x)), 2 <= x < 2^13
841 (p13) ldfpd FR_A10,FR_A9 = [GR_ad_Co],16
842 fcvt.xf FR_N = FR_int_N
843 cmp.eq.unc p6,p7 = GR_ExpOf2,GR_SignExp
844 }
845 { .mib
846 (p13) ldfpd FR_A8,FR_A7 = [GR_ad_Ce]
847 (p14) cmp.le.unc p9,p0 = GR_StirlBound,GR_Exp
848 // jump if x is less or equal to -2^13
849 (p9) br.cond.spnt lgammaf_negstirling
850 };;
851 .pred.rel "mutex",p6,p7
852 { .mfi
853 (p13) ldfpd FR_A6,FR_A5 = [GR_ad_C650],16
854 (p6) fma.s1 FR_x = f0,f0,FR_NormX
855 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data
856 }
857 { .mfi
858 (p13) ldfpd FR_A4,FR_A3 = [GR_ad_C43]
859 (p7) fms.s1 FR_x = FR_x,f1,f1
860 (p14) mov GR_ReqBound = 0x20005
861 };;
862 { .mfi
863 (p13) ldfpd FR_A2,FR_A1 = [GR_ad_Co],16
864 fms.s1 FR_xm2 = FR_xm2,f1,f1
865 (p14) extr.u GR_Arg = GR_Sig,60,4
866 }
867 { .mfi
868 mov GR_SignOfGamma = 1 // set sign of gamma(x) to 1
869 fcvt.xf FR_Ntrunc = FR_int_Ntrunc
870 nop.i 0
871 };;
872 { .mfi
873 ldfd FR_T = [GR_ad_T]
874 fma.s1 FR_r2 = FR_r,FR_r,f0
875 shl GR_ReqBound = GR_ReqBound,3
876 }
877 { .mfi
878 add GR_ad_Co = 0xCA0,GR_ad_Data
879 fnma.s1 FR_Req = FR_Xp1,FR_NormX,f0 // -x*(x+1)
880 (p14) shladd GR_Arg = GR_Exp,4,GR_Arg
881 };;
882 { .mfi
883 (p13) ldfd FR_A0 = [GR_ad_C650]
884 fma.s1 FR_Xp3 = FR_2,f1,FR_Xp1 // (x+3)
885 (p14) cmp.le.unc p9,p0 = GR_Arg,GR_ReqBound
886 }
887 { .mfi
888 (p14) add GR_ad_Ce = 0x20,GR_ad_Co
889 fma.s1 FR_Xp4 = FR_2,FR_2,FR_NormX // (x+4)
890 (p15) add GR_ad_OvfBound = 0xBB8,GR_ad_Data
891 };;
892 { .mfi
893 // load coefficients of polynomial approximation
894 // of ln(sin(Pi*xf)/(Pi*xf)), |xf| <= 0.5
895 (p14) ldfpd FR_S16,FR_S14 = [GR_ad_Co],16
896 (p14) fms.s1 FR_Xf = FR_NormX,f1,FR_N // xf = x - [x]
897 (p14) sub GR_SignOfGamma = r0,GR_SignOfGamma // set sign of
898 // gamma(x) to -1
899 }
900 { .mfb
901 (p14) ldfpd FR_S12,FR_S10 = [GR_ad_Ce],16
902 fma.s1 FR_Xp5 = FR_2,FR_2,FR_Xp1 // (x+5)
903 // jump if x is from the interval (-9; 0)
904 (p9) br.cond.spnt lgammaf_negrecursion
905 };;
906 { .mfi
907 (p14) ldfpd FR_S8,FR_S6 = [GR_ad_Co],16
908 fma.s1 FR_P32 = FR_P3,FR_r,FR_P2
909 nop.i 0
910 }
911 { .mfb
912 (p14) ldfpd FR_S4,FR_S2 = [GR_ad_Ce],16
913 fma.s1 FR_x2 = FR_x,FR_x,f0
914 // jump if x is from the interval (-2^13; -9)
915 (p14) br.cond.spnt lgammaf_negpoly
916 };;
917 { .mfi
918 ldfd FR_OverflowBound = [GR_ad_OvfBound]
919 (p12) fcvt.xf FR_N = FR_int_Ln
920 // set p9 if signgum is 32-bit int
921 // set p10 if signgum is 64-bit int
922 cmp.eq p10,p9 = 8,r34
923 }
924 { .mfi
925 nop.m 0
926 (p12) fma.s1 FR_P10 = FR_P1,FR_r,f1
927 nop.i 0
928 };;
929 .pred.rel "mutex",p6,p7
930 .pred.rel "mutex",p9,p10
931 { .mfi
932 // store sign of gamma(x) as 32-bit int
933 (p9) st4 [r33] = GR_SignOfGamma
934 (p6) fma.s1 FR_xx = FR_x,FR_xm2,f0
935 nop.i 0
936 }
937 { .mfi
938 // store sign of gamma(x) as 64-bit int
939 (p10) st8 [r33] = GR_SignOfGamma
940 (p7) fma.s1 FR_xx = f0,f0,FR_x
941 nop.i 0
942 };;
943 { .mfi
944 nop.m 0
945 (p13) fma.s1 FR_A9 = FR_A10,FR_x,FR_A9
946 nop.i 0
947 }
948 { .mfi
949 nop.m 0
950 (p13) fma.s1 FR_A7 = FR_A8,FR_x,FR_A7
951 nop.i 0
952 };;
953 { .mfi
954 nop.m 0
955 (p13) fma.s1 FR_A5 = FR_A6,FR_x,FR_A5
956 nop.i 0
957 }
958 { .mfi
959 nop.m 0
960 (p13) fma.s1 FR_A3 = FR_A4,FR_x,FR_A3
961 nop.i 0
962 };;
963 { .mfi
964 nop.m 0
965 (p15) fcmp.eq.unc.s1 p8,p0 = FR_NormX,FR_2 // is input argument 2.0?
966 nop.i 0
967 }
968 { .mfi
969 nop.m 0
970 (p13) fma.s1 FR_A1 = FR_A2,FR_x,FR_A1
971 nop.i 0
972 };;
973 { .mfi
974 nop.m 0
975 (p12) fma.s1 FR_T = FR_N,FR_Ln2,FR_T
976 nop.i 0
977 }
978 { .mfi
979 nop.m 0
980 (p12) fma.s1 FR_P32 = FR_P32,FR_r2,FR_P10
981 nop.i 0
982 };;
983 { .mfi
984 nop.m 0
985 (p13) fma.s1 FR_x4 = FR_x2,FR_x2,f0
986 nop.i 0
987 }
988 { .mfi
989 nop.m 0
990 (p13) fma.s1 FR_x3 = FR_x2,FR_xx,f0
991 nop.i 0
992 };;
993 { .mfi
994 nop.m 0
995 (p13) fma.s1 FR_A7 = FR_A9,FR_x2,FR_A7
996 nop.i 0
997 }
998 { .mfb
999 nop.m 0
1000 (p8) fma.s.s0 f8 = f0,f0,f0
1001 (p8) br.ret.spnt b0 // fast exit for 2.0
1002 };;
1003 { .mfi
1004 nop.m 0
1005 (p6) fma.s1 FR_A0 = FR_A0,FR_xm2,f0
1006 nop.i 0
1007 }
1008 { .mfi
1009 nop.m 0
1010 (p13) fma.s1 FR_A3 = FR_A5,FR_x2,FR_A3
1011 nop.i 0
1012 };;
1013 { .mfi
1014 nop.m 0
1015 (p15) fcmp.le.unc.s1 p8,p0 = FR_OverflowBound,FR_NormX // overflow test
1016 nop.i 0
1017 }
1018 { .mfi
1019 nop.m 0
1020 (p12) fms.s1 FR_xm05 = FR_NormX,f1,FR_05
1021 nop.i 0
1022 };;
1023 { .mfi
1024 nop.m 0
1025 (p12) fma.s1 FR_Ln = FR_P32,FR_r,FR_T
1026 nop.i 0
1027 }
1028 { .mfi
1029 nop.m 0
1030 (p12) fms.s1 FR_LnSqrt2Pi = FR_LnSqrt2Pi,f1,FR_NormX
1031 nop.i 0
1032 };;
1033 { .mfi
1034 nop.m 0
1035 (p13) fma.s1 FR_A0 = FR_A1,FR_xx,FR_A0
1036 nop.i 0
1037 }
1038 { .mfb
1039 nop.m 0
1040 (p13) fma.s1 FR_A3 = FR_A7,FR_x4,FR_A3
1041 // jump if result overflows
1042 (p8) br.cond.spnt lgammaf_overflow
1043 };;
1044 .pred.rel "mutex",p12,p13
1045 { .mfi
1046 nop.m 0
1047 (p12) fma.s.s0 f8 = FR_Ln,FR_xm05,FR_LnSqrt2Pi
1048 nop.i 0
1049 }
1050 { .mfb
1051 nop.m 0
1052 (p13) fma.s.s0 f8 = FR_A3,FR_x3,FR_A0
1053 br.ret.sptk b0
1054 };;
1055 // branch for calculating of ln(GAMMA(x)) for 0 < x < 1
1056 //---------------------------------------------------------------------
1057 .align 32
1058 lgammaf_0_1:
1059 { .mfi
1060 getf.sig GR_Ind = FR_Xp1
1061 fma.s1 FR_r2 = FR_r,FR_r,f0
1062 mov GR_fff7 = 0xFFF7
1063 }
1064 { .mfi
1065 ldfpd FR_Ln2,FR_05 = [GR_ad_Data],16
1066 fma.s1 FR_P32 = FR_P3,FR_r,FR_P2
1067 // input argument cann't be equal to 1.0
1068 cmp.eq p0,p14 = r0,r0
1069 };;
1070 { .mfi
1071 getf.exp GR_Exp = FR_w
1072 fcvt.xf FR_N = FR_int_Ln
1073 add GR_ad_Co = 0xCE0,GR_ad_Data
1074 }
1075 { .mfi
1076 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data
1077 fma.s1 FR_P10 = FR_P1,FR_r,f1
1078 add GR_ad_Ce = 0xD00,GR_ad_Data
1079 };;
1080 { .mfi
1081 ldfd FR_T = [GR_ad_T]
1082 fma.s1 FR_w2 = FR_w,FR_w,f0
1083 extr.u GR_Ind = GR_Ind,61,2
1084 }
1085 { .mfi
1086 nop.m 0
1087 fma.s1 FR_Q32 = FR_P3,FR_w,FR_P2
1088 //// add GR_ad_C0 = 0xB30,GR_ad_Data
1089 add GR_ad_C0 = 0xB38,GR_ad_Data
1090 };;
1091 { .mfi
1092 and GR_Exp = GR_Exp,GR_ExpMask
1093 nop.f 0
1094 shladd GR_IndX8 = GR_Ind,3,r0
1095 }
1096 { .mfi
1097 shladd GR_IndX2 = GR_Ind,1,r0
1098 fma.s1 FR_Q10 = FR_P1,FR_w,f1
1099 cmp.eq p6,p15 = 0,GR_Ind
1100 };;
1101 { .mfi
1102 shladd GR_ad_Co = GR_IndX8,3,GR_ad_Co
1103 (p6) fma.s1 FR_x = f0,f0,FR_NormX
1104 shladd GR_ad_C0 = GR_IndX2,4,GR_ad_C0
1105 }
1106 { .mfi
1107 shladd GR_ad_Ce = GR_IndX8,3,GR_ad_Ce
1108 nop.f 0
1109 (p15) cmp.eq.unc p7,p8 = 1,GR_Ind
1110 };;
1111 .pred.rel "mutex",p7,p8
1112 { .mfi
1113 ldfpd FR_A8,FR_A7 = [GR_ad_Co],16
1114 (p7) fms.s1 FR_x = FR_NormX,f1,FR_LocalMin
1115 cmp.ge p10,p11 = GR_Exp,GR_fff7
1116 }
1117 { .mfb
1118 ldfpd FR_A6,FR_A5 = [GR_ad_Ce],16
1119 (p8) fma.s1 FR_x = f1,f1,FR_NormX
1120 br.cond.sptk lgamma_0_2_core
1121 };;
1122 // branch for calculating of ln(GAMMA(x)) for 1 <= x < 2
1123 //---------------------------------------------------------------------
1124 .align 32
1125 lgammaf_1_2:
1126 { .mfi
1127 add GR_ad_Co = 0xCF0,GR_ad_Data
1128 fcmp.eq.s1 p14,p0 = f1,FR_NormX // is input argument 1.0?
1129 extr.u GR_Ind = GR_Sig,61,2
1130 }
1131 { .mfi
1132 add GR_ad_Ce = 0xD10,GR_ad_Data
1133 nop.f 0
1134 //// add GR_ad_C0 = 0xB40,GR_ad_Data
1135 add GR_ad_C0 = 0xB48,GR_ad_Data
1136 };;
1137 { .mfi
1138 shladd GR_IndX8 = GR_Ind,3,r0
1139 nop.f 0
1140 shladd GR_IndX2 = GR_Ind,1,r0
1141 }
1142 { .mfi
1143 cmp.eq p6,p15 = 0,GR_Ind // p6 <- x from [1;1.25)
1144 nop.f 0
1145 cmp.ne p9,p0 = r0,r0
1146 };;
1147 { .mfi
1148 shladd GR_ad_Co = GR_IndX8,3,GR_ad_Co
1149 (p6) fms.s1 FR_x = FR_NormX,f1,f1 // reduced x for [1;1.25)
1150 shladd GR_ad_C0 = GR_IndX2,4,GR_ad_C0
1151 }
1152 { .mfi
1153 shladd GR_ad_Ce = GR_IndX8,3,GR_ad_Ce
1154 (p14) fma.s.s0 f8 = f0,f0,f0
1155 (p15) cmp.eq.unc p7,p8 = 1,GR_Ind // p7 <- x from [1.25;1.5)
1156 };;
1157 .pred.rel "mutex",p7,p8
1158 { .mfi
1159 ldfpd FR_A8,FR_A7 = [GR_ad_Co],16
1160 (p7) fms.s1 FR_x = FR_xm2,f1,FR_LocalMin
1161 nop.i 0
1162 }
1163 { .mfi
1164 ldfpd FR_A6,FR_A5 = [GR_ad_Ce],16
1165 (p8) fma.s1 FR_x = f0,f0,FR_NormX
1166 (p9) cmp.eq.unc p10,p11 = r0,r0
1167 };;
1168 lgamma_0_2_core:
1169 { .mmi
1170 ldfpd FR_A4,FR_A3 = [GR_ad_Co],16
1171 ldfpd FR_A2,FR_A1 = [GR_ad_Ce],16
1172 mov GR_SignOfGamma = 1 // set sign of gamma(x) to 1
1173 };;
1174 { .mfi
1175 // add GR_ad_C0 = 8,GR_ad_C0
1176 ldfd FR_A0 = [GR_ad_C0]
1177 nop.f 0
1178 // set p13 if signgum is 32-bit int
1179 // set p15 if signgum is 64-bit int
1180 cmp.eq p15,p13 = 8,r34
1181 };;
1182 .pred.rel "mutex",p13,p15
1183 { .mmf
1184 // store sign of gamma(x)
1185 (p13) st4 [r33] = GR_SignOfGamma // as 32-bit int
1186 (p15) st8 [r33] = GR_SignOfGamma // as 64-bit int
1187 (p11) fma.s1 FR_Q32 = FR_Q32,FR_w2,FR_Q10
1188 };;
1189 { .mfb
1190 nop.m 0
1191 (p10) fma.s1 FR_P32 = FR_P32,FR_r2,FR_P10
1192 (p14) br.ret.spnt b0 // fast exit for 1.0
1193 };;
1194 { .mfi
1195 nop.m 0
1196 (p10) fma.s1 FR_T = FR_N,FR_Ln2,FR_T
1197 cmp.eq p6,p7 = 0,GR_Ind // p6 <- x from [1;1.25)
1198 }
1199 { .mfi
1200 nop.m 0
1201 fma.s1 FR_x2 = FR_x,FR_x,f0
1202 cmp.eq p8,p0 = r0,r0 // set p8 to 1 that means we on [1;2]
1203 };;
1204 { .mfi
1205 nop.m 0
1206 (p11) fma.s1 FR_Ln = FR_Q32,FR_w,f0
1207 nop.i 0
1208 }
1209 { .mfi
1210 nop.m 0
1211 nop.f 0
1212 nop.i 0
1213 };;
1214 .pred.rel "mutex",p6,p7
1215 { .mfi
1216 nop.m 0
1217 (p6) fma.s1 FR_xx = f0,f0,FR_x
1218 nop.i 0
1219 }
1220 { .mfi
1221 nop.m 0
1222 (p7) fma.s1 FR_xx = f0,f0,f1
1223 nop.i 0
1224 };;
1225 { .mfi
1226 nop.m 0
1227 fma.s1 FR_A7 = FR_A8,FR_x,FR_A7
1228 nop.i 0
1229 }
1230 { .mfi
1231 nop.m 0
1232 fma.s1 FR_A5 = FR_A6,FR_x,FR_A5
1233 (p9) cmp.ne p8,p0 = r0,r0 // set p8 to 0 that means we on [0;1]
1234 };;
1235 { .mfi
1236 nop.m 0
1237 fma.s1 FR_A3 = FR_A4,FR_x,FR_A3
1238 nop.i 0
1239 }
1240 { .mfi
1241 nop.m 0
1242 fma.s1 FR_A1 = FR_A2,FR_x,FR_A1
1243 nop.i 0
1244 };;
1245 { .mfi
1246 nop.m 0
1247 fma.s1 FR_x4 = FR_x2,FR_x2,f0
1248 nop.i 0
1249 }
1250 { .mfi
1251 nop.m 0
1252 (p10) fma.s1 FR_Ln = FR_P32,FR_r,FR_T
1253 nop.i 0
1254 };;
1255 { .mfi
1256 nop.m 0
1257 fma.s1 FR_A5 = FR_A7,FR_x2,FR_A5
1258 nop.i 0
1259 }
1260 { .mfi
1261 nop.m 0
1262 fma.s1 FR_A1 = FR_A3,FR_x2,FR_A1
1263 nop.i 0
1264 };;
1265 .pred.rel "mutex",p9,p8
1266 { .mfi
1267 nop.m 0
1268 (p9) fms.d.s1 FR_A0 = FR_A0,FR_xx,FR_Ln
1269 nop.i 0
1270 }
1271 { .mfi
1272 nop.m 0
1273 (p8) fms.s1 FR_A0 = FR_A0,FR_xx,f0
1274 nop.i 0
1275 };;
1276 { .mfi
1277 nop.m 0
1278 fma.d.s1 FR_A1 = FR_A5,FR_x4,FR_A1
1279 nop.i 0
1280 }
1281 { .mfi
1282 nop.m 0
1283 nop.f 0
1284 nop.i 0
1285 };;
1286 .pred.rel "mutex",p6,p7
1287 { .mfi
1288 nop.m 0
1289 (p6) fma.s.s0 f8 = FR_A1,FR_x2,FR_A0
1290 nop.i 0
1291 }
1292 { .mfb
1293 nop.m 0
1294 (p7) fma.s.s0 f8 = FR_A1,FR_x,FR_A0
1295 br.ret.sptk b0
1296 };;
1297 // branch for calculating of ln(GAMMA(x)) for -9 < x < 1
1298 //---------------------------------------------------------------------
1299 .align 32
1300 lgammaf_negrecursion:
1301 { .mfi
1302 getf.sig GR_N = FR_int_Ntrunc
1303 fms.s1 FR_1pXf = FR_Xp2,f1,FR_Ntrunc // 1 + (x+1) - [x]
1304 mov GR_Neg2 = 2
1305 }
1306 { .mfi
1307 add GR_ad_Co = 0xCE0,GR_ad_Data
1308 fms.s1 FR_Xf = FR_Xp1,f1,FR_Ntrunc // (x+1) - [x]
1309 mov GR_Neg4 = 4
1310 };;
1311 { .mfi
1312 add GR_ad_Ce = 0xD00,GR_ad_Data
1313 fma.s1 FR_Xp6 = FR_2,FR_2,FR_Xp2 // (x+6)
1314 add GR_ad_C0 = 0xB30,GR_ad_Data
1315 }
1316 { .mfi
1317 sub GR_Neg2 = r0,GR_Neg2
1318 fma.s1 FR_Xp7 = FR_2,FR_3,FR_Xp1 // (x+7)
1319 sub GR_Neg4 = r0,GR_Neg4
1320 };;
1321 { .mfi
1322 cmp.ne p8,p0 = r0,GR_N
1323 fcmp.eq.s1 p13,p0 = FR_NormX,FR_Ntrunc
1324 and GR_IntNum = 0xF,GR_N
1325 }
1326 { .mfi
1327 cmp.lt p6,p0 = GR_N,GR_Neg2
1328 fma.s1 FR_Xp8 = FR_2,FR_3,FR_Xp2 // (x+8)
1329 cmp.lt p7,p0 = GR_N,GR_Neg4
1330 };;
1331 { .mfi
1332 getf.d GR_Arg = FR_NormX
1333 (p6) fma.s1 FR_Xp2 = FR_Xp2,FR_Xp3,f0
1334 (p8) tbit.z.unc p14,p15 = GR_IntNum,0
1335 }
1336 { .mfi
1337 sub GR_RootInd = 0xE,GR_IntNum
1338 (p7) fma.s1 FR_Xp4 = FR_Xp4,FR_Xp5,f0
1339 add GR_ad_Root = 0xDE0,GR_ad_Data
1340 };;
1341 { .mfi
1342 shladd GR_ad_Root = GR_RootInd,3,GR_ad_Root
1343 fms.s1 FR_x = FR_Xp1,f1,FR_Ntrunc // (x+1) - [x]
1344 nop.i 0
1345 }
1346 { .mfb
1347 nop.m 0
1348 nop.f 0
1349 (p13) br.cond.spnt lgammaf_singularity
1350 };;
1351 .pred.rel "mutex",p14,p15
1352 { .mfi
1353 cmp.gt p6,p0 = 0xA,GR_IntNum
1354 (p14) fma.s1 FR_Req = FR_Req,FR_Xf,f0
1355 cmp.gt p7,p0 = 0xD,GR_IntNum
1356 }
1357 { .mfi
1358 (p15) mov GR_SignOfGamma = 1 // set sign of gamma(x) to 1
1359 (p15) fnma.s1 FR_Req = FR_Req,FR_Xf,f0
1360 cmp.leu p0,p13 = 2,GR_RootInd
1361 };;
1362 { .mfi
1363 nop.m 0
1364 (p6) fma.s1 FR_Xp6 = FR_Xp6,FR_Xp7,f0
1365 (p13) add GR_ad_RootCo = 0xE00,GR_ad_Data
1366 };;
1367 { .mfi
1368 nop.m 0
1369 fcmp.eq.s1 p12,p11 = FR_1pXf,FR_2
1370 nop.i 0
1371 };;
1372 { .mfi
1373 getf.sig GR_Sig = FR_1pXf
1374 fcmp.le.s1 p9,p0 = FR_05,FR_Xf
1375 nop.i 0
1376 }
1377 { .mfi
1378 (p13) shladd GR_RootInd = GR_RootInd,4,r0
1379 (p7) fma.s1 FR_Xp2 = FR_Xp2,FR_Xp4,f0
1380 (p8) cmp.gt.unc p10,p0 = 0x9,GR_IntNum
1381 };;
1382 .pred.rel "mutex",p11,p12
1383 { .mfi
1384 nop.m 0
1385 (p10) fma.s1 FR_Req = FR_Req,FR_Xp8,f0
1386 (p11) extr.u GR_Ind = GR_Sig,61,2
1387 }
1388 { .mfi
1389 (p13) add GR_RootInd = GR_RootInd,GR_RootInd
1390 nop.f 0
1391 (p12) mov GR_Ind = 3
1392 };;
1393 { .mfi
1394 shladd GR_IndX2 = GR_Ind,1,r0
1395 nop.f 0
1396 cmp.gt p14,p0 = 2,GR_Ind
1397 }
1398 { .mfi
1399 shladd GR_IndX8 = GR_Ind,3,r0
1400 nop.f 0
1401 cmp.eq p6,p0 = 1,GR_Ind
1402 };;
1403 .pred.rel "mutex",p6,p9
1404 { .mfi
1405 shladd GR_ad_Co = GR_IndX8,3,GR_ad_Co
1406 (p6) fms.s1 FR_x = FR_Xf,f1,FR_LocalMin
1407 cmp.gt p10,p0 = 0xB,GR_IntNum
1408 }
1409 { .mfi
1410 shladd GR_ad_Ce = GR_IndX8,3,GR_ad_Ce
1411 (p9) fma.s1 FR_x = f0,f0,FR_1pXf
1412 shladd GR_ad_C0 = GR_IndX2,4,GR_ad_C0
1413 };;
1414 { .mfi
1415 // load coefficients of polynomial approximation
1416 // of ln(GAMMA(x)), 1 <= x < 2
1417 ldfpd FR_A8,FR_A7 = [GR_ad_Co],16
1418 (p10) fma.s1 FR_Xp2 = FR_Xp2,FR_Xp6,f0
1419 add GR_ad_C0 = 8,GR_ad_C0
1420 }
1421 { .mfi
1422 ldfpd FR_A6,FR_A5 = [GR_ad_Ce],16
1423 nop.f 0
1424 (p14) add GR_ad_Root = 0x10,GR_ad_Root
1425 };;
1426 { .mfi
1427 ldfpd FR_A4,FR_A3 = [GR_ad_Co],16
1428 nop.f 0
1429 add GR_ad_RootCe = 0xE10,GR_ad_Data
1430 }
1431 { .mfi
1432 ldfpd FR_A2,FR_A1 = [GR_ad_Ce],16
1433 nop.f 0
1434 (p14) add GR_RootInd = 0x40,GR_RootInd
1435 };;
1436 { .mmi
1437 ldfd FR_A0 = [GR_ad_C0]
1438 (p13) add GR_ad_RootCo = GR_ad_RootCo,GR_RootInd
1439 (p13) add GR_ad_RootCe = GR_ad_RootCe,GR_RootInd
1440 };;
1441 { .mmi
1442 (p13) ld8 GR_Root = [GR_ad_Root]
1443 (p13) ldfd FR_Root = [GR_ad_Root]
1444 mov GR_ExpBias = 0xffff
1445 };;
1446 { .mfi
1447 nop.m 0
1448 fma.s1 FR_x2 = FR_x,FR_x,f0
1449 nop.i 0
1450 }
1451 { .mlx
1452 (p8) cmp.gt.unc p10,p0 = 0xF,GR_IntNum
1453 movl GR_Dx = 0x000000014F8B588E
1454 };;
1455 { .mfi
1456 // load coefficients of polynomial approximation
1457 // of ln(GAMMA(x)), x is close to one of negative roots
1458 (p13) ldfpd FR_R3,FR_R2 = [GR_ad_RootCo]
1459 // argumenth for logarithm
1460 (p10) fma.s1 FR_Req = FR_Req,FR_Xp2,f0
1461 mov GR_ExpMask = 0x1ffff
1462 }
1463 { .mfi
1464 (p13) ldfpd FR_R1,FR_R0 = [GR_ad_RootCe]
1465 nop.f 0
1466 // set p9 if signgum is 32-bit int
1467 // set p8 if signgum is 64-bit int
1468 cmp.eq p8,p9 = 8,r34
1469 };;
1470 .pred.rel "mutex",p9,p8
1471 { .mfi
1472 (p9) st4 [r33] = GR_SignOfGamma // as 32-bit int
1473 fma.s1 FR_A7 = FR_A8,FR_x,FR_A7
1474 (p13) sub GR_Root = GR_Arg,GR_Root
1475 }
1476 { .mfi
1477 (p8) st8 [r33] = GR_SignOfGamma // as 64-bit int
1478 fma.s1 FR_A5 = FR_A6,FR_x,FR_A5
1479 nop.i 0
1480 };;
1481 { .mfi
1482 nop.m 0
1483 fms.s1 FR_w = FR_Req,f1,f1
1484 (p13) add GR_Root = GR_Root,GR_Dx
1485 }
1486 { .mfi
1487 nop.m 0
1488 nop.f 0
1489 (p13) add GR_2xDx = GR_Dx,GR_Dx
1490 };;
1491 { .mfi
1492 nop.m 0
1493 fma.s1 FR_A3 = FR_A4,FR_x,FR_A3
1494 nop.i 0
1495 }
1496 { .mfi
1497 nop.m 0
1498 fma.s1 FR_A1 = FR_A2,FR_x,FR_A1
1499 (p13) cmp.leu.unc p10,p0 = GR_Root,GR_2xDx
1500 };;
1501 { .mfi
1502 nop.m 0
1503 frcpa.s1 FR_InvX,p0 = f1,FR_Req
1504 nop.i 0
1505 }
1506 { .mfi
1507 nop.m 0
1508 (p10) fms.s1 FR_rx = FR_NormX,f1,FR_Root
1509 nop.i 0
1510 };;
1511 { .mfi
1512 getf.exp GR_SignExp = FR_Req
1513 fma.s1 FR_x4 = FR_x2,FR_x2,f0
1514 nop.i 0
1515 };;
1516 { .mfi
1517 getf.sig GR_Sig = FR_Req
1518 fma.s1 FR_A5 = FR_A7,FR_x2,FR_A5
1519 nop.i 0
1520 };;
1521 { .mfi
1522 sub GR_PureExp = GR_SignExp,GR_ExpBias
1523 fma.s1 FR_w2 = FR_w,FR_w,f0
1524 nop.i 0
1525 }
1526 { .mfi
1527 nop.m 0
1528 fma.s1 FR_Q32 = FR_P3,FR_w,FR_P2
1529 nop.i 0
1530 };;
1531 { .mfi
1532 setf.sig FR_int_Ln = GR_PureExp
1533 fma.s1 FR_A1 = FR_A3,FR_x2,FR_A1
1534 extr.u GR_Ind4T = GR_Sig,55,8
1535 }
1536 { .mfi
1537 nop.m 0
1538 fma.s1 FR_Q10 = FR_P1,FR_w,f1
1539 nop.i 0
1540 };;
1541 { .mfi
1542 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data
1543 fms.s1 FR_r = FR_InvX,FR_Req,f1
1544 nop.i 0
1545 }
1546 { .mfi
1547 nop.m 0
1548 (p10) fms.s1 FR_rx2 = FR_rx,FR_rx,f0
1549 nop.i 0
1550 };;
1551 { .mfi
1552 ldfd FR_T = [GR_ad_T]
1553 (p10) fma.s1 FR_R2 = FR_R3,FR_rx,FR_R2
1554 nop.i 0
1555 }
1556 { .mfi
1557 nop.m 0
1558 (p10) fma.s1 FR_R0 = FR_R1,FR_rx,FR_R0
1559 nop.i 0
1560 };;
1561 { .mfi
1562 getf.exp GR_Exp = FR_w
1563 fma.s1 FR_A1 = FR_A5,FR_x4,FR_A1
1564 mov GR_ExpMask = 0x1ffff
1565 }
1566 { .mfi
1567 nop.m 0
1568 fma.s1 FR_Q32 = FR_Q32, FR_w2,FR_Q10
1569 nop.i 0
1570 };;
1571 { .mfi
1572 nop.m 0
1573 fma.s1 FR_r2 = FR_r,FR_r,f0
1574 mov GR_fff7 = 0xFFF7
1575 }
1576 { .mfi
1577 nop.m 0
1578 fma.s1 FR_P32 = FR_P3,FR_r,FR_P2
1579 nop.i 0
1580 };;
1581 { .mfi
1582 nop.m 0
1583 fma.s1 FR_P10 = FR_P1,FR_r,f1
1584 and GR_Exp = GR_ExpMask,GR_Exp
1585 }
1586 { .mfb
1587 nop.m 0
1588 (p10) fma.s.s0 f8 = FR_R2,FR_rx2,FR_R0
1589 (p10) br.ret.spnt b0 // exit for arguments close to negative roots
1590 };;
1591 { .mfi
1592 nop.m 0
1593 fcvt.xf FR_N = FR_int_Ln
1594 nop.i 0
1595 }
1596 { .mfi
1597 cmp.ge p14,p15 = GR_Exp,GR_fff7
1598 nop.f 0
1599 nop.i 0
1600 };;
1601 { .mfi
1602 nop.m 0
1603 fma.s1 FR_A0 = FR_A1,FR_x,FR_A0
1604 nop.i 0
1605 }
1606 { .mfi
1607 nop.m 0
1608 (p15) fma.s1 FR_Ln = FR_Q32,FR_w,f0
1609 nop.i 0
1610 };;
1611 { .mfi
1612 nop.m 0
1613 (p14) fma.s1 FR_P32 = FR_P32,FR_r2,FR_P10
1614 cmp.eq p6,p7 = 0,GR_Ind
1615 };;
1616 { .mfi
1617 nop.m 0
1618 (p14) fma.s1 FR_T = FR_N,FR_Ln2,FR_T
1619 nop.i 0
1620 };;
1621 { .mfi
1622 nop.m 0
1623 (p14) fma.s1 FR_Ln = FR_P32,FR_r,FR_T
1624 nop.i 0
1625 };;
1626 .pred.rel "mutex",p6,p7
1627 { .mfi
1628 nop.m 0
1629 (p6) fms.s.s0 f8 = FR_A0,FR_x,FR_Ln
1630 nop.i 0
1631 }
1632 { .mfb
1633 nop.m 0
1634 (p7) fms.s.s0 f8 = FR_A0,f1,FR_Ln
1635 br.ret.sptk b0
1636 };;
1638 // branch for calculating of ln(GAMMA(x)) for x < -2^13
1639 //---------------------------------------------------------------------
1640 .align 32
1641 lgammaf_negstirling:
1642 { .mfi
1643 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data
1644 fms.s1 FR_Xf = FR_NormX,f1,FR_N // xf = x - [x]
1645 mov GR_SingBound = 0x10016
1646 }
1647 { .mfi
1648 add GR_ad_Co = 0xCA0,GR_ad_Data
1649 fma.s1 FR_P32 = FR_P3,FR_r,FR_P2
1650 nop.i 0
1651 };;
1652 { .mfi
1653 ldfd FR_T = [GR_ad_T]
1654 fcvt.xf FR_int_Ln = FR_int_Ln
1655 cmp.le p6,p0 = GR_SingBound,GR_Exp
1656 }
1657 { .mfb
1658 add GR_ad_Ce = 0x20,GR_ad_Co
1659 fma.s1 FR_r2 = FR_r,FR_r,f0
1660 (p6) br.cond.spnt lgammaf_singularity
1661 };;
1662 { .mfi
1663 // load coefficients of polynomial approximation
1664 // of ln(sin(Pi*xf)/(Pi*xf)), |xf| <= 0.5
1665 ldfpd FR_S16,FR_S14 = [GR_ad_Co],16
1666 fma.s1 FR_P10 = FR_P1,FR_r,f1
1667 nop.i 0
1668 }
1669 { .mfi
1670 ldfpd FR_S12,FR_S10 = [GR_ad_Ce],16
1671 fms.s1 FR_xm05 = FR_NormX,f1,FR_05
1672 nop.i 0
1673 };;
1674 { .mmi
1675 ldfpd FR_S8,FR_S6 = [GR_ad_Co],16
1676 ldfpd FR_S4,FR_S2 = [GR_ad_Ce],16
1677 nop.i 0
1678 };;
1679 { .mfi
1680 getf.sig GR_N = FR_int_Ntrunc // signgam calculation
1681 fma.s1 FR_Xf2 = FR_Xf,FR_Xf,f0
1682 nop.i 0
1683 };;
1684 { .mfi
1685 nop.m 0
1686 frcpa.s1 FR_InvXf,p0 = f1,FR_Xf
1687 nop.i 0
1688 };;
1689 { .mfi
1690 getf.d GR_Arg = FR_Xf
1691 fcmp.eq.s1 p6,p0 = FR_NormX,FR_N
1692 mov GR_ExpBias = 0x3FF
1693 };;
1694 { .mfi
1695 nop.m 0
1696 fma.s1 FR_T = FR_int_Ln,FR_Ln2,FR_T
1697 extr.u GR_Exp = GR_Arg,52,11
1698 }
1699 { .mfi
1700 nop.m 0
1701 fma.s1 FR_P32 = FR_P32,FR_r2,FR_P10
1702 nop.i 0
1703 };;
1704 { .mfi
1705 sub GR_PureExp = GR_Exp,GR_ExpBias
1706 fma.s1 FR_S14 = FR_S16,FR_Xf2,FR_S14
1707 extr.u GR_Ind4T = GR_Arg,44,8
1708 }
1709 { .mfb
1710 mov GR_SignOfGamma = 1 // set signgam to -1
1711 fma.s1 FR_S10 = FR_S12,FR_Xf2,FR_S10
1712 (p6) br.cond.spnt lgammaf_singularity
1713 };;
1714 { .mfi
1715 setf.sig FR_int_Ln = GR_PureExp
1716 fms.s1 FR_rf = FR_InvXf,FR_Xf,f1
1717 // set p14 if GR_N is even
1718 tbit.z p14,p0 = GR_N,0
1719 }
1720 { .mfi
1721 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data
1722 fma.s1 FR_Xf4 = FR_Xf2,FR_Xf2,f0
1723 nop.i 0
1724 };;
1725 { .mfi
1726 (p14) sub GR_SignOfGamma = r0,GR_SignOfGamma // set signgam to -1
1727 fma.s1 FR_S6 = FR_S8,FR_Xf2,FR_S6
1728 nop.i 0
1729 }
1730 { .mfi
1731 // set p9 if signgum is 32-bit int
1732 // set p10 if signgum is 64-bit int
1733 cmp.eq p10,p9 = 8,r34
1734 fma.s1 FR_S2 = FR_S4,FR_Xf2,FR_S2
1735 nop.i 0
1736 };;
1737 { .mfi
1738 ldfd FR_Tf = [GR_ad_T]
1739 fma.s1 FR_Ln = FR_P32,FR_r,FR_T
1740 nop.i 0
1741 }
1742 { .mfi
1743 nop.m 0
1744 fma.s1 FR_LnSqrt2Pi = FR_LnSqrt2Pi,f1,FR_NormX
1745 nop.i 0
1746 };;
1747 .pred.rel "mutex",p9,p10
1748 { .mfi
1749 (p9) st4 [r33] = GR_SignOfGamma // as 32-bit int
1750 fma.s1 FR_rf2 = FR_rf,FR_rf,f0
1751 nop.i 0
1752 }
1753 { .mfi
1754 (p10) st8 [r33] = GR_SignOfGamma // as 64-bit int
1755 fma.s1 FR_S10 = FR_S14,FR_Xf4,FR_S10
1756 nop.i 0
1757 };;
1758 { .mfi
1759 nop.m 0
1760 fma.s1 FR_P32f = FR_P3,FR_rf,FR_P2
1761 nop.i 0
1762 }
1763 { .mfi
1764 nop.m 0
1765 fma.s1 FR_Xf8 = FR_Xf4,FR_Xf4,f0
1766 nop.i 0
1767 };;
1768 { .mfi
1769 nop.m 0
1770 fma.s1 FR_P10f = FR_P1,FR_rf,f1
1771 nop.i 0
1772 }
1773 { .mfi
1774 nop.m 0
1775 fma.s1 FR_S2 = FR_S6,FR_Xf4,FR_S2
1776 nop.i 0
1777 };;
1778 { .mfi
1779 nop.m 0
1780 fms.s1 FR_Ln = FR_Ln,FR_xm05,FR_LnSqrt2Pi
1781 nop.i 0
1782 };;
1783 { .mfi
1784 nop.m 0
1785 fcvt.xf FR_Nf = FR_int_Ln
1786 nop.i 0
1787 };;
1788 { .mfi
1789 nop.m 0
1790 fma.s1 FR_S2 = FR_S10,FR_Xf8,FR_S2
1791 nop.i 0
1792 };;
1793 { .mfi
1794 nop.m 0
1795 fma.s1 FR_Tf = FR_Nf,FR_Ln2,FR_Tf
1796 nop.i 0
1797 }
1798 { .mfi
1799 nop.m 0
1800 fma.s1 FR_P32f = FR_P32f,FR_rf2,FR_P10f // ??????
1801 nop.i 0
1802 };;
1803 { .mfi
1804 nop.m 0
1805 fnma.s1 FR_Ln = FR_S2,FR_Xf2,FR_Ln
1806 nop.i 0
1807 };;
1808 { .mfi
1809 nop.m 0
1810 fma.s1 FR_Lnf = FR_P32f,FR_rf,FR_Tf
1811 nop.i 0
1812 };;
1813 { .mfb
1814 nop.m 0
1815 fms.s.s0 f8 = FR_Ln,f1,FR_Lnf
1816 br.ret.sptk b0
1817 };;
1818 // branch for calculating of ln(GAMMA(x)) for -2^13 < x < -9
1819 //---------------------------------------------------------------------
1820 .align 32
1821 lgammaf_negpoly:
1822 { .mfi
1823 getf.d GR_Arg = FR_Xf
1824 frcpa.s1 FR_InvXf,p0 = f1,FR_Xf
1825 mov GR_ExpBias = 0x3FF
1826 }
1827 { .mfi
1828 nop.m 0
1829 fma.s1 FR_Xf2 = FR_Xf,FR_Xf,f0
1830 nop.i 0
1831 };;
1832 { .mfi
1833 getf.sig GR_N = FR_int_Ntrunc
1834 fcvt.xf FR_N = FR_int_Ln
1835 mov GR_SignOfGamma = 1
1836 }
1837 { .mfi
1838 nop.m 0
1839 fma.s1 FR_A9 = FR_A10,FR_x,FR_A9
1840 nop.i 0
1841 };;
1842 { .mfi
1843 nop.m 0
1844 fma.s1 FR_P10 = FR_P1,FR_r,f1
1845 extr.u GR_Exp = GR_Arg,52,11
1846 }
1847 { .mfi
1848 nop.m 0
1849 fma.s1 FR_x4 = FR_x2,FR_x2,f0
1850 nop.i 0
1851 };;
1852 { .mfi
1853 sub GR_PureExp = GR_Exp,GR_ExpBias
1854 fma.s1 FR_A7 = FR_A8,FR_x,FR_A7
1855 tbit.z p14,p0 = GR_N,0
1856 }
1857 { .mfi
1858 nop.m 0
1859 fma.s1 FR_A5 = FR_A6,FR_x,FR_A5
1860 nop.i 0
1861 };;
1862 { .mfi
1863 setf.sig FR_int_Ln = GR_PureExp
1864 fma.s1 FR_A3 = FR_A4,FR_x,FR_A3
1865 nop.i 0
1866 }
1867 { .mfi
1868 nop.m 0
1869 fma.s1 FR_A1 = FR_A2,FR_x,FR_A1
1870 (p14) sub GR_SignOfGamma = r0,GR_SignOfGamma
1871 };;
1872 { .mfi
1873 nop.m 0
1874 fms.s1 FR_rf = FR_InvXf,FR_Xf,f1
1875 nop.i 0
1876 }
1877 { .mfi
1878 nop.m 0
1879 fma.s1 FR_Xf4 = FR_Xf2,FR_Xf2,f0
1880 nop.i 0
1881 };;
1882 { .mfi
1883 nop.m 0
1884 fma.s1 FR_S14 = FR_S16,FR_Xf2,FR_S14
1885 nop.i 0
1886 }
1887 { .mfi
1888 nop.m 0
1889 fma.s1 FR_S10 = FR_S12,FR_Xf2,FR_S10
1890 nop.i 0
1891 };;
1892 { .mfi
1893 nop.m 0
1894 fma.s1 FR_T = FR_N,FR_Ln2,FR_T
1895 nop.i 0
1896 }
1897 { .mfi
1898 nop.m 0
1899 fma.s1 FR_P32 = FR_P32,FR_r2,FR_P10
1900 nop.i 0
1901 };;
1902 { .mfi
1903 nop.m 0
1904 fma.s1 FR_S6 = FR_S8,FR_Xf2,FR_S6
1905 extr.u GR_Ind4T = GR_Arg,44,8
1906 }
1907 { .mfi
1908 nop.m 0
1909 fma.s1 FR_S2 = FR_S4,FR_Xf2,FR_S2
1910 nop.i 0
1911 };;
1912 { .mfi
1913 nop.m 0
1914 fma.s1 FR_A7 = FR_A9,FR_x2,FR_A7
1915 nop.i 0
1916 }
1917 { .mfi
1918 shladd GR_ad_T = GR_Ind4T,3,GR_ad_Data
1919 fma.s1 FR_A3 = FR_A5,FR_x2,FR_A3
1920 nop.i 0
1921 };;
1922 { .mfi
1923 nop.m 0
1924 fma.s1 FR_Xf8 = FR_Xf4,FR_Xf4,f0
1925 nop.i 0
1926 }
1927 { .mfi
1928 nop.m 0
1929 fma.s1 FR_rf2 = FR_rf,FR_rf,f0
1930 nop.i 0
1931 };;
1932 { .mfi
1933 nop.m 0
1934 fma.s1 FR_P32f = FR_P3,FR_rf,FR_P2
1935 nop.i 0
1936 }
1937 { .mfi
1938 nop.m 0
1939 fma.s1 FR_P10f = FR_P1,FR_rf,f1
1940 nop.i 0
1941 };;
1942 { .mfi
1943 ldfd FR_Tf = [GR_ad_T]
1944 fma.s1 FR_Ln = FR_P32,FR_r,FR_T
1945 nop.i 0
1946 }
1947 { .mfi
1948 nop.m 0
1949 fma.s1 FR_A0 = FR_A1,FR_x,FR_A0
1950 nop.i 0
1951 };;
1952 { .mfi
1953 nop.m 0
1954 fma.s1 FR_S10 = FR_S14,FR_Xf4,FR_S10
1955 nop.i 0
1956 }
1957 { .mfi
1958 nop.m 0
1959 fma.s1 FR_S2 = FR_S6,FR_Xf4,FR_S2
1960 nop.i 0
1961 };;
1962 { .mfi
1963 nop.m 0
1964 fcvt.xf FR_Nf = FR_int_Ln
1965 nop.i 0
1966 }
1967 { .mfi
1968 nop.m 0
1969 fma.s1 FR_A3 = FR_A7,FR_x4,FR_A3
1970 nop.i 0
1971 };;
1972 { .mfi
1973 nop.m 0
1974 fcmp.eq.s1 p13,p0 = FR_NormX,FR_Ntrunc
1975 nop.i 0
1976 }
1977 { .mfi
1978 nop.m 0
1979 fnma.s1 FR_x3 = FR_x2,FR_x,f0 // -x^3
1980 nop.i 0
1981 };;
1982 { .mfi
1983 nop.m 0
1984 fma.s1 FR_P32f = FR_P32f,FR_rf2,FR_P10f
1985 nop.i 0
1986 };;
1987 { .mfb
1988 // set p9 if signgum is 32-bit int
1989 // set p10 if signgum is 64-bit int
1990 cmp.eq p10,p9 = 8,r34
1991 fma.s1 FR_S2 = FR_S10,FR_Xf8,FR_S2
1992 (p13) br.cond.spnt lgammaf_singularity
1993 };;
1994 .pred.rel "mutex",p9,p10
1995 { .mmf
1996 (p9) st4 [r33] = GR_SignOfGamma // as 32-bit int
1997 (p10) st8 [r33] = GR_SignOfGamma // as 64-bit int
1998 fms.s1 FR_A0 = FR_A3,FR_x3,FR_A0 // -A3*x^3-A0
1999 };;
2000 { .mfi
2001 nop.m 0
2002 fma.s1 FR_Tf = FR_Nf,FR_Ln2,FR_Tf
2003 nop.i 0
2004 };;
2005 { .mfi
2006 nop.m 0
2007 fma.s1 FR_Ln = FR_S2,FR_Xf2,FR_Ln // S2*Xf^2+Ln
2008 nop.i 0
2009 };;
2010 { .mfi
2011 nop.m 0
2012 fma.s1 FR_Lnf = FR_P32f,FR_rf,FR_Tf
2013 nop.i 0
2014 };;
2015 { .mfi
2016 nop.m 0
2017 fms.s1 FR_Ln = FR_A0,f1,FR_Ln
2018 nop.i 0
2019 };;
2020 { .mfb
2021 nop.m 0
2022 fms.s.s0 f8 = FR_Ln,f1,FR_Lnf
2023 br.ret.sptk b0
2024 };;
2025 // branch for handling +/-0, NaT, QNaN, +/-INF and denormalised numbers
2026 //---------------------------------------------------------------------
2027 .align 32
2028 lgammaf_spec:
2029 { .mfi
2030 getf.exp GR_SignExp = FR_NormX
2031 fclass.m p6,p0 = f8,0x21 // is arg +INF?
2032 mov GR_SignOfGamma = 1 // set signgam to 1
2033 };;
2034 { .mfi
2035 getf.sig GR_Sig = FR_NormX
2036 fclass.m p7,p0 = f8,0xB // is x deno?
2037 // set p11 if signgum is 32-bit int
2038 // set p12 if signgum is 64-bit int
2039 cmp.eq p12,p11 = 8,r34
2040 };;
2041 .pred.rel "mutex",p11,p12
2042 { .mfi
2043 // store sign of gamma(x) as 32-bit int
2044 (p11) st4 [r33] = GR_SignOfGamma
2045 fclass.m p8,p0 = f8,0x1C0 // is arg NaT or NaN?
2046 dep.z GR_Ind = GR_SignExp,3,4
2047 }
2048 { .mib
2049 // store sign of gamma(x) as 64-bit int
2050 (p12) st8 [r33] = GR_SignOfGamma
2051 and GR_Exp = GR_ExpMask,GR_SignExp
2052 (p6) br.ret.spnt b0 // exit for +INF
2053 };;
2054 { .mfi
2055 sub GR_PureExp = GR_Exp,GR_ExpBias
2056 fclass.m p9,p0 = f8,0x22 // is arg -INF?
2057 extr.u GR_Ind4T = GR_Sig,55,8
2058 }
2059 { .mfb
2060 nop.m 0
2061 (p7) fma.s0 FR_tmp = f1,f1,f8
2062 (p7) br.cond.sptk lgammaf_core
2063 };;
2064 { .mfb
2065 nop.m 0
2066 (p8) fms.s.s0 f8 = f8,f1,f8
2067 (p8) br.ret.spnt b0 // exit for NaT and NaN
2068 };;
2069 { .mfb
2070 nop.m 0
2071 (p9) fmerge.s f8 = f1,f8
2072 (p9) br.ret.spnt b0 // exit -INF
2073 };;
2074 // branch for handling negative integers and +/-0
2075 //---------------------------------------------------------------------
2076 .align 32
2077 lgammaf_singularity:
2078 { .mfi
2079 mov GR_SignOfGamma = 1 // set signgam to 1
2080 fclass.m p6,p0 = f8,0x6 // is x -0?
2081 mov GR_TAG = 109 // negative
2082 }
2083 { .mfi
2084 mov GR_ad_SignGam = r33
2085 fma.s1 FR_X = f0,f0,f8
2086 nop.i 0
2087 };;
2088 { .mfi
2089 nop.m 0
2090 frcpa.s0 f8,p0 = f1,f0
2091 // set p9 if signgum is 32-bit int
2092 // set p10 if signgum is 64-bit int
2093 cmp.eq p10,p9 = 8,r34
2094 }
2095 { .mib
2096 nop.m 0
2097 (p6) sub GR_SignOfGamma = r0,GR_SignOfGamma
2098 br.cond.sptk lgammaf_libm_err
2099 };;
2100 // overflow (x > OVERFLOV_BOUNDARY)
2101 //---------------------------------------------------------------------
2102 .align 32
2103 lgammaf_overflow:
2104 { .mfi
2105 nop.m 0
2106 nop.f 0
2107 mov r8 = 0x1FFFE
2108 };;
2109 { .mfi
2110 setf.exp f9 = r8
2111 fmerge.s FR_X = f8,f8
2112 mov GR_TAG = 108 // overflow
2113 };;
2114 { .mfi
2115 mov GR_ad_SignGam = r33
2116 nop.f 0
2117 // set p9 if signgum is 32-bit int
2118 // set p10 if signgum is 64-bit int
2119 cmp.eq p10,p9 = 8,r34
2120 }
2121 { .mfi
2122 nop.m 0
2123 fma.s.s0 f8 = f9,f9,f0 // Set I,O and +INF result
2124 nop.i 0
2125 };;
2126 // gate to __libm_error_support#
2127 //---------------------------------------------------------------------
2128 .align 32
2129 lgammaf_libm_err:
2130 { .mmi
2131 alloc r32 = ar.pfs,1,4,4,0
2132 mov GR_Parameter_TAG = GR_TAG
2133 nop.i 0
2134 };;
2135 .pred.rel "mutex",p9,p10
2136 { .mmi
2137 // store sign of gamma(x) as 32-bit int
2138 (p9) st4 [GR_ad_SignGam] = GR_SignOfGamma
2139 // store sign of gamma(x) as 64-bit int
2140 (p10) st8 [GR_ad_SignGam] = GR_SignOfGamma
2141 nop.i 0
2142 };;
2143 GLOBAL_LIBM_END(__libm_lgammaf)
2146 LOCAL_LIBM_ENTRY(__libm_error_region)
2147 .prologue
2148 { .mfi
2149 add GR_Parameter_Y=-32,sp // Parameter 2 value
2150 nop.f 0
2151 .save ar.pfs,GR_SAVE_PFS
2152 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
2153 }
2154 { .mfi
2155 .fframe 64
2156 add sp=-64,sp // Create new stack
2157 nop.f 0
2158 mov GR_SAVE_GP=gp // Save gp
2159 };;
2160 { .mmi
2161 stfs [GR_Parameter_Y] = FR_Y,16 // STORE Parameter 2 on stack
2162 add GR_Parameter_X = 16,sp // Parameter 1 address
2163 .save b0, GR_SAVE_B0
2164 mov GR_SAVE_B0=b0 // Save b0
2165 };;
2166 .body
2167 { .mib
2168 stfs [GR_Parameter_X] = FR_X // STORE Parameter 1
2169 // on stack
2170 add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address
2171 nop.b 0
2172 }
2173 { .mib
2174 stfs [GR_Parameter_Y] = FR_RESULT // STORE Parameter 3
2175 // on stack
2176 add GR_Parameter_Y = -16,GR_Parameter_Y
2177 br.call.sptk b0=__libm_error_support# // Call error handling
2178 // function
2179 };;
2180 { .mmi
2181 nop.m 0
2182 nop.m 0
2183 add GR_Parameter_RESULT = 48,sp
2184 };;
2185 { .mmi
2186 ldfs f8 = [GR_Parameter_RESULT] // Get return result off stack
2187 .restore sp
2188 add sp = 64,sp // Restore stack pointer
2189 mov b0 = GR_SAVE_B0 // Restore return address
2190 };;
2191 { .mib
2192 mov gp = GR_SAVE_GP // Restore gp
2193 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
2194 br.ret.sptk b0 // Return
2195 };;
2197 LOCAL_LIBM_END(__libm_error_region)
2198 .type __libm_error_support#,@function
2199 .global __libm_error_support#