开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 通讯/手机编程 > 3G开发 > 查看源码
laguerre.cpp
资源名称：simulation_wireless_c++.rar　[点击查看]
上传用户：jtjnyq9001
上传日期：2014-11-21
资源大小：3974k
文件大小：2k
源码类别：

3G开发

开发平台：

Visual C++

laguerre.cpp：源码内容
  1. //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2. //
  3. //  File = laguerre.cpp
  4. //
  5. //  Laguerre method for finding polynomial roots
  6. //
  8. #include <math.h>
  9. #include <stdlib.h>
  10. #include <iostream>
  11. #include "laguerre.h"
  12. using namespace std;
  14. #ifdef _DEBUG
  15.   extern ofstream *DebugFile;
  16. #endif
  18. int LaguerreMethod( CmplxPolynomial *poly,
  19.                     std::complex<double> *root_ptr,
  20.                     double epsilon,
  21.                     double epsilon2,
  22.                     int max_iter)
  23. {
  24. int iter, j, order;
  25. complex<double> p_eval, p_p, p_p_p;
  26. complex<double> root, f, f_sq, g, rad;
  27. complex<double> f_p_r, f_m_r, delta_root;
  28. double old_delta_mag, root_mag, error;
  29. complex<double> *coeff;
  30. order = poly->GetDegree();
  31. coeff = new complex<double>[order+1];
  32. poly->CopyCoeffs(coeff);
  33. root = *root_ptr;
  34. old_delta_mag = 1000.0;
  36. for(iter=1; iter<=max_iter; iter++)
  37.   {
  38.   p_p_p = complex<double>(0.0,0.0);
  39.   p_p = complex<double>(0.0,0.0);
  40.   p_eval = coeff[order];
  41.   error = std::abs<double>(p_eval);
  42.   root_mag = std::abs<double>(root);
  44.   for( j=order-1; j>=0; j--)
  45.     {
  46.     p_p_p = p_p + root * p_p_p;
  47.     p_p = p_eval + root * p_p;
  48.     p_eval = coeff[j] + root * p_eval;
  49.     error = std::abs<double>(p_eval) + root_mag * error;
  50.     }
  51.   error = epsilon2 * error;
  52.   p_p_p = 2.0 * p_p_p;
  54.   if(std::abs<double>(p_eval) < error)
  55.     {
  56.     cout << "mag(p_eval) = " << std::abs<double>(p_eval) << "  error = "
  57.          << error << endl;
  58.     *root_ptr = root;
  59.     delete [] coeff;
  60.     return(1);
  61.     }
  62.   f = p_p/p_eval;
  63.   f_sq = f * f;
  64.   g = f_sq - p_p_p/p_eval;
  65.   rad = double(order-1)*(double(order) * g - f_sq);
  66.   rad = std::sqrt<double>(rad);
  67.   f_p_r = f + rad;
  68.   f_m_r = f - rad;
  69.   if( std::abs<double>(f_p_r) > std::abs<double>(f_m_r) )
  70.     {
  71.     delta_root = std::complex<double>(double(order), 0.0) / f_p_r;
  72.     }
  73.   else
  74.     {
  75.     delta_root = std::complex<double>(double(order), 0.0) / f_m_r;
  76.     }
  77.   root = root - delta_root;
  78.   if( (iter > 6) && (std::abs<double>(delta_root) > old_delta_mag) )
  79.     {
  80.     *root_ptr = root;
  81.     delete [] coeff;
  82.     return(2);
  83.     }
  84.   if( std::abs<double>(delta_root) < (epsilon * std::abs<double>(root)))
  85.     {
  86.     *root_ptr = root;
  87.     delete [] coeff;
  88.     return(3);
  89.     }
  90.   old_delta_mag = std::abs<double>(delta_root);
  91.   }
  92. #ifdef _DEBUG
  93.   *DebugFile << "Laguerre method failed to converge" << endl;
  94. #endif
  95. delete [] coeff;
  96. return(-1);
  97. }