1 | */* Implementation of gamma function according to ISO C.* |

2 | * Copyright (C) 1997, 1999, 2002, 2004 Free Software Foundation, Inc.* |

3 | * This file is part of the GNU C Library.* |

4 | * Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and* |

5 | * Jakub Jelinek <jj@ultra.linux.cz, 1999.* |

6 | |

7 | * The GNU C Library is free software; you can redistribute it and/or* |

8 | * modify it under the terms of the GNU Lesser General Public* |

9 | * License as published by the Free Software Foundation; either* |

10 | * version 2.1 of the License, or (at your option) any later version.* |

11 | |

12 | * The GNU C Library is distributed in the hope that it will be useful,* |

13 | * but WITHOUT ANY WARRANTY; without even the implied warranty of* |

14 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU* |

15 | * Lesser General Public License for more details.* |

16 | |

17 | * You should have received a copy of the GNU Lesser General Public* |

18 | * License along with the GNU C Library; if not, write to the Free* |

19 | * Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA* |

20 | * 02111-1307 USA. */* |

21 | |

22 | __#include "quadmath-imp.h"__ |

23 | |

24 | |

25 | **__float128** |

26 | tgammaq (**__float128** x) |

27 | { |

28 | */* We don't have a real gamma implementation now. We'll use lgamma* |

29 | * and the exp function. But due to the required boundary* |

30 | * conditions we must check some values separately. */* |

31 | int64_t hx; |

32 | uint64_t lx; |

33 | **__float128** res; |

34 | *int* sign; |

35 | |

36 | GET_FLT128_WORDS64 (hx, lx, x); |

37 | |

38 | **if** (((hx & `0x7fffffffffffffffLL`) | lx) == `0`) |

39 | */* Return value for x == 0 is Inf with divide by zero exception. */* |

40 | **return** `1.0` / x; |

41 | |

42 | **if** (hx < `0` && (uint64_t) hx < `0xffff000000000000ULL` && rintq (x) == x) |

43 | */* Return value for integer x < 0 is NaN with invalid exception. */* |

44 | **return** (x - x) / (x - x); |

45 | |

46 | **if** (hx == `0xffff000000000000ULL` && lx == `0`) |

47 | */* x == -Inf. According to ISO this is NaN. */* |

48 | **return** x - x; |

49 | |

50 | */* XXX FIXME. */* |

51 | res = expq (lgammaq (x)); |

52 | **return** signbitq (x) ? -res : res; |

53 | } |

54 | |