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?

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!

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.

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?

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

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> 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.