895 lines
20 KiB
C
895 lines
20 KiB
C
|
|
#ifdef SALONEINSTLIB
|
|
|
|
/*
|
|
* A very small subset of icclib, copied to here.
|
|
* This is just enough to support the standalone instruments
|
|
*/
|
|
|
|
/*
|
|
* Argyll Color Management 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];
|
|
}
|
|
|
|
/* Add two 3 vectors */
|
|
void sa_Add3(double out[3], double in1[3], double in2[3]) {
|
|
out[0] = in1[0] + in2[0];
|
|
out[1] = in1[1] + in2[1];
|
|
out[2] = in1[2] + in2[2];
|
|
}
|
|
|
|
/* 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;
|
|
}
|
|
}
|
|
|
|
/* - - - - - - - - - - - - - - - - - - - - - - - - */
|
|
|
|
/* 32 bit linear congruent generator */
|
|
/* generates number between 0 and 4294967295 */
|
|
/* (From Knuth & H.W.Lewis) */
|
|
#define PSRAND32L(S) ((S) * 1664525L + 1013904223L)
|
|
|
|
/* Return a 32 bit number between 0 and 4294967295 */
|
|
unsigned int
|
|
sa_rand32( /* Return 32 bit random number */
|
|
unsigned int seed /* Optional seed. Non-zero re-initialized with that seed */
|
|
) {
|
|
static unsigned int pval = 12345678;
|
|
|
|
if (seed != 0)
|
|
pval = seed;
|
|
|
|
pval = PSRAND32L(pval);
|
|
|
|
return pval;
|
|
}
|
|
|
|
/* - - - - - - - - - - - - - - - - - - - - - - - - */
|
|
|
|
/* 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 */
|
|
|
|
|