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
More information about the Gsl-discuss
mailing list