![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
Changeset 514 in Sophya for trunk/SophyaLib/NTools/poly.cc
- Timestamp:
- Oct 25, 1999, 6:43:04 PM (26 years ago)
- File:
-
- 1 edited
-
trunk/SophyaLib/NTools/poly.cc (modified) (40 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/NTools/poly.cc
r508 r514 13 13 //++ 14 14 // Links Parents 15 // OVector15 // Vector 16 16 //-- 17 17 … … 21 21 22 22 23 ////////////////////////////////////////////////////////////////////////// 23 24 //++ 24 25 // Poly::Poly(int degre=0) … … 29 30 30 31 Poly::Poly(int degre) 31 : OVector(degre+1), dirty(0), deg(0)32 : Vector(degre+1), dirty(0), deg(0) 32 33 { 33 34 END_CONSTRUCTOR 34 35 } 35 36 37 //++ 38 Poly::Poly(Poly const& a) 39 // 40 // Constructeur par copie. 41 //-- 42 :Vector(a), dirty(a.dirty), deg(a.deg) 43 { 44 END_CONSTRUCTOR 45 } 46 36 47 37 48 void Poly::UpdateDeg() const 38 49 { 39 int i = nalloc-1;40 while ( data[i]== 0 && i>0) i--;50 int i = NElts()-1; 51 while (Element(i) == 0 && i>0) i--; 41 52 42 53 (int&) deg = i; // bientot mutable dans ANSI C++ … … 61 72 { 62 73 UpdateDegIfDirty(); 63 double res = data[deg];74 double res = Element(deg); 64 75 for (int i=deg-1; i>=0; i--) { 65 76 res *= x; 66 res += data[i];77 res += Element(i); 67 78 } 68 79 return res; … … 76 87 { 77 88 UpdateDegIfDirty(); 78 if (deg == 0) { data[0]= 0; return;}89 if (deg == 0) { Element(0) = 0; return;} 79 90 for (int i=1; i<=deg; i++) 80 data[i-1] = data[i]*i;81 data[deg]= 0;91 Element(i-1) = Element(i)*i; 92 Element(deg) = 0; 82 93 deg--; 83 94 } … … 94 105 // der.Zero(); // on sait que Realloc met a zero le reste. 95 106 for (int i=1; i<=deg; i++) 96 der. data[i-1] = data[i]*i;97 } 98 99 100 //++ 101 int Poly::Roots( OVector& roots) const107 der.Element(i-1) = Element(i)*i; 108 } 109 110 111 //++ 112 int Poly::Roots(Vector& roots) const 102 113 // 103 114 // Retourne dans roots les racines réelles, si on sait … … 133 144 if (deg != 1) THROW(sizeMismatchErr); 134 145 135 if ( data[1]== 0) return 0;136 r = - data[0]/data[1];146 if (Element(1) == 0) return 0; 147 r = -Element(0)/Element(1); 137 148 return 1; 138 149 } … … 149 160 if (deg != 2) THROW(sizeMismatchErr); 150 161 151 double delta = data[1]*data[1] - 4*data[0]*data[2];162 double delta = Element(1)*Element(1) - 4*Element(0)*Element(2); 152 163 if (delta < 0) return 0; 153 164 if (delta == 0) { 154 r1 = r2 = - data[1]/2/data[0];165 r1 = r2 = -Element(1)/2/Element(0); 155 166 return 1; 156 167 } 157 r1 = (- data[1] - sqrt(delta))/2/data[0];158 r2 = (- data[1] + sqrt(delta))/2/data[0];168 r1 = (-Element(1) - sqrt(delta))/2/Element(0); 169 r2 = (-Element(1) + sqrt(delta))/2/Element(0); 159 170 return 2; 160 171 } … … 167 178 { 168 179 if (this == &a) return *this; 169 OVector::operator=(a);180 Vector::operator=(a); 170 181 171 182 UpdateDeg(); … … 185 196 b.UpdateDegIfDirty(); 186 197 187 if (b.deg > deg) 188 Realloc(b.deg); 198 if (b.deg > deg) Realloc(b.deg); 189 199 190 200 int n = (deg > b.deg) ? deg : b.deg; 191 for (int i=0; i<=n; i++) 192 data[i] += b.data[i]; 201 for (int i=0; i<=n; i++) Element(i) += b.Element(i); 193 202 194 203 UpdateDeg(); … … 204 213 b.UpdateDegIfDirty(); 205 214 206 if (b.deg > deg) 207 Realloc(b.deg); 215 if (b.deg > deg) Realloc(b.deg); 208 216 209 217 int n = (deg > b.deg) ? deg : b.deg; 210 for (int i=0; i<=n; i++) 211 data[i] -= b.data[i]; 218 for (int i=0; i<=n; i++) Element(i) -= b.Element(i); 212 219 213 220 UpdateDeg(); … … 215 222 } 216 223 217 218 //++ 219 Poly operator+ (Poly const& a, Poly const& b) 220 // 221 //-- 222 #if HAS_NAMED_RETURN 223 return c((a.deg > b.deg)?(a.deg+1):(b.deg+1)) 224 #endif 225 { 226 #if !HAS_NAMED_RETURN 227 Poly c((a.deg > b.deg)?(a.deg+1):(b.deg+1)); 228 #endif 229 c = a; 230 c += b; 231 #if !HAS_NAMED_RETURN 232 return c; 233 #endif 234 } 235 236 //++ 237 Poly operator- (Poly const& a, Poly const& b) 238 // 239 //-- 240 #if HAS_NAMED_RETURN 241 return c((a.deg > b.deg)?(a.deg+1):(b.deg+1)) 242 #endif 243 { 244 #if !HAS_NAMED_RETURN 245 Poly c((a.deg > b.deg)?(a.deg+1):(b.deg+1)); 246 #endif 247 c = a; 248 c -= b; 249 #if !HAS_NAMED_RETURN 250 return c; 251 #endif 252 } 253 254 //++ 255 Poly operator* (Poly const& a, Poly const& b) 256 // 257 //-- 258 //#if HAS_NAMED_RETURN 259 // return c(a.deg + b.deg) 260 //#endif 261 { 262 //#if !HAS_NAMED_RETURN 263 Poly c(a.deg + b.deg); 264 //#endif 265 a.UpdateDegIfDirty(); 224 //++ 225 Poly& Poly::operator *= (double a) 226 // 227 //-- 228 { 229 UpdateDegIfDirty(); 230 for (int i=0; i<=deg; i++) Element(i) *= a; 231 return *this; 232 } 233 234 //++ 235 Poly Poly::Mult(Poly const& b) const 236 // 237 //-- 238 { 239 Poly c(deg + b.deg); 240 UpdateDegIfDirty(); 266 241 b.UpdateDegIfDirty(); 267 242 268 c.deg = a.deg + b.deg;243 c.deg = deg + b.deg; 269 244 270 245 for (int i=0; i<=c.deg; i++) { 271 246 c[i] = 0; 272 int kmin = (i <= a.deg) ? 0 : i - a.deg;273 int kmax = (i <= a.deg) ? i : a.deg;247 int kmin = (i <= deg) ? 0 : i - deg; 248 int kmax = (i <= deg) ? i : deg; 274 249 for (int k=kmin; k<=kmax; k++) 275 c[i] += a[k] * b[i-k]; 276 } 277 //#if !HAS_NAMED_RETURN 250 c[i] += (*this)[k] * b[i-k]; 251 } 278 252 return c; 279 //#endif280 }281 282 //++283 Poly& Poly::operator *= (double a)284 //285 //--286 {287 UpdateDegIfDirty();288 289 for (int i=0; i<=deg; i++)290 data[i] *= a;291 292 return *this;293 }294 295 296 //++297 Poly operator* (double a, Poly const& b)298 //299 //--300 #if HAS_NAMED_RETURN301 return c(b)302 #endif303 {304 #if !HAS_NAMED_RETURN305 Poly c(b);306 #endif307 308 c *= a;309 310 #if !HAS_NAMED_RETURN311 return c;312 #endif313 }314 315 316 void Poly::ReadSelf(PInPersist& s)317 {318 int_4 dg;319 s >> dg;320 Realloc(dg, true);321 322 #ifdef DEBUG323 int r = nr; int c = nc;324 #endif325 326 OVector::ReadSelf(s);327 328 #ifdef DEBUG329 DBASSERT(r == nr && c == nc);330 #endif331 UpdateDeg();332 333 #ifdef DEBUG334 DBASSERT(dg == deg);335 #endif336 }337 338 void Poly::WriteSelf(POutPersist& s) const339 {340 UpdateDegIfDirty();341 ((Poly*)(this))->Realloc(deg, true);342 s << deg;343 OVector::WriteSelf(s);344 253 } 345 254 … … 365 274 366 275 //++ 367 ostream& operator << (ostream& s, const Poly& a) 368 // 369 // Impressions sur le stream "s". 370 //-- 371 { 372 a.Print(s); 373 return s; 374 } 375 376 //++ 377 double Poly::Fit(OVector const& x, OVector const& y, int degre) 276 double Poly::Fit(Vector const& x, Vector const& y, int degre) 378 277 // 379 278 // Ajustement polynomial par moindre carrés. Un polynôme de … … 387 286 Realloc(degre); 388 287 389 OMatrix a(degre+1, n);288 Matrix a(degre+1, n); 390 289 391 290 for (int c=0; c<n; c++) { … … 397 296 } 398 297 399 double rc = LinFit(a,y,( OVector&)*this);298 double rc = LinFit(a,y,(Vector&)*this); 400 299 UpdateDeg(); 401 300 return rc; … … 403 302 404 303 //++ 405 double Poly::Fit( OVector const& x, OVector const& y, OVector const& erry2, int degre,406 OVector& errCoef)304 double Poly::Fit(Vector const& x, Vector const& y, Vector const& erry2, int degre, 305 Vector& errCoef) 407 306 // 408 307 // Ajustement polynomial par moindre carrés. Un polynôme de … … 419 318 errCoef.Realloc(degre+1); 420 319 421 OMatrix a(degre+1, n);320 Matrix a(degre+1, n); 422 321 423 322 for (int c=0; c<n; c++) { … … 429 328 } 430 329 431 double rc = LinFit(a,y,erry2,( OVector&)*this,errCoef);330 double rc = LinFit(a,y,erry2,(Vector&)*this,errCoef); 432 331 UpdateDeg(); 433 332 return rc; 434 333 } 435 334 335 336 //++ 337 // Poly Poly::power(int n) const 338 // 339 // Retourne le polynôme à la puissance n 340 //-- 341 342 Poly Poly::power(int n) const // a accelerer !!! 343 { 344 if (n < 0) THROW(rangeCheckErr); 345 if (n == 0) { Poly r(0); r[0] = 1; return r;} 346 if (n == 1) { return *this; } 347 return *this * power(n-1); 348 } 349 350 //++ 351 Poly Poly::operator() (Poly const& b) const 352 // 353 // Substitution d'un polynôme dans un autre. 354 //-- 355 { 356 Poly c(b.Degre()*Degre()); 357 for (int i=0; i<= Degre(); i++) 358 c += (*this)[i] * b.power(i); 359 return c; 360 } 361 362 363 ////////////////////////////////////////////////////////////////////////// 364 void ObjFileIO<Poly>::ReadSelf(PInPersist& is) 365 { 366 if(dobj==NULL) dobj=new Poly; 367 int_4 dg; 368 is >> dg; 369 dobj->Realloc(dg,true); 370 is >> *((Vector *) dobj); 371 dobj->UpdateDeg(); 372 } 373 374 void ObjFileIO<Poly>::WriteSelf(POutPersist& os) const 375 { 376 if(dobj == NULL) return; 377 dobj->UpdateDegIfDirty(); 378 dobj->Realloc(dobj->deg,true); 379 os << dobj->deg; 380 os << *((Vector *) dobj); 381 } 382 383 ////////////////////////////////////////////////////////////////////////// 436 384 int binomial(int n, int p) 437 385 { … … 442 390 } 443 391 444 445 //++ 446 // Poly Poly::power(int n) const 447 // 448 // Retourne le polynôme à la puissance n 449 //-- 450 451 Poly Poly::power(int n) const // a accelerer !!! 452 { 453 if (n < 0) THROW(rangeCheckErr); 454 if (n == 0) { Poly r(0); r[0] = 1; return r;} 455 if (n == 1) { return *this; } 456 return *this * power(n-1); 457 } 458 459 //++ 460 Poly Poly::operator() (Poly const& b) const 461 // 462 // Substitution d'un polynôme dans un autre. 463 //-- 464 { 465 Poly c(b.Degre()*Degre()); 466 for (int i=0; i<= Degre(); i++) 467 c += (*this)[i] * b.power(i); 468 return c; 469 } 470 471 392 ////////////////////////////////////////////////////////////////////////// 472 393 // ******************* POLY 2 VARIABLES ****************** 473 474 394 //++ 475 395 // Class Poly2 … … 482 402 //++ 483 403 // Links Parents 484 // OVector404 // Vector 485 405 //-- 486 406 … … 494 414 // Crée un polynôme de degrés partiels degreX et degreY. 495 415 //-- 496 : OVector((degreX+1)*(degreY+1)), dirty(0),416 :Vector((degreX+1)*(degreY+1)), dirty(0), 497 417 maxDegX(degreX), maxDegY(degreY), degX(0), degY(0), deg(0) 498 418 { … … 506 426 // de deux polynômes à une variable, p2(x,y)=px(x)py(y) 507 427 //-- 508 : OVector((polX.Degre()+1)*(polY.Degre()+1)), dirty(0),428 :Vector((polX.Degre()+1)*(polY.Degre()+1)), dirty(0), 509 429 maxDegX(polX.Degre()), maxDegY(polY.Degre()), 510 430 degX(polX.Degre()), degY(polY.Degre()), deg(degX+degY) … … 521 441 // Constructeur par copie. 522 442 //-- 523 : OVector(a), dirty(a.dirty),443 :Vector(a), dirty(a.dirty), 524 444 maxDegX(a.maxDegX), maxDegY(a.maxDegY), 525 445 degX(a.degX), degY(a.degY), deg(a.deg) … … 561 481 UpdateDegIfDirty(); 562 482 Poly2 tmp(*this); 563 OVector::Realloc((degreX+1)*(degreY+1));564 Zero();483 Vector::Realloc((degreX+1)*(degreY+1)); 484 Reset(); 565 485 maxDegX = degreX; 566 486 maxDegY = degreY; … … 631 551 632 552 //++ 633 double Poly2::Fit( OVector const& x, OVector const& y, OVector const& z,553 double Poly2::Fit(Vector const& x, Vector const& y, Vector const& z, 634 554 int degreX, int degreY) 635 555 // … … 643 563 Realloc(degreX, degreY); 644 564 645 OMatrix a((degreX+1)*(degreY+1), n);565 Matrix a((degreX+1)*(degreY+1), n); 646 566 647 567 for (int c=0; c<n; c++) { … … 657 577 } 658 578 659 double rc = LinFit(a,z,( OVector&)*this);579 double rc = LinFit(a,z,(Vector&)*this); 660 580 UpdateDeg(); 661 581 return rc; … … 664 584 665 585 //++ 666 double Poly2::Fit( OVector const& x, OVector const& y, OVector const& z,667 OVector const& errz2, int degreX, int degreY,668 OVector& errCoef)586 double Poly2::Fit(Vector const& x, Vector const& y, Vector const& z, 587 Vector const& errz2, int degreX, int degreY, 588 Vector& errCoef) 669 589 // 670 590 // Ajustement par moindre carrés z = P(x,y), degrés partiels imposés, … … 680 600 errCoef.Realloc((degreX+1)*(degreY+1)); 681 601 682 OMatrix a((degreX+1)*(degreY+1), n);602 Matrix a((degreX+1)*(degreY+1), n); 683 603 684 604 for (int c=0; c<n; c++) { … … 694 614 } 695 615 696 double rc = LinFit(a,z,errz2,( OVector&)*this,errCoef);616 double rc = LinFit(a,z,errz2,(Vector&)*this,errCoef); 697 617 UpdateDeg(); 698 618 return rc; … … 700 620 701 621 //++ 702 double Poly2::Fit( OVector const& x, OVector const& y, OVector const& z,622 double Poly2::Fit(Vector const& x, Vector const& y, Vector const& z, 703 623 int degre) 704 624 // … … 712 632 Realloc(degre, degre); // certains vaudront 0, impose. 713 633 714 OMatrix a((degre+1)*(degre+2)/2, n);715 OVector cf((degre+1)*(degre+2)/2);634 Matrix a((degre+1)*(degre+2)/2, n); 635 Vector cf((degre+1)*(degre+2)/2); 716 636 #define C_INDEX(i,j) ((i) + (j)*(2*degre+3-(j))/2) 717 637 … … 744 664 745 665 //++ 746 double Poly2::Fit( OVector const& x, OVector const& y, OVector const& z,747 OVector const& errz2, int degre,748 OVector& errCoef)666 double Poly2::Fit(Vector const& x, Vector const& y, Vector const& z, 667 Vector const& errz2, int degre, 668 Vector& errCoef) 749 669 // 750 670 // Ajustement par moindre carrés z = P(x,y), degré total imposé, … … 761 681 #define C_INDEX(i,j) ((i) + (j)*(2*degre+3-(j))/2) 762 682 763 OMatrix a((degre+1)*(degre+2)/2, n);764 OVector cf((degre+1)*(degre+2)/2);765 OVector ecf((degre+1)*(degre+2)/2);683 Matrix a((degre+1)*(degre+2)/2, n); 684 Vector cf((degre+1)*(degre+2)/2); 685 Vector ecf((degre+1)*(degre+2)/2); 766 686 767 687 for (int c=0; c<n; c++) { … … 792 712 UpdateDeg(); 793 713 return rc; 794 }795 796 797 void Poly2::ReadSelf(PInPersist& s)798 {799 int_4 dgx, dgy;800 s >> dgx >> dgy;801 802 Realloc(dgx, dgy);803 804 #ifdef DEBUG805 int r = nr; int c = nc;806 #endif807 808 OVector::ReadSelf(s);809 810 #ifdef DEBUG811 DBASSERT(r == nr && c == nc);812 #endif813 UpdateDeg();814 }815 816 void Poly2::WriteSelf(POutPersist& s) const817 {818 s << maxDegX << maxDegY;819 OVector::WriteSelf(s);820 714 } 821 715 … … 850 744 851 745 //++ 852 ostream& operator << (ostream& s, const Poly2& a)853 //854 // Impression sur stream s.855 //--856 {857 a.Print(s);858 return s;859 }860 861 862 //++863 746 // Titre Opérations 864 747 //-- … … 905 788 906 789 //++ 907 Poly2 operator+ (Poly2 const& a, Poly2 const& b)908 //909 //--910 #if HAS_NAMED_RETURN911 return c(a)912 #endif913 {914 #if !HAS_NAMED_RETURN915 Poly2 c(a);916 #endif917 c += b;918 #if !HAS_NAMED_RETURN919 return c;920 #endif921 }922 923 //++924 Poly2 operator- (Poly2 const& a, Poly2 const& b)925 //926 //--927 #if HAS_NAMED_RETURN928 return c(a)929 #endif930 {931 #if !HAS_NAMED_RETURN932 Poly2 c(a);933 #endif934 c -= b;935 #if !HAS_NAMED_RETURN936 return c;937 #endif938 }939 940 //++941 790 Poly2& Poly2::operator *= (double a) 942 791 // … … 944 793 { 945 794 for (int i=0; i<NElts(); i++) 946 data[i]*= a;795 Element(i) *= a; 947 796 948 797 return *this; … … 950 799 951 800 //++ 952 Poly2 operator * (double a, Poly2 const& b) 953 // 954 //-- 955 #if HAS_NAMED_RETURN 956 return c(b) 957 #endif 958 { 959 #if !HAS_NAMED_RETURN 960 Poly2 c(b); 961 #endif 962 c *= a; 963 #if !HAS_NAMED_RETURN 964 return c; 965 #endif 966 } 967 968 //++ 969 Poly2 operator* (Poly2 const& a, Poly2 const& b) 970 // 971 //-- 972 { 973 Poly2 c(a.DegX() + b.DegX(), a.DegY() + b.DegY()); 974 a.UpdateDegIfDirty(); 801 Poly2 Poly2::Mult(Poly2 const& b) const 802 // 803 //-- 804 { 805 Poly2 c(DegX() + b.DegX(), DegY() + b.DegY()); 806 UpdateDegIfDirty(); 975 807 b.UpdateDegIfDirty(); 976 808 977 for (int i=0; i<= a.DegX(); i++)978 for (int j=0; j<= a.DegY(); j++)809 for (int i=0; i<=DegX(); i++) 810 for (int j=0; j<=DegY(); j++) 979 811 for (int k=0; k<=b.DegX(); k++) 980 812 for (int l=0; l<=b.DegY(); l++) 981 c.Coef(i+k,j+l) += a.Coef(i,j)*b.Coef(k,l);813 c.Coef(i+k,j+l) += Coef(i,j)*b.Coef(k,l); 982 814 return c; 983 815 } … … 1029 861 } 1030 862 1031 863 ////////////////////////////////////////////////////////////////////////// 864 void ObjFileIO<Poly2>::ReadSelf(PInPersist& is) 865 { 866 if(dobj==NULL) dobj=new Poly2; 867 int_4 dgx, dgy; 868 is >> dgx >> dgy; 869 dobj->Realloc(dgx,dgy); 870 is >> *((Vector *) dobj); 871 dobj->UpdateDeg(); 872 } 873 874 void ObjFileIO<Poly2>::WriteSelf(POutPersist& os) const 875 { 876 if(dobj == NULL) return; 877 os << dobj->maxDegX << dobj->maxDegY; 878 os << *((Vector *) dobj); 879 } 880 881 882 ////////////////////////////////////////////////////////////////////////// 883 #ifdef __CXX_PRAGMA_TEMPLATES__ 884 #pragma define_template ObjFileIO<Poly> 885 #pragma define_template ObjFileIO<Poly2> 886 #endif 887 888 #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES) 889 template class ObjFileIO<Poly>; 890 template class ObjFileIO<Poly2>; 891 #endif
Note:
See TracChangeset
for help on using the changeset viewer.
