binaries/src/mafft/extensions/mxscarna_src/nrutil.h

   1 /*
   2   McCaskill's Algorithm -- The algorithm calculates a base paring probability matrix from the input of one sequence.
   3
   4   $Id: nrutil.h,v 1.0 2005/10/20 14:22 $;
   5
   6   Copyright (C) 2005 Yasuo Tabei <tabei@cb.k.u-tokyo.ac.jp>
   7
   8   This is free software with ABSOLUTELY NO WARRANTY.
   9
  10   This library is free software; you can redistribute it and/or
  11   modify it under the terms of the GNU Lesser General Public
  12   License as published by the Free Software Foundation; either
  13   version 2.1 of the License, or (at your option) any later version.
  14
  15   This library is distributed in the hope that it will be useful,
  16   but WITHOUT ANY WARRANTY; without even the implied warranty of
  17   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18   Lesser General Public License for more details.
  19
  20   You should have received a copy of the GNU Lesser General Public
  21   License along with this library; if not, write to the Free Software
  22   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  23 */
  24
  25 #ifndef _NR_UTIL_H_
  26 #define _NR_UTIL_H_
  27 #include <string>
  28 #include <cmath>
  29 #include <complex>
  30 #include <iostream>
  31 #include <cstdlib> // by katoh
  32
  33 using namespace std;
  34
  35 typedef double DP;
  36
  37 template<class T>
  38 inline const T SQR(const T a) {return a*a;}
  39
  40 template<class T>
  41 inline const T MAX(const T &a, const T &b)
  42 {return b > a ? (b) : (a);}
  43
  44 inline float MAX(const double &a, const float &b)
  45 {return b > a ? (b) : float(a);}
  46
  47 inline float MAX(const float &a, const double &b)
  48 {return b > a ? float(b) : (a);}
  49
  50 template<class T>
  51 inline const T MIN(const T &a, const T &b)
  52 {return b < a ? (b) : (a);}
  53
  54 inline float MIN(const double &a, const float &b)
  55 {return b < a ? (b) : float(a);}
  56
  57 inline float MIN(const float &a, const double &b)
  58 {return b < a ? float(b) : (a);}
  59
  60 template<class T>
  61 inline const T SIGN(const T &a, const T &b)
  62 {return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}
  63
  64 inline float SIGN(const float &a, const double &b)
  65 {return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}
  66
  67 inline float SIGN(const double &a, const float &b)
  68 {return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}
  69
  70 template<class T>
  71 inline void SWAP(T &a, T &b)
  72 {T dum=a; a=b; b=dum;}
  73 namespace NR {
  74     inline void nrerror(const string error_text)
  75 // Numerical Recipes standard error handler
  76         {
  77             cerr << "Numerical Recipes run-time error..." << endl;
  78             cerr << error_text << endl;
  79             cerr << "...now exiting to system..." << endl;
  80             exit(1);
  81         }
  82 }
  83
  84 template <class T>
  85 class NRVec {
  86  private:
  87     int nn; // size of array. upper index is nn-1
  88     T *v;
  89  public:
  90     NRVec();
  91     explicit NRVec(int n);    // Zero-based array
  92     NRVec(const T &a, int n); //initialize to constant value
  93     NRVec(const T *a, int n); // Initialize to array
  94     NRVec(const NRVec &rhs);  // Copy constructor
  95     NRVec & operator=(const NRVec &rhs); //assignment
  96     NRVec & operator=(const T &a); //assign a to every element
  97     inline T & operator[](const int i); //i¡Çth element
  98     inline const T & operator[](const int i) const;
  99     void Allocator(int i);
 100     inline int size() const;
 101     ~NRVec();
 102 };
 103
 104 template <class T>
 105 NRVec<T>::NRVec() : nn(0), v(0) {}
 106
 107 template <class T>
 108 NRVec<T>::NRVec(int n) : nn(n), v(new T[n]) {}
 109
 110 template <class T>
 111 NRVec<T>::NRVec(const T& a, int n) : nn(n), v(new T[n])
 112 {
 113     for(int i=0; i<n; i++)
 114         v[i] = a;
 115 }
 116
 117 template <class T>
 118 NRVec<T>::NRVec(const T *a, int n) : nn(n), v(new T[n])
 119 {
 120 for(int i=0; i<n; i++)
 121     v[i] = *a++;
 122 }
 123
 124 template <class T>
 125 void NRVec<T>::Allocator(int n = 0)
 126 {
 127    v = new T[n];
 128 }
 129
 130 template <class T>
 131 NRVec<T>::NRVec(const NRVec<T> &rhs) : nn(rhs.nn), v(new T[nn])
 132 {
 133     for(int i=0; i<nn; i++)
 134         v[i] = rhs[i];
 135 }
 136
 137 template <class T>
 138 NRVec<T> & NRVec<T>::operator=(const NRVec<T> &rhs)
 139 // postcondition: normal assignment via copying has been performed;
 140 // if vector and rhs were different sizes, vector
 141 // has been resized to match the size of rhs
 142 {
 143     if (this != &rhs)
 144 {
 145     if (nn != rhs.nn) {
 146         if (v != 0) delete [] (v);
 147         nn=rhs.nn;
 148         v= new T[nn];
 149     }
 150     for (int i=0; i<nn; i++)
 151         v[i]=rhs[i];
 152 }
 153     return *this;
 154 }
 155
 156 template <class T>
 157 NRVec<T> & NRVec<T>::operator=(const T &a) //assign a to every element
 158 {
 159     for (int i=0; i<nn; i++)
 160         v[i]=a;
 161     return *this;
 162 }
 163
 164 template <class T>
 165 inline T & NRVec<T>::operator[](const int i) //subscripting
 166 {
 167     return v[i];
 168 }
 169
 170 template <class T>
 171 inline const T & NRVec<T>::operator[](const int i) const //subscripting
 172 {
 173     return v[i];
 174 }
 175
 176 template <class T>
 177 inline int NRVec<T>::size() const
 178 {
 179     return nn;
 180 }
 181
 182 template <class T>
 183 NRVec<T>::~NRVec()
 184 {
 185     if (v != 0)
 186         delete[] (v);
 187 }
 188
 189 template <class T>
 190 class NRMat {
 191  private:
 192     int nn;
 193     int mm;
 194     T **v;
 195  public:
 196     NRMat();
 197     NRMat(int n, int m); // Zero-based array
 198     NRMat(const T &a, int n, int m); //Initialize to constant
 199     NRMat(const T *a, int n, int m); // Initialize to array
 200     NRMat(const NRMat &rhs); // Copy constructor
 201     void Allocator(int n, int m);
 202     void Allocator(const T &a, int n, int m);
 203     void Allocator(const T *a, int n, int m);
 204     NRMat & operator=(const NRMat &rhs); //assignment
 205     NRMat & operator=(const T &a); //assign a to every element
 206     inline T* operator[](const int i); //subscripting: pointer to row i
 207     inline const T* operator[](const int i) const;
 208     inline T &  ref(const int i, const int j);
 209     inline const T ref(const int i, const int j) const;
 210     inline int nrows() const;
 211     inline int ncols() const;
 212     ~NRMat();
 213 };
 214
 215 template <class T>
 216 NRMat<T>::NRMat() : nn(0), mm(0), v(0) {}
 217
 218 template <class T>
 219 NRMat<T>::NRMat(int n, int m) : nn(n), mm(m), v(new T*[n])
 220 {
 221     v[0] = new T[m*n];
 222     for (int i=1; i< n; i++)
 223         v[i] = v[i-1] + m;
 224 }
 225
 226 template <class T>
 227 NRMat<T>::NRMat(const T &a, int n, int m) : nn(n), mm(m), v(new T*[n])
 228 {
 229     int i,j;
 230     v[0] = new T[m*n];
 231     for (i=1; i< n; i++)
 232         v[i] = v[i-1] + m;
 233     for (i=0; i< n; i++)
 234         for (j=0; j<m; j++)
 235             v[i][j] = a;
 236 }
 237
 238 template <class T>
 239 NRMat<T>::NRMat(const T *a, int n, int m) : nn(n), mm(m), v(new T*[n])
 240 {
 241     int i,j;
 242     v[0] = new T[m*n];
 243     for (i=1; i< n; i++)
 244         v[i] = v[i-1] + m;
 245     for (i=0; i< n; i++)
 246         for (j=0; j<m; j++)
 247             v[i][j] = *a++;
 248 }
 249
 250 template <class T>
 251 void NRMat<T>::Allocator(int n, int m)
 252 {
 253     if( v != 0 ) {
 254         delete[] (v[0]); delete (v);
 255     }
 256
 257     nn = n; mm = m; v = new T*[n];
 258
 259     v[0] = new T[m*n];
 260     for (int i=1; i< n; i++)
 261         v[i] = v[i-1] + m;
 262 }
 263
 264 template <class T>
 265 void NRMat<T>::Allocator(const T &a, int n, int m)
 266 {
 267     if( v != 0 ) {
 268         delete[] (v[0]); delete (v);
 269     }
 270
 271     int i,j;
 272
 273     nn = n; mm = m; v = new T*[n];
 274
 275     v[0] = new T[m*n];
 276     for (i=1; i< n; i++)
 277         v[i] = v[i-1] + m;
 278     for (i=0; i< n; i++)
 279         for (j=0; j<m; j++)
 280             v[i][j] = a;
 281 }
 282
 283 template <class T>
 284 void NRMat<T>::Allocator(const T *a, int n, int m)
 285 {
 286     if( v != 0 ) {
 287         delete[] (v[0]); delete (v);
 288     }
 289
 290     int i,j;
 291
 292     nn = n; mm = m; v = new T*[n];
 293
 294     v[0] = new T[m*n];
 295     for (i=1; i< n; i++)
 296         v[i] = v[i-1] + m;
 297     for (i=0; i< n; i++)
 298         for (j=0; j<m; j++)
 299             v[i][j] = *a++;
 300 }
 301
 302 template <class T>
 303 NRMat<T>::NRMat(const NRMat &rhs) : nn(rhs.nn), mm(rhs.mm), v(new T*[nn])
 304 {
 305     int i,j;
 306     v[0] = new T[mm*nn];
 307     for (i=1; i< nn; i++)
 308         v[i] = v[i-1] + mm;
 309     for (i=0; i< nn; i++)
 310         for (j=0; j<mm; j++)
 311             v[i][j] = rhs[i][j];
 312 }
 313 template <class T>
 314 NRMat<T> & NRMat<T>::operator=(const NRMat<T> &rhs)
 315 // postcondition: normal assignment via copying has been performed;
 316 // if matrix and rhs were different sizes, matrix
 317 // has been resized to match the size of rhs
 318 {
 319     if (this != &rhs) {
 320         int i,j;
 321         if (nn != rhs.nn || mm != rhs.mm) {
 322             if (v != 0) {
 323                 delete[] (v[0]);
 324                 delete[] (v);
 325             }
 326             nn=rhs.nn;
 327             mm=rhs.mm;
 328             v = new T*[nn];
 329             v[0] = new T[mm*nn];
 330         }
 331         for (i=1; i< nn; i++)
 332             v[i] = v[i-1] + mm;
 333         for (i=0; i< nn; i++)
 334             for (j=0; j<mm; j++)
 335                 v[i][j] = rhs[i][j];
 336     }
 337     return *this;
 338 }
 339
 340 template <class T>
 341 NRMat<T> & NRMat<T>::operator=(const T &a) //assign a to every element
 342 {
 343     for (int i=0; i< nn; i++)
 344         for (int j=0; j<mm; j++)
 345             v[i][j] = a;
 346     return *this;
 347 }
 348
 349 template <class T>
 350 inline T* NRMat<T>::operator[](const int i) //subscripting: pointer to row i
 351 {
 352     return v[i];
 353 }
 354
 355 template <class T>
 356 inline const T* NRMat<T>::operator[](const int i) const
 357 {
 358     return v[i];
 359 }
 360
 361 template <class T>
 362 inline T &  NRMat<T>::ref(const int i, const int j)
 363 {
 364     return v[i][j];
 365 }
 366
 367 template <class T>
 368 inline const T NRMat<T>::ref(const int i, const int j) const
 369 {
 370     return v[i][j];
 371 }
 372
 373 template <class T>
 374 inline int NRMat<T>::nrows() const
 375 {
 376     return nn;
 377 }
 378
 379 template <class T>
 380 inline int NRMat<T>::ncols() const
 381 {
 382     return mm;
 383 }
 384
 385 template <class T>
 386 NRMat<T>::~NRMat()
 387 {
 388     if (v != 0) {
 389         delete[] (v[0]);
 390         delete[] (v);
 391     }
 392 }
 393
 394 template <class T>
 395 class NRMat3d {
 396  private:
 397     int nn;
 398     int mm;
 399     int kk;
 400     T ***v;
 401  public:
 402     NRMat3d();
 403     NRMat3d(int n, int m, int k);
 404     inline void Allocator(int n, int m, int k);
 405     inline T** operator[](const int i); //subscripting: pointer to row i
 406     inline const T* const * operator[](const int i) const;
 407     inline int dim1() const;
 408     inline int dim2() const;
 409     inline int dim3() const;
 410     ~NRMat3d();
 411 };
 412
 413 template <class T>
 414 NRMat3d<T>::NRMat3d(): nn(0), mm(0), kk(0), v(0) {}
 415 template <class T>
 416 NRMat3d<T>::NRMat3d(int n, int m, int k) : nn(n), mm(m), kk(k), v(new T**[n])
 417 {
 418     int i,j;
 419     v[0] = new T*[n*m];
 420     v[0][0] = new T[n*m*k];
 421     for(j=1; j<m; j++)
 422         v[0][j] = v[0][j-1] + k;
 423     for(i=1; i<n; i++) {
 424         v[i] = v[i-1] + m;
 425         v[i][0] = v[i-1][0] + m*k;
 426         for(j=1; j<m; j++)
 427             v[i][j] = v[i][j-1] + k;
 428     }
 429 }
 430
 431 template <class T>
 432 inline void NRMat3d<T>::Allocator(int n, int m, int k)
 433 {
 434     int i,j;
 435     v[0] = new T*[n*m];
 436     v[0][0] = new T[n*m*k];
 437     for(j=1; j<m; j++)
 438         v[0][j] = v[0][j-1] + k;
 439     for(i=1; i<n; i++) {
 440         v[i] = v[i-1] + m;
 441         v[i][0] = v[i-1][0] + m*k;
 442         for(j=1; j<m; j++)
 443             v[i][j] = v[i][j-1] + k;
 444     }
 445 }
 446
 447 template <class T>
 448 inline T** NRMat3d<T>::operator[](const int i) //subscripting: pointer to row i
 449 {
 450     return v[i];
 451 }
 452
 453 template <class T>
 454 inline const T* const * NRMat3d<T>::operator[](const int i) const
 455 {
 456     return v[i];
 457 }
 458
 459 template <class T>
 460 inline int NRMat3d<T>::dim1() const
 461 {
 462     return nn;
 463 }
 464
 465 template <class T>
 466 inline int NRMat3d<T>::dim2() const
 467 {
 468     return mm;
 469 }
 470
 471 template <class T>
 472 inline int NRMat3d<T>::dim3() const
 473 {
 474     return kk;
 475 }
 476
 477 template <class T>
 478 NRMat3d<T>::~NRMat3d()
 479 {
 480     if (v != 0) {
 481         delete[] (v[0][0]);
 482         delete[] (v[0]);
 483         delete[] (v);
 484     }
 485 }
 486
 487 //The next 3 classes are used in artihmetic coding, Huffman coding, and
 488 //wavelet transforms respectively. This is as good a place as any to put them!
 489 class arithcode {
 490  private:
 491     NRVec<unsigned long> *ilob_p,*iupb_p,*ncumfq_p;
 492  public:
 493     NRVec<unsigned long> &ilob,&iupb,&ncumfq;
 494     unsigned long jdif,nc,minint,nch,ncum,nrad;
 495     arithcode(unsigned long n1, unsigned long n2, unsigned long n3)
 496         : ilob_p(new NRVec<unsigned long>(n1)),
 497         iupb_p(new NRVec<unsigned long>(n2)),
 498         ncumfq_p(new NRVec<unsigned long>(n3)),
 499         ilob(*ilob_p),iupb(*iupb_p),ncumfq(*ncumfq_p) {}
 500     ~arithcode() {
 501         if (ilob_p != 0) delete ilob_p;
 502         if (iupb_p != 0) delete iupb_p;
 503         if (ncumfq_p != 0) delete ncumfq_p;
 504     }
 505 };
 506
 507 class huffcode {
 508  private:
 509     NRVec<unsigned long> *icod_p,*ncod_p,*left_p,*right_p;
 510  public:
 511     NRVec<unsigned long> &icod,&ncod,&left,&right;
 512     int nch,nodemax;
 513     huffcode(unsigned long n1, unsigned long n2, unsigned long n3,
 514              unsigned long n4) :
 515         icod_p(new NRVec<unsigned long>(n1)),
 516         ncod_p(new NRVec<unsigned long>(n2)),
 517         left_p(new NRVec<unsigned long>(n3)),
 518         right_p(new NRVec<unsigned long>(n4)),
 519         icod(*icod_p),ncod(*ncod_p),left(*left_p),right(*right_p) {}
 520     ~huffcode() {
 521         if (icod_p != 0) delete icod_p;
 522         if (ncod_p != 0) delete ncod_p;
 523         if (left_p != 0) delete left_p;
 524         if (right_p != 0) delete right_p;
 525     }
 526 };
 527
 528 class wavefilt {
 529  private:
 530     NRVec<DP> *cc_p,*cr_p;
 531  public:
 532     int ncof,ioff,joff;
 533     NRVec<DP> &cc,&cr;
 534     wavefilt() : cc(*cc_p),cr(*cr_p) {}
 535     wavefilt(const DP *a, const int n) : //initialize to array
 536         cc_p(new NRVec<DP>(n)),cr_p(new NRVec<DP>(n)),
 537         ncof(n),ioff(-(n >> 1)),joff(-(n >> 1)),cc(*cc_p),cr(*cr_p) {
 538         int i;
 539         for (i=0; i<n; i++)
 540             cc[i] = *a++;
 541         DP sig = -1.0;
 542         for (i=0; i<n; i++) {
 543             cr[n-1-i]=sig*cc[i];
 544             sig = -sig;
 545         }
 546     }
 547     ~wavefilt() {
 548         if (cc_p != 0) delete cc_p;
 549         if (cr_p != 0) delete cr_p;
 550     }
 551 };
 552
 553
 554 /* Triangle Matrix Class
 555    ---------------------------------------------------------
 556    |v[0][0]|v[0][1]|      |            |         |v[0][n-1]|
 557    |-------|-------|------|------------|---------|---------|
 558    |v[1][1]|v[1][2]|      |            |v[1][n-2]|         |
 559    |-------|-------|------|------------|---------|---------|
 560    |       |       |      |            |         |         |
 561    |-------|-------|------|------------|---------|---------|
 562    |       |                                               |
 563    |       |                                               |
 564    |       |                                               |
 565    |-------|-----------------------------------------------|
 566    |       |                                               |
 567    |-------|-----------------------------------------------|
 568    |v[n-2][0]|v[n-2][1]|                                   |
 569    |-------|-----------------------------------------------|
 570    |v[n-1][0]|                                             |
 571    |-------------------------------------------------------|
 572  */
 573 template <class T>
 574 class Trimat {
 575  private:
 576     int nn;
 577     T **v;
 578     inline T* operator[](const int i); //subscripting: pointer to row i
 579     inline const T* operator[](const int i) const;
 580  public:
 581     Trimat();
 582     Trimat(int n); // Zero-based array
 583     Trimat(const T &a, int n); //Initialize to constant
 584     Trimat(const T *a, int n); // Initialize to array
 585     Trimat(const Trimat &rhs); // Copy constructor
 586     void Allocator(int n);
 587     void Allocator(const T &a, int n);
 588     void Allocator(const T *a, int n);
 589     Trimat & operator=(const Trimat &rhs); //assignment
 590     Trimat & operator=(const T &a); //assign a to every element
 591     inline T & ref(const int i, const int j);
 592     inline T * getPointer(const int i, const int j);
 593     inline T * begin() const;
 594     inline T * end() const;
 595     inline const T ref(const int i, const int j) const;
 596     inline int nrows() const;
 597     ~Trimat();
 598 };
 599
 600 template <class T>
 601 Trimat<T>::Trimat() : nn(0), v(0) {}
 602
 603 template <class T>
 604 Trimat<T>::Trimat(int n) : nn(n), v(new T*[n])
 605 {
 606     v[0] = new T[n*(n+1)/2];
 607     for (int i=1; i< n; i++)
 608         v[i] = v[i-1] + (n-i+1);
 609
 610     for (int i=0; i< n; i++)
 611         for (int j=0; j<(n-i); j++)
 612             v[i][j] = 0;
 613 }
 614 template <class T>
 615 Trimat<T>::Trimat(const T &a, int n) : nn(n), v(new T*[n])
 616 {
 617     int i,j;
 618     v[0] = new T[n*(n+1)/2];
 619     for (i=1; i< n; i++)
 620         v[i] = v[i-1] + (n-i+1);
 621     for (i=0; i< n; i++)
 622         for (j=0; j<(n-i); j++)
 623             v[i][j] = a;
 624 }
 625
 626 template <class T>
 627 Trimat<T>::Trimat(const T *a, int n) : nn(n), v(new T*[n])
 628 {
 629     int i,j;
 630     v[0] = new T[n*(n+1)/2];
 631     for (i=1; i< n; i++)
 632         v[i] = v[i-1] + (n-i+1);
 633     for (i=0; i< n; i++)
 634         for (j=0; j<(n-i); j++)
 635             v[i][j] = *a++;
 636 }
 637
 638
 639 template <class T>
 640 void Trimat<T>::Allocator(int n)
 641 {
 642     nn = n; v = new T*[n];
 643
 644     v[0] = new T[n*(n+1)/2];
 645     for (int i=1; i< n; i++)
 646         v[i] = v[i-1] + (n-i+1);
 647 }
 648
 649 template <class T>
 650 void Trimat<T>::Allocator(const T &a, int n)
 651 {
 652     nn = n; v = new T*[n];
 653
 654     int i,j;
 655     v[0] = new T[n*(n+1)/2];
 656     for (i=1; i < n; i++)
 657         v[i] = v[i-1] + (n-i+1);
 658     for (i=0; i < n; i++)
 659         for (j=0; j < (n-i); j++)
 660             v[i][j] = a;
 661 }
 662
 663 template <class T>
 664 void Trimat<T>::Allocator(const T *a, int n)
 665 {
 666     nn = n; v = new T*[n];
 667     int i,j;
 668     v[0] = new T[n*(n+1)/2];
 669     for (i=1; i< n; i++)
 670         v[i] = v[i-1] + (n-i+1);
 671     for (i=0; i< n; i++)
 672         for (j=0; j<(n-i); j++)
 673             v[i][j] = *a++;
 674 }
 675
 676
 677 template <class T>
 678 Trimat<T>::Trimat(const Trimat &rhs) : nn(rhs.nn), v(new T*[nn])
 679 {
 680     int i,j;
 681     v[0] = new T[nn*(nn+1)/2];
 682     for (i=1; i< nn; i++)
 683         v[i] = v[i-1] + (nn-i+1);
 684     for (i=0; i< nn; i++)
 685         for (j=0; j<(nn-i); j++)
 686             v[i][j] = rhs[i][j];
 687 }
 688 template <class T>
 689 Trimat<T> & Trimat<T>::operator=(const Trimat<T> &rhs)
 690 // postcondition: normal assignment via copying has been performed;
 691 // if matrix and rhs were different sizes, matrix
 692 // has been resized to match the size of rhs
 693 {
 694     if (this != &rhs) {
 695         int i,j;
 696         if (nn != rhs.nn) {
 697             if (v != 0) {
 698                 delete[] (v[0]);
 699                 delete[] (v);
 700             }
 701             nn=rhs.nn;
 702             v = new T*[nn];
 703             v[0] = new T[nn*(nn+1)/2];
 704         }
 705         for (i=1; i< nn; i++)
 706             v[i] = v[i-1] + (nn-i+1);
 707         for (i=0; i< nn; i++)
 708             for (j=0; j<(nn-i); j++)
 709                 v[i][j] = rhs[i][j];
 710     }
 711     return *this;
 712 }
 713
 714 template <class T>
 715 Trimat<T> & Trimat<T>::operator=(const T &a) //assign a to every element
 716 {
 717     for (int i=0; i< nn; i++)
 718         for (int j=0; j<nn-i; j++)
 719             v[i][j] = a;
 720     return *this;
 721 }
 722
 723 template <class T>
 724 inline T &  Trimat<T>::ref(const int i, const int j)
 725 {
 726     return v[i][j-i];
 727 }
 728
 729 template <class T>
 730 inline const T Trimat<T>::ref(const int i, const int j) const
 731 {
 732     return v[i][j-i];
 733 }
 734
 735 template <class T>
 736 inline T * Trimat<T>::getPointer(const int i, const int j)
 737 {
 738     return &v[i][j-i];
 739 }
 740
 741 template <class T>
 742 inline T * Trimat<T>::begin() const
 743 {
 744     return &v[0][0];
 745 }
 746
 747 template <class T>
 748 inline T * Trimat<T>::end() const
 749 {
 750     return (&v[nn-1][0] + 1);
 751 }
 752
 753 template <class T>
 754 inline int Trimat<T>::nrows() const
 755 {
 756     return nn;
 757 }
 758
 759 template <class T>
 760 Trimat<T>::~Trimat()
 761 {
 762     if (v != 0) {
 763         delete[] (v[0]);
 764         delete[] (v);
 765     }
 766 }
 767
 768
 769 /* Triangle Vertical Matrix Class
 770    ---------------------------------------------------------
 771    |v[0][0]|v[1][0]|      |            |         |v[n-1][0]|
 772    |-------|-------|------|------------|---------|---------|
 773    |v[0][1]|v[1][1]|      |            |v[n-2][1]|         |
 774    |-------|-------|------|------------|---------|---------|
 775    |       |       |      |            |         |         |
 776    |-------|-------|------|------------|---------|---------|
 777    |       |                                               |
 778    |       |                                               |
 779    |       |                                               |
 780    |-------|-----------------------------------------------|
 781    |       |                                               |
 782    |-------|-----------------------------------------------|
 783    |v[0][n-2]|v[n-2][n-2]|                                   |
 784    |-------|-----------------------------------------------|
 785    |v[0][n-1]|                                             |
 786    |-------------------------------------------------------|
 787  */
 788 template <class T>
 789 class TriVertMat {
 790  private:
 791     int nn;
 792     T **v;
 793     inline T* operator[](const int i); //subscripting: pointer to row i
 794     inline const T* operator[](const int i) const;
 795  public:
 796     TriVertMat();
 797     TriVertMat(int n); // Zero-based array
 798     TriVertMat(const T &a, int n); //Initialize to constant
 799     TriVertMat(const T *a, int n); // Initialize to array
 800     TriVertMat(const TriVertMat &rhs); // Copy constructor
 801     void Allocator(int n);
 802     void Allocator(const T &a, int n);
 803     void Allocator(const T *a, int n);
 804     TriVertMat & operator=(const TriVertMat &rhs); //assignment
 805     TriVertMat & operator=(const T &a); //assign a to every element
 806     inline T & ref(const int i, const int j);
 807     inline T * getPointer(const int i, const int j);
 808     inline const T ref(const int i, const int j) const;
 809     inline int nrows() const;
 810     ~TriVertMat();
 811 };
 812
 813 template <class T>
 814 TriVertMat<T>::TriVertMat() : nn(0), v(0) {}
 815
 816 template <class T>
 817 TriVertMat<T>::TriVertMat(int n) : nn(n), v(new T*[n])
 818 {
 819     v[0] = new T[n*(n+1)/2];
 820     for (int i=1; i< n; i++)
 821         v[i] = v[i-1] + i;
 822 }
 823
 824 template <class T>
 825 TriVertMat<T>::TriVertMat(const T &a, int n) : nn(n), v(new T*[n])
 826 {
 827     int i,j;
 828     v[0] = new T[n*(n+1)/2];
 829     for (i=1; i< n; i++)
 830         v[i] = v[i-1] + i;
 831     for (i=0; i< n; i++)
 832         for (j=0; j<(n-i); j++)
 833             v[i][j] = a;
 834 }
 835
 836 template <class T>
 837 TriVertMat<T>::TriVertMat(const T *a, int n) : nn(n), v(new T*[n])
 838 {
 839     int i,j;
 840     v[0] = new T[n*(n+1)/2];
 841     for (i=1; i< n; i++)
 842         v[i] = v[i-1] + i;
 843     for (i=0; i< n; i++)
 844         for (j=0; j<(n-i); j++)
 845             v[i][j] = *a++;
 846 }
 847
 848
 849 template <class T>
 850 void TriVertMat<T>::Allocator(int n)
 851 {
 852     nn = n; v = new T*[n];
 853
 854     v[0] = new T[n*(n+1)/2];
 855     for (int i=1; i< n; i++)
 856         v[i] = v[i-1] + i;
 857 }
 858
 859 template <class T>
 860 void TriVertMat<T>::Allocator(const T &a, int n)
 861 {
 862     nn = n; v = new T*[n];
 863
 864     int i,j;
 865     v[0] = new T[n*(n+1)/2];
 866     for (i=1; i< n; i++)
 867         v[i] = v[i-1] + i;
 868     for (i=0; i< n; i++)
 869         for (j=0; j<(n-i); j++)
 870             v[i][j] = a;
 871 }
 872
 873 template <class T>
 874 void TriVertMat<T>::Allocator(const T *a, int n)
 875 {
 876     nn = n; v = new T*[n];
 877     int i,j;
 878     v[0] = new T[n*(n+1)/2];
 879     for (i=1; i< n; i++)
 880         v[i] = v[i-1] + i;
 881     for (i=0; i< n; i++)
 882         for (j=0; j<(n-i); j++)
 883             v[i][j] = *a++;
 884 }
 885
 886
 887 template <class T>
 888 TriVertMat<T>::TriVertMat(const TriVertMat &rhs) : nn(rhs.nn), v(new T*[nn])
 889 {
 890     int i,j;
 891     v[0] = new T[nn*(nn+1)/2];
 892     for (i=1; i< nn; i++)
 893         v[i] = v[i-1] + i;
 894     for (i=0; i< nn; i++)
 895         for (j=0; j<(nn-i); j++)
 896             v[i][j] = rhs[i][j];
 897 }
 898 template <class T>
 899 TriVertMat<T> & TriVertMat<T>::operator=(const TriVertMat<T> &rhs)
 900 // postcondition: normal assignment via copying has been performed;
 901 // if matrix and rhs were different sizes, matrix
 902 // has been resized to match the size of rhs
 903 {
 904     if (this != &rhs) {
 905         int i,j;
 906         if (nn != rhs.nn) {
 907             if (v != 0) {
 908                 delete[] (v[0]);
 909                 delete[] (v);
 910             }
 911             nn=rhs.nn;
 912             v = new T*[nn];
 913             v[0] = new T[nn*(nn+1)/2];
 914         }
 915         for (i=1; i< nn; i++)
 916             v[i] = v[i-1] + i;
 917         for (i=0; i< nn; i++)
 918             for (j=0; j<(nn-i); j++)
 919                 v[i][j] = rhs[i][j];
 920     }
 921     return *this;
 922 }
 923
 924 template <class T>
 925 TriVertMat<T> & TriVertMat<T>::operator=(const T &a) //assign a to every element
 926 {
 927     for (int i=0; i< nn; i++)
 928         for (int j=0; j<nn-i; j++)
 929             v[i][j] = a;
 930     return *this;
 931 }
 932
 933 template <class T>
 934 inline T &  TriVertMat<T>::ref(const int i, const int j)
 935 {
 936     return v[j][i];
 937 }
 938
 939 template <class T>
 940 inline const T TriVertMat<T>::ref(const int i, const int j) const
 941 {
 942     return v[j][i];
 943 }
 944
 945 template <class T>
 946 inline T * TriVertMat<T>::getPointer(const int i, const int j)
 947 {
 948     return &v[j][i];
 949 }
 950
 951 template <class T>
 952 inline int TriVertMat<T>::nrows() const
 953 {
 954     return nn;
 955 }
 956
 957 template <class T>
 958 TriVertMat<T>::~TriVertMat()
 959 {
 960     if (v != 0) {
 961         delete[] (v[0]);
 962         delete[] (v);
 963     }
 964 }
 965
 966
 967 //Overloaded complex operations to handle mixed float and double
 968 //This takes care of e.g. 1.0/z, z complex<float>
 969 inline const complex<float> operator+(const double &a,
 970                                       const complex<float> &b) { return float(a)+b; }
 971 inline const complex<float> operator+(const complex<float> &a,
 972                                       const double &b) { return a+float(b); }
 973 inline const complex<float> operator-(const double &a,
 974                                       const complex<float> &b) { return float(a)-b; }
 975 inline const complex<float> operator-(const complex<float> &a,
 976                                       const double &b) { return a-float(b); }
 977 inline const complex<float> operator*(const double &a,
 978                                       const complex<float> &b) { return float(a)*b; }
 979 inline const complex<float> operator*(const complex<float> &a,
 980                                       const double &b) { return a*float(b); }
 981 inline const complex<float> operator/(const double &a,
 982                                       const complex<float> &b) { return float(a)/b; }
 983 inline const complex<float> operator/(const complex<float> &a,
 984                                       const double &b) { return a/float(b); }
 985 //some compilers choke on pow(float,double) in single precision. also atan2
 986 inline float pow (float x, double y) {return pow(double(x),y);}
 987 inline float pow (double x, float y) {return pow(x,double(y));}
 988 inline float atan2 (float x, double y) {return atan2(double(x),y);}
 989 inline float atan2 (double x, float y) {return atan2(x,double(y));}
 990
 991 #endif /* _NR_UTIL_H_ */