Updated version 1.9.1+repack from 'upstream/1.9.1+repack'

with Debian dir 98a996367aa69ae41accf9c6d369f600bc94de80

1 files changed, 865 insertions, 0 deletions

diff --git a/spectro/sa_conv.c b/spectro/sa_conv.c
new file mode 100644
index 0000000..84286ba
--- /dev/null
+++ b/spectro/sa_conv.c
@@ -0,0 +1,865 @@
+
+#ifdef SALONEINSTLIB
+
+/*
+ * A very small subset of icclib, copied to here.
+ * This is just enough to support the standalone instruments
+ */
+
+/* 
+ * Argyll Color Correction System
+ *
+ * Author: Graeme W. Gill
+ * Date:   28/9/97
+ *
+ * Copyright 1997 - 2013 Graeme W. Gill
+ * All rights reserved.
+ *
+ * This material is licenced under the GNU GENERAL PUBLIC LICENSE Version 2 or later :-
+ * see the License2.txt file for licencing details.
+ */
+
+#include "sa_config.h"
+#include "numsup.h"
+#include "sa_conv.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+
+sa_XYZNumber sa_D50 = {
+    0.9642, 1.0000, 0.8249
+};
+
+sa_XYZNumber sa_D65 = {
+    0.9505, 1.0000, 1.0890
+};
+
+sa_XYZNumber sa_D50_100 = {
+    96.42, 100.00, 82.49
+};
+
+sa_XYZNumber sa_D65_100 = {
+    95.05, 100.00, 108.90
+};
+
+unsigned int sa_CSSig2nchan(icColorSpaceSignature sig) {
+	switch(sig) {
+		case icSigXYZData:
+			return 3;
+		case icSigLabData:
+			return 3;
+		case icSigLuvData:
+			return 3;
+		case icSigYCbCrData:
+			return 3;
+		case icSigYxyData:
+			return 3;
+		case icSigRgbData:
+			return 3;
+		case icSigGrayData:
+			return 1;
+		case icSigHsvData:
+			return 3;
+		case icSigHlsData:
+			return 3;
+		case icSigCmykData:
+			return 4;
+		case icSigCmyData:
+			return 3;
+		case icSig2colorData:
+			return 2;
+		case icSig3colorData:
+			return 3;
+		case icSig4colorData:
+			return 4;
+		case icSig5colorData:
+		case icSigMch5Data:
+			return 5;
+		case icSig6colorData:
+		case icSigMch6Data:
+			return 6;
+		case icSig7colorData:
+		case icSigMch7Data:
+			return 7;
+		case icSig8colorData:
+		case icSigMch8Data:
+			return 8;
+		case icSig9colorData:
+			return 9;
+		case icSig10colorData:
+			return 10;
+		case icSig11colorData:
+			return 11;
+		case icSig12colorData:
+			return 12;
+		case icSig13colorData:
+			return 13;
+		case icSig14colorData:
+			return 14;
+		case icSig15colorData:
+			return 15;
+
+#ifdef NEVER
+		/* Non-standard and Pseudo spaces */
+		case icmSigYData:
+			return 1;
+		case icmSigLData:
+			return 1;
+		case icmSigL8Data:
+			return 1;
+		case icmSigLV2Data:
+			return 1;
+		case icmSigLV4Data:
+			return 1;
+		case icmSigPCSData:
+			return 3;
+		case icmSigLab8Data:
+			return 3;
+		case icmSigLabV2Data:
+			return 3;
+		case icmSigLabV4Data:
+			return 3;
+#endif /* NEVER */
+
+		default:
+			break;
+	}
+	return 0;
+}
+
+void sa_SetUnity3x3(double mat[3][3]) {
+	int i, j;
+	for (j = 0; j < 3; j++) {
+		for (i = 0; i < 3; i++) {
+			if (i == j)
+				mat[j][i] = 1.0;
+			else
+				mat[j][i] = 0.0;
+		}
+	}
+	
+}
+
+void sa_Cpy3x3(double dst[3][3], double src[3][3]) {
+	int i, j;
+
+	for (j = 0; j < 3; j++) {
+		for (i = 0; i < 3; i++)
+			dst[j][i] = src[j][i];
+	}
+}
+
+void sa_MulBy3x3(double out[3], double mat[3][3], double in[3]) {
+	double tt[3];
+
+	tt[0] = mat[0][0] * in[0] + mat[0][1] * in[1] + mat[0][2] * in[2];
+	tt[1] = mat[1][0] * in[0] + mat[1][1] * in[1] + mat[1][2] * in[2];
+	tt[2] = mat[2][0] * in[0] + mat[2][1] * in[1] + mat[2][2] * in[2];
+
+	out[0] = tt[0];
+	out[1] = tt[1];
+	out[2] = tt[2];
+}
+
+void sa_Mul3x3_2(double dst[3][3], double src1[3][3], double src2[3][3]) {
+	int i, j, k;
+	double td[3][3];		/* Temporary dest */
+
+	for (j = 0; j < 3; j++) {
+		for (i = 0; i < 3; i++) {
+			double tt = 0.0;
+			for (k = 0; k < 3; k++)
+				tt += src1[j][k] * src2[k][i];
+			td[j][i] = tt;
+		}
+	}
+
+	/* Copy result out */
+	for (j = 0; j < 3; j++)
+		for (i = 0; i < 3; i++)
+			dst[j][i] = td[j][i];
+}
+
+
+/* Matrix Inversion by Richard Carling from "Graphics Gems", Academic Press, 1990 */
+#define det2x2(a, b, c, d) (a * d - b * c)
+
+static void adjoint(
+double out[3][3],
+double in[3][3]
+) {
+    double a1, a2, a3, b1, b2, b3, c1, c2, c3;
+
+    /* assign to individual variable names to aid  */
+    /* selecting correct values  */
+
+	a1 = in[0][0]; b1 = in[0][1]; c1 = in[0][2];
+	a2 = in[1][0]; b2 = in[1][1]; c2 = in[1][2];
+	a3 = in[2][0]; b3 = in[2][1]; c3 = in[2][2];
+
+    /* row column labeling reversed since we transpose rows & columns */
+
+    out[0][0]  =   det2x2(b2, b3, c2, c3);
+    out[1][0]  = - det2x2(a2, a3, c2, c3);
+    out[2][0]  =   det2x2(a2, a3, b2, b3);
+        
+    out[0][1]  = - det2x2(b1, b3, c1, c3);
+    out[1][1]  =   det2x2(a1, a3, c1, c3);
+    out[2][1]  = - det2x2(a1, a3, b1, b3);
+        
+    out[0][2]  =   det2x2(b1, b2, c1, c2);
+    out[1][2]  = - det2x2(a1, a2, c1, c2);
+    out[2][2]  =   det2x2(a1, a2, b1, b2);
+}
+
+static double sa_Det3x3(double in[3][3]) {
+    double a1, a2, a3, b1, b2, b3, c1, c2, c3;
+    double ans;
+
+	a1 = in[0][0]; b1 = in[0][1]; c1 = in[0][2];
+	a2 = in[1][0]; b2 = in[1][1]; c2 = in[1][2];
+	a3 = in[2][0]; b3 = in[2][1]; c3 = in[2][2];
+
+    ans = a1 * det2x2(b2, b3, c2, c3)
+        - b1 * det2x2(a2, a3, c2, c3)
+        + c1 * det2x2(a2, a3, b2, b3);
+    return ans;
+}
+
+#define SA__SMALL_NUMBER	1.e-8
+
+int sa_Inverse3x3(double out[3][3], double in[3][3]) {
+    int i, j;
+    double det;
+
+    /*  calculate the 3x3 determinant
+     *  if the determinant is zero, 
+     *  then the inverse matrix is not unique.
+     */
+    det = sa_Det3x3(in);
+
+    if ( fabs(det) < SA__SMALL_NUMBER)
+        return 1;
+
+    /* calculate the adjoint matrix */
+    adjoint(out, in);
+
+    /* scale the adjoint matrix to get the inverse */
+    for (i = 0; i < 3; i++)
+        for(j = 0; j < 3; j++)
+		    out[i][j] /= det;
+	return 0;
+}
+
+#undef SA__SMALL_NUMBER	
+#undef det2x2
+
+/* - - - - - - - - - - - - - - - - - - - - - - - - */
+/* Transpose a 3x3 matrix */
+void sa_Transpose3x3(double out[3][3], double in[3][3]) {
+	int i, j;
+	if (out != in) {
+		for (i = 0; i < 3; i++)
+			for (j = 0; j < 3; j++)
+				out[i][j] = in[j][i];
+	} else {
+		double tt[3][3];
+		for (i = 0; i < 3; i++)
+			for (j = 0; j < 3; j++)
+				tt[i][j] = in[j][i];
+		for (i = 0; i < 3; i++)
+			for (j = 0; j < 3; j++)
+				out[i][j] = tt[i][j];
+	}
+}
+
+/* Scale a 3 vector by the given ratio */
+void sa_Scale3(double out[3], double in[3], double rat) {
+	out[0] = in[0] * rat;
+	out[1] = in[1] * rat;
+	out[2] = in[2] * rat;
+}
+
+/* Clamp a 3 vector to be +ve */
+void sa_Clamp3(double out[3], double in[3]) {
+	int i;
+	for (i = 0; i < 3; i++)
+		out[i] = in[i] < 0.0 ? 0.0 : in[i];
+}
+
+/* Return the normal Delta E given two Lab values */
+double sa_LabDE(double *Lab0, double *Lab1) {
+	double rv = 0.0, tt;
+
+	tt = Lab0[0] - Lab1[0];
+	rv += tt * tt;
+	tt = Lab0[1] - Lab1[1];
+	rv += tt * tt;
+	tt = Lab0[2] - Lab1[2];
+	rv += tt * tt;
+
+	return sqrt(rv);
+}
+
+/* Return the CIE94 Delta E color difference measure, squared */
+double sa_CIE94sq(double Lab0[3], double Lab1[3]) {
+	double desq, dhsq;
+	double dlsq, dcsq;
+	double c12;
+
+	{
+		double dl, da, db;
+		dl = Lab0[0] - Lab1[0];
+		dlsq = dl * dl;		/* dl squared */
+		da = Lab0[1] - Lab1[1];
+		db = Lab0[2] - Lab1[2];
+
+		/* Compute normal Lab delta E squared */
+		desq = dlsq + da * da + db * db;
+	}
+
+	{
+		double c1, c2, dc;
+
+		/* Compute chromanance for the two colors */
+		c1 = sqrt(Lab0[1] * Lab0[1] + Lab0[2] * Lab0[2]);
+		c2 = sqrt(Lab1[1] * Lab1[1] + Lab1[2] * Lab1[2]);
+		c12 = sqrt(c1 * c2);	/* Symetric chromanance */
+
+		/* delta chromanance squared */
+		dc = c1 - c2;
+		dcsq = dc * dc;
+	}
+
+	/* Compute delta hue squared */
+	if ((dhsq = desq - dlsq - dcsq) < 0.0)
+		dhsq = 0.0;
+	{
+		double sc, sh;
+
+		/* Weighting factors for delta chromanance & delta hue */
+		sc = 1.0 + 0.045 * c12;
+		sh = 1.0 + 0.015 * c12;
+		return dlsq + dcsq/(sc * sc) + dhsq/(sh * sh);
+	}
+}
+
+/* Return the CIE94 Delta E color difference measure */
+double sa_CIE94(double Lab0[3], double Lab1[3]) {
+	return sqrt(sa_CIE94sq(Lab0, Lab1));
+}
+
+/* Return the CIE94 Delta E color difference measure for two XYZ values */
+double sa_XYZCIE94(sa_XYZNumber *w, double *in0, double *in1) {
+	double lab0[3], lab1[3];
+
+	sa_XYZ2Lab(w, lab0, in0);
+	sa_XYZ2Lab(w, lab1, in1);
+	return sqrt(sa_CIE94sq(lab0, lab1));
+}
+
+/* CIE XYZ to perceptual CIE 1976 L*a*b* */
+void
+sa_XYZ2Lab(sa_XYZNumber *w, double *out, double *in) {
+	double X = in[0], Y = in[1], Z = in[2];
+	double x,y,z,fx,fy,fz;
+
+	x = X/w->X;
+	y = Y/w->Y;
+	z = Z/w->Z;
+
+	if (x > 0.008856451586)
+		fx = pow(x,1.0/3.0);
+	else
+		fx = 7.787036979 * x + 16.0/116.0;
+
+	if (y > 0.008856451586)
+		fy = pow(y,1.0/3.0);
+	else
+		fy = 7.787036979 * y + 16.0/116.0;
+
+	if (z > 0.008856451586)
+		fz = pow(z,1.0/3.0);
+	else
+		fz = 7.787036979 * z + 16.0/116.0;
+
+	out[0] = 116.0 * fy - 16.0;
+	out[1] = 500.0 * (fx - fy);
+	out[2] = 200.0 * (fy - fz);
+}
+
+void sa_Lab2XYZ(sa_XYZNumber *w, double *out, double *in) {
+	double L = in[0], a = in[1], b = in[2];
+	double x,y,z,fx,fy,fz;
+
+	fy = (L + 16.0)/116.0;
+	fx = a/500.0 + fy;
+	fz = fy - b/200.0;
+
+	if (fy > 24.0/116.0)
+		y = pow(fy,3.0);
+	else
+		y = (fy - 16.0/116.0)/7.787036979;
+
+	if (fx > 24.0/116.0)
+		x = pow(fx,3.0);
+	else
+		x = (fx - 16.0/116.0)/7.787036979;
+
+	if (fz > 24.0/116.0)
+		z = pow(fz,3.0);
+	else
+		z = (fz - 16.0/116.0)/7.787036979;
+
+	out[0] = x * w->X;
+	out[1] = y * w->Y;
+	out[2] = z * w->Z;
+}
+
+
+void sa_Yxy2XYZ(double *out, double *in) {
+	double Y = in[0];
+	double x = in[1];
+	double y = in[2];
+	double z = 1.0 - x - y;
+	double sum;
+	if (y < 1e-9) {
+		out[0] = out[1] = out[2] = 0.0;
+	} else {
+		sum = Y/y;
+		out[0] = x * sum;
+		out[1] = Y;
+		out[2] = z * sum;
+	}
+}
+
+/* - - - - - - - - - - - - - - - - - - - - - - - - */
+
+/* Object for computing RFC 1321 MD5 checksums. */
+/* Derived from Colin Plumb's 1993 public domain code. */
+
+/* Reset the checksum */
+static void sa_MD5_reset(sa_MD5 *p) {
+	p->tlen = 0;
+
+	p->sum[0] = 0x67452301;
+	p->sum[1] = 0xefcdab89;
+	p->sum[2] = 0x98badcfe;
+	p->sum[3] = 0x10325476;
+
+	p->fin = 0;
+}
+
+#define F1(x, y, z) (z ^ (x & (y ^ z)))
+#define F2(x, y, z) F1(z, x, y)
+#define F3(x, y, z) (x ^ y ^ z)
+#define F4(x, y, z) (y ^ (x | ~z))
+
+#define MD5STEP(f, w, x, y, z, pp, xtra, s) \
+		data = (pp)[0] + ((pp)[3] << 24) + ((pp)[2] << 16) + ((pp)[1] << 8); \
+		w += f(x, y, z) + data + xtra; \
+		w = (w << s) | (w >> (32-s)); \
+		w += x;
+
+/* Add another 64 bytes to the checksum */
+static void sa_MD5_accume(sa_MD5 *p, ORD8 *in) {
+	ORD32 data, a, b, c, d;
+
+	a = p->sum[0];
+	b = p->sum[1];
+	c = p->sum[2];
+	d = p->sum[3];
+
+	MD5STEP(F1, a, b, c, d, in + (4 * 0),  0xd76aa478, 7);
+	MD5STEP(F1, d, a, b, c, in + (4 * 1),  0xe8c7b756, 12);
+	MD5STEP(F1, c, d, a, b, in + (4 * 2),  0x242070db, 17);
+	MD5STEP(F1, b, c, d, a, in + (4 * 3),  0xc1bdceee, 22);
+	MD5STEP(F1, a, b, c, d, in + (4 * 4),  0xf57c0faf, 7);
+	MD5STEP(F1, d, a, b, c, in + (4 * 5),  0x4787c62a, 12);
+	MD5STEP(F1, c, d, a, b, in + (4 * 6),  0xa8304613, 17);
+	MD5STEP(F1, b, c, d, a, in + (4 * 7),  0xfd469501, 22);
+	MD5STEP(F1, a, b, c, d, in + (4 * 8),  0x698098d8, 7);
+	MD5STEP(F1, d, a, b, c, in + (4 * 9),  0x8b44f7af, 12);
+	MD5STEP(F1, c, d, a, b, in + (4 * 10), 0xffff5bb1, 17);
+	MD5STEP(F1, b, c, d, a, in + (4 * 11), 0x895cd7be, 22);
+	MD5STEP(F1, a, b, c, d, in + (4 * 12), 0x6b901122, 7);
+	MD5STEP(F1, d, a, b, c, in + (4 * 13), 0xfd987193, 12);
+	MD5STEP(F1, c, d, a, b, in + (4 * 14), 0xa679438e, 17);
+	MD5STEP(F1, b, c, d, a, in + (4 * 15), 0x49b40821, 22);
+
+	MD5STEP(F2, a, b, c, d, in + (4 * 1),  0xf61e2562, 5);
+	MD5STEP(F2, d, a, b, c, in + (4 * 6),  0xc040b340, 9);
+	MD5STEP(F2, c, d, a, b, in + (4 * 11), 0x265e5a51, 14);
+	MD5STEP(F2, b, c, d, a, in + (4 * 0),  0xe9b6c7aa, 20);
+	MD5STEP(F2, a, b, c, d, in + (4 * 5),  0xd62f105d, 5);
+	MD5STEP(F2, d, a, b, c, in + (4 * 10), 0x02441453, 9);
+	MD5STEP(F2, c, d, a, b, in + (4 * 15), 0xd8a1e681, 14);
+	MD5STEP(F2, b, c, d, a, in + (4 * 4),  0xe7d3fbc8, 20);
+	MD5STEP(F2, a, b, c, d, in + (4 * 9),  0x21e1cde6, 5);
+	MD5STEP(F2, d, a, b, c, in + (4 * 14), 0xc33707d6, 9);
+	MD5STEP(F2, c, d, a, b, in + (4 * 3),  0xf4d50d87, 14);
+	MD5STEP(F2, b, c, d, a, in + (4 * 8),  0x455a14ed, 20);
+	MD5STEP(F2, a, b, c, d, in + (4 * 13), 0xa9e3e905, 5);
+	MD5STEP(F2, d, a, b, c, in + (4 * 2),  0xfcefa3f8, 9);
+	MD5STEP(F2, c, d, a, b, in + (4 * 7),  0x676f02d9, 14);
+	MD5STEP(F2, b, c, d, a, in + (4 * 12), 0x8d2a4c8a, 20);
+
+	MD5STEP(F3, a, b, c, d, in + (4 * 5),  0xfffa3942, 4);
+	MD5STEP(F3, d, a, b, c, in + (4 * 8),  0x8771f681, 11);
+	MD5STEP(F3, c, d, a, b, in + (4 * 11), 0x6d9d6122, 16);
+	MD5STEP(F3, b, c, d, a, in + (4 * 14), 0xfde5380c, 23);
+	MD5STEP(F3, a, b, c, d, in + (4 * 1),  0xa4beea44, 4);
+	MD5STEP(F3, d, a, b, c, in + (4 * 4),  0x4bdecfa9, 11);
+	MD5STEP(F3, c, d, a, b, in + (4 * 7),  0xf6bb4b60, 16);
+	MD5STEP(F3, b, c, d, a, in + (4 * 10), 0xbebfbc70, 23);
+	MD5STEP(F3, a, b, c, d, in + (4 * 13), 0x289b7ec6, 4);
+	MD5STEP(F3, d, a, b, c, in + (4 * 0),  0xeaa127fa, 11);
+	MD5STEP(F3, c, d, a, b, in + (4 * 3),  0xd4ef3085, 16);
+	MD5STEP(F3, b, c, d, a, in + (4 * 6),  0x04881d05, 23);
+	MD5STEP(F3, a, b, c, d, in + (4 * 9),  0xd9d4d039, 4);
+	MD5STEP(F3, d, a, b, c, in + (4 * 12), 0xe6db99e5, 11);
+	MD5STEP(F3, c, d, a, b, in + (4 * 15), 0x1fa27cf8, 16);
+	MD5STEP(F3, b, c, d, a, in + (4 * 2),  0xc4ac5665, 23);
+
+	MD5STEP(F4, a, b, c, d, in + (4 * 0),  0xf4292244, 6);
+	MD5STEP(F4, d, a, b, c, in + (4 * 7),  0x432aff97, 10);
+	MD5STEP(F4, c, d, a, b, in + (4 * 14), 0xab9423a7, 15);
+	MD5STEP(F4, b, c, d, a, in + (4 * 5),  0xfc93a039, 21);
+	MD5STEP(F4, a, b, c, d, in + (4 * 12), 0x655b59c3, 6);
+	MD5STEP(F4, d, a, b, c, in + (4 * 3),  0x8f0ccc92, 10);
+	MD5STEP(F4, c, d, a, b, in + (4 * 10), 0xffeff47d, 15);
+	MD5STEP(F4, b, c, d, a, in + (4 * 1),  0x85845dd1, 21);
+	MD5STEP(F4, a, b, c, d, in + (4 * 8),  0x6fa87e4f, 6);
+	MD5STEP(F4, d, a, b, c, in + (4 * 15), 0xfe2ce6e0, 10);
+	MD5STEP(F4, c, d, a, b, in + (4 * 6),  0xa3014314, 15);
+	MD5STEP(F4, b, c, d, a, in + (4 * 13), 0x4e0811a1, 21);
+	MD5STEP(F4, a, b, c, d, in + (4 * 4),  0xf7537e82, 6);
+	MD5STEP(F4, d, a, b, c, in + (4 * 11), 0xbd3af235, 10);
+	MD5STEP(F4, c, d, a, b, in + (4 * 2),  0x2ad7d2bb, 15);
+	MD5STEP(F4, b, c, d, a, in + (4 * 9),  0xeb86d391, 21);
+
+	p->sum[0] += a;
+	p->sum[1] += b;
+	p->sum[2] += c;
+	p->sum[3] += d;
+}
+
+#undef F1
+#undef F2
+#undef F3
+#undef F4
+#undef MD5STEP
+
+/* Add some bytes */
+static void sa_MD5_add(sa_MD5 *p, ORD8 *ibuf, unsigned int len) {
+	unsigned int bs;
+
+	if (p->fin)
+		return;				/* This is actually an error */
+
+	bs = p->tlen;			/* Current bytes added */
+	p->tlen = bs + len;		/* Update length after adding this buffer */
+	bs &= 0x3f;				/* Bytes already in buffer */
+
+	/* Deal with any existing partial bytes in p->buf */
+	if (bs) {
+		ORD8 *np = (ORD8 *)p->buf + bs;	/* Next free location in partial buffer */
+
+		bs = 64 - bs;		/* Free space in partial buffer */
+
+		if (len < bs) {		/* Not enought new to make a full buffer */
+			memmove(np, ibuf, len);
+			return;
+		}
+
+		memmove(np, ibuf, bs);	/* Now got one full buffer */
+		sa_MD5_accume(p, np);
+		ibuf += bs;
+		len -= bs;
+	}
+
+	/* Deal with input data 64 bytes at a time */
+	while (len >= 64) {
+		sa_MD5_accume(p, ibuf);
+		ibuf += 64;
+		len -= 64;
+	}
+
+	/* Deal with any remaining bytes */
+	memmove(p->buf, ibuf, len);
+}
+
+/* Finalise the checksum and return the result. */
+static void sa_MD5_get(sa_MD5 *p, ORD8 chsum[16]) {
+	int i;
+	unsigned count;
+	ORD32 bits1, bits0;
+	ORD8 *pp;
+
+	if (p->fin == 0) {
+	
+		/* Compute number of bytes processed mod 64 */
+		count = p->tlen & 0x3f;
+
+		/* Set the first char of padding to 0x80.  This is safe since there is
+		   always at least one byte free */
+		pp = p->buf + count;
+		*pp++ = 0x80;
+
+		/* Bytes of padding needed to make 64 bytes */
+		count = 64 - 1 - count;
+
+		/* Pad out to 56 mod 64, allowing 8 bytes for length in bits. */
+		if (count < 8) {	/* Not enough space for padding and length */
+
+			memset(pp, 0, count);
+			sa_MD5_accume(p, p->buf);
+
+			/* Now fill the next block with 56 bytes */
+			memset(p->buf, 0, 56);
+		} else {
+			/* Pad block to 56 bytes */
+			memset(pp, 0, count - 8);
+		}
+
+		/* Compute number of bits */
+		bits1 = 0x7 & (p->tlen >> (32 - 3)); 
+		bits0 = p->tlen << 3;
+
+		/* Append number of bits */
+		p->buf[64 - 8] = bits0 & 0xff; 
+		p->buf[64 - 7] = (bits0 >> 8) & 0xff; 
+		p->buf[64 - 6] = (bits0 >> 16) & 0xff; 
+		p->buf[64 - 5] = (bits0 >> 24) & 0xff; 
+		p->buf[64 - 4] = bits1 & 0xff; 
+		p->buf[64 - 3] = (bits1 >> 8) & 0xff; 
+		p->buf[64 - 2] = (bits1 >> 16) & 0xff; 
+		p->buf[64 - 1] = (bits1 >> 24) & 0xff; 
+
+		sa_MD5_accume(p, p->buf);
+
+		p->fin = 1;
+	}
+
+	/* Return the result, lsb to msb */
+	pp = chsum;
+	for (i = 0; i < 4; i++) {
+		*pp++ = p->sum[i] & 0xff; 
+		*pp++ = (p->sum[i] >> 8) & 0xff; 
+		*pp++ = (p->sum[i] >> 16) & 0xff; 
+		*pp++ = (p->sum[i] >> 24) & 0xff; 
+	}
+}
+
+
+/* Delete the instance */
+static void sa_MD5_del(sa_MD5 *p) {
+
+	/* This object */
+	if (p != NULL)
+		free(p);
+}
+
+/* Create a new MD5 checksumming object, with a reset checksum value */
+/* Return it or NULL if there is an error */
+sa_MD5 *new_sa_MD5() {
+	sa_MD5 *p;
+
+	if ((p = (sa_MD5 *)calloc(1,sizeof(sa_MD5))) == NULL)
+		return NULL;
+
+	p->reset    = sa_MD5_reset;
+	p->add      = sa_MD5_add;
+	p->get      = sa_MD5_get;
+	p->del      = sa_MD5_del;
+
+	p->reset(p);
+
+	return p;
+}
+
+/* - - - - - - - - - - - - - - - - - - - - - - - - */
+/* A sub-set of ludecomp code from numlib          */
+
+int sa_lu_decomp(double **a, int n, int *pivx, double *rip) {
+	int    i, j;
+	double *rscale, RSCALE[10];		
+
+	if (n <= 10)
+		rscale = RSCALE;
+	else
+		rscale = dvector(0, n-1);
+
+	for (i = 0; i < n; i++) {
+		double big;
+		for (big = 0.0, j=0; j < n; j++) {
+			double temp;
+			temp = fabs(a[i][j]);
+			if (temp > big)
+				big = temp;
+		}
+		if (fabs(big) <= DBL_MIN) {
+			if (rscale != RSCALE)
+				free_dvector(rscale, 0, n-1);
+			return 1;		
+		}
+		rscale[i] = 1.0/big;	
+	}
+
+	for (*rip = 1.0, j = 0; j < n; j++) {
+		double big;
+		int k, bigi = 0;
+
+		for (i = 0; i < j; i++) {
+			double sum;
+			sum = a[i][j];
+			for (k = 0; k < i; k++)
+				sum -= a[i][k] * a[k][j];
+			a[i][j] = sum;
+		}
+
+		for (big = 0.0, i = j; i < n; i++) {
+			double sum, temp;
+			sum = a[i][j];
+			for (k = 0; k < j; k++)
+				sum -= a[i][k] * a[k][j];
+			a[i][j] = sum;
+			temp = rscale[i] * fabs(sum);	
+			if (temp >= big) {
+				big = temp;		
+				bigi = i;		
+			}
+		}
+		
+		if (j != bigi) {
+			{	
+				double *temp;
+				temp = a[bigi];
+				a[bigi] = a[j];
+				a[j] = temp;
+			}
+			*rip = -(*rip);				
+			rscale[bigi] = rscale[j];	
+		}
+		
+		pivx[j] = bigi;					
+		if (fabs(a[j][j]) <= DBL_MIN) {
+			if (rscale != RSCALE)
+				free_dvector(rscale, 0, n-1);
+			return 1; 					
+		}
+
+		if (j != (n-1)) {
+			double temp;
+			temp = 1.0/a[j][j];
+			for (i = j+1; i < n; i++)
+				a[i][j] *= temp;
+		}
+	}
+	if (rscale != RSCALE)
+		free_dvector(rscale, 0, n-1);
+	return 0;
+}
+
+void sa_lu_backsub(double **a, int n, int *pivx, double  *b) {
+	int i, j;
+	int nvi;		
+	
+	for (nvi = -1, i = 0; i < n; i++) {
+		int px;
+		double sum;
+
+		px = pivx[i];
+		sum = b[px];
+		b[px] = b[i];
+		if (nvi >= 0) {
+			for (j = nvi; j < i; j++)
+				sum -= a[i][j] * b[j];
+		} else {
+			if (sum != 0.0)
+				nvi = i;			
+		}
+		b[i] = sum;
+	}
+
+	for (i = (n-1); i >= 0; i--) {
+		double sum;
+		sum = b[i];
+		for (j = i+1; j < n; j++)
+			sum -= a[i][j] * b[j];
+		b[i] = sum/a[i][i];
+	}
+}
+
+int sa_lu_invert(double **a, int n) {
+	int i, j;
+	double rip;		
+	int *pivx, PIVX[10];
+	double **y;
+
+	if (n <= 10)
+		pivx = PIVX;
+	else
+		pivx = ivector(0, n-1);
+
+	if (sa_lu_decomp(a, n, pivx, &rip)) {
+		if (pivx != PIVX)
+			free_ivector(pivx, 0, n-1);
+		return 1;
+	}
+
+	y = dmatrix(0, n-1, 0, n-1);
+	for (i = 0; i < n; i++) {
+		for (j = 0; j < n; j++) {
+			y[i][j] = a[i][j];
+		}
+	}
+
+	for (i = 0; i < n; i++) {
+		for (j = 0; j < n; j++)
+			a[i][j] = 0.0;
+		a[i][i] = 1.0;
+		sa_lu_backsub(y, n, pivx, a[i]);
+	}
+
+	free_dmatrix(y, 0, n-1, 0, n-1);
+	if (pivx != PIVX)
+		free_ivector(pivx, 0, n-1);
+
+	return 0;
+}
+
+int sa_lu_psinvert(double **out, double **in, int m, int n) {
+	int rv = 0;
+	double **tr;		
+	double **sq;		
+
+	tr = dmatrix(0, n-1, 0, m-1);
+	matrix_trans(tr, in, m,  n);
+
+	if (m > n) {
+		sq = dmatrix(0, n-1, 0, n-1);
+		if ((rv = matrix_mult(sq, n, n, tr, n, m, in, m, n)) == 0) {
+			if ((rv = sa_lu_invert(sq, n)) == 0) {
+				rv = matrix_mult(out, n, m, sq, n, n, tr, n, m);
+			}
+		}
+		free_dmatrix(sq, 0, n-1, 0, n-1);
+	} else {
+		sq = dmatrix(0, m-1, 0, m-1);
+		if ((rv = matrix_mult(sq, m, m, in, m, n, tr, n, m)) == 0) {
+			if ((rv = sa_lu_invert(sq, m)) == 0) {
+				rv = matrix_mult(out, n, m, tr, n, m, sq, m, m);
+			}
+		}
+		free_dmatrix(sq, 0, m-1, 0, m-1);
+	}
+
+	free_dmatrix(tr, 0, n-1, 0, m-1);
+	return rv;
+}
+
+
+#endif /* SALONEINSTLIB */
+
+