/**********************************************************************
 * File:        quadlsq.cpp  (Formerly qlsq.c)
 * Description: Code for least squares approximation of quadratics.
 * Author:		Ray Smith
 * Created:		Wed Oct  6 15:14:23 BST 1993
 *
 * (C) Copyright 1993, Hewlett-Packard Ltd.
 ** Licensed under the Apache License, Version 2.0 (the "License");
 ** you may not use this file except in compliance with the License.
 ** You may obtain a copy of the License at
 ** http://www.apache.org/licenses/LICENSE-2.0
 ** Unless required by applicable law or agreed to in writing, software
 ** distributed under the License is distributed on an "AS IS" BASIS,
 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 ** See the License for the specific language governing permissions and
 ** limitations under the License.
 *
 **********************************************************************/

#include "mfcpch.h"
#include          <stdio.h>
#include          <math.h>
#include          "errcode.h"
#include          "quadlsq.h"

const ERRCODE EMPTY_QLSQ = "Can't delete from an empty QLSQ";

#define EXTERN

/**********************************************************************
 * QLSQ::clear
 *
 * Function to initialize a QLSQ.
 **********************************************************************/

void QLSQ::clear() {  //initialize
  a = 0;
  b = 0;
  c = 0;
  n = 0;                         //no elements
  sigx = 0;                      //update accumulators
  sigy = 0;
  sigxx = 0;
  sigxy = 0;
  sigyy = 0;
  sigxxx = 0;
  sigxxy = 0;
  sigxxxx = 0;
}


/**********************************************************************
 * QLSQ::add
 *
 * Add an element to the accumulator.
 **********************************************************************/

void QLSQ::add(           //add an element
               double x,  //xcoord
               double y   //ycoord
              ) {
  n++;                           //count elements
  sigx += x;                     //update accumulators
  sigy += y;
  sigxx += x * x;
  sigxy += x * y;
  sigyy += y * y;
  sigxxx += (long double) x *x * x;
  sigxxy += (long double) x *x * y;
  sigxxxx += (long double) x *x * x * x;
}


/**********************************************************************
 * QLSQ::remove
 *
 * Delete an element from the acculuator.
 **********************************************************************/

void QLSQ::remove(           //delete an element
                  double x,  //xcoord
                  double y   //ycoord
                 ) {
  if (n <= 0)
                                 //illegal
    EMPTY_QLSQ.error ("QLSQ::remove", ABORT, NULL);
  n--;                           //count elements
  sigx -= x;                     //update accumulators
  sigy -= y;
  sigxx -= x * x;
  sigxy -= x * y;
  sigyy -= y * y;
  sigxxx -= (long double) x *x * x;
  sigxxy -= (long double) x *x * y;
  sigxxxx -= (long double) x *x * x * x;
}


/**********************************************************************
 * QLSQ::fit
 *
 * Fit the given degree of polynomial and store the result.
 **********************************************************************/

void QLSQ::fit(            //fit polynomial
               int degree  //degree to fit
              ) {
  long double cubetemp;          //intermediates
  long double squaretemp;
  long double top96, bottom96;   /*accurate top & bottom */

  if (n >= 4 && degree >= 2) {
    cubetemp = sigxxx * n - (long double) sigxx *sigx;

    top96 =
      cubetemp * ((long double) sigxy * n - (long double) sigx * sigy);

    squaretemp = (long double) sigxx *n - (long double) sigx *sigx;

    top96 += squaretemp * ((long double) sigxx * sigy - sigxxy * n);

    bottom96 = cubetemp * cubetemp;

    bottom96 -= squaretemp * (sigxxxx * n - (long double) sigxx * sigxx);

    a = top96 / bottom96;

    top96 = ((long double) sigxx * sigx - sigxxx * n) * a
      + (long double) sigxy *n - (long double) sigx *sigy;
    bottom96 = (long double) sigxx *n - (long double) sigx *sigx;
    b = top96 / bottom96;

    c = (sigy - a * sigxx - b * sigx) / n;
  }
  else if (n == 0 || degree < 0) {
    a = b = c = 0;
  }
  else {
    a = 0;
    if (n > 1 && degree > 0) {
      b = (sigxy * n - sigx * sigy) / (sigxx * n - sigx * sigx);
    }
    else
      b = 0;
    c = (sigy - b * sigx) / n;
  }
}