• Don Clugston (42/42) Sep 22 2005 I have ported the Cephes gamma and lgamma functions to D.
• Walter Bright (22/22) Sep 22 2005 This is good news! Thanks!
• Don Clugston (30/52) Sep 23 2005 I'm certainly no expert on them. Actually, the functions I'm most
• Georg Wrede (5/45) Sep 23 2005 I wish the above were in comments next to the code. This way it would be...
• Don Clugston (5/45) Sep 25 2005 Yes, OK.
• Regan Heath (3/9) Sep 24 2005 So lets use an alias, call it gamma and alias the other name(s).
• ElfQT (7/17) Sep 26 2005 wrote:

Don Clugston <dac nospam.com.au> writes:

I have ported the Cephes gamma and lgamma functions to D.

I've changed gamma() so that it behaves like the C99 tgamma().
eg, for +-0, it returns +-infinity instead of nan.

But, I have a few questions about some of the details:

(a) What should the names be?

C calls it tgamma() to avoid naming conflicts with an older gamma().
D does not have that problem, so I suggest it should be called
real gamma(real x)

Log Gamma is trickier. For reals, it returns log(abs(gamma(x)).
Possible signatures include:
  real lgamma(real x)
  real loggamma(real x)

Personally, I prefer the latter, as it's rather more informative.
(Is the short name just because of identifier length limitations in
some C compilers?)


(b) Behaviour at poles.

gamma(x) has poles at x= 0, -1, -2, -3, ....
At each of these points, the result flips from +inf to -inf.

In C, lgamma(pole) = +inf.
also  tgamma(+-0) = +-inf
but tgamma(-1,-2,...) = nan.

so log(abs(tgamma(x)) != lgamma(x) for x=-1,-2,...
I've read that C is considering changing the behaviour so that
the invariant is valid at all x; tgamma(-2) would return +-inf.

If a complex gamma(z) was provided, it could choose + or -inf based
on the sign of the +-0.0i, so I think it would make sense to return
+inf instead of nan.


(c) Returning the sign of log gamma.

In C, lgamma(x) returns the sign of gamma in a global variable!
(an int called sgamma or sgammal).
The sign can actually be calculated easily, so this could be a seperate 
function. But what would it be called?
Possibly:

real sgamma(real x)
real sgngamma(real x)
real signgamma(real x)

would return +1 or -1, nan if x is a pole.

Another option would be to drop this feature completely.

In C, the contents of that variable don't seem be defined for poles.

I just think some aspects of the gamma function are a mess in C, and 
shouldn't be directly duplicated.

Does anyone have an opinion on any of this?

"Walter Bright" <newshound digitalmars.com> writes:

This is good news! Thanks!

a) I don't know squat about the uses of lgamma or tgamma.

b) D needs an lgamma and a tgamma that match C's behavior, if for no other
reason than to support porting C code to D. They should be named lgamma and
tgamma in D, and should be simple shells to call the C versions.

c) The C99 spec doesn't say what should happen in the special cases. Is
there an authoritative last word reference on how these special cases should
be handled?

d) C99 also says nothing about returning the sign of lgamma. I suspect
that's a kludge in the implementation you're using, a kludge to maintain
conformance to the C99 spec. Clearly, that is a bad idea for D, and I
suggest overloading lgamma with a version that has an out parameter for the
sign.

e) I want to backport your version into C, and put it in the DMC library,
and shell it for D. The D version will stay using a version declaration to
enable it for those platforms that don't have a gamma, or that have a broken
one. This is a general strategy D should follow for functions that should
exist in a C99 library. Note that I need to do this for pow() on Linux
because gcc's pow() doesn't handle special cases correctly.

f) In the future I'm going to have to split std.math into two modules, one
with the implementation and one without (the "header" version). It's getting
too big!

Don Clugston <dac nospam.com.au> writes:

Walter Bright wrote:
 This is good news! Thanks!

 a) I don't know squat about the uses of lgamma or tgamma.

I'm certainly no expert on them. Actually, the functions I'm most 
interested in having from Cephes are the statistical distribution ones 
-- student t and its ilk.

The doc below is a good one, it proposes all of the Cephes functions 
with better names.

www.open-std.org/jtc1/sc22/wg14/www/docs/n1069.pdf

I believe it is correct when it argues that some of these statistical 
functions are more widely used than even some of the trig functions.
Maybe D could have a "std.statistical" or similar. The unit tests, IEEE 
math, and inbuilt complex types make it much easier to implement such 
things in D than C. (eg, gamma looks much better in D than in C!).

 b) D needs an lgamma and a tgamma that match C's behavior, if for no other
 reason than to support porting C code to D. They should be named lgamma and
 tgamma in D, and should be simple shells to call the C versions.

OK. (Sigh) It's a shame when poor names get perpetuated, though.

 c) The C99 spec doesn't say what should happen in the special cases. Is
 there an authoritative last word reference on how these special cases should
 be handled?

No idea. It probably doesn't matter, the only reason I care is that it 
shows up in the unit tests as an inconsistency between tgamma and lgamma.

 d) C99 also says nothing about returning the sign of lgamma. I suspect
 that's a kludge in the implementation you're using, a kludge to maintain
 conformance to the C99 spec. Clearly, that is a bad idea for D, and I
 suggest overloading lgamma with a version that has an out parameter for the
 sign.

I've just discovered that it's an extension, albeit a commonly 
implemented one:

http://www.opengroup.org/onlinepubs/009695399/functions/lgamma.html

I have found one reference that suggests it's a Cephes thing that was 
dumped when it was ported to C99. So I will just abandon the sign.

 e) I want to backport your version into C, and put it in the DMC library,
 and shell it for D. The D version will stay using a version declaration to
 enable it for those platforms that don't have a gamma, or that have a broken
 one. This is a general strategy D should follow for functions that should
 exist in a C99 library. Note that I need to do this for pow() on Linux
 because gcc's pow() doesn't handle special cases correctly.

OK, that makes sense.

 f) In the future I'm going to have to split std.math into two modules, one
 with the implementation and one without (the "header" version). It's getting
 too big!

That seems to be a general issue with D now. The 'strip function bodies' 
  tool seems to be a necessity. And <math.h> is becoming a monster in 
the C/C++ world. IMHO, the category of "mathematics" is far too general.

See also my reply about floating point literals. I think implicit 
conversions from real to creal need to be removed.
In general, a type should only have one possible implicit conversion,
otherwise you need C++ rules for matching.

I've also made complex versions of all the trig functions, but they 
won't work properly while those conversions still exist.

Georg Wrede <georg.wrede nospam.org> writes:

Don Clugston wrote:
 Walter Bright wrote:

 Maybe D could have a "std.statistical" or similar. The unit tests,
 IEEE math, and inbuilt complex types make it much easier to implement
 such things in D than C. (eg, gamma looks much better in D than in
 C!).

Those ought to be put up somewhere side by side, so folks see this.

 b) D needs an lgamma and a tgamma that match C's behavior, if for
 no other reason than to support porting C code to D. They should be
 named lgamma and tgamma in D, and should be simple shells to call
 the C versions.

 
 OK. (Sigh) It's a shame when poor names get perpetuated, though.
 
 c) The C99 spec doesn't say what should happen in the special
 cases. Is there an authoritative last word reference on how these
 special cases should be handled?

 
 No idea. It probably doesn't matter, the only reason I care is that
 it shows up in the unit tests as an inconsistency between tgamma and
 lgamma.
 
 d) C99 also says nothing about returning the sign of lgamma. I
 suspect that's a kludge in the implementation you're using, a
 kludge to maintain conformance to the C99 spec. Clearly, that is a
 bad idea for D, and I suggest overloading lgamma with a version
 that has an out parameter for the sign.

 
 I've just discovered that it's an extension, albeit a commonly 
 implemented one:
 
 http://www.opengroup.org/onlinepubs/009695399/functions/lgamma.html
 
 I have found one reference that suggests it's a Cephes thing that was
  dumped when it was ported to C99. So I will just abandon the sign.

I wish the above were in comments next to the code. This way it would be 
clear to anybody reading the code why it is like it is.

 See also my reply about floating point literals. I think implicit 
 conversions from real to creal need to be removed. In general, a type
 should only have one possible implicit conversion, otherwise you need
 C++ rules for matching.
 
 I've also made complex versions of all the trig functions, but they 
 won't work properly while those conversions still exist.

Should they be put in the code already, as commented-out code?

Don Clugston <dac nospam.com.au> writes:

Georg Wrede wrote:
 Don Clugston wrote:
 
 Walter Bright wrote:

 
 
 Maybe D could have a "std.statistical" or similar. The unit tests,
 IEEE math, and inbuilt complex types make it much easier to implement
 such things in D than C. (eg, gamma looks much better in D than in
 C!).

 
 
 Those ought to be put up somewhere side by side, so folks see this.

They're too long for that. But I'm sure I'll find a good example elsewhere.

 d) C99 also says nothing about returning the sign of lgamma. I
 suspect that's a kludge in the implementation you're using, a
 kludge to maintain conformance to the C99 spec. Clearly, that is a
 bad idea for D, and I suggest overloading lgamma with a version
 that has an out parameter for the sign.


 I've just discovered that it's an extension, albeit a commonly 
 implemented one:

 http://www.opengroup.org/onlinepubs/009695399/functions/lgamma.html

 I have found one reference that suggests it's a Cephes thing that was
  dumped when it was ported to C99. So I will just abandon the sign.

 
 
 I wish the above were in comments next to the code. This way it would be 
 clear to anybody reading the code why it is like it is.

Yes, OK.

 See also my reply about floating point literals. I think implicit 
 conversions from real to creal need to be removed. In general, a type
 should only have one possible implicit conversion, otherwise you need
 C++ rules for matching.

 I've also made complex versions of all the trig functions, but they 
 won't work properly while those conversions still exist.

 
 Should they be put in the code already, as commented-out code?

Maybe if I submit them all, Walter will be more likely to do something 
about the situation :-)

"Regan Heath" <regan netwin.co.nz> writes:

On Fri, 23 Sep 2005 09:30:57 +0200, Don Clugston <dac nospam.com.au> wrote:
 b) D needs an lgamma and a tgamma that match C's behavior, if for no  
 other
 reason than to support porting C code to D. They should be named lgamma  
 and
 tgamma in D, and should be simple shells to call the C versions.

 OK. (Sigh) It's a shame when poor names get perpetuated, though.

So lets use an alias, call it gamma and alias the other name(s).

Regan

"ElfQT" <dethjunk yahoo.com> writes:

"Regan Heath" <regan netwin.co.nz> wrote in message
news:opsxmi98dr23k2f5 nrage.netwin.co.nz...
 On Fri, 23 Sep 2005 09:30:57 +0200, Don Clugston <dac nospam.com.au>

wrote:
 b) D needs an lgamma and a tgamma that match C's behavior, if for no
 other
 reason than to support porting C code to D. They should be named lgamma
 and
 tgamma in D, and should be simple shells to call the C versions.

 OK. (Sigh) It's a shame when poor names get perpetuated, though.

 So lets use an alias, call it gamma and alias the other name(s).

 Regan

And, of course, with the great new DDoc, clearly indentify the purpose and
usage in DDoc comments.
Hmm, does aliases have the right DDoc feature for this? This case shows
there should be.