/***************************************************************************
           HypbSolSrch1.cpp  -  search for solutions on H^5 x S^2
                             -------------------
    begin                : Fri 4 March 2011
    copyright            : (C) 2011 by Chris Austin
    email                : chris@chrisaustin.info
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <iostream>
#include <cmath>

using namespace std;

struct doubpow
{
  double value;
  doubpow(double v)
  {
    value = v;
  }
  operator double()
  {
    return value;
  }
  double operator ^(double c)
  {
    return pow(value, c);
  }
};

int main(int argc, char *argv[])
{
  double Xa, X2, X3, qua, qu2, qu3, qfa, qf2, qf3;
  double Cdb, Fdb, Hdb, M, Lc, Jc, minval, maxval, nval, multval;
  int nC, nF, nH;
  cout << "Enter minval:" << endl;
  cin >> minval;
  cout << "Enter maxval:" << endl;
  cin >> maxval;
  cout << "Enter number of values:" << endl;
  cin >> nval;
  multval = pow(maxval/minval,1/(nval - 1));
  cout << "minval = " << minval << " maxval = " << maxval;
  cout << " nval = " << nval << " multval = " << multval << endl;
  nC = 0;
  Cdb = minval;
  while ( nC < nval )
  {
    nF = 0;
    Fdb = minval;
    while ( nF < nval )
    {
      nH = 0;
      Hdb = minval;
      while ( nH < nval )
      {
        doubpow C(Cdb), F(Fdb), H(Hdb);
        Xa = (20)*(C^6)*(H) + (30)*(C^8)*(F) + (8)*(C^7);
        X2 = - (10)*(C^6)*(H) + (15)*(C^8)*(F) - (120)*(C^8) + (20)*(C^7);
        X3 = (10)*(C^6)*(H) - (5)*(C^8)*(F) - (40)*(C^8);
        qfa = (66465)*(H^4) + (88290)*(C^2)*(F)*(H^3) + (78480)*(C^2)*(H^3) - (26700)*(C)*(H^3) +
        (88605)*(C^4)*(F^2)*(H^2) + (171840)*(C^4)*(F)*(H^2) - (17640)*(C^3)*(F)*(H^2) +
        (71600)*(C^4)*(H^2) - (13440)*(C^3)*(H^2) + (6420)*(C^2)*(H^2) + (21060)*(C^6)*(F^3)*(H) +
        (64800)*(C^6)*(F^2)*(H) - (8100)*(C^5)*(F^2)*(H) + (60480)*(C^6)*(F)*(H) -
        (14400)*(C^5)*(F)*(H) + (600)*(C^4)*(F)*(H) + (15360)*(C^6)*(H) - (4800)*(C^5)*(H) +
        (320)*(C^4)*(H) - (480)*(C^3)*(H) + (18225)*(C^8)*(F^4) + (77760)*(C^8)*(F^3) +
        (116640)*(C^8)*(F^2) + (1260)*(C^6)*(F^2) + (69120)*(C^8)*(F) + (1920)*(C^6)*(F) +
        (11520)*(C^8) + (320)*(C^6) - (48)*(C^4);
        qf2 = (20889)*(H^4) + (29430)*(C^2)*(F)*(H^3) + (23544)*(C^2)*(H^3) - (7476)*(C)*(H^3) +
        (30609)*(C^4)*(F^2)*(H^2) + (55848)*(C^4)*(F)*(H^2) - (5544)*(C^3)*(F)*(H^2) +
        (20048)*(C^4)*(H^2) - (3360)*(C^3)*(H^2) + (1284)*(C^2)*(H^2) + (7452)*(C^6)*(F^3)*(H) +
        (22032)*(C^6)*(F^2)*(H) - (2700)*(C^5)*(F^2)*(H) + (19008)*(C^6)*(F)*(H) -
        (4320)*(C^5)*(F)*(H)  + (168)*(C^4)*(F)*(H) + (3840)*(C^6)*(H) - (960)*(C^5)*(H) +
        (32)*(C^4)*(H) + (96)*(C^3)*(H) +  (6561)*(C^8)*(F^4) + (27216)*(C^8)*(F^3) +
        (38880)*(C^8)*(F^2) + (396)*(C^6)*(F^2) +  (20736)*(C^8)*(F) + (480)*(C^6)*(F) +
        (2304)*(C^8) - (64)*(C^6) - (48)*(C^4);
        qf3 = (66465)*(H^4) + (49050)*(C^2)*(F)*(H^3) + (65400)*(C^2)*(H^3) - (32040)*(C)*(H^3) +
        (24165)*(C^4)*(F^2)*(H^2) + (64440)*(C^4)*(F)*(H^2) - (10080)*(C^3)*(F)*(H^2) +
        (42960)*(C^4)*(H^2) - (13440)*(C^3)*(H^2) + (10700)*(C^2)*(H^2) + (1620)*(C^6)*(F^3)*(H) +
        (6480)*(C^6)*(F^2)*(H) - (1800)*(C^5)*(F^2)*(H) + (8640)*(C^6)*(F)*(H) -
        (4800)*(C^5)*(F)*(H) +
        (360)*(C^4)*(F)*(H) + (3840)*(C^6)*(H) - (3200)*(C^5)*(H) + (480)*(C^4)*(H) -
        (1920)*(C^3)*(H)
        - (1215)*(C^8)*(F^4) - (6480)*(C^8)*(F^3) - (12960)*(C^8)*(F^2) + (180)*(C^6)*(F^2) -
        (11520)*(C^8)*(F) + (480)*(C^6)*(F) - (3840)*(C^8) + (320)*(C^6) + (144)*(C^4);
      //  cout << "Xa = " << Xa << " X2 = " << X2 << " X3 = " << X3 << endl;
      //  cout << "qfa = " << qfa << " qf2 = " << qf2 << " qf3 = " << qf3 << endl;
        qua = 48*qfa;
        qu2 = 120*qf2;
        qu3 = 24*qf3;
      //  cout << "qua = " << qua << " qu2 = " << qu2 << " qu3 = " << qu3 << endl;
        M = (X2/2 + 3*X3 + Xa)/(2*qu2/4 + 6*qu3/2 + qua);
        Lc = (X2 - M*qu2)/4;
        Jc = (X3 - M*qu3)/2;
        if ( ( M >= 0 ) && ( Lc >= 0 ) && ( Jc >= 0 ) )
        {
          cout << " *** ALL NON-NEGATIVE *** " << endl;
          cout << "C = " << Cdb << " F = " << Fdb << " H = " << Hdb << endl;
          cout << "M = " << M << " Lc = " << Lc << " Jc = " << Jc << endl;
        }
        nH++;
        Hdb *= multval;
      }
      nF++;
      Fdb *= multval;
    }
    nC++;
    Cdb *= multval;
  }
  
  return 0;
}