double integrate(function f, double a, double b, double epsilon, int rdepth) {
  /* Berechnet die erste Naeherung nach Simpson fuer das Intervall [a,b] */
  double h = b - a;
  double mid = 0.5 * (a + b);
  
  double fa = f.func(a, f.args);
  double fmid = f.func(mid, f.args);
  double fb = f.func(b, f.args);
  
  double simp = (fa + 4*fmid + fb) * h / 6;
  
  /* Leitet die Rekursion ein. Dabei werden alle schon berechneten Funktions-
   * werte weitergegeben, damit die Anzahl der Funktionsaufrufe minimiert wird.
   * Es wird auch der aktuelle Wert des Integrals ueberreicht, damit verglichen
   * werden kann, welchen Einfluss die Verfeinerung auf das Ergebnis hat */
  return recursive_simpson(f, simp, a, mid, b, fa, fmid, fb, epsilon, rdepth);
}

bool QGraph::addFunction(const function& func){
	bool exist=false;
	QTime t;
	
	for (QValueList<function>::iterator it = funclist.begin(); it != funclist.end() && !exist; ++it ){
		if((*it).expression() == func.expression()){
			exist=true;
			(*it)=func;
		}
	}
	
	if(!exist)
		funclist << func;
	update_points();
	this->repaint(false);
	sendStatus(i18n("%1 function added").arg(func.expression()));
	return exist;
}

bool remember_table_entry::is_equal(function const & f) const
{
	GINAC_ASSERT(f.seq.size()==seq.size());
	if (f.gethash()!=hashvalue) return false;
	size_t num = seq.size();
	for (size_t i=0; i<num; ++i)
		if (!seq[i].is_equal(f.seq[i])) return false;
	++last_access = access_counter;
	++successful_hits;
	return true;
}

/**
 * 简单的枚举文件夹内容，返回内容数量
 */
int SimpleEnumFolder(LPCTSTR lpszPath		// 文件夹路径
	, CShellManager* pShellManager			// Shell管理器
	, function<void(LPITEMIDLIST)> filter)	// 过滤器函数
{
	ENSURE(lpszPath != nullptr);
	ASSERT_VALID(pShellManager);

	AFX_SHELLITEMINFO info;
	HRESULT hr = pShellManager->ItemFromPath(lpszPath, info.pidlRel);
	if (FAILED(hr))	{
		return 0;
	}

	int nFolderCount = 0;

	LPSHELLFOLDER pDesktopFolder;
	hr = SHGetDesktopFolder(&pDesktopFolder);

	if (SUCCEEDED(hr)) {

		IShellFolder* psfCurFolder = nullptr;
		hr = pDesktopFolder->BindToObject(info.pidlRel, nullptr, IID_IShellFolder, (LPVOID*)&psfCurFolder);

		LPENUMIDLIST pEnum = nullptr;
		HRESULT hRes = psfCurFolder->EnumObjects(nullptr, (SHCONTF)(SHCONTF_FOLDERS), &pEnum);
		if (SUCCEEDED(hRes) && pEnum != nullptr) {

			DWORD dwFetched = 1;
			LPITEMIDLIST pidlTemp;
			while (pEnum->Next(1, &pidlTemp, &dwFetched) == S_OK && dwFetched) {

				if (!filter._Empty()) {
					LPITEMIDLIST itemID = pShellManager->ConcatenateItem(info.pidlRel, pidlTemp);
					filter(itemID);
					pShellManager->FreeItem(itemID);
				}

				pShellManager->FreeItem(pidlTemp);

				nFolderCount++;
				dwFetched = 0;
			}

			pEnum->Release();
		}

		psfCurFolder->Release();
		pDesktopFolder->Release();
	}

	pShellManager->FreeItem(info.pidlRel);
	return nFolderCount;
}

int find_root(function f, double x1, double x2, double epsilon, double *root) {
  /* Es wird davon ausgegangen, dass eine Nullstelle mit Vorzeichenwechsel
   * im Intervall [x1,x2] vorliegt und die Funktion stetig ist. Es werden die
   * Funktionswerte am Rand und in der Mitte des Intervalls berechnet */
  double f1, f2;
  double mid, fmid;
  
  f1 = f.func(x1, f.args);
  f2 = f.func(x2, f.args);
  
  /* Ueberpruefung ob das Intervall sinnvoll gewaehlt wurde */
  if (f1 * f2 > 0) {
    return 1;
  }
  
  /* Bisektion solange bis die Nullstelle mit einem absoluten Fehler <= epsilon
   * bestimmt wurde. */
  while (fabs(f1 - f2) > epsilon) {
    mid = 0.5 * (x1 + x2);
    fmid = f.func(mid, f.args);
    
    /* Bisektion: das Intervall wurde halbiert und der Funktionswert der Mitte
     * berechnet. Einige Faelle (Rest analog):
     * Fall (+,+,-): Die Nullstelle liegt im rechten Teilintervall 
     * Fall (+,-,-): Die Nullstelle liegt im linken Teilintervall */
    if (f1 * fmid > 0) {
      x1 = mid;
      f1 = fmid;
    } else {
      x2 = mid;
      f2 = fmid;
    }
  }
  
  *root = 0.5 * (x1 + x2);
  return 0;
}

double recursive_simpson(function f, double prev_simp,
                         double x0, double x1, double x2,
                         double f0, double f1, double f2,
                         double epsilon, int rdepth) {
  /* Wir haben zwei Intervalle [x0,x1] und [x1,x2] mit den Funktionswerten f0, 
   * f1 und f2. Es wird die naechste Verfeinerung nach Simpson berechnet. */
  double h = x1 - x0;
  /* Berechne den Beitrag des Intervalls [x0,x1] zum Integral, dazu muss die
   * Stuetzstelle in der Haelfte des Intervalls berechnet werden: */
  double x01 = 0.5 * (x0 + x1);
  double f01 = f.func(x01, f.args);
  
  double simp01 = (f0 + 4*f01 + f1) * h / 6;
  
  /* Berechne den Beitrag des Intervalls [x1,x2] zum Integral: */
  double x12 = x01 + h;
  double f12 = f.func(x12, f.args);
  
  double simp12 = (f1 + 4*f12 + f2) * h / 6;
  
  /* Neue Abschaetzung des Integrals: */
  double simp = simp01 + simp12;
  
  /* Rekursionszaehler */
  rdepth--;
  
  /* Vergleich des neuen Schaetzwertes mit dem Alten. Sollte die Differenz
   * groesser als epsilon sein, werden die Intervalle halbiert und der Prozess
   * erneut auf die halbierten Intervalle durchgefuehrt. */
  if (rdepth > 0 && fabs(prev_simp - simp) > epsilon) {
    return recursive_simpson(f, simp01, x0, x01, x1, f0, f01, f1, 0.5*epsilon, rdepth)
           + recursive_simpson(f, simp12, x1, x12, x2, f1, f12, f2, 0.5*epsilon, rdepth);
  } else {
    return simp;
  }
}

  dnnPSO_optimized_function<T>::dnnPSO_optimized_function(neuralnetwork<T>& nn,
							  const dataset<T>* input,
							  const function<std::vector< math::vertex<T> >,T>& dist)
  {
    this->testnet = new neuralnetwork<T>(nn);
    this->input   = input;
    this->pseudodist = dist.clone();
    
    {
      math::vertex<T> v;
      
      if(this->testnet->exportdata(v) == false)
	throw std::logic_error("bad neural network parameter/copy");
      
      this->fvector_dimension = v.size();
    }
    
  }

 inline void swap(function<Signature>& f1, function<Signature>& f2)
 {
   f1.swap(f2);
 }

inline void swap(function<Signature, Allocator>& f1,
                 function<Signature, Allocator>& f2)
{
    f1.swap(f2);
}

bool remember_table::lookup_entry(function const & f, ex & result) const
{
	long entry = f.gethash() & (table_size-1);
	GINAC_ASSERT(entry<size());
	return operator[](entry).lookup_entry(f,result);
}

  void TaskManager :: CreateJob (const function<void(TaskInfo&)> & afunc,
                                 int antasks)
  {
    if (num_threads == 1 || !task_manager || func)
      {
        if (startup_function) (*startup_function)();
        
        TaskInfo ti;
        ti.ntasks = antasks;
        ti.thread_nr = 0; ti.nthreads = 1;
        // ti.node_nr = 0; ti.nnodes = 1;
        for (ti.task_nr = 0; ti.task_nr < antasks; ti.task_nr++)
          afunc(ti);

        if (cleanup_function) (*cleanup_function)();        
        return;
      }

    
    trace->StartJob(jobnr, afunc.target_type());

    func = &afunc;

    ntasks.store (antasks); // , memory_order_relaxed);
    ex = nullptr;


    nodedata[0]->start_cnt.store (0, memory_order_relaxed);

    jobnr++;
    
    for (int j = 0; j < num_nodes; j++)
      nodedata[j]->participate |= 1;

    if (startup_function) (*startup_function)();
    
    int thd = 0;
    int thds = GetNumThreads();
    int mynode = num_nodes * thd/thds;

    IntRange mytasks = Range(int(ntasks)).Split (mynode, num_nodes);
    NodeData & mynode_data = *(nodedata[mynode]);

    TaskInfo ti;
    ti.nthreads = thds;
    ti.thread_nr = thd;
    // ti.nnodes = num_nodes;
    // ti.node_nr = mynode;

    try
      {
        while (1)
          {
            int mytask = mynode_data.start_cnt++;
            if (mytask >= mytasks.Size()) break;
            
            ti.task_nr = mytasks.First()+mytask;
            ti.ntasks = ntasks;

            {
              RegionTracer t(ti.thread_nr, jobnr, RegionTracer::ID_JOB, ti.task_nr);
              (*func)(ti); 
            }
          }

      }
    catch (Exception e)
      {
        {
          lock_guard<mutex> guard(copyex_mutex);
          delete ex;
          ex = new Exception (e);
          mynode_data.start_cnt = mytasks.Size();
        }
      }

    if (cleanup_function) (*cleanup_function)();
    
    for (int j = 0; j < num_nodes; j++)
      if (workers_on_node[j])
        {
          while (complete[j] != jobnr)
            _mm_pause();
        }

    func = nullptr;
    if (ex)
      throw Exception (*ex);

    trace->StopJob();
  }

remember_table_entry::remember_table_entry(function const & f, ex  r)
  : hashvalue(f.gethash()), seq(f.seq), result(std::move(r))
{
	++last_access = access_counter;
	successful_hits = 0;
}

void remember_table::add_entry(function const & f, ex const & result)
{
	long entry = f.gethash() & (table_size-1);
	GINAC_ASSERT(entry<size());
	operator[](entry).add_entry(f,result);
}        

 bool is_empty_function(function<Sig, IArchive, OArchive> const& f)
 {
     return f.empty();
 }

	bool check() {
		if(active) {
			if(!time_left()) {
				if(!callback.empty()) {
					callback();
					callback.clear();
				}
				active = false;
				return true;
			}
		}
		return false;
	}

void GaussianPDF::set_parameter_values(ConstVecRef values)
{
  verify_set_parameter_values(values, "GaussianPDF");
  assign_parameters(values, x0_, s_);
}

	void cancel() {
		active = false;
		callback.clear();
	}