[ChangeSet] Beta function...

[John W. Eaton]

On 18-Mar-2008, David Bateman wrote:

| Ben Abbott wrote:
| >
| > On Mar 18, 2008, at 6:37 AM, David Bateman wrote:
| >
| >> Ben Abbott wrote:
| >>>
| >>> On Mar 17, 2008, at 3:40 PM, Ben Abbott wrote:
| >>>
| >>>> On Monday, March 17, 2008, at 02:56PM, "Jorge Londoño"
| >>>> <jmlon at yahoo.com> wrote:
| >>>>> I just happen to jump into a problem that require me to compute the
| >>>>> Beta function with negative arguments.
| >>>>> If you look at the definition of the beta function, for example here:
| >>>>
| >>>
| >>> I replaced gammaln(x) with log(gamma(x)), modified the help text to
| >>> reflect that inputs must be real, and added a few extra tests.
| >>>
| >>> Ben
| >>
| >> Is that the right (TM) thing to do, as Marco says it increases the risk
| >> of overflows? Shouldn't we rather get lgamma to work correctly for
| >> negative arguments. What about something like the attached patch
| >> instead?
| >>
| >> D.
| >>
| >
| > I agree that changing lgamma is more proper!
| >
| > Thanks
| > Ben
| The attached version might be better as the reason I used the function
| name xclgamma was that I was trying to have a function like
| 
| extern Complex xclgamma (double);
| 
| but that is not possible to use within the current mapper
| framework..

I think it would be fairly easy to make it work.  We currently have
things like this (in dNDArray.h):


  typedef double (*dmapper) (double);
  typedef Complex (*cmapper) (const Complex&);
  typedef bool (*bmapper) (double);

  NDArray map (dmapper fcn) const;
  ComplexNDArray map (cmapper fcn) const;
  boolNDArray map (bmapper fcn) const;

Because functions can't be overloaded on return type, we can't simply
add

  typedef Complex (*cdmapper) (double);

and

  ComplexNDArray map (cdmapper fcn) const;

but we could add the new typedef and

  ComplexNDArray cdmap (cdmapper fcn) const;

However, it might make more sense if all the functions with the simple
name "map" were reserved for functions that take arguments of the same
type as the elements of their class.  So we could use

  typedef double (*dmapper) (double);
  typedef Complex (*cmapper) (double);
  typedef bool (*bmapper) (double);

  NDArray map (dmapper fcn) const;
  ComplexNDArray map (cmapper fcn) const;
  boolNDArray map (bmapper fcn) const;

  typedef Complex (*ccmapper) (const Complex&);

  ComplexNDArray ccmap (ccmapper fcn) const;

instead.  The implementation of these functions is all template based,
so not much work is needed to make them work.  We currently have

  NDArray
  NDArray::map (dmapper fcn) const
  {
    return MArrayN<double>::map<double> (func_ptr (fcn));
  }

  ComplexNDArray
  NDArray::map (cmapper fcn) const
  {
    return MArrayN<double>::map<Complex> (func_ptr (fcn));
  }

  boolNDArray
  NDArray::map (bmapper fcn) const
  {
    return MArrayN<double>::map<bool> (func_ptr (fcn));
  }

With the above changes to the declarations, we would need

  NDArray
  NDArray::map (dmapper fcn) const
  {
    return MArrayN<double>::map<double> (func_ptr (fcn));
  }

  ComplexNDArray
  NDArray::map (cmapper fcn) const
  {
    return MArrayN<double>::map<Complex> (func_ptr (fcn));
  }

  boolNDArray
  NDArray::map (bmapper fcn) const
  {
    return MArrayN<double>::map<bool> (func_ptr (fcn));
  }

  ComplexNDArray
  NDArray::ccmap (ccmapper fcn) const
  {
    return MArrayN<double>::map<Complex> (func_ptr (fcn));
  }


What do you think?  Is this change worth making?

Also, it would be good to have the same mechanism working in the
Sparse classes.  If you don't have time, I could take a look at it and
propose a patch if you could review it.

Thanks,

jwe