The R Project SVN R

Rev

Rev 68947 | Rev 79203 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 68947 Rev 75839
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
- 
 
 3 
 *  Copyright (C) 2000-2018 The R Core Team
3 
 *  Copyright (C) 1998 Ross Ihaka
 4 
 *  Copyright (C) 1998 Ross Ihaka
4 
 *  Copyright (C) 2000-2012 The R Core Team
 - 
 
5 
 *
 5 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 6 
 *  This program is free software; you can redistribute it and/or modify
7 
 *  it under the terms of the GNU General Public License as published by
 7 
 *  it under the terms of the GNU General Public License as published by
8 
 *  the Free Software Foundation; either version 2 of the License, or
 8 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 9 
 *  (at your option) any later version.
Line 70... Line 70...
70 
 
 70 
 
71 
    if (sgn != NULL && x < 0 && fmod(floor(-x), 2.) == 0)
 71 
    if (sgn != NULL && x < 0 && fmod(floor(-x), 2.) == 0)
72 
	*sgn = -1;
 72 
	*sgn = -1;
73 
 
 73 
 
74 
    if (x <= 0 && x == trunc(x)) { /* Negative integer argument */
 74 
    if (x <= 0 && x == trunc(x)) { /* Negative integer argument */
75 
	ML_ERROR(ME_RANGE, "lgamma");
 75 
	// No warning: this is the best answer; was  ML_ERROR(ME_RANGE, "lgamma");
76 
	return ML_POSINF;/* +Inf, since lgamma(x) = log|gamma(x)| */
 76 
	return ML_POSINF;/* +Inf, since lgamma(x) = log|gamma(x)| */
77 
    }
 77 
    }
78 
 
 78 
 
79 
    y = fabs(x);
 79 
    y = fabs(x);
80 
 
 80 
 
Line 82... Line 82...
82 
    if (y <= 10) return log(fabs(gammafn(x)));
 82 
    if (y <= 10) return log(fabs(gammafn(x)));
83 
    /*
 83 
    /*
84 
      ELSE  y = |x| > 10 ---------------------- */
 84 
      ELSE  y = |x| > 10 ---------------------- */
85 
 
 85 
 
86 
    if (y > xmax) {
 86 
    if (y > xmax) {
87 
	ML_ERROR(ME_RANGE, "lgamma");
 87 
	// No warning: +Inf is the best answer
88 
	return ML_POSINF;
 88 
	return ML_POSINF;
89 
    }
 89 
    }
90 
 
 90 
 
91 
    if (x > 0) { /* i.e. y = x > 10 */
 91 
    if (x > 0) { /* i.e. y = x > 10 */
92 
#ifdef IEEE_754
 92 
#ifdef IEEE_754