00001
00020 #include "mfcpch.h"
00021 #include <stdio.h>
00022 #include <math.h>
00023 #include "errcode.h"
00024 #include "quadlsq.h"
00025
00026 const ERRCODE EMPTY_QLSQ = "Can't delete from an empty QLSQ";
00027
00028 #define EXTERN
00029
00033 void QLSQ::clear() {
00034 a = 0;
00035 b = 0;
00036 c = 0;
00037 n = 0;
00038 sigx = 0;
00039 sigy = 0;
00040 sigxx = 0;
00041 sigxy = 0;
00042 sigyy = 0;
00043 sigxxx = 0;
00044 sigxxy = 0;
00045 sigxxxx = 0;
00046 }
00047
00048
00052 void QLSQ::add(
00053 double x,
00054 double y
00055 ) {
00056 n++;
00057 sigx += x;
00058 sigy += y;
00059 sigxx += x * x;
00060 sigxy += x * y;
00061 sigyy += y * y;
00062 sigxxx += (long double) x *x * x;
00063 sigxxy += (long double) x *x * y;
00064 sigxxxx += (long double) x *x * x * x;
00065 }
00066
00067
00071 void QLSQ::remove(
00072 double x,
00073 double y
00074 ) {
00075 if (n <= 0)
00076
00077 EMPTY_QLSQ.error ("QLSQ::remove", ABORT, NULL);
00078 n--;
00079 sigx -= x;
00080 sigy -= y;
00081 sigxx -= x * x;
00082 sigxy -= x * y;
00083 sigyy -= y * y;
00084 sigxxx -= (long double) x *x * x;
00085 sigxxy -= (long double) x *x * y;
00086 sigxxxx -= (long double) x *x * x * x;
00087 }
00088
00089
00093 void QLSQ::fit(
00094 int degree
00095 ) {
00096 long double cubetemp;
00097 long double squaretemp;
00098 long double top96, bottom96;
00099
00100 if (n >= 4 && degree >= 2) {
00101 cubetemp = sigxxx * n - (long double) sigxx *sigx;
00102
00103 top96 =
00104 cubetemp * ((long double) sigxy * n - (long double) sigx * sigy);
00105
00106 squaretemp = (long double) sigxx *n - (long double) sigx *sigx;
00107
00108 top96 += squaretemp * ((long double) sigxx * sigy - sigxxy * n);
00109
00110 bottom96 = cubetemp * cubetemp;
00111
00112 bottom96 -= squaretemp * (sigxxxx * n - (long double) sigxx * sigxx);
00113
00114 a = top96 / bottom96;
00115
00116 top96 = ((long double) sigxx * sigx - sigxxx * n) * a
00117 + (long double) sigxy *n - (long double) sigx *sigy;
00118 bottom96 = (long double) sigxx *n - (long double) sigx *sigx;
00119 b = top96 / bottom96;
00120
00121 c = (sigy - a * sigxx - b * sigx) / n;
00122 }
00123 else if (n == 0 || degree < 0) {
00124 a = b = c = 0;
00125 }
00126 else {
00127 a = 0;
00128 if (n > 1 && degree > 0) {
00129 b = (sigxy * n - sigx * sigy) / (sigxx * n - sigx * sigx);
00130 }
00131 else
00132 b = 0;
00133 c = (sigy - b * sigx) / n;
00134 }
00135 }