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

源码开发语言/平台
当前位置：首页 > 源码/资料 > 多媒体 > 语音合成与识别 > 查看源码
FE_polynomial.h
资源名称：speech_analysis.rar　[点击查看]
上传用户：italyroyal
上传日期：2013-05-06
资源大小：473k
文件大小：21k
源码类别：

语音合成与识别

开发平台：

Visual C++

FE_polynomial.h：源码内容
  1. ///////////////////////////////////////////////////////////////////////////////
  2. // This is a part of the Feature program.
  3. // Version: 1.0
  4. // Date: February 22, 2003
  5. // Programmer: Oh-Wook Kwon
  6. // Copyright(c) 2003 Oh-Wook Kwon. All rights reserved. owkwon@ucsd.edu
  7. ///////////////////////////////////////////////////////////////////////////////
  9. #ifndef _FE_POLYNOMIAL_H_
  10. #define _FE_POLYNOMIAL_H_
  12. #include <cstdio>
  13. #include <cstdarg>
  14. #include <cmath>
  15. #include <cstring>
  16. #include <cstdlib>
  17. #include <typeinfo>
  18. #include <vector>
  19. using namespace std;
  20. #pragma warning(disable:4786)
  22. #include "FE_complex.h"
  24. template <typename T> class CPolynomial;
  26. // friend template functions
  27. template <typename T> CPolynomial<T> poly_plus(const CPolynomial<T>& A, const CPolynomial<T>& B);
  28. template <typename T> CPolynomial<T> poly_plus(const CPolynomial<T>& A, const double b);
  29. template <typename T> CPolynomial<T> poly_plus(const double b, const CPolynomial<T>& A);
  31. template <typename T> CPolynomial<T> poly_minus(const CPolynomial<T>& A, const CPolynomial<T>& B);
  32. template <typename T> CPolynomial<T> poly_minus(const CPolynomial<T>& A, const double b);
  33. template <typename T> CPolynomial<T> poly_minus(const double b, const CPolynomial<T>& A);
  35. template <typename T> CPolynomial<T> poly_times(const CPolynomial<T>& A, const CPolynomial<T>& B);
  36. template <typename T> CPolynomial<T> poly_times(const CPolynomial<T>& A, const double b);
  37. template <typename T> CPolynomial<T> poly_times(const double b, const CPolynomial<T>& A);
  39. template <typename T> CPolynomial<T> poly_divide(const CPolynomial<T>& A, const CPolynomial<T>& B);
  40. template <typename T> CPolynomial<T> poly_divide(const CPolynomial<T>& A, const double b);
  41. template <typename T> CPolynomial<T> poly_divide(const double b, const CPolynomial<T>& A);
  43. template <typename T> CPolynomial<T> poly_modulo(const CPolynomial<T>& A, const double b);
  44. template <typename T> CPolynomial<T> poly_modulo(const CPolynomial<T>& A, const double b);
  45. template <typename T> CPolynomial<T> poly_modulo(const double b, const CPolynomial<T>& A);
  47. // GNU compiler gcc-2.95.2 cannot handle operator overloading with user-defined class type.
  48. #if 0
  49. template <typename T> CPolynomial<T> poly_plus(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  50. template <typename T> CPolynomial<T> poly_minus(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  51. template <typename T> CPolynomial<T> poly_times(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  52. template <typename T> CPolynomial<T> poly_divide(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  53. #endif
  56. template<typename T> 
  57. class CPolynomial
  58. {
  59. typedef Complex<T> CComplex;
  60. public:
  61. CPolynomial();
  62. CPolynomial(const int v0) { v.resize(1); v[0]=(T)v0; }
  63. CPolynomial(const double v0) { v.resize(1); v[0]=(T)v0; }
  64. CPolynomial(int m, const T* a);
  65. CPolynomial(const char* a);
  66. CPolynomial(const CPolynomial<T>& a) : v(a.v) {}
  67. virtual ~CPolynomial();
  68. void Print() const;
  69. T Eval(T x, T* pd = NULL, int nd = 0) const;
  70. void Roots(CComplex* roots) const;
  71. void zroots(CComplex* a, int m, CComplex* roots, int polish) const;
  72. void laguer(CComplex* a, int m, CComplex* x, int *its) const;
  74. public:
  75. vector<T> v;
  77. CPolynomial<T> operator+();
  78. CPolynomial<T> operator-();
  79. CPolynomial<T>& operator=(const CPolynomial<T>& a);
  80. CPolynomial<T>& operator+=(const CPolynomial<T>& a);
  81. CPolynomial<T>& operator-=(const CPolynomial<T>& a);
  82. CPolynomial<T>& operator*=(const CPolynomial<T>& a);
  83. CPolynomial<T>& operator/=(const CPolynomial<T>& a);
  85. // GCC-2.95-2 did not allow friend operator overloading definition outside class declaration.
  86.         friend CPolynomial<T> operator+(const CPolynomial<T>& A, const CPolynomial<T>& B) {
  87.                 return poly_plus(A,B);
  88.         }
  89.         friend CPolynomial<T> operator+(const CPolynomial<T>& A, const double b) {
  90.                 return poly_plus(A,b);
  91.         }
  92.         friend CPolynomial<T> operator+(const double b, const CPolynomial<T>& A) {
  93.                 return poly_plus(b,A);
  94.         }
  95.         friend CPolynomial<T> operator-(const CPolynomial<T>& A, const CPolynomial<T>& B) {
  96.                 return poly_minus(A,B);
  97.         }
  98.         friend CPolynomial<T> operator-(const CPolynomial<T>& A, const double b) {
  99.                 return poly_minus(A,b);
  100.         }
  101.         friend CPolynomial<T> operator-(const double b, const CPolynomial<T>& A) {
  102.                 return poly_minus(b,A);
  103.         }
  104.         friend CPolynomial<T> operator*(const CPolynomial<T>& A, const CPolynomial<T>& B) {
  105.                 return poly_times(A,B);
  106.         }
  107.         friend CPolynomial<T> operator*(const CPolynomial<T>& A, const double b) {
  108.                 return poly_times(A,b);
  109.         }
  110.         friend CPolynomial<T> operator*(const double b, const CPolynomial<T>& A) {
  111.                 return poly_times(b,A);
  112.         }
  113.         friend CPolynomial<T> operator/(const CPolynomial<T>& A, const CPolynomial<T>& B) {
  114.                 return poly_divide(A,B);
  115.         }
  116.         friend CPolynomial<T> operator/(const CPolynomial<T>& A, const double b) {
  117.                 return poly_divide(A,b);
  118.         }
  119.         friend CPolynomial<T> operator/(const double b, const CPolynomial<T>& A) {
  120.                 return poly_divide(b,A);
  121.         }
  122.         friend CPolynomial<T> operator%(const CPolynomial<T>& A, const CPolynomial<T>& B) {
  123.                 return poly_modulo(A,B);
  124.         }
  125.         friend CPolynomial<T> operator%(const CPolynomial<T>& A, const double b) {
  126.                 return poly_modulo(A,b);
  127.         }
  128.         friend CPolynomial<T> operator%(const double b, const CPolynomial<T>& A) {
  129.                 return poly_modulo(b,A);
  130.         }
  132. #if 0
  133.         friend CPolynomial<T> operator+(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B) {
  134.                 return poly_plus(A,B);
  135.         }
  136.         friend CPolynomial<T> operator+(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B) {
  137.                 return poly_minus(A,B);
  138.         }
  139.         friend CPolynomial<T> operator+(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B) {
  140.                 return poly_times(A,B);
  141.         }
  142.         friend CPolynomial<T> operator+(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B) {
  143.                 return poly_divide(A,B);
  144.         }
  145. #endif
  147. #ifdef __GNUC__
  148.         friend CPolynomial<T> poly_plus<T>(const CPolynomial<T>& A, const CPolynomial<T>& B);
  149.         friend CPolynomial<T> poly_plus<T>(const CPolynomial<T>& A, const double b);
  150.         friend CPolynomial<T> poly_plus<T>(const double b, const CPolynomial<T>& A);
  151.         friend CPolynomial<T> poly_minus<T>(const CPolynomial<T>& A, const CPolynomial<T>& B);
  152.         friend CPolynomial<T> poly_minus<T>(const CPolynomial<T>& A, const double b);
  153.         friend CPolynomial<T> poly_minus<T>(const double b, const CPolynomial<T>& A);
  154.         friend CPolynomial<T> poly_times<T>(const CPolynomial<T>& A, const CPolynomial<T>& B);
  155.         friend CPolynomial<T> poly_times<T>(const CPolynomial<T>& A, const double b);
  156.         friend CPolynomial<T> poly_times<T>(const double b, const CPolynomial<T>& A);
  157.         friend CPolynomial<T> poly_divide<T>(const CPolynomial<T>& A, const CPolynomial<T>& B);
  158.         friend CPolynomial<T> poly_divide<T>(const CPolynomial<T>& A, const double b);
  159.         friend CPolynomial<T> poly_divide<T>(const double b, const CPolynomial<T>& A);
  160.         friend CPolynomial<T> poly_modulo<T>(const CPolynomial<T>& A, const CPolynomial<T>& B);
  161.         friend CPolynomial<T> poly_modulo<T>(const CPolynomial<T>& A, const double b);
  162.         friend CPolynomial<T> poly_modulo<T>(const double b, const CPolynomial<T>& A);
  163. #if 0
  164. friend CPolynomial<T> poly_plus<T>(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  165. friend CPolynomial<T> poly_minus<T>(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  166. friend CPolynomial<T> poly_times<T>(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  167. friend CPolynomial<T> poly_divide<T>(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  168. #endif
  169. #else
  170. friend CPolynomial<T> poly_plus(const CPolynomial<T>& A, const CPolynomial<T>& B);
  171.         friend CPolynomial<T> poly_plus(const CPolynomial<T>& A, const double b);
  172.         friend CPolynomial<T> poly_plus(const double b, const CPolynomial<T>& A);
  173.         friend CPolynomial<T> poly_minus(const CPolynomial<T>& A, const CPolynomial<T>& B);
  174.         friend CPolynomial<T> poly_minus(const CPolynomial<T>& A, const double b);
  175.         friend CPolynomial<T> poly_minus(const double b, const CPolynomial<T>& A);
  176.         friend CPolynomial<T> poly_times(const CPolynomial<T>& A, const CPolynomial<T>& B);
  177.         friend CPolynomial<T> poly_times(const CPolynomial<T>& A, const double b);
  178.         friend CPolynomial<T> poly_times(const double b, const CPolynomial<T>& A);
  179.         friend CPolynomial<T> poly_divide(const CPolynomial<T>& A, const CPolynomial<T>& B);
  180.         friend CPolynomial<T> poly_divide(const CPolynomial<T>& A, const double b);
  181.         friend CPolynomial<T> poly_divide(const double b, const CPolynomial<T>& A);
  182.         friend CPolynomial<T> poly_modulo(const CPolynomial<T>& A, const CPolynomial<T>& B);
  183.         friend CPolynomial<T> poly_modulo(const CPolynomial<T>& A, const double b);
  184.         friend CPolynomial<T> poly_modulo(const double b, const CPolynomial<T>& A);
  185. #if 0
  186. friend CPolynomial<T> poly_plus(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  187. friend CPolynomial<T> poly_minus(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  188. friend CPolynomial<T> poly_times(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  189. friend CPolynomial<T> poly_divide(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
  190. #endif
  191. #endif
  193. #if 0
  194. friend CPolynomial<T> operator+(const CPolynomial<T>& a, const CPolynomial<T>& b);
  195. friend CPolynomial<T> operator+(const CPolynomial<T>& a, const double& b);
  196. friend CPolynomial<T> operator+(const double& a, const CPolynomial<T>& b);
  197. friend CPolynomial<T> operator+(const CPolynomial<T>& a);
  198. friend CPolynomial<T> operator-(const CPolynomial<T>& a, const CPolynomial<T>& b);
  199. friend CPolynomial<T> operator-(const CPolynomial<T>& a, const double& b);
  200. friend CPolynomial<T> operator-(const double& a, const CPolynomial<T>& b);
  201. friend CPolynomial<T> operator-(const CPolynomial<T>& a);
  202. friend CPolynomial<T> operator*(const CPolynomial<T>& a, const CPolynomial<T>& b);
  203. friend CPolynomial<T> operator*(const CPolynomial<T>& a, const double& b);
  204. friend CPolynomial<T> operator*(const double& a, const CPolynomial<T>& b);
  205. friend CPolynomial<T> operator/(const CPolynomial<T>& a, const CPolynomial<T>& b);
  206. friend CPolynomial<T> operator/(const CPolynomial<T>& a, const double& b);
  207. friend CPolynomial<T> operator/(const double& a, const CPolynomial<T>& b);
  208. friend CPolynomial<T> operator%(const CPolynomial<T>& a, const CPolynomial<T>& b);
  209. friend CPolynomial<T> operator%(const CPolynomial<T>& a, const double& b);
  210. friend CPolynomial<T> operator%(const double& a, const CPolynomial<T>& b);
  211. #endif
  213. };
  215. template<typename T>
  216. inline void PolynomialDivision(const CPolynomial<T>& a, const CPolynomial<T>& b, CPolynomial<T>& q, CPolynomial<T>& r)
  217. {
  218. int an = a.v.size() - 1;
  219. int bn = a.v.size() - 1;
  220. int qn = ((an >= bn) ? an-bn : 0);
  221. int rn = ((bn >= 1) ? ((an <bn-1) ? an : bn-1) : 0);
  222. q.v.resize(qn+1);
  223. r.v.resize(rn+1);
  225. int k,j;
  226. for(j=0;j<=an;j++){
  227. r.v[j] = a.v[j];
  228. q.v[j] = 0;
  229. }
  230. for(k=an-bn;k>=0;k--){
  231. q.v[k] = r.v[bn+k]/b.v[bn];
  232. for(j=bn+k-1;j>=k;j--) r.v[j] -= q.v[k]*b.v[j-k];
  233. }
  234. for(j=bn;j<=an;j++) r.v[j] = 0;
  235. }
  237. template<typename T>
  238. inline CPolynomial<T> poly_plus(const CPolynomial<T>& a, const CPolynomial<T>& b)
  239. {
  240. CPolynomial<T> c = a;
  241. c += b;
  242. return c;
  243. }
  245. template<typename T>
  246. inline CPolynomial<T> poly_plus(const CPolynomial<T>& a, const double b)
  247. {
  248. CPolynomial<T> c = a;
  249. c += CPolynomial<T>((T)b);
  250. return c;
  251. }
  253. template<typename T>
  254. inline CPolynomial<T> poly_plus(const double a, const CPolynomial<T>& b)
  255. {
  256. CPolynomial<T> c((T)a);
  257. c += b;
  258. return c;
  259. }
  261. template<typename T>
  262. inline CPolynomial<T> operator+(const CPolynomial<T>& a)
  263. {
  264. CPolynomial<T> c = a;
  265. return c;
  266. }
  268. template<typename T>
  269. inline CPolynomial<T> poly_minus(const CPolynomial<T>& a, const CPolynomial<T>& b)
  270. {
  271. CPolynomial<T> c = a;
  272. c -= b;
  273. return c;
  274. }
  276. template<typename T>
  277. inline CPolynomial<T> poly_minus(const CPolynomial<T>& a, const double b)
  278. {
  279. CPolynomial<T> c = a;
  280. c -= CPolynomial<T>((T)b);
  281. return c;
  282. }
  284. template<typename T>
  285. inline CPolynomial<T> poly_minus(const double a, const CPolynomial<T>& b)
  286. {
  287. CPolynomial<T> c((T)a);
  288. c -= b;
  289. return c;
  290. }
  292. template<typename T>
  293. inline CPolynomial<T> operator-(const CPolynomial<T>& a)
  294. {
  295. CPolynomial<T> c = 0;
  296. c -= a;
  297. return c;
  298. }
  300. template<typename T>
  301. inline CPolynomial<T> poly_times(const CPolynomial<T>& a, const CPolynomial<T>& b)
  302. {
  303. CPolynomial<T> c = a;
  304. c *= b;
  305. return c;
  306. }
  308. template<typename T>
  309. inline CPolynomial<T> poly_times(const CPolynomial<T>& a, const double b)
  310. {
  311. CPolynomial<T> c = a;
  312. c *= CPolynomial<T>((T)b);
  313. return c;
  314. }
  316. template<typename T>
  317. inline CPolynomial<T> poly_times(const double a, const CPolynomial<T>& b)
  318. {
  319. CPolynomial<T> c((T)a);
  320. c *= b;
  321. return c;
  322. }
  324. template<typename T>
  325. inline CPolynomial<T> poly_divide(const CPolynomial<T>& a, const CPolynomial<T>& b)
  326. {
  327. CPolynomial<T> q,r;
  328. PolynomialDivision(a,b,q,r);
  329. return q;
  330. }
  332. template<typename T>
  333. inline CPolynomial<T> poly_divide(const CPolynomial<T>& a, const double b)
  334. {
  335. CPolynomial<T> q,r;
  336. PolynomialDivision(a,CPolynomial<T>((T)b),q,r);
  337. return q;
  338. }
  340. template<typename T>
  341. inline CPolynomial<T> poly_divide(const double a, const CPolynomial<T>& b)
  342. {
  343. CPolynomial<T> q,r;
  344. PolynomialDivision(CPolynomial<T>((T)a),b,q,r);
  345. return q;
  346. }
  348. template<typename T>
  349. inline CPolynomial<T> poly_modulo(const CPolynomial<T>& a, const CPolynomial<T>& b)
  350. {
  351. CPolynomial<T> q,r;
  352. PolynomialDivision(a,b,q,r);
  353. return r;
  354. }
  356. template<typename T>
  357. inline CPolynomial<T> poly_modulo(const CPolynomial<T>& a, const double b)
  358. {
  359. CPolynomial<T> q,r;
  360. PolynomialDivision(a,CPolynomial<T>((T)b),q,r);
  361. return r;
  362. }
  364. template<typename T>
  365. inline CPolynomial<T> poly_modulo(const double a, const CPolynomial<T>& b)
  366. {
  367. CPolynomial<T> q,r;
  368. PolynomialDivision(CPolynomial<T>((T)a),b,q,r);
  369. return r;
  370. }
  373. #if 0
  374. template<typename T> inline CPolynomial<T> poly_plus(const CPolynomial<Complex<T> >& a, const CPolynomial<Complex<T> >& b)
  375. {
  376. CPolynomial<T> c = a; c += b; return c;
  377. }
  379. template<typename T> inline CPolynomial<T> poly_minus(const CPolynomial<Complex<T> >& a, const CPolynomial<Complex<T> >& b)
  380. {
  381. CPolynomial<T> c = a; c -= b; return c;
  382. }
  384. template<typename T> inline CPolynomial<T> poly_times(const CPolynomial<Complex<T> >& a, const CPolynomial<Complex<T> >& b)
  385. {
  386. CPolynomial<T> c = a; c *= b; return c;
  387. }
  389. template<typename T> inline CPolynomial<T> poly_divide(const CPolynomial<Complex<T> >& a, const CPolynomial<Complex<T> >& b)
  390. {
  391. CPolynomial<T> c = a; c /= b; return c;
  392. }
  393. #endif
  395. template<typename T>
  396. CPolynomial<T>::CPolynomial(int m, const T* a)
  397. {
  398. for(int i=0; i<=m; i++){
  399. v.push_back(a[i]);
  400. }
  401. }
  403. #if 0
  404. // do not use <cstdarg> and <typeinfo> because they are very fragile
  405. // and depend on compilers
  406. template<typename T>
  407. CPolynomial<T>::CPolynomial(int m, const T v0, ...)
  408. {
  409. va_list args;
  410. va_start (args, m); // Although GNU compiler warns here, we can ignore it.
  411. const type_info& a = typeid(v0);
  412. if(strcmp(a.name(),"float")==0 || strcmp(a.name(),"double")==0){
  413. for(int i=1; i<=m; i++){
  414. v.push_back((T)va_arg(args, double));
  415. }
  416. }
  417. else{
  418. printf("Unsupported Polynomial type (%s)!n", a.name());
  419. exit(1);
  420. }
  421.     va_end (args);
  422. }
  423. #endif
  426. template<typename T>
  427. CPolynomial<T>::CPolynomial(const char* a)
  428. {
  429. char *t = strdup(a);
  430. char *p = strtok(t, " trn");
  431. v.push_back((T)atof(p));
  432. while((p = strtok(NULL, " trn"))){
  433. v.push_back((T)atof(p));
  434. }
  435. free(t);
  436. }
  438. template<typename T>
  439. CPolynomial<T>::CPolynomial()
  440. {
  441. v.push_back((T)0);
  442. }
  444. template<typename T>
  445. CPolynomial<T>::~CPolynomial()
  446. {
  447. }
  449. #define EPS 2.0e-6
  450. #define MAXM 100
  452. /* Roots() yields n roots into roots[1]...roots[n] */
  453. template<typename T>
  454. void CPolynomial<T>::Roots(CComplex* roots) const
  455. {
  456. int n = v.size()-1;
  457. vector<CComplex> a(n+1);
  458. for(int i=0;i<=n;i++) a[i] = CComplex(v[i],0);
  459. zroots(&a[0], n, roots, 1);
  460. }
  462. #define EPSS 1.0e-7
  463. #define MR 8
  464. #define MT 10
  465. #define MAXIT (MT*MR)
  467. #ifndef my_max
  468. #define my_max(a, b) (((a) > (b)) ? (a) : (b)) 
  469. #endif
  470. #ifndef my_min
  471. #define my_min(a, b) (((a) < (b)) ? (a) : (b)) 
  472. #endif
  475. template<typename T>
  476. void CPolynomial<T>::zroots(CComplex* a, int m, CComplex* roots, int polish) const
  477. {
  478. int i,its,j,jj;
  479. CComplex x,b,c,ad[MAXM];
  481. for(j=0;j<=m;j++) ad[j] = a[j];
  482. for(j=m;j>=1;j--){
  483. //x = 0;
  484. x = roots[j];
  485. laguer(ad,j,&x,&its);
  486. if(fabs(imag(x)) <= 2.0*EPS*fabs(real(x))) x=CComplex(real(x), 0);
  487. roots[j] = x;
  488. b = ad[j];
  489. for(jj=j-1;jj>=0;jj--){
  490. c = ad[jj];
  491. ad[jj] = b;
  492. b = x*b+c;
  493. }
  494. }
  495. if(polish)
  496. for(j=1;j<=m;j++)
  497. laguer(a,m,&roots[j],&its);
  498. for(j=2;j<=m;j++){
  499. x = roots[j];
  500. for(i=j-1;i>=1;i--){
  501. if(real(roots[i]) <= real(x)) break;
  502. roots[i+1] = roots[i];
  503. }
  504. roots[i+1] = x;
  505. }
  506. }
  508. template<typename T>
  509. void CPolynomial<T>::laguer(CComplex* a, int m, CComplex* x, int *its) const
  510. {
  511. int iter,j;
  512. T abx,abp,abm,err;
  513. CComplex dx,x1,b,d,f,g,h,sq,gp,gm,g2;
  514. static T frac[MR+1]={
  515. (T)0.00,(T)0.50,(T)0.25,
  516. (T)0.75,(T)0.13,(T)0.38,
  517. (T)0.62,(T)0.88,(T)1.00
  518. };
  520. for(iter=1;iter<=MAXIT;iter++){
  521. *its = iter;
  522. b = a[m];
  523. err = (T)abs(b);
  524. d = f = 0;
  525. abx = (T)abs(*x);
  526. for(j=m-1;j>=0;j--){
  527. f = (*x)*f+d;
  528. d = (*x)*d+b;
  529. b = (*x)*b+a[j];
  530. err = (T)abs(b)+abx*err;
  531. }
  532. err *= (T)EPSS;
  533. if(abs(b) <= err) return;
  534. g = d/b;
  535. g2 = g*g;
  536. h = g2 - (T)2*(f/b);
  537. sq = sqrt((T)(m-1)*(((T)m*h)-g2));
  538. gp = g+sq;
  539. gm = g-sq;
  540. abp = (T)abs(gp);
  541. abm = (T)abs(gm);
  542. if(abp < abm) gp = gm;
  543. dx = ((my_max(abp,abm) > 0.0 ? CComplex(m,0)/gp
  544. : ((T)exp(log(1+abx))*CComplex(cos(iter),sin(iter)))));
  545. x1 = (*x)-dx;
  546. if(real(*x) == real(x1) && imag(*x) == imag(x1)) return;
  547. if(iter%MT) *x = x1;
  548. else *x = (*x)-frac[iter/MT]*dx;
  549. }
  550. printf("too many iteration in laguer");
  551. return;
  552. }
  554. template<typename T>
  555. T CPolynomial<T>::Eval(T x, T* pd, int nd) const
  556. {
  557. int nnd, j, i;
  558. T cnst = 1;
  559. T ftmp;
  561. if(!pd){
  562. nd = 0;
  563. }
  564. int n = v.size()-1;
  565. ftmp = v[n];
  566. if(pd) pd[0] = v[n];
  567. for(j=1;j<=nd;j++) pd[j] = 0;
  568. for(i=n-1;i>=0;i--){
  569. nnd = (nd < (n-i) ? nd : n-i);
  570. for(j=nnd;j>=1;j--)
  571. pd[j] = (T)(pd[j]*x+pd[j-1]);
  572. ftmp = (T)(ftmp*x+v[i]);
  573. if(pd) pd[0] = (T)(pd[0]*x+v[i]);
  574. }
  575. for(i=2;i<=nd;i++){
  576. cnst *= i;
  577. pd[i] *= cnst;
  578. }
  579. return ftmp;
  580. }
  582. template<typename T>
  583. void CPolynomial<T>::Print() const
  584. {
  585. int n=v.size()-1;
  586. int bAllZero = 1;
  587. for(int i=n;i>=0;i--){
  588. if(fabs(v[i]) < 1e-37){
  589. if(i==0 && bAllZero) printf("0");
  590. continue;
  591. }
  592. bAllZero = 0;
  593. if(v[i] > 0 && i != n) printf(" + ");
  594. else if(v[i] < 0) printf(" - ");
  596. if(fabs(fabs(v[i])-1) > 1e-37){
  597. printf("%gx(%d)",fabs(v[i]),i);
  598. }
  599. else{
  600. printf("x(%d)",i);
  601. }
  602. }
  603. printf("n");
  604. }
  606. template<typename T>
  607. CPolynomial<T>& CPolynomial<T>::operator=(const CPolynomial<T>& a)
  608. {
  609. int an = a.v.size()-1;
  610. int n = an;
  611. v.resize(n+1);
  612. int m = 0;
  613. for(int i=0; i<=n; i++){
  614. v[i] = a.v[i];
  615. if(v[i] != 0.0) 
  616. m = i;
  617. }
  618. if(m != n) v.resize(m+1);
  619. return *this;
  620. }
  622. template<typename T>
  623. CPolynomial<T>& CPolynomial<T>::operator+=(const CPolynomial<T>& a)
  624. {
  625. int m = v.size()-1;
  626. int an = a.v.size()-1;
  627. int n = my_max(m,an);
  628. v.resize(n+1);
  629. for(int i=0; i<=n; i++)
  630. v[i] = ((i<=m)? v[i]: 0) + ((i<=an)? a.v[i]: 0);
  631. return *this;
  632. }
  634. template<typename T>
  635. CPolynomial<T>& CPolynomial<T>::operator-=(const CPolynomial<T>& a)
  636. {
  637. int m = v.size()-1;
  638. int an = a.v.size()-1;
  639. int n = my_max(m,an);
  640. v.resize(n+1);
  641. for(int i=0; i<=n; i++)
  642. v[i] = ((i<=m)? v[i]: 0) - ((i<=an)? a.v[i]: 0);
  643. return *this;
  644. }
  646. template<typename T>
  647. CPolynomial<T>& CPolynomial<T>::operator*=(const CPolynomial<T>& a)
  648. {
  649. int p,i;
  650. int m = v.size()-1;
  651. int an = a.v.size()-1;
  652. int n = m + an;
  653. v.resize(n+1);
  654. for(p=n;p>m;p--) v[p] = 0;
  655. for(p=n; p>=0; p--){
  656. T sum = 0;
  657. for(i=my_min(p,m); i>=0 && p-i<=an; i--){
  658. sum += v[i] * a.v[p-i];
  659. }
  660. v[p] = sum;
  661. }
  662. return *this;
  663. }
  665. template<typename T>
  666. CPolynomial<T>& CPolynomial<T>::operator/=(const CPolynomial<T>& a)
  667. {
  668. CPolynomial<T> q,r;
  669. PolynomialDivision(*this,a,q,r);
  670. return q;
  671. }
  673. #endif