summaryrefslogtreecommitdiff
path: root/spectro/sa_conv.c
diff options
context:
space:
mode:
Diffstat (limited to 'spectro/sa_conv.c')
-rw-r--r--spectro/sa_conv.c865
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 */
+
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-01 14:01:53 +0000