# Making some minor contributions to GSL

linas@austin.ibm.com linas@austin.ibm.com
Sat Jan 17 00:24:00 GMT 2004

```On Thu, Jan 15, 2004 at 05:57:31PM -0700, Gerard Jungman wrote:
> On Tue, 2004-01-13 at 17:57, linas@austin.ibm.com wrote:
>
> > The only "deficiency" is that the GSL implementation generates
> > a "domain error" when gsl_sf_psi_n (1,x) is called with
> > negative x.  The psi's should be perfectly well defined
> > for negative, non-integer x's.
>
> It might take me a while to get around to that. I don't think
> I ever found a good way, or at least not one that I understood.

?
I'm not sure what you mean.  For small negative values,
abram & steg 6.4.6 can be applied a few times and for larger
values, 6.4.7.  I'm not sure, I guess you must be worried about
precision/rounding or performance or something like that?
Or maybe large values of n, for which this isn't practical?

>     gsl_sf_zetam1_e(double s, gsl_sf_result * r);

excellent!

--linas

p.s. I don't know if its of interest, but just in case:
the attachment contains some code to compute the determinant
of a matrix.  It wasn't written to fit into gsl, but could
probably be modified to do so easily enough.

-------------- next part --------------

/* This is a quick-n-dirty implementation for computing the determinant
* of a matrix.
*/

#define MS 100

typedef long double matrix[MS][MS];
typedef long double vector[MS];

typedef long double my_dbl;

/* alternate() returns the determinant of a submatrix of matrix m.
*   vlen: input, size of submatrix within m
*   ir, ic: vectors containing row numbers and column numbers
*         from which the alternating product should be composed.
*
*  The implementation is recursive: the alternating product for
*  larger matrices is expressed in terms of smaller ones.  The
*  cases for the first 5 orders are computed explicitly, to minimize
*  the overhead of a recursive call for these smaller sizes.
*
*  The code for the first 5 cases was generated automatically
*  to ensure correctness, i.e. to avoid typing mistakes.  See
*  the ifdef'd portion below for the generator.
*
*  Note that execution time grows like N factorial. Thus n=14
*  might take less than an hour, n=15 might take more, and n=16
*  might take days.
*/
my_dbl
alternate (matrix *m, vector *ir, vector *ic, int vlen)
{
my_dbl det;
if (2 == vlen)
{
int ra, rb, ca, cb;
ra = (*ir)[0];
rb = (*ir)[1];

ca = (*ic)[0];
cb = (*ic)[1];

det = (*m)[ra][ca] * (*m)[rb][cb];
det -= (*m)[ra][cb] * (*m)[rb][ca];
return det;
}
else if (3 == vlen)
{
int r0, r1, r2, c0, c1, c2;
r0 = (*ir)[0];
r1 = (*ir)[1];
r2 = (*ir)[2];

c0 = (*ic)[0];
c1 = (*ic)[1];
c2 = (*ic)[2];

det =
(  ((*m)[r0][c0]) * (((*m)[r1][c1]) * ((*m)[r2][c2]) - ((*m)[r1][c2]) * ((*m)[r2][c1]))
- ((*m)[r0][c1]) * (((*m)[r1][c0]) * ((*m)[r2][c2]) - ((*m)[r1][c2]) * ((*m)[r2][c0]))
+ ((*m)[r0][c2]) * (((*m)[r1][c0]) * ((*m)[r2][c1]) - ((*m)[r1][c1]) * ((*m)[r2][c0]))
);
return det;
}
else if (4 == vlen)
{
int r0, r1, r2, r3, c0, c1, c2, c3;
r0 = (*ir)[0];
r1 = (*ir)[1];
r2 = (*ir)[2];
r3 = (*ir)[3];

c0 = (*ic)[0];
c1 = (*ic)[1];
c2 = (*ic)[2];
c3 = (*ic)[3];

det =
(  ((*m)[r0][c0]) * (  ((*m)[r1][c1]) * (((*m)[r2][c2]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c2]))
- ((*m)[r1][c2]) * (((*m)[r2][c1]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c1]))
+ ((*m)[r1][c3]) * (((*m)[r2][c1]) * ((*m)[r3][c2]) - ((*m)[r2][c2]) * ((*m)[r3][c1]))
)
- ((*m)[r0][c1]) * (  ((*m)[r1][c0]) * (((*m)[r2][c2]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c2]))
- ((*m)[r1][c2]) * (((*m)[r2][c0]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c0]))
+ ((*m)[r1][c3]) * (((*m)[r2][c0]) * ((*m)[r3][c2]) - ((*m)[r2][c2]) * ((*m)[r3][c0]))
)
+ ((*m)[r0][c2]) * (  ((*m)[r1][c0]) * (((*m)[r2][c1]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c1]))
- ((*m)[r1][c1]) * (((*m)[r2][c0]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c0]))
+ ((*m)[r1][c3]) * (((*m)[r2][c0]) * ((*m)[r3][c1]) - ((*m)[r2][c1]) * ((*m)[r3][c0]))
)
- ((*m)[r0][c3]) * (  ((*m)[r1][c0]) * (((*m)[r2][c1]) * ((*m)[r3][c2]) - ((*m)[r2][c2]) * ((*m)[r3][c1]))
- ((*m)[r1][c1]) * (((*m)[r2][c0]) * ((*m)[r3][c2]) - ((*m)[r2][c2]) * ((*m)[r3][c0]))
+ ((*m)[r1][c2]) * (((*m)[r2][c0]) * ((*m)[r3][c1]) - ((*m)[r2][c1]) * ((*m)[r3][c0]))
));
return det;
}
else if (5 == vlen)
{
int r0, r1, r2, r3, r4, c0, c1, c2, c3, c4;
r0 = (*ir)[0];
r1 = (*ir)[1];
r2 = (*ir)[2];
r3 = (*ir)[3];
r4 = (*ir)[4];

c0 = (*ic)[0];
c1 = (*ic)[1];
c2 = (*ic)[2];
c3 = (*ic)[3];
c4 = (*ic)[4];

det =
(  ((*m)[r0][c0]) * (  ((*m)[r1][c1]) * (  ((*m)[r2][c2]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
- ((*m)[r2][c3]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
+ ((*m)[r2][c4]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
)
- ((*m)[r1][c2]) * (  ((*m)[r2][c1]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
- ((*m)[r2][c3]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
+ ((*m)[r2][c4]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
)
+ ((*m)[r1][c3]) * (  ((*m)[r2][c1]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
- ((*m)[r2][c2]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
+ ((*m)[r2][c4]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
)
- ((*m)[r1][c4]) * (  ((*m)[r2][c1]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
- ((*m)[r2][c2]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
+ ((*m)[r2][c3]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
)
)
- ((*m)[r0][c1]) * (  ((*m)[r1][c0]) * (  ((*m)[r2][c2]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
- ((*m)[r2][c3]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
+ ((*m)[r2][c4]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
)
- ((*m)[r1][c2]) * (  ((*m)[r2][c0]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
- ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
+ ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
)
+ ((*m)[r1][c3]) * (  ((*m)[r2][c0]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
- ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
+ ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
)
- ((*m)[r1][c4]) * (  ((*m)[r2][c0]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
- ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
+ ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
)
)
+ ((*m)[r0][c2]) * (  ((*m)[r1][c0]) * (  ((*m)[r2][c1]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
- ((*m)[r2][c3]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
+ ((*m)[r2][c4]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
)
- ((*m)[r1][c1]) * (  ((*m)[r2][c0]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
- ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
+ ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
)
+ ((*m)[r1][c3]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
- ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
+ ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
)
- ((*m)[r1][c4]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
- ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
+ ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
)
)
- ((*m)[r0][c3]) * (  ((*m)[r1][c0]) * (  ((*m)[r2][c1]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
- ((*m)[r2][c2]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
+ ((*m)[r2][c4]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
)
- ((*m)[r1][c1]) * (  ((*m)[r2][c0]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
- ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
+ ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
)
+ ((*m)[r1][c2]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
- ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
+ ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
)
- ((*m)[r1][c4]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
- ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
+ ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
)
)
+ ((*m)[r0][c4]) * (  ((*m)[r1][c0]) * (  ((*m)[r2][c1]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
- ((*m)[r2][c2]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
+ ((*m)[r2][c3]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
)
- ((*m)[r1][c1]) * (  ((*m)[r2][c0]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
- ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
+ ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
)
+ ((*m)[r1][c2]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
- ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
+ ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
)
- ((*m)[r1][c3]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
- ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
+ ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
)));
return det;
}
else
{
vector row, col;

int rrr = (*ir)[0];

// skip the top row when working submatrices
int i,j;
for (i=0; i<vlen-1; i++)
{
row[i] = (*ir)[i+1];
}

my_dbl sign = 1.0L;
my_dbl acc = 0.0L;
for (i=0; i<vlen; i++)
{
int k=0;
for (j=0; j<vlen-1; j++)
{
if (k==i) k++;  // skip this column
col[j] = (*ic)[k];
k++;
}

int ccc = (*ic)[i];
acc += sign * (*m)[rrr][ccc] * alt (m, &row, &col, vlen-1);
sign = -sign;
}
return acc;
}
}

/* Return determinant of matrix m, of dimension dim */

my_dbl
determinant (matrix *m, int dim)
{
int i;
vector row, col;

for (i=0; i<dim; i++)
{
row[i] = i;
col[i] = i;
}

return alternate (m, &row, &col, dim);
}

#ifdef TEST
main ()
{
matrix m;
int  i, j;
for (i=0; i<5; i++)
{
for (j=0; j<5; j++)
{
// m[i][j] = j*i;
m[i][j] = 0.0;
}
m[i][i] = 1.0;
}

printf ("its %Lg\n", determinant (&m, 5));
return 0;
}
#endif

/* ===================================================================== */
#ifdef DET_GENERATOR

typedef int vector[30];

/* This code generates the source code for the lower orders.
* Its straightforward to audit
*/
void
alt (vector *ir, vector *ic, int vlen)
{
if (2 == vlen)
{
int ra, rb, ca, cb;
ra = (*ir)[0];
rb = (*ir)[1];

ca = (*ic)[0];
cb = (*ic)[1];

printf ("(((*m)[r%d][c%d]) * ((*m)[r%d][c%d]) - ((*m)[r%d][c%d]) * ((*m)[r%d][c%d]))\n",
ra,ca,rb,cb, ra,cb,rb,ca);
}
else
{
vector row, col;

int rrr = (*ir)[0];
int i,j;
for (i=0; i<vlen-1; i++)
{
row[i] = (*ir)[i+1];
}

int sign = 1;
printf ("(  ");
for (i=0; i<vlen; i++)
{
int k=0;
for (j=0; j<vlen-1; j++)
{
if (k==i) k++;  // skip this column
col[j] = (*ic)[k];
k++;
}

int ccc = (*ic)[i];
if (i != 0) {
if (sign > 0) printf (" + ");
else printf (" - ");
}
printf ("((*m)[r%d][c%d]) * ", rrr, ccc);
alt (&row, &col, vlen-1);
sign = -sign;
}
printf (")\n");
}
}

void
det_gen (int dim)
{
int i;
vector row, col;

for (i=0; i<dim; i++)
{
row[i] = i;
col[i] = i;
}

alt (&row, &col, dim);
}

main ()
{
det_gen (5);
return 0;
}
#endif

```