72 complex(
const double &_x,
const double &_y):
x(_x),
y(_y){};
75 complex& operator= (
const double&
v){
x =
v;
y = 0.0;
return *
this; };
90 if( fabs(z.
y)<fabs(z.
x) )
94 result.
x = (z.
x+z.
y*e)/f;
95 result.
y = (z.
y-z.
x*e)/f;
101 result.
x = (z.
y+z.
x*e)/f;
102 result.
y = (-z.
x+z.
y*e)/f;
148 pData(const_cast<
T*>(Data)),iLength(Length),iStep(Step){};
202 for(i=imax; i!=0; i--)
204 r += (*p1)*(*p2) + p1[1]*p2[1] + p1[2]*p2[2] + p1[3]*p2[3];
209 r += (*(p1++))*(*(p2++));
217 int offset11 = v1.
GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
218 int offset21 = v2.
GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
224 for(i=0; i<imax; i++)
226 r += (*p1)*(*p2) + p1[offset11]*p2[offset21] + p1[offset12]*p2[offset22] + p1[offset13]*p2[offset23];
257 for(i=imax; i!=0; i--)
273 int offset11 = vdst.
GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
274 int offset21 = vsrc.
GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
279 for(i=0; i<imax; i++)
282 p1[offset11] = p2[offset21];
283 p1[offset12] = p2[offset22];
284 p1[offset13] = p2[offset23];
315 for(i=0; i<imax; i++)
331 int offset11 = vdst.
GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
332 int offset21 = vsrc.
GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
337 for(i=imax; i!=0; i--)
340 p1[offset11] = -p2[offset21];
341 p1[offset12] = -p2[offset22];
342 p1[offset13] = -p2[offset23];
360 template<
class T,
class T2>
373 for(i=imax; i!=0; i--)
383 *(p1++) = alpha*(*(p2++));
391 int offset11 = vdst.
GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
392 int offset21 = vsrc.
GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
397 for(i=0; i<imax; i++)
400 p1[offset11] = alpha*p2[offset21];
401 p1[offset12] = alpha*p2[offset22];
402 p1[offset13] = alpha*p2[offset23];
433 for(i=imax; i!=0; i--)
451 int offset11 = vdst.
GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
452 int offset21 = vsrc.
GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
457 for(i=0; i<imax; i++)
460 p1[offset11] += p2[offset21];
461 p1[offset12] += p2[offset22];
462 p1[offset13] += p2[offset23];
480 template<
class T,
class T2>
493 for(i=imax; i!=0; i--)
496 p1[1] += alpha*p2[1];
497 p1[2] += alpha*p2[2];
498 p1[3] += alpha*p2[3];
503 *(p1++) += alpha*(*(p2++));
511 int offset11 = vdst.
GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
512 int offset21 = vsrc.
GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
517 for(i=0; i<imax; i++)
520 p1[offset11] += alpha*p2[offset21];
521 p1[offset12] += alpha*p2[offset22];
522 p1[offset13] += alpha*p2[offset23];
553 for(i=imax; i!=0; i--)
571 int offset11 = vdst.
GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
572 int offset21 = vsrc.
GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
577 for(i=0; i<imax; i++)
580 p1[offset11] -= p2[offset21];
581 p1[offset12] -= p2[offset22];
582 p1[offset13] -= p2[offset23];
600 template<
class T,
class T2>
603 vadd(vdst, vsrc, -alpha);
610 template<
class T,
class T2>
621 for(i=imax; i!=0; i--)
638 int offset11 = vdst.
GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
642 for(i=0; i<imax; i++)
645 p1[offset11] *=
alpha;
646 p1[offset12] *=
alpha;
647 p1[offset13] *=
alpha;
686 m_Vec =
new T[m_iVecSize];
687 #ifndef UNSAFE_MEM_COPY 688 for(
int i=0;
i<m_iVecSize;
i++)
691 memcpy(m_Vec, rhs.
m_Vec, m_iVecSize*
sizeof(
T));
711 m_Vec =
new T[m_iVecSize];
712 #ifndef UNSAFE_MEM_COPY 713 for(
int i=0;
i<m_iVecSize;
i++)
716 memcpy(m_Vec, rhs.
m_Vec, m_iVecSize*
sizeof(
T));
730 return m_Vec[ i-m_iLow ];
739 return m_Vec[ i-m_iLow ];
749 m_iVecSize = iHigh-iLow+1;
750 m_Vec =
new T[m_iVecSize];
756 setbounds(iLow, iHigh);
757 for(
int i=iLow;
i<=iHigh;
i++)
758 (*
this)(
i) = pContent[
i-iLow];
786 if( iStart>iEnd || wrongIdx(iStart) || wrongIdx(iEnd) )
795 if( iStart>iEnd || wrongIdx(iStart) || wrongIdx(iEnd) )
801 bool wrongIdx(
int i)
const {
return i<m_iLow || i>m_iHigh; };
805 long m_iLow, m_iHigh;
840 m_Vec =
new T[m_iVecSize];
841 #ifndef UNSAFE_MEM_COPY 842 for(
int i=0;
i<m_iVecSize;
i++)
845 memcpy(m_Vec, rhs.
m_Vec, m_iVecSize*
sizeof(
T));
867 m_Vec =
new T[m_iVecSize];
868 #ifndef UNSAFE_MEM_COPY 869 for(
int i=0;
i<m_iVecSize;
i++)
872 memcpy(m_Vec, rhs.
m_Vec, m_iVecSize*
sizeof(
T));
886 return m_Vec[ m_iConstOffset + i2 +i1*m_iLinearMember];
895 return m_Vec[ m_iConstOffset + i2 +i1*m_iLinearMember];
898 void setbounds(
int iLow1,
int iHigh1,
int iLow2,
int iHigh2 )
902 m_iVecSize = (iHigh1-iLow1+1)*(iHigh2-iLow2+1);
903 m_Vec =
new T[m_iVecSize];
908 m_iConstOffset = -m_iLow2-m_iLow1*(m_iHigh2-m_iLow2+1);
909 m_iLinearMember = (m_iHigh2-m_iLow2+1);
912 void setcontent(
int iLow1,
int iHigh1,
int iLow2,
int iHigh2,
const T *pContent )
914 setbounds(iLow1, iHigh1, iLow2, iHigh2);
915 for(
int i=0;
i<m_iVecSize;
i++)
916 m_Vec[
i]=pContent[
i];
931 return iBoundNum==1 ? m_iLow1 : m_iLow2;
936 return iBoundNum==1 ? m_iHigh1 : m_iHigh2;
941 if( (iRowStart>iRowEnd) || wrongColumn(iColumn) || wrongRow(iRowStart) ||wrongRow(iRowEnd) )
944 return raw_vector<T>(&((*this)(iRowStart, iColumn)), iRowEnd-iRowStart+1, m_iLinearMember);
947 raw_vector<T>
getrow(
int iRow,
int iColumnStart,
int iColumnEnd)
949 if( (iColumnStart>iColumnEnd) || wrongRow(iRow) || wrongColumn(iColumnStart) || wrongColumn(iColumnEnd))
950 return raw_vector<T>(0, 0, 1);
952 return raw_vector<T>(&((*this)(iRow, iColumnStart)), iColumnEnd-iColumnStart+1, 1);
957 if( (iRowStart>iRowEnd) || wrongColumn(iColumn) || wrongRow(iRowStart) ||wrongRow(iRowEnd) )
960 return const_raw_vector<T>(&((*this)(iRowStart, iColumn)), iRowEnd-iRowStart+1, m_iLinearMember);
963 const_raw_vector<T>
getrow(
int iRow,
int iColumnStart,
int iColumnEnd)
const 965 if( (iColumnStart>iColumnEnd) || wrongRow(iRow) || wrongColumn(iColumnStart) || wrongColumn(iColumnEnd))
966 return const_raw_vector<T>(0, 0, 1);
968 return const_raw_vector<T>(&((*this)(iRow, iColumnStart)), iColumnEnd-iColumnStart+1, 1);
971 bool wrongRow(
int i)
const {
return i<m_iLow1 || i>m_iHigh1; };
976 long m_iLow1, m_iLow2, m_iHigh1, m_iHigh2;
977 long m_iConstOffset, m_iLinearMember;
1006 double sqr(
double x);
1007 int maxint(
int m1,
int m2);
1008 int minint(
int m1,
int m2);
1009 double maxreal(
double m1,
double m2);
1010 double minreal(
double m1,
double m2);
1020 #include <stdexcept> 1034 class incorrectPrecision :
public exception {};
1035 class overflow :
public exception {};
1036 class divisionByZero :
public exception {};
1037 class sqrtOfNegativeNumber :
public exception {};
1038 class invalidConversion :
public exception {};
1039 class invalidString :
public exception {};
1040 class internalError :
public exception {};
1041 class domainError :
public exception {};
1048 unsigned int refCount;
1049 unsigned int Precision;
1062 static mpfr_record* newMpfr(
unsigned int Precision);
1063 static void deleteMpfr(mpfr_record* ref);
1065 static gmp_randstate_t* getRandState();
1067 static mpfr_record_ptr&
getList(
unsigned int Precision);
1073 class mpfr_reference
1077 mpfr_reference(
const mpfr_reference& r);
1078 mpfr_reference& operator= (
const mpfr_reference &r);
1084 mpfr_srcptr getReadPtr()
const;
1085 mpfr_ptr getWritePtr();
1093 template<
unsigned int Precision>
1103 if( rval->refCount==0 )
1104 mpfr_storage::deleteMpfr(rval);
1113 ampf (
long double v) { InitializeAsDouble(v); }
1114 ampf (
double v) { InitializeAsDouble(v); }
1115 ampf (
float v) { InitializeAsDouble(v); }
1116 ampf (
signed long v) { InitializeAsSLong(v); }
1117 ampf (
unsigned long v) { InitializeAsULong(v); }
1118 ampf (
signed int v) { InitializeAsSLong(v); }
1119 ampf (
unsigned int v) { InitializeAsULong(v); }
1120 ampf (
signed short v) { InitializeAsSLong(v); }
1121 ampf (
unsigned short v) { InitializeAsULong(v); }
1122 ampf (
signed char v) { InitializeAsSLong(v); }
1123 ampf (
unsigned char v) { InitializeAsULong(v); }
1130 ampf (
const char *
s) { InitializeAsString(s); }
1140 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS 1141 template<
unsigned int Precision2>
1145 rval = mpfr_storage::newMpfr(Precision);
1146 mpfr_set(getWritePtr(), r.getReadPtr(), GMP_RNDN);
1153 ampf& operator= (
long double v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
return *
this; }
1154 ampf& operator= (
double v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
return *
this; }
1155 ampf& operator= (
float v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
return *
this; }
1156 ampf& operator= (
signed long v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN);
return *
this; }
1157 ampf& operator= (
unsigned long v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
return *
this; }
1158 ampf& operator= (
signed int v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN);
return *
this; }
1159 ampf& operator= (
unsigned int v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
return *
this; }
1160 ampf& operator= (
signed short v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN);
return *
this; }
1161 ampf& operator= (
unsigned short v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
return *
this; }
1162 ampf& operator= (
signed char v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN);
return *
this; }
1163 ampf& operator= (
unsigned char v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
return *
this; }
1164 ampf& operator= (
const char *
s) { mpfr_strtofr(getWritePtr(), s,
NULL, 0, GMP_RNDN);
return *
this; }
1165 ampf& operator= (
const std::string &
s) { mpfr_strtofr(getWritePtr(), s.c_str(),
NULL, 0, GMP_RNDN);
return *
this; }
1174 if( rval->refCount==0 )
1175 mpfr_storage::deleteMpfr(rval);
1181 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS 1182 template<
unsigned int Precision2>
1185 if( (
void*)
this==(
void*)(&r) )
1187 mpfr_set(getWritePtr(), r.getReadPtr(), GMP_RNDN);
1204 mpfr_srcptr getReadPtr()
const;
1205 mpfr_ptr getWritePtr();
1210 bool isFiniteNumber()
const;
1211 bool isPositiveNumber()
const;
1213 bool isNegativeNumber()
const;
1214 const ampf getUlpOf();
1219 double toDouble()
const;
1228 static const ampf getUlpOf(
const ampf &
x);
1229 static const ampf getUlp();
1230 static const ampf getUlp256();
1231 static const ampf getUlp512();
1232 static const ampf getMaxNumber();
1233 static const ampf getMinNumber();
1234 static const ampf getAlgoPascalEpsilon();
1235 static const ampf getAlgoPascalMaxNumber();
1236 static const ampf getAlgoPascalMinNumber();
1237 static const ampf getRandom();
1239 void CheckPrecision();
1240 void InitializeAsZero();
1241 void InitializeAsSLong(
signed long v);
1242 void InitializeAsULong(
unsigned long v);
1243 void InitializeAsDouble(
long double v);
1244 void InitializeAsString(
const char *
s);
1256 template<
unsigned int Precision>
1261 WerrorS(
"incorrectPrecision");
1264 template<
unsigned int Precision>
1268 rval = mpfr_storage::newMpfr(Precision);
1269 mpfr_set_ui(getWritePtr(), 0, GMP_RNDN);
1272 template<
unsigned int Precision>
1276 rval = mpfr_storage::newMpfr(Precision);
1277 mpfr_set_si(getWritePtr(), sv, GMP_RNDN);
1280 template<
unsigned int Precision>
1284 rval = mpfr_storage::newMpfr(Precision);
1285 mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
1288 template<
unsigned int Precision>
1292 rval = mpfr_storage::newMpfr(Precision);
1293 mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
1296 template<
unsigned int Precision>
1300 rval = mpfr_storage::newMpfr(Precision);
1301 mpfr_strtofr(getWritePtr(), s,
NULL, 0, GMP_RNDN);
1304 template<
unsigned int Precision>
1310 template<
unsigned int Precision>
1313 if( rval->refCount==1 )
1315 mpfr_record *newrval = mpfr_storage::newMpfr(Precision);
1316 mpfr_set(newrval->value, rval->value, GMP_RNDN);
1322 template<
unsigned int Precision>
1325 return mpfr_number_p(getReadPtr())!=0;
1328 template<
unsigned int Precision>
1331 if( !isFiniteNumber() )
1333 return mpfr_sgn(getReadPtr())>0;
1336 template<
unsigned int Precision>
1339 return mpfr_zero_p(getReadPtr())!=0;
1342 template<
unsigned int Precision>
1345 if( !isFiniteNumber() )
1347 return mpfr_sgn(getReadPtr())<0;
1350 template<
unsigned int Precision>
1353 return getUlpOf(*
this);
1356 template<
unsigned int Precision>
1359 return mpfr_get_d(getReadPtr(), GMP_RNDN);
1362 template<
unsigned int Precision>
1368 if( !isFiniteNumber() )
1373 ptr = mpfr_get_str(
NULL, &_e, 16, 0, getReadPtr(), GMP_RNDN);
1384 signed long iexpval;
1388 ptr = mpfr_get_str(
NULL, &expval, 16, 0, getReadPtr(), GMP_RNDN);
1391 if( iexpval!=expval )
1394 sprintf(buf_e,
"%ld",
long(iexpval));
1404 mpfr_free_str(ptr2);
1408 template<
unsigned int Precision>
1416 if( !isFiniteNumber() )
1421 ptr = mpfr_get_str(
NULL, &_e, 10, 0, getReadPtr(), GMP_RNDN);
1432 signed long iexpval;
1436 ptr = mpfr_get_str(
NULL, &expval, 10, 0, getReadPtr(), GMP_RNDN);
1439 if( iexpval!=expval )
1442 sprintf(buf_e,
"%ld",
long(iexpval));
1452 mpfr_free_str(ptr2);
1455 template<
unsigned int Precision>
1458 char *toString_Block=(
char *)
omAlloc(256);
1462 if( !isFiniteNumber() )
1466 ptr = mpfr_get_str(
NULL, &_e, 10, 0, getReadPtr(), GMP_RNDN);
1467 strcpy(toString_Block, ptr);
1469 return toString_Block;
1477 signed long iexpval;
1481 ptr = mpfr_get_str(
NULL, &expval, 10, 0, getReadPtr(), GMP_RNDN);
1484 if( iexpval!=expval )
1487 sprintf(buf_e,
"%ld",
long(iexpval));
1491 sprintf(toString_Block,
"-0.%sE%s",ptr,buf_e);
1494 sprintf(toString_Block,
"0.%sE%s",ptr,buf_e);
1495 mpfr_free_str(ptr2);
1496 return toString_Block;
1499 template<
unsigned int Precision>
1502 if( !x.isFiniteNumber() )
1507 mpfr_nextabove(r.getWritePtr());
1508 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1512 mpfr_get_exp(x.getReadPtr()),
1522 template<
unsigned int Precision>
1526 mpfr_nextabove(r.getWritePtr());
1527 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1531 template<
unsigned int Precision>
1535 mpfr_nextabove(r.getWritePtr());
1536 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1545 template<
unsigned int Precision>
1549 mpfr_nextabove(r.getWritePtr());
1550 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1559 template<
unsigned int Precision>
1563 mpfr_nextbelow(r.getWritePtr());
1564 mpfr_set_exp(r.getWritePtr(),mpfr_get_emax());
1568 template<
unsigned int Precision>
1572 mpfr_set_exp(r.getWritePtr(),mpfr_get_emin());
1576 template<
unsigned int Precision>
1582 template<
unsigned int Precision>
1586 mp_exp_t e1 = mpfr_get_emax();
1587 mp_exp_t e2 = -mpfr_get_emin();
1588 mp_exp_t e = e1>e2 ? e1 : e2;
1589 mpfr_set_exp(r.getWritePtr(), e-5);
1593 template<
unsigned int Precision>
1597 mp_exp_t e1 = mpfr_get_emax();
1598 mp_exp_t e2 = -mpfr_get_emin();
1599 mp_exp_t e = e1>e2 ? e1 : e2;
1600 mpfr_set_exp(r.getWritePtr(), 2-(e-5));
1604 template<
unsigned int Precision>
1615 template<
unsigned int Precision>
1618 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())==0;
1621 template<
unsigned int Precision>
1624 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())!=0;
1627 template<
unsigned int Precision>
1630 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())<0;
1633 template<
unsigned int Precision>
1636 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())>0;
1639 template<
unsigned int Precision>
1642 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())<=0;
1645 template<
unsigned int Precision>
1648 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())>=0;
1654 template<
unsigned int Precision>
1660 template<
unsigned int Precision>
1663 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1664 mpfr_neg(v->value, op1.getReadPtr(), GMP_RNDN);
1668 template<
unsigned int Precision>
1671 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1672 mpfr_add(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1676 template<
unsigned int Precision>
1679 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1680 mpfr_sub(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1685 template<
unsigned int Precision>
1688 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1689 mpfr_mul(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1693 template<
unsigned int Precision>
1696 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1697 mpfr_div(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1704 template<
unsigned int Precision>
1709 mpfr_sqr(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
1713 template<
unsigned int Precision>
1716 int s = mpfr_sgn(x.getReadPtr());
1724 template<
unsigned int Precision>
1729 mpfr_abs(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
1733 template<
unsigned int Precision>
1738 mpfr_max(res.getWritePtr(), x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
1742 template<
unsigned int Precision>
1747 mpfr_min(res.getWritePtr(), x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
1751 template<
unsigned int Precision>
1756 mpfr_sqrt(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
1760 template<
unsigned int Precision>
1765 mpfr_trunc(tmp.getWritePtr(), x.getReadPtr());
1766 if( mpfr_integer_p(tmp.getReadPtr())==0 )
1769 mpfr_clear_erangeflag();
1770 r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
1771 if( mpfr_erangeflag_p()!=0 )
1777 template<
unsigned int Precision>
1782 mpfr_frac(r.getWritePtr(), x.getReadPtr(), GMP_RNDN);
1786 template<
unsigned int Precision>
1791 mpfr_floor(tmp.getWritePtr(), x.getReadPtr());
1792 if( mpfr_integer_p(tmp.getReadPtr())==0 )
1795 mpfr_clear_erangeflag();
1796 r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
1797 if( mpfr_erangeflag_p()!=0 )
1803 template<
unsigned int Precision>
1808 mpfr_ceil(tmp.getWritePtr(), x.getReadPtr());
1809 if( mpfr_integer_p(tmp.getReadPtr())==0 )
1812 mpfr_clear_erangeflag();
1813 r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
1814 if( mpfr_erangeflag_p()!=0 )
1820 template<
unsigned int Precision>
1825 mpfr_round(tmp.getWritePtr(), x.getReadPtr());
1826 if( mpfr_integer_p(tmp.getReadPtr())==0 )
1829 mpfr_clear_erangeflag();
1830 r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
1831 if( mpfr_erangeflag_p()!=0 )
1837 template<
unsigned int Precision>
1842 if( !x.isFiniteNumber() )
1852 *exponent = mpfr_get_exp(r.getReadPtr());
1853 mpfr_set_exp(r.getWritePtr(),0);
1857 template<
unsigned int Precision>
1862 mpfr_mul_2si(r.getWritePtr(), x.getReadPtr(),
exponent, GMP_RNDN);
1869 #define __AMP_BINARY_OPI(type) \ 1870 template<unsigned int Precision> const ampf<Precision> operator+(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \ 1871 template<unsigned int Precision> const ampf<Precision> operator+(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \ 1872 template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const signed type& op2) { return op1+ampf<Precision>(op2); } \ 1873 template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const unsigned type& op2) { return op1+ampf<Precision>(op2); } \ 1874 template<unsigned int Precision> const ampf<Precision> operator-(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \ 1875 template<unsigned int Precision> const ampf<Precision> operator-(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \ 1876 template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const signed type& op2) { return op1-ampf<Precision>(op2); } \ 1877 template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const unsigned type& op2) { return op1-ampf<Precision>(op2); } \ 1878 template<unsigned int Precision> const ampf<Precision> operator*(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \ 1879 template<unsigned int Precision> const ampf<Precision> operator*(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \ 1880 template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const signed type& op2) { return op1*ampf<Precision>(op2); } \ 1881 template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const unsigned type& op2) { return op1*ampf<Precision>(op2); } \ 1882 template<unsigned int Precision> const ampf<Precision> operator/(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \ 1883 template<unsigned int Precision> const ampf<Precision> operator/(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \ 1884 template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const signed type& op2) { return op1/ampf<Precision>(op2); } \ 1885 template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const unsigned type& op2) { return op1/ampf<Precision>(op2); } \ 1886 template<unsigned int Precision> bool operator==(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \ 1887 template<unsigned int Precision> bool operator==(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \ 1888 template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const signed type& op2) { return op1==ampf<Precision>(op2); } \ 1889 template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const unsigned type& op2) { return op1==ampf<Precision>(op2); } \ 1890 template<unsigned int Precision> bool operator!=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \ 1891 template<unsigned int Precision> bool operator!=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \ 1892 template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const signed type& op2) { return op1!=ampf<Precision>(op2); } \ 1893 template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const unsigned type& op2) { return op1!=ampf<Precision>(op2); } \ 1894 template<unsigned int Precision> bool operator<=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \ 1895 template<unsigned int Precision> bool operator<=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \ 1896 template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const signed type& op2) { return op1<=ampf<Precision>(op2); } \ 1897 template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const unsigned type& op2) { return op1<=ampf<Precision>(op2); } \ 1898 template<unsigned int Precision> bool operator>=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \ 1899 template<unsigned int Precision> bool operator>=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \ 1900 template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const signed type& op2) { return op1>=ampf<Precision>(op2); } \ 1901 template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const unsigned type& op2) { return op1>=ampf<Precision>(op2); } \ 1902 template<unsigned int Precision> bool operator<(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \ 1903 template<unsigned int Precision> bool operator<(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \ 1904 template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const signed type& op2) { return op1<ampf<Precision>(op2); } \ 1905 template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const unsigned type& op2) { return op1<ampf<Precision>(op2); } \ 1906 template<unsigned int Precision> bool operator>(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \ 1907 template<unsigned int Precision> bool operator>(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \ 1908 template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const signed type& op2) { return op1>ampf<Precision>(op2); } \ 1909 template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const unsigned type& op2) { return op1>ampf<Precision>(op2); } 1914 #undef __AMP_BINARY_OPI 1915 #define __AMP_BINARY_OPF(type) \ 1916 template<unsigned int Precision> const ampf<Precision> operator+(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \ 1917 template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const type& op2) { return op1+ampf<Precision>(op2); } \ 1918 template<unsigned int Precision> const ampf<Precision> operator-(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \ 1919 template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const type& op2) { return op1-ampf<Precision>(op2); } \ 1920 template<unsigned int Precision> const ampf<Precision> operator*(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \ 1921 template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const type& op2) { return op1*ampf<Precision>(op2); } \ 1922 template<unsigned int Precision> const ampf<Precision> operator/(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \ 1923 template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const type& op2) { return op1/ampf<Precision>(op2); } \ 1924 template<unsigned int Precision> bool operator==(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \ 1925 template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const type& op2) { return op1==ampf<Precision>(op2); } \ 1926 template<unsigned int Precision> bool operator!=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \ 1927 template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const type& op2) { return op1!=ampf<Precision>(op2); } \ 1928 template<unsigned int Precision> bool operator<=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \ 1929 template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const type& op2) { return op1<=ampf<Precision>(op2); } \ 1930 template<unsigned int Precision> bool operator>=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \ 1931 template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const type& op2) { return op1>=ampf<Precision>(op2); } \ 1932 template<unsigned int Precision> bool operator<(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \ 1933 template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const type& op2) { return op1<ampf<Precision>(op2); } \ 1934 template<unsigned int Precision> bool operator>(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \ 1935 template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const type& op2) { return op1>ampf<Precision>(op2); } 1939 #undef __AMP_BINARY_OPF 1944 template<
unsigned int Precision>
1947 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1948 mpfr_const_pi(v->value, GMP_RNDN);
1952 template<
unsigned int Precision>
1955 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1956 mpfr_const_pi(v->value, GMP_RNDN);
1957 mpfr_mul_2si(v->value, v->value, -1, GMP_RNDN);
1961 template<
unsigned int Precision>
1964 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1965 mpfr_const_pi(v->value, GMP_RNDN);
1966 mpfr_mul_2si(v->value, v->value, +1, GMP_RNDN);
1970 template<
unsigned int Precision>
1973 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1974 mpfr_sin(v->value, x.getReadPtr(), GMP_RNDN);
1978 template<
unsigned int Precision>
1981 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1982 mpfr_cos(v->value, x.getReadPtr(), GMP_RNDN);
1986 template<
unsigned int Precision>
1989 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1990 mpfr_tan(v->value, x.getReadPtr(), GMP_RNDN);
1994 template<
unsigned int Precision>
1997 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1998 mpfr_asin(v->value, x.getReadPtr(), GMP_RNDN);
2002 template<
unsigned int Precision>
2005 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2006 mpfr_acos(v->value, x.getReadPtr(), GMP_RNDN);
2010 template<
unsigned int Precision>
2013 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2014 mpfr_atan(v->value, x.getReadPtr(), GMP_RNDN);
2018 template<
unsigned int Precision>
2021 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2022 mpfr_atan2(v->value, y.getReadPtr(), x.getReadPtr(), GMP_RNDN);
2026 template<
unsigned int Precision>
2029 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2030 mpfr_log(v->value, x.getReadPtr(), GMP_RNDN);
2034 template<
unsigned int Precision>
2037 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2038 mpfr_log2(v->value, x.getReadPtr(), GMP_RNDN);
2042 template<
unsigned int Precision>
2045 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2046 mpfr_log10(v->value, x.getReadPtr(), GMP_RNDN);
2050 template<
unsigned int Precision>
2053 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2054 mpfr_exp(v->value, x.getReadPtr(), GMP_RNDN);
2058 template<
unsigned int Precision>
2061 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2062 mpfr_sinh(v->value, x.getReadPtr(), GMP_RNDN);
2066 template<
unsigned int Precision>
2069 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2070 mpfr_cosh(v->value, x.getReadPtr(), GMP_RNDN);
2074 template<
unsigned int Precision>
2077 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2078 mpfr_tanh(v->value, x.getReadPtr(), GMP_RNDN);
2082 template<
unsigned int Precision>
2085 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2086 mpfr_pow(v->value, x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
2093 template<
unsigned int Precision>
2112 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS 2113 template<
unsigned int Prec2>
2117 campf& operator= (
long double v) { x=
v; y=0;
return *
this; }
2118 campf& operator= (
double v) { x=
v; y=0;
return *
this; }
2119 campf& operator= (
float v) { x=
v; y=0;
return *
this; }
2120 campf& operator= (
signed long v) { x=
v; y=0;
return *
this; }
2121 campf& operator= (
unsigned long v) { x=
v; y=0;
return *
this; }
2122 campf& operator= (
signed int v) { x=
v; y=0;
return *
this; }
2123 campf& operator= (
unsigned int v) { x=
v; y=0;
return *
this; }
2124 campf& operator= (
signed short v) { x=
v; y=0;
return *
this; }
2125 campf& operator= (
unsigned short v) { x=
v; y=0;
return *
this; }
2126 campf& operator= (
signed char v) { x=
v; y=0;
return *
this; }
2127 campf& operator= (
unsigned char v) { x=
v; y=0;
return *
this; }
2128 campf& operator= (
const char *s) { x=
s; y=0;
return *
this; }
2136 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS 2137 template<
unsigned int Precision2>
2152 template<
unsigned int Precision>
2154 {
return lhs.
x==rhs.
x && lhs.
y==rhs.
y; }
2156 template<
unsigned int Precision>
2158 {
return lhs.
x!=rhs.
x || lhs.
y!=rhs.
y; }
2160 template<
unsigned int Precision>
2164 template<
unsigned int Precision>
2166 { lhs.
x += rhs.
x; lhs.
y += rhs.
y;
return lhs; }
2168 template<
unsigned int Precision>
2172 template<
unsigned int Precision>
2176 template<
unsigned int Precision>
2178 { lhs.
x -= rhs.
x; lhs.
y -= rhs.
y;
return lhs; }
2180 template<
unsigned int Precision>
2184 template<
unsigned int Precision>
2193 template<
unsigned int Precision>
2197 template<
unsigned int Precision>
2207 result.
x = (lhs.
x+lhs.
y*e)/f;
2208 result.
y = (lhs.
y-lhs.
x*e)/f;
2214 result.
x = (lhs.
y+lhs.
x*e)/f;
2215 result.
y = (-lhs.
x+lhs.
y*e)/f;
2220 template<
unsigned int Precision>
2227 template<
unsigned int Precision>
2234 w = xabs>yabs ? xabs : yabs;
2235 v = xabs<yabs ? xabs : yabs;
2245 template<
unsigned int Precision>
2251 template<
unsigned int Precision>
2261 #define __AMP_BINARY_OPI(type) \ 2262 template<unsigned int Precision> const campf<Precision> operator+ (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \ 2263 template<unsigned int Precision> const campf<Precision> operator+ (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \ 2264 template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \ 2265 template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \ 2266 template<unsigned int Precision> const campf<Precision> operator- (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \ 2267 template<unsigned int Precision> const campf<Precision> operator- (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \ 2268 template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \ 2269 template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \ 2270 template<unsigned int Precision> const campf<Precision> operator* (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \ 2271 template<unsigned int Precision> const campf<Precision> operator* (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \ 2272 template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \ 2273 template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \ 2274 template<unsigned int Precision> const campf<Precision> operator/ (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(ampf<Precision>(op1),ampf<Precision>(0))/op2; } \ 2275 template<unsigned int Precision> const campf<Precision> operator/ (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(ampf<Precision>(op1),ampf<Precision>(0))/op2; } \ 2276 template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \ 2277 template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \ 2278 template<unsigned int Precision> bool operator==(const signed type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \ 2279 template<unsigned int Precision> bool operator==(const unsigned type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \ 2280 template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const signed type& op2) { return op1.x==op2 && op1.y==0; } \ 2281 template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const unsigned type& op2) { return op1.x==op2 && op1.y==0; } \ 2282 template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const signed type& op2) { return op1.x!=op2 || op1.y!=0; } \ 2283 template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const unsigned type& op2) { return op1.x!=op2 || op1.y!=0; } \ 2284 template<unsigned int Precision> bool operator!=(const signed type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } \ 2285 template<unsigned int Precision> bool operator!=(const unsigned type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } 2290 #undef __AMP_BINARY_OPI 2291 #define __AMP_BINARY_OPF(type) \ 2292 template<unsigned int Precision> const campf<Precision> operator+ (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \ 2293 template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \ 2294 template<unsigned int Precision> const campf<Precision> operator- (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \ 2295 template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \ 2296 template<unsigned int Precision> const campf<Precision> operator* (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \ 2297 template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \ 2298 template<unsigned int Precision> const campf<Precision> operator/ (const type& op1, const campf<Precision>& op2) { return campf<Precision>(ampf<Precision>(op1),ampf<Precision>(0))/op2; } \ 2299 template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \ 2300 template<unsigned int Precision> bool operator==(const type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \ 2301 template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const type& op2) { return op1.x==op2 && op1.y==0; } \ 2302 template<unsigned int Precision> bool operator!=(const type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } \ 2303 template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const type& op2) { return op1.x!=op2 || op1.y!=0; } 2308 #undef __AMP_BINARY_OPF 2313 template<
unsigned int Precision>
2317 int i, cnt = v1.GetLength();
2318 const ampf<Precision> *p1 = v1.GetData();
2319 const ampf<Precision> *p2 = v2.GetData();
2320 mpfr_record *r =
NULL;
2321 mpfr_record *t =
NULL;
2324 r = mpfr_storage::newMpfr(Precision);
2325 t = mpfr_storage::newMpfr(Precision);
2326 mpfr_set_ui(r->value, 0, GMP_RNDN);
2327 for(i=0; i<cnt; i++)
2329 mpfr_mul(t->value, p1->getReadPtr(), p2->getReadPtr(), GMP_RNDN);
2330 mpfr_add(r->value, r->value, t->value, GMP_RNDN);
2334 mpfr_storage::deleteMpfr(t);
2347 template<
unsigned int Precision>
2351 int i, cnt = vDst.GetLength();
2352 ampf<Precision> *pDst = vDst.GetData();
2353 const ampf<Precision> *pSrc = vSrc.GetData();
2356 for(i=0; i<cnt; i++)
2359 pDst += vDst.GetStep();
2360 pSrc += vSrc.GetStep();
2364 template<
unsigned int Precision>
2368 int i, cnt = vDst.GetLength();
2369 ampf<Precision> *pDst = vDst.GetData();
2370 const ampf<Precision> *pSrc = vSrc.GetData();
2371 for(i=0; i<cnt; i++)
2374 mpfr_ptr v = pDst->getWritePtr();
2375 mpfr_neg(v, v, GMP_RNDN);
2376 pDst += vDst.GetStep();
2377 pSrc += vSrc.GetStep();
2381 template<
unsigned int Precision,
class T2>
2385 int i, cnt = vDst.GetLength();
2386 ampf<Precision> *pDst = vDst.GetData();
2387 const ampf<Precision> *pSrc = vSrc.GetData();
2388 ampf<Precision> a(alpha);
2389 for(i=0; i<cnt; i++)
2392 mpfr_ptr v = pDst->getWritePtr();
2393 mpfr_mul(v, v, a.getReadPtr(), GMP_RNDN);
2394 pDst += vDst.GetStep();
2395 pSrc += vSrc.GetStep();
2399 template<
unsigned int Precision>
2403 int i, cnt = vDst.GetLength();
2404 ampf<Precision> *pDst = vDst.GetData();
2405 const ampf<Precision> *pSrc = vSrc.GetData();
2406 for(i=0; i<cnt; i++)
2408 mpfr_ptr v = pDst->getWritePtr();
2409 mpfr_srcptr vs = pSrc->getReadPtr();
2410 mpfr_add(v, v, vs, GMP_RNDN);
2411 pDst += vDst.GetStep();
2412 pSrc += vSrc.GetStep();
2416 template<
unsigned int Precision,
class T2>
2420 int i, cnt = vDst.GetLength();
2421 ampf<Precision> *pDst = vDst.GetData();
2422 const ampf<Precision> *pSrc = vSrc.GetData();
2423 ampf<Precision> a(alpha), tmp;
2424 for(i=0; i<cnt; i++)
2426 mpfr_ptr v = pDst->getWritePtr();
2427 mpfr_srcptr vs = pSrc->getReadPtr();
2428 mpfr_mul(tmp.getWritePtr(), a.getReadPtr(), vs, GMP_RNDN);
2429 mpfr_add(v, v, tmp.getWritePtr(), GMP_RNDN);
2430 pDst += vDst.GetStep();
2431 pSrc += vSrc.GetStep();
2435 template<
unsigned int Precision>
2439 int i, cnt = vDst.GetLength();
2440 ampf<Precision> *pDst = vDst.GetData();
2441 const ampf<Precision> *pSrc = vSrc.GetData();
2442 for(i=0; i<cnt; i++)
2444 mpfr_ptr v = pDst->getWritePtr();
2445 mpfr_srcptr vs = pSrc->getReadPtr();
2446 mpfr_sub(v, v, vs, GMP_RNDN);
2447 pDst += vDst.GetStep();
2448 pSrc += vSrc.GetStep();
2452 template<
unsigned int Precision,
class T2>
2455 vAdd(vDst, vSrc, -alpha);
2458 template<
unsigned int Precision,
class T2>
2461 int i, cnt = vDst.GetLength();
2462 ampf<Precision> *pDst = vDst.GetData();
2463 ampf<Precision> a(alpha);
2464 for(i=0; i<cnt; i++)
2466 mpfr_ptr v = pDst->getWritePtr();
2467 mpfr_mul(v, a.getReadPtr(),
v, GMP_RNDN);
2468 pDst += vDst.GetStep();
2515 template<
unsigned int Precision>
2519 template<
unsigned int Precision>
2528 template<
unsigned int Precision>
2579 template<
unsigned int Precision>
2609 mx = amp::maximum<Precision>(amp::abs<Precision>(
x(j)), mx);
2616 xnorm = xnorm+amp::sqr<Precision>(
x(j)/mx);
2618 xnorm = amp::sqrt<Precision>(xnorm)*mx;
2633 mx = amp::maximum<Precision>(amp::abs<Precision>(
alpha), amp::abs<Precision>(xnorm));
2634 beta = -mx*amp::sqrt<Precision>(amp::sqr<Precision>(alpha/mx)+amp::sqr<Precision>(xnorm/mx));
2639 tau = (beta-
alpha)/beta;
2674 template<
unsigned int Precision>
2687 if( tau==0 || n1>n2 || m1>m2 )
2695 for(i=n1; i<=n2; i++)
2699 for(i=m1; i<=m2; i++)
2702 ap::vadd(work.getvector(n1, n2), c.getrow(i, n1, n2), t);
2708 for(i=m1; i<=m2; i++)
2711 ap::vsub(c.getrow(i, n1, n2), work.getvector(n1, n2), t);
2744 template<
unsigned int Precision>
2759 if( tau==0 || n1>n2 || m1>m2 )
2768 for(i=m1; i<=m2; i++)
2777 for(i=m1; i<=m2; i++)
2780 ap::vsub(c.getrow(i, n1, n2), v.getvector(1, vm), t);
2827 template<
unsigned int Precision>
2833 template<
unsigned int Precision>
2840 template<
unsigned int Precision>
2850 template<
unsigned int Precision>
2857 template<
unsigned int Precision>
2867 template<
unsigned int Precision>
2874 template<
unsigned int Precision>
2880 template<
unsigned int Precision>
2887 template<
unsigned int Precision>
2897 template<
unsigned int Precision>
2904 template<
unsigned int Precision>
2914 template<
unsigned int Precision>
2973 template<
unsigned int Precision>
3003 tauq.setbounds(0, n-1);
3004 taup.setbounds(0, n-1);
3008 tauq.setbounds(0, m-1);
3009 taup.setbounds(0, m-1);
3017 for(i=0; i<=n-1; i++)
3024 reflections::generatereflection<Precision>(t, m-
i, ltau);
3026 ap::vmove(a.getcolumn(i, i, m-1), t.getvector(1, m-i));
3032 reflections::applyreflectionfromtheleft<Precision>(a, ltau, t,
i, m-1, i+1, n-1, work);
3040 ap::vmove(t.getvector(1, n-i-1), a.getrow(i, i+1, n-1));
3041 reflections::generatereflection<Precision>(t, n-1-
i, ltau);
3043 ap::vmove(a.getrow(i, i+1, n-1), t.getvector(1, n-1-i));
3049 reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1, m-1, i+1, n-1, work);
3063 for(i=0; i<=m-1; i++)
3070 reflections::generatereflection<Precision>(t, n-
i, ltau);
3072 ap::vmove(a.getrow(i, i, n-1), t.getvector(1, n-i));
3078 reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1, m-1,
i, n-1, work);
3086 ap::vmove(t.getvector(1, m-1-i), a.getcolumn(i, i+1, m-1));
3087 reflections::generatereflection<Precision>(t, m-1-
i, ltau);
3089 ap::vmove(a.getcolumn(i, i+1, m-1), t.getvector(1, m-1-i));
3095 reflections::applyreflectionfromtheleft<Precision>(a, ltau, t, i+1, m-1, i+1, n-1, work);
3127 template<
unsigned int Precision>
3141 if( m==0 || n==0 || qcolumns==0 )
3149 q.setbounds(0, m-1, 0, qcolumns-1);
3150 for(i=0; i<=m-1; i++)
3152 for(j=0; j<=qcolumns-1; j++)
3168 rmatrixbdmultiplybyq<Precision>(qp,
m, n, tauq, q,
m, qcolumns,
false,
false);
3201 template<
unsigned int Precision>
3221 if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
3271 reflections::applyreflectionfromtheright<Precision>(z, tauq(i),
v, 0, zrows-1,
i, m-1, work);
3275 reflections::applyreflectionfromtheleft<Precision>(z, tauq(i),
v,
i, m-1, 0, zcolumns-1, work);
3279 while( i!=i2+istep );
3319 reflections::applyreflectionfromtheright<Precision>(z, tauq(i),
v, 0, zrows-1, i+1, m-1, work);
3323 reflections::applyreflectionfromtheleft<Precision>(z, tauq(i),
v, i+1, m-1, 0, zcolumns-1, work);
3327 while( i!=i2+istep );
3354 template<
unsigned int Precision>
3368 if( m==0 || n==0 || ptrows==0 )
3376 pt.setbounds(0, ptrows-1, 0, n-1);
3377 for(i=0; i<=ptrows-1; i++)
3379 for(j=0; j<=n-1; j++)
3395 rmatrixbdmultiplybyp<Precision>(qp,
m, n, taup, pt, ptrows, n,
true,
true);
3428 template<
unsigned int Precision>
3448 if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
3502 reflections::applyreflectionfromtheright<Precision>(z, taup(i),
v, 0, zrows-1, i+1, n-1, work);
3506 reflections::applyreflectionfromtheleft<Precision>(z, taup(i),
v, i+1, n-1, 0, zcolumns-1, work);
3510 while( i!=i2+istep );
3549 reflections::applyreflectionfromtheright<Precision>(z, taup(i),
v, 0, zrows-1,
i, n-1, work);
3553 reflections::applyreflectionfromtheleft<Precision>(z, taup(i),
v,
i, n-1, 0, zcolumns-1, work);
3557 while( i!=i2+istep );
3584 template<
unsigned int Precision>
3602 d.setbounds(0, n-1);
3603 e.setbounds(0, n-1);
3604 for(i=0; i<=n-2; i++)
3609 d(n-1) =
b(n-1,n-1);
3613 d.setbounds(0, m-1);
3614 e.setbounds(0, m-1);
3615 for(i=0; i<=m-2; i++)
3620 d(m-1) =
b(m-1,m-1);
3629 template<
unsigned int Precision>
3653 taup.setbounds(1, minmn);
3654 tauq.setbounds(1, minmn);
3669 reflections::generatereflection<Precision>(t, mmip1, ltau);
3671 ap::vmove(a.getcolumn(i, i, m), t.getvector(1, mmip1));
3677 reflections::applyreflectionfromtheleft<Precision>(a, ltau, t,
i,
m, i+1, n, work);
3687 ap::vmove(t.getvector(1, nmi), a.getrow(i, ip1, n));
3688 reflections::generatereflection<Precision>(t, nmi, ltau);
3690 ap::vmove(a.getrow(i, ip1, n), t.getvector(1, nmi));
3696 reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1,
m, i+1, n, work);
3718 reflections::generatereflection<Precision>(t, nmip1, ltau);
3720 ap::vmove(a.getrow(i, i, n), t.getvector(1, nmip1));
3726 reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1,
m,
i, n, work);
3736 ap::vmove(t.getvector(1, mmi), a.getcolumn(i, ip1, m));
3737 reflections::generatereflection<Precision>(t, mmi, ltau);
3739 ap::vmove(a.getcolumn(i, ip1, m), t.getvector(1, mmi));
3745 reflections::applyreflectionfromtheleft<Precision>(a, ltau, t, i+1,
m, i+1, n, work);
3760 template<
unsigned int Precision>
3777 if( m==0 || n==0 || qcolumns==0 )
3785 q.setbounds(1, m, 1, qcolumns);
3794 for(j=1; j<=qcolumns; j++)
3813 reflections::applyreflectionfromtheleft<Precision>(q, tauq(i),
v,
i,
m, 1, qcolumns, work);
3818 for(i=
ap::minint(m-1, qcolumns-1); i>=1; i--)
3824 reflections::applyreflectionfromtheleft<Precision>(q, tauq(i),
v, i+1,
m, 1, qcolumns, work);
3834 template<
unsigned int Precision>
3856 if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
3907 reflections::applyreflectionfromtheright<Precision>(z, tauq(i),
v, 1, zrows,
i,
m, work);
3911 reflections::applyreflectionfromtheleft<Precision>(z, tauq(i),
v,
i,
m, 1, zcolumns, work);
3915 while( i!=i2+istep );
3957 reflections::applyreflectionfromtheright<Precision>(z, tauq(i),
v, 1, zrows, i+1,
m, work);
3961 reflections::applyreflectionfromtheleft<Precision>(z, tauq(i),
v, i+1,
m, 1, zcolumns, work);
3965 while( i!=i2+istep );
3975 template<
unsigned int Precision>
3992 if( m==0 || n==0 || ptrows==0 )
4000 pt.setbounds(1, ptrows, 1, n);
4007 for(i=1; i<=ptrows; i++)
4029 reflections::applyreflectionfromtheright<Precision>(pt, taup(i),
v, 1, ptrows, i+1, n, work);
4039 reflections::applyreflectionfromtheright<Precision>(pt, taup(i),
v, 1, ptrows,
i, n, work);
4049 template<
unsigned int Precision>
4071 if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
4127 reflections::applyreflectionfromtheright<Precision>(z, taup(i),
v, 1, zrows, i+1, n, work);
4131 reflections::applyreflectionfromtheleft<Precision>(z, taup(i),
v, i+1, n, 1, zcolumns, work);
4135 while( i!=i2+istep );
4175 reflections::applyreflectionfromtheright<Precision>(z, taup(i),
v, 1, zrows,
i, n, work);
4179 reflections::applyreflectionfromtheleft<Precision>(z, taup(i),
v,
i, n, 1, zcolumns, work);
4183 while( i!=i2+istep );
4192 template<
unsigned int Precision>
4212 for(i=1; i<=n-1; i++)
4223 for(i=1; i<=m-1; i++)
4275 template<
unsigned int Precision>
4280 template<
unsigned int Precision>
4287 template<
unsigned int Precision>
4292 template<
unsigned int Precision>
4297 template<
unsigned int Precision>
4304 template<
unsigned int Precision>
4350 template<
unsigned int Precision>
4371 tau.setbounds(0, minmn-1);
4377 for(i=0; i<=k-1; i++)
4384 reflections::generatereflection<Precision>(t, m-
i, tmp);
4386 ap::vmove(a.getcolumn(i, i, m-1), t.getvector(1, m-i));
4394 reflections::applyreflectionfromtheleft<Precision>(a,
tau(i), t,
i, m-1, i+1, n-1, work);
4420 template<
unsigned int Precision>
4437 if( m<=0 || n<=0 || qcolumns<=0 )
4447 q.setbounds(0, m-1, 0, qcolumns-1);
4450 for(i=0; i<=m-1; i++)
4452 for(j=0; j<=qcolumns-1; j++)
4468 for(i=k-1; i>=0; i--)
4476 reflections::applyreflectionfromtheleft<Precision>(q,
tau(i),
v,
i, m-1, 0, qcolumns-1, work);
4496 template<
unsigned int Precision>
4511 r.setbounds(0, m-1, 0, n-1);
4512 for(i=0; i<=n-1; i++)
4516 for(i=1; i<=m-1; i++)
4518 ap::vmove(r.getrow(i, 0, n-1), r.getrow(0, 0, n-1));
4520 for(i=0; i<=k-1; i++)
4522 ap::vmove(r.getrow(i, i, n-1), a.getrow(i, i, n-1));
4530 template<
unsigned int Precision>
4548 tau.setbounds(1, minmn);
4562 reflections::generatereflection<Precision>(t, mmip1, tmp);
4564 ap::vmove(a.getcolumn(i, i, m), t.getvector(1, mmip1));
4572 reflections::applyreflectionfromtheleft<Precision>(a,
tau(i), t,
i,
m, i+1, n, work);
4581 template<
unsigned int Precision>
4599 if( m==0 || n==0 || qcolumns==0 )
4609 q.setbounds(1, m, 1, qcolumns);
4614 for(j=1; j<=qcolumns; j++)
4639 reflections::applyreflectionfromtheleft<Precision>(q,
tau(i),
v,
i,
m, 1, qcolumns, work);
4647 template<
unsigned int Precision>
4668 q.setbounds(1, m, 1, m);
4669 r.setbounds(1, m, 1, n);
4674 qrdecomposition<Precision>(a,
m, n,
tau);
4685 ap::vmove(r.getrow(i, 1, n), r.getrow(1, 1, n));
4689 ap::vmove(r.getrow(i, i, n), a.getrow(i, i, n));
4695 unpackqfromqr<Precision>(a,
m, n,
tau,
m, q);
4735 template<
unsigned int Precision>
4740 template<
unsigned int Precision>
4747 template<
unsigned int Precision>
4752 template<
unsigned int Precision>
4757 template<
unsigned int Precision>
4764 template<
unsigned int Precision>
4806 template<
unsigned int Precision>
4825 tau.setbounds(0, minmn-1);
4827 for(i=0; i<=k-1; i++)
4834 reflections::generatereflection<Precision>(t, n-
i, tmp);
4836 ap::vmove(a.getrow(i, i, n-1), t.getvector(1, n-i));
4844 reflections::applyreflectionfromtheright<Precision>(a,
tau(i), t, i+1, m-1,
i, n-1, work);
4870 template<
unsigned int Precision>
4887 if( m<=0 || n<=0 || qrows<=0 )
4897 q.setbounds(0, qrows-1, 0, n-1);
4900 for(i=0; i<=qrows-1; i++)
4902 for(j=0; j<=n-1; j++)
4918 for(i=k-1; i>=0; i--)
4926 reflections::applyreflectionfromtheright<Precision>(q,
tau(i),
v, 0, qrows-1,
i, n-1, work);
4946 template<
unsigned int Precision>
4960 l.setbounds(0, m-1, 0, n-1);
4961 for(i=0; i<=n-1; i++)
4965 for(i=1; i<=m-1; i++)
4967 ap::vmove(
l.getrow(i, 0, n-1),
l.getrow(0, 0, n-1));
4969 for(i=0; i<=m-1; i++)
4972 ap::vmove(
l.getrow(i, 0, k), a.getrow(i, 0, k));
4981 template<
unsigned int Precision>
5001 tau.setbounds(1, minmn);
5015 reflections::generatereflection<Precision>(t, nmip1, tmp);
5017 ap::vmove(a.getrow(i, i, n), t.getvector(1, nmip1));
5025 reflections::applyreflectionfromtheright<Precision>(a,
tau(i), t, i+1,
m,
i, n, work);
5035 template<
unsigned int Precision>
5053 if( m==0 || n==0 || qrows==0 )
5063 q.setbounds(1, qrows, 1, n);
5066 for(i=1; i<=qrows; i++)
5093 reflections::applyreflectionfromtheright<Precision>(q,
tau(i),
v, 1, qrows,
i, n, work);
5101 template<
unsigned int Precision>
5117 q.setbounds(1, n, 1, n);
5118 l.setbounds(1, m, 1, n);
5123 lqdecomposition<Precision>(a,
m, n,
tau);
5146 unpackqfromlq<Precision>(a,
m, n,
tau, n, q);
5186 template<
unsigned int Precision>
5190 template<
unsigned int Precision>
5194 template<
unsigned int Precision>
5199 template<
unsigned int Precision>
5204 template<
unsigned int Precision>
5211 template<
unsigned int Precision>
5222 template<
unsigned int Precision>
5229 template<
unsigned int Precision>
5240 template<
unsigned int Precision>
5255 template<
unsigned int Precision>
5258 template<
unsigned int Precision>
5281 template<
unsigned int Precision>
5302 result = amp::abs<Precision>(
x(i1));
5307 for(ix=i1; ix<=i2; ix++)
5311 absxi = amp::abs<Precision>(
x(ix));
5314 ssq = 1+ssq*amp::sqr<Precision>(scl/absxi);
5319 ssq = ssq+amp::sqr<Precision>(absxi/scl);
5323 result = scl*amp::sqrt<Precision>(ssq);
5328 template<
unsigned int Precision>
5339 a = amp::abs<Precision>(
x(result));
5340 for(i=i1+1; i<=i2; i++)
5342 if( amp::abs<Precision>(
x(i))>amp::abs<Precision>(
x(result)) )
5351 template<
unsigned int Precision>
5363 a = amp::abs<Precision>(
x(result,j));
5364 for(i=i1+1; i<=i2; i++)
5366 if( amp::abs<Precision>(
x(i,j))>amp::abs<Precision>(
x(result,j)) )
5375 template<
unsigned int Precision>
5387 a = amp::abs<Precision>(
x(i,result));
5388 for(j=j1+1; j<=j2; j++)
5390 if( amp::abs<Precision>(
x(i,j))>amp::abs<Precision>(
x(i,result)) )
5399 template<
unsigned int Precision>
5413 for(j=j1; j<=j2; j++)
5417 for(i=i1; i<=i2; i++)
5421 work(j) = work(j)+amp::abs<Precision>(a(i,j));
5425 for(j=j1; j<=j2; j++)
5427 result = amp::maximum<Precision>(
result, work(j));
5433 template<
unsigned int Precision>
5449 if( is1>is2 || js1>js2 )
5455 for(isrc=is1; isrc<=is2; isrc++)
5457 idst = isrc-is1+id1;
5458 ap::vmove(
b.getrow(idst, jd1, jd2), a.getrow(isrc, js1, js2));
5463 template<
unsigned int Precision>
5478 if( i1>i2 || j1>j2 )
5483 for(i=i1; i<=i2-1; i++)
5489 ap::vmove(work.getvector(1, l), a.getcolumn(j, ips, i2));
5490 ap::vmove(a.getcolumn(j, ips, i2), a.getrow(i, jps, j2));
5491 ap::vmove(a.getrow(i, jps, j2), work.getvector(1, l));
5496 template<
unsigned int Precision>
5512 if( is1>is2 || js1>js2 )
5518 for(isrc=is1; isrc<=is2; isrc++)
5520 jdst = isrc-is1+jd1;
5521 ap::vmove(
b.getcolumn(jdst, id1, id2), a.getrow(isrc, js1, js2));
5526 template<
unsigned int Precision>
5552 if( i1>i2 || j1>j2 )
5564 for(i=iy1; i<=iy2; i++)
5577 for(i=i1; i<=i2; i++)
5580 y(iy1+i-i1) =
y(iy1+i-i1)+alpha*
v;
5589 if( i1>i2 || j1>j2 )
5601 for(i=iy1; i<=iy2; i++)
5614 for(i=i1; i<=i2; i++)
5616 v = alpha*
x(ix1+i-i1);
5617 ap::vadd(y.getvector(iy1, iy2), a.getrow(i, j1, j2),
v);
5623 template<
unsigned int Precision>
5634 xabs = amp::abs<Precision>(
x);
5635 yabs = amp::abs<Precision>(
y);
5636 w = amp::maximum<Precision>(xabs, yabs);
5637 z = amp::minimum<Precision>(xabs, yabs);
5644 result = w*amp::sqrt<Precision>(1+amp::sqr<Precision>(z/
w));
5650 template<
unsigned int Precision>
5711 if( arows<=0 || acols<=0 || brows<=0 || bcols<=0 )
5732 for(i=ci1; i<=ci2; i++)
5734 for(j=cj1; j<=cj2; j++)
5742 for(i=ci1; i<=ci2; i++)
5751 if( !transa && !transb )
5753 for(l=ai1; l<=ai2; l++)
5755 for(r=bi1; r<=bi2; r++)
5757 v = alpha*a(l,aj1+r-bi1);
5759 ap::vadd(c.getrow(k, cj1, cj2),
b.getrow(r, bj1, bj2),
v);
5768 if( !transa && transb )
5770 if( arows*acols<brows*bcols )
5772 for(r=bi1; r<=bi2; r++)
5774 for(l=ai1; l<=ai2; l++)
5777 c(ci1+l-ai1,cj1+r-bi1) = c(ci1+l-ai1,cj1+r-bi1)+alpha*
v;
5784 for(l=ai1; l<=ai2; l++)
5786 for(r=bi1; r<=bi2; r++)
5789 c(ci1+l-ai1,cj1+r-bi1) = c(ci1+l-ai1,cj1+r-bi1)+alpha*
v;
5799 if( transa && !transb )
5801 for(l=aj1; l<=aj2; l++)
5803 for(r=bi1; r<=bi2; r++)
5805 v = alpha*a(ai1+r-bi1,l);
5807 ap::vadd(c.getrow(k, cj1, cj2),
b.getrow(r, bj1, bj2),
v);
5816 if( transa && transb )
5818 if( arows*acols<brows*bcols )
5820 for(r=bi1; r<=bi2; r++)
5822 for(i=1; i<=crows; i++)
5826 for(l=ai1; l<=ai2; l++)
5828 v = alpha*
b(r,bj1+l-ai1);
5830 ap::vadd(work.getvector(1, crows), a.getrow(l, aj1, aj2),
v);
5832 ap::vadd(c.getcolumn(k, ci1, ci2), work.getvector(1, crows));
5838 for(l=aj1; l<=aj2; l++)
5841 ap::vmove(work.getvector(1, k), a.getcolumn(l, ai1, ai2));
5842 for(r=bi1; r<=bi2; r++)
5845 c(ci1+l-aj1,cj1+r-bi1) = c(ci1+l-aj1,cj1+r-bi1)+alpha*
v;
5896 template<
unsigned int Precision>
5906 template<
unsigned int Precision>
5916 template<
unsigned int Precision>
5949 template<
unsigned int Precision>
5967 if( m1>m2 || n1>n2 )
5983 for(j=m1; j<=m2-1; j++)
5987 if( ctemp!=1 || stemp!=0 )
5990 ap::vmove(work.getvector(n1, n2), a.getrow(jp1, n1, n2), ctemp);
5991 ap::vsub(work.getvector(n1, n2), a.getrow(j, n1, n2), stemp);
5992 ap::vmul(a.getrow(j, n1, n2), ctemp);
5993 ap::vadd(a.getrow(j, n1, n2), a.getrow(jp1, n1, n2), stemp);
5994 ap::vmove(a.getrow(jp1, n1, n2), work.getvector(n1, n2));
6004 for(j=m1; j<=m2-1; j++)
6008 if( ctemp!=1 || stemp!=0 )
6011 a(j+1,n1) = ctemp*temp-stemp*a(j,n1);
6012 a(j,n1) = stemp*temp+ctemp*a(j,n1);
6025 for(j=m2-1; j>=m1; j--)
6029 if( ctemp!=1 || stemp!=0 )
6032 ap::vmove(work.getvector(n1, n2), a.getrow(jp1, n1, n2), ctemp);
6033 ap::vsub(work.getvector(n1, n2), a.getrow(j, n1, n2), stemp);
6034 ap::vmul(a.getrow(j, n1, n2), ctemp);
6035 ap::vadd(a.getrow(j, n1, n2), a.getrow(jp1, n1, n2), stemp);
6036 ap::vmove(a.getrow(jp1, n1, n2), work.getvector(n1, n2));
6046 for(j=m2-1; j>=m1; j--)
6050 if( ctemp!=1 || stemp!=0 )
6053 a(j+1,n1) = ctemp*temp-stemp*a(j,n1);
6054 a(j,n1) = stemp*temp+ctemp*a(j,n1);
6087 template<
unsigned int Precision>
6117 for(j=n1; j<=n2-1; j++)
6121 if( ctemp!=1 || stemp!=0 )
6124 ap::vmove(work.getvector(m1, m2), a.getcolumn(jp1, m1, m2), ctemp);
6125 ap::vsub(work.getvector(m1, m2), a.getcolumn(j, m1, m2), stemp);
6126 ap::vmul(a.getcolumn(j, m1, m2), ctemp);
6127 ap::vadd(a.getcolumn(j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
6128 ap::vmove(a.getcolumn(jp1, m1, m2), work.getvector(m1, m2));
6138 for(j=n1; j<=n2-1; j++)
6142 if( ctemp!=1 || stemp!=0 )
6145 a(m1,j+1) = ctemp*temp-stemp*a(m1,j);
6146 a(m1,j) = stemp*temp+ctemp*a(m1,j);
6159 for(j=n2-1; j>=n1; j--)
6163 if( ctemp!=1 || stemp!=0 )
6166 ap::vmove(work.getvector(m1, m2), a.getcolumn(jp1, m1, m2), ctemp);
6167 ap::vsub(work.getvector(m1, m2), a.getcolumn(j, m1, m2), stemp);
6168 ap::vmul(a.getcolumn(j, m1, m2), ctemp);
6169 ap::vadd(a.getcolumn(j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
6170 ap::vmove(a.getcolumn(jp1, m1, m2), work.getvector(m1, m2));
6180 for(j=n2-1; j>=n1; j--)
6184 if( ctemp!=1 || stemp!=0 )
6187 a(m1,j+1) = ctemp*temp-stemp*a(m1,j);
6188 a(m1,j) = stemp*temp+ctemp*a(m1,j);
6204 template<
unsigned int Precision>
6233 r = amp::sqrt<Precision>(amp::sqr<Precision>(f1)+amp::sqr<Precision>(g1));
6236 if( amp::abs<Precision>(f)>amp::abs<Precision>(
g) && cs<0 )
6289 template<
unsigned int Precision>
6294 bool isfractionalaccuracyrequired,
6301 template<
unsigned int Precision>
6306 bool isfractionalaccuracyrequired,
6313 template<
unsigned int Precision>
6318 bool isfractionalaccuracyrequired,
6328 template<
unsigned int Precision>
6331 template<
unsigned int Precision>
6337 template<
unsigned int Precision>
6427 template<
unsigned int Precision>
6432 bool isfractionalaccuracyrequired,
6452 result = bidiagonalsvddecompositioninternal<Precision>(d1, e1, n, isupper, isfractionalaccuracyrequired, u, 0, nru, c, 0, ncc, vt, 0, ncvt);
6453 ap::vmove(d.getvector(0, n-1), d1.getvector(1, n));
6470 template<
unsigned int Precision>
6475 bool isfractionalaccuracyrequired,
6486 result = bidiagonalsvddecompositioninternal<Precision>(d, e, n, isupper, isfractionalaccuracyrequired, u, 1, nru, c, 1, ncc, vt, 1, ncvt);
6494 template<
unsigned int Precision>
6499 bool isfractionalaccuracyrequired,
6556 bool matrixsplitflag;
6584 ap::vmul(vt.getrow(vstart, vstart, vstart+ncvt-1), -1);
6610 for(i=1; i<=n-1; i++)
6615 for(i=1; i<=n-1; i++)
6634 for(i=1; i<=n-1; i++)
6636 rotations::generaterotation<Precision>(d(i), e(i), cs, sn, r);
6649 rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, 1+ustart-1, n+ustart-1, work0, work1, u, utemp);
6653 rotations::applyrotationsfromtheleft<Precision>(fwddir, 1+cstart-1, n+cstart-1, cstart, cend, work0, work1, c, ctemp);
6662 tolmul = amp::maximum<Precision>(10, amp::minimum<Precision>(100, amp::pow<Precision>(eps, -
amp::ampf<Precision>(
"0.125"))));
6664 if( !isfractionalaccuracyrequired )
6675 smax = amp::maximum<Precision>(smax, amp::abs<Precision>(d(i)));
6677 for(i=1; i<=n-1; i++)
6679 smax = amp::maximum<Precision>(smax, amp::abs<Precision>(e(i)));
6688 sminoa = amp::abs<Precision>(d(1));
6694 mu = amp::abs<Precision>(d(i))*(mu/(mu+amp::abs<Precision>(e(i-1))));
6695 sminoa = amp::minimum<Precision>(sminoa,
mu);
6702 sminoa = sminoa/amp::sqrt<Precision>(n);
6703 thresh = amp::maximum<Precision>(tol*sminoa, maxitr*n*n*unfl);
6711 thresh = amp::maximum<Precision>(amp::abs<Precision>(tol)*smax, maxitr*n*n*unfl);
6751 if( tol<0 && amp::abs<Precision>(d(m))<=thresh )
6755 smax = amp::abs<Precision>(d(m));
6757 matrixsplitflag =
false;
6758 for(lll=1; lll<=m-1; lll++)
6761 abss = amp::abs<Precision>(d(ll));
6762 abse = amp::abs<Precision>(e(ll));
6763 if( tol<0 && abss<=thresh )
6769 matrixsplitflag =
true;
6772 smin = amp::minimum<Precision>(smin, abss);
6773 smax = amp::maximum<Precision>(smax, amp::maximum<Precision>(abss, abse));
6775 if( !matrixsplitflag )
6807 svdv2x2<Precision>(d(m-1), e(m-1), d(m), sigmn, sigmx, sinr, cosr, sinl, cosl);
6818 mm1 = m-1+(vstart-1);
6821 ap::vmul(vt.getrow(mm0, vstart, vend), cosr);
6822 ap::vsub(vt.getrow(mm0, vstart, vend), vt.getrow(mm1, vstart, vend), sinr);
6831 ap::vmul(u.getcolumn(mm0, ustart, uend), cosl);
6832 ap::vsub(u.getcolumn(mm0, ustart, uend), u.getcolumn(mm1, ustart, uend), sinl);
6841 ap::vmul(c.getrow(mm0, cstart, cend), cosl);
6842 ap::vsub(c.getrow(mm0, cstart, cend), c.getrow(mm1, cstart, cend), sinl);
6859 if( idir==1 && amp::abs<Precision>(d(ll))<
amp::ampf<Precision>(
"1.0E-3")*amp::abs<Precision>(d(m)) )
6863 if( idir==2 && amp::abs<Precision>(d(m))<
amp::ampf<Precision>(
"1.0E-3")*amp::abs<Precision>(d(ll)) )
6867 if( ll!=oldll || m!=oldm || bchangedir )
6869 if( amp::abs<Precision>(d(ll))>=amp::abs<Precision>(d(m)) )
6897 if( amp::abs<Precision>(e(m-1))<=amp::abs<Precision>(tol)*amp::abs<Precision>(d(m)) || tol<0 && amp::abs<Precision>(e(m-1))<=thresh )
6909 mu = amp::abs<Precision>(d(ll));
6912 for(lll=ll; lll<=m-1; lll++)
6914 if( amp::abs<Precision>(e(lll))<=tol*
mu )
6921 mu = amp::abs<Precision>(d(lll+1))*(mu/(mu+amp::abs<Precision>(e(lll))));
6922 sminl = amp::minimum<Precision>(sminl,
mu);
6937 if( amp::abs<Precision>(e(ll))<=amp::abs<Precision>(tol)*amp::abs<Precision>(d(ll)) || tol<0 && amp::abs<Precision>(e(ll))<=thresh )
6949 mu = amp::abs<Precision>(d(m));
6952 for(lll=m-1; lll>=ll; lll--)
6954 if( amp::abs<Precision>(e(lll))<=tol*
mu )
6961 mu = amp::abs<Precision>(d(lll))*(mu/(mu+amp::abs<Precision>(e(lll))));
6962 sminl = amp::minimum<Precision>(sminl,
mu);
6977 if( tol>=0 && n*tol*(sminl/smax)<=amp::maximum<Precision>(eps,
amp::ampf<Precision>(
"0.01")*tol) )
6993 sll = amp::abs<Precision>(d(ll));
6994 svd2x2<Precision>(d(m-1), e(m-1), d(m), shift, r);
6998 sll = amp::abs<Precision>(d(m));
6999 svd2x2<Precision>(d(ll), e(ll), d(ll+1), shift, r);
7007 if( amp::sqr<Precision>(shift/sll)<eps )
7033 for(i=ll; i<=m-1; i++)
7035 rotations::generaterotation<Precision>(d(i)*cs, e(i), cs, sn, r);
7040 rotations::generaterotation<Precision>(oldcs*r, d(i+1)*sn, oldcs, oldsn, tmp);
7044 work2(i-ll+1) = oldcs;
7045 work3(i-ll+1) = oldsn;
7056 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work0, work1, vt, vttemp);
7060 rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work2, work3, u, utemp);
7064 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work2, work3, c, ctemp);
7070 if( amp::abs<Precision>(e(m-1))<=thresh )
7084 for(i=m; i>=ll+1; i--)
7086 rotations::generaterotation<Precision>(d(i)*cs, e(i-1), cs, sn, r);
7091 rotations::generaterotation<Precision>(oldcs*r, d(i-1)*sn, oldcs, oldsn, tmp);
7095 work2(i-ll) = oldcs;
7096 work3(i-ll) = -oldsn;
7107 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work2, work3, vt, vttemp);
7111 rotations::applyrotationsfromtheright<Precision>(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work0, work1, u, utemp);
7115 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work0, work1, c, ctemp);
7121 if( amp::abs<Precision>(e(ll))<=thresh )
7140 f = (amp::abs<Precision>(d(ll))-shift)*(extsignbdsqr<Precision>(1, d(ll))+shift/d(ll));
7142 for(i=ll; i<=m-1; i++)
7144 rotations::generaterotation<Precision>(
f,
g, cosr, sinr, r);
7149 f = cosr*d(i)+sinr*e(i);
7150 e(i) = cosr*e(i)-sinr*d(i);
7152 d(i+1) = cosr*d(i+1);
7153 rotations::generaterotation<Precision>(
f,
g, cosl, sinl, r);
7155 f = cosl*e(i)+sinl*d(i+1);
7156 d(i+1) = cosl*d(i+1)-sinl*e(i);
7160 e(i+1) = cosl*e(i+1);
7162 work0(i-ll+1) = cosr;
7163 work1(i-ll+1) = sinr;
7164 work2(i-ll+1) = cosl;
7165 work3(i-ll+1) = sinl;
7174 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work0, work1, vt, vttemp);
7178 rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work2, work3, u, utemp);
7182 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work2, work3, c, ctemp);
7188 if( amp::abs<Precision>(e(m-1))<=thresh )
7200 f = (amp::abs<Precision>(d(m))-shift)*(extsignbdsqr<Precision>(1, d(m))+shift/d(m));
7202 for(i=m; i>=ll+1; i--)
7204 rotations::generaterotation<Precision>(
f,
g, cosr, sinr, r);
7209 f = cosr*d(i)+sinr*e(i-1);
7210 e(i-1) = cosr*e(i-1)-sinr*d(i);
7212 d(i-1) = cosr*d(i-1);
7213 rotations::generaterotation<Precision>(
f,
g, cosl, sinl, r);
7215 f = cosl*e(i-1)+sinl*d(i-1);
7216 d(i-1) = cosl*d(i-1)-sinl*e(i-1);
7220 e(i-2) = cosl*e(i-2);
7223 work1(i-ll) = -sinr;
7225 work3(i-ll) = -sinl;
7232 if( amp::abs<Precision>(e(ll))<=thresh )
7242 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work2, work3, vt, vttemp);
7246 rotations::applyrotationsfromtheright<Precision>(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work0, work1, u, utemp);
7250 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work0, work1, c, ctemp);
7275 ap::vmul(vt.getrow(i+vstart-1, vstart, vend), -1);
7284 for(i=1; i<=n-1; i++)
7292 for(j=2; j<=n+1-
i; j++)
7312 ap::vmove(vt.getrow(isub+vstart-1, vstart, vend), vt.getrow(j+vstart-1, vstart, vend));
7319 ap::vmove(u.getcolumn(isub+ustart-1, ustart, uend), u.getcolumn(j+ustart-1, ustart, uend));
7326 ap::vmove(c.getrow(isub+cstart-1, cstart, cend), c.getrow(j+cstart-1, cstart, cend));
7335 template<
unsigned int Precision>
7344 result = amp::abs<Precision>(a);
7348 result = -amp::abs<Precision>(a);
7354 template<
unsigned int Precision>
7372 fa = amp::abs<Precision>(
f);
7373 ga = amp::abs<Precision>(
g);
7374 ha = amp::abs<Precision>(
h);
7375 fhmn = amp::minimum<Precision>(
fa, ha);
7376 fhmx = amp::maximum<Precision>(
fa, ha);
7386 ssmax = amp::maximum<Precision>(fhmx, ga)*amp::sqrt<Precision>(1+amp::sqr<Precision>(amp::minimum<Precision>(fhmx, ga)/amp::maximum<Precision>(fhmx, ga)));
7394 at = (fhmx-fhmn)/fhmx;
7395 au = amp::sqr<Precision>(ga/fhmx);
7396 c = 2/(amp::sqrt<Precision>(aas*aas+au)+amp::sqrt<Precision>(at*at+au));
7411 ssmin = fhmn*fhmx/ga;
7417 at = (fhmx-fhmn)/fhmx;
7418 c = 1/(amp::sqrt<Precision>(1+amp::sqr<Precision>(aas*au))+amp::sqrt<Precision>(1+amp::sqr<Precision>(at*au)));
7420 ssmin = ssmin+ssmin;
7428 template<
unsigned int Precision>
7467 fa = amp::abs<Precision>(ft);
7469 ha = amp::abs<Precision>(
h);
7494 ga = amp::abs<Precision>(gt);
7557 s = amp::sqrt<Precision>(tt+mm);
7560 r = amp::abs<Precision>(
m);
7564 r = amp::sqrt<Precision>(l*l+mm);
7577 t = extsignbdsqr<Precision>(2, ft)*extsignbdsqr<Precision>(1, gt);
7581 t = gt/extsignbdsqr<Precision>(d, ft)+m/t;
7586 t = (m/(s+t)+m/(r+l))*(1+a);
7588 l = amp::sqrt<Precision>(t*t+4);
7591 clt = (crt+srt*
m)/a;
7616 tsign = extsignbdsqr<Precision>(1, csr)*extsignbdsqr<Precision>(1, csl)*extsignbdsqr<Precision>(1,
f);
7620 tsign = extsignbdsqr<Precision>(1, snr)*extsignbdsqr<Precision>(1, csl)*extsignbdsqr<Precision>(1,
g);
7624 tsign = extsignbdsqr<Precision>(1, snr)*extsignbdsqr<Precision>(1, snl)*extsignbdsqr<Precision>(1,
h);
7626 ssmax = extsignbdsqr<Precision>(ssmax, tsign);
7627 ssmin = extsignbdsqr<Precision>(ssmin, tsign*extsignbdsqr<Precision>(1,
f)*extsignbdsqr<Precision>(1, h));
7685 template<
unsigned int Precision>
7691 int additionalmemory,
7695 template<
unsigned int Precision>
7701 int additionalmemory,
7757 template<
unsigned int Precision>
7763 int additionalmemory,
7800 w.setbounds(1, minmn);
7807 u.setbounds(0, nru-1, 0, ncu-1);
7813 u.setbounds(0, nru-1, 0, ncu-1);
7821 vt.setbounds(0, nrvt-1, 0, ncvt-1);
7827 vt.setbounds(0, nrvt-1, 0, ncvt-1);
7842 qr::rmatrixqr<Precision>(a,
m, n,
tau);
7843 for(i=0; i<=n-1; i++)
7845 for(j=0; j<=i-1; j++)
7850 bidiagonal::rmatrixbd<Precision>(a, n, n, tauq, taup);
7851 bidiagonal::rmatrixbdunpackpt<Precision>(a, n, n, taup, nrvt, vt);
7852 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, n, n, isupper,
w, e);
7853 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, n, isupper,
false, u, 0, a, 0, vt, ncvt);
7862 qr::rmatrixqr<Precision>(a,
m, n,
tau);
7863 qr::rmatrixqrunpackq<Precision>(a,
m, n,
tau, ncu, u);
7864 for(i=0; i<=n-1; i++)
7866 for(j=0; j<=i-1; j++)
7871 bidiagonal::rmatrixbd<Precision>(a, n, n, tauq, taup);
7872 bidiagonal::rmatrixbdunpackpt<Precision>(a, n, n, taup, nrvt, vt);
7873 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, n, n, isupper,
w, e);
7874 if( additionalmemory<1 )
7880 bidiagonal::rmatrixbdmultiplybyq<Precision>(a, n, n, tauq, u,
m, n,
true,
false);
7881 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, n, isupper,
false, u,
m, a, 0, vt, ncvt);
7890 bidiagonal::rmatrixbdunpackq<Precision>(a, n, n, tauq, n, t2);
7891 blas::copymatrix<Precision>(u, 0, m-1, 0, n-1, a, 0, m-1, 0, n-1);
7892 blas::inplacetranspose<Precision>(t2, 0, n-1, 0, n-1, work);
7893 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, n, isupper,
false, u, 0, t2, n, vt, ncvt);
7894 blas::matrixmatrixmultiply<Precision>(a, 0, m-1, 0, n-1,
false, t2, 0, n-1, 0, n-1,
true,
amp::ampf<Precision>(
"1.0"), u, 0, m-1, 0, n-1,
amp::ampf<Precision>(
"0.0"), work);
7912 lq::rmatrixlq<Precision>(a,
m, n,
tau);
7913 for(i=0; i<=m-1; i++)
7915 for(j=i+1; j<=m-1; j++)
7920 bidiagonal::rmatrixbd<Precision>(a,
m,
m, tauq, taup);
7921 bidiagonal::rmatrixbdunpackq<Precision>(a,
m,
m, tauq, ncu, u);
7922 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a,
m,
m, isupper,
w, e);
7924 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7925 result = bdsvd::rmatrixbdsvd<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, 0);
7926 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7935 lq::rmatrixlq<Precision>(a,
m, n,
tau);
7936 lq::rmatrixlqunpackq<Precision>(a,
m, n,
tau, nrvt, vt);
7937 for(i=0; i<=m-1; i++)
7939 for(j=i+1; j<=m-1; j++)
7944 bidiagonal::rmatrixbd<Precision>(a,
m,
m, tauq, taup);
7945 bidiagonal::rmatrixbdunpackq<Precision>(a,
m,
m, tauq, ncu, u);
7946 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a,
m,
m, isupper,
w, e);
7948 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7949 if( additionalmemory<1 )
7955 bidiagonal::rmatrixbdmultiplybyp<Precision>(a,
m,
m, taup, vt,
m, n,
false,
true);
7956 result = bdsvd::rmatrixbdsvd<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, n);
7964 bidiagonal::rmatrixbdunpackpt<Precision>(a,
m,
m, taup,
m, t2);
7965 result = bdsvd::rmatrixbdsvd<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, t2,
m);
7966 blas::copymatrix<Precision>(vt, 0, m-1, 0, n-1, a, 0, m-1, 0, n-1);
7967 blas::matrixmatrixmultiply<Precision>(t2, 0, m-1, 0, m-1,
false, a, 0, m-1, 0, n-1,
false,
amp::ampf<Precision>(
"1.0"), vt, 0, m-1, 0, n-1,
amp::ampf<Precision>(
"0.0"), work);
7969 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7980 bidiagonal::rmatrixbd<Precision>(a,
m, n, tauq, taup);
7981 bidiagonal::rmatrixbdunpackq<Precision>(a,
m, n, tauq, ncu, u);
7982 bidiagonal::rmatrixbdunpackpt<Precision>(a,
m, n, taup, nrvt, vt);
7983 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a,
m, n, isupper,
w, e);
7985 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7986 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, minmn, isupper,
false, a, 0, u, nru, vt, ncvt);
7987 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7994 bidiagonal::rmatrixbd<Precision>(a,
m, n, tauq, taup);
7995 bidiagonal::rmatrixbdunpackq<Precision>(a,
m, n, tauq, ncu, u);
7996 bidiagonal::rmatrixbdunpackpt<Precision>(a,
m, n, taup, nrvt, vt);
7997 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a,
m, n, isupper,
w, e);
7998 if( additionalmemory<2 || uneeded==0 )
8004 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, minmn, isupper,
false, u, nru, a, 0, vt, ncvt);
8013 blas::copyandtranspose<Precision>(u, 0, m-1, 0, minmn-1, t2, 0, minmn-1, 0, m-1);
8014 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, minmn, isupper,
false, u, 0, t2,
m, vt, ncvt);
8015 blas::copyandtranspose<Precision>(t2, 0, minmn-1, 0, m-1, u, 0, m-1, 0, minmn-1);
8025 template<
unsigned int Precision>
8031 int additionalmemory,
8068 w.setbounds(1, minmn);
8075 u.setbounds(1, nru, 1, ncu);
8081 u.setbounds(1, nru, 1, ncu);
8089 vt.setbounds(1, nrvt, 1, ncvt);
8095 vt.setbounds(1, nrvt, 1, ncvt);
8110 qr::qrdecomposition<Precision>(a,
m, n,
tau);
8113 for(j=1; j<=i-1; j++)
8118 bidiagonal::tobidiagonal<Precision>(a, n, n, tauq, taup);
8119 bidiagonal::unpackptfrombidiagonal<Precision>(a, n, n, taup, nrvt, vt);
8120 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, n, n, isupper,
w, e);
8121 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, n, isupper,
false, u, 0, a, 0, vt, ncvt);
8130 qr::qrdecomposition<Precision>(a,
m, n,
tau);
8131 qr::unpackqfromqr<Precision>(a,
m, n,
tau, ncu, u);
8134 for(j=1; j<=i-1; j++)
8139 bidiagonal::tobidiagonal<Precision>(a, n, n, tauq, taup);
8140 bidiagonal::unpackptfrombidiagonal<Precision>(a, n, n, taup, nrvt, vt);
8141 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, n, n, isupper,
w, e);
8142 if( additionalmemory<1 )
8148 bidiagonal::multiplybyqfrombidiagonal<Precision>(a, n, n, tauq, u,
m, n,
true,
false);
8149 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, n, isupper,
false, u,
m, a, 0, vt, ncvt);
8158 bidiagonal::unpackqfrombidiagonal<Precision>(a, n, n, tauq, n, t2);
8159 blas::copymatrix<Precision>(u, 1,
m, 1, n, a, 1,
m, 1, n);
8160 blas::inplacetranspose<Precision>(t2, 1, n, 1, n, work);
8161 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, n, isupper,
false, u, 0, t2, n, vt, ncvt);
8162 blas::matrixmatrixmultiply<Precision>(a, 1,
m, 1, n,
false, t2, 1, n, 1, n,
true,
amp::ampf<Precision>(
"1.0"), u, 1, m, 1, n,
amp::ampf<Precision>(
"0.0"), work);
8180 lq::lqdecomposition<Precision>(a,
m, n,
tau);
8181 for(i=1; i<=m-1; i++)
8183 for(j=i+1; j<=
m; j++)
8188 bidiagonal::tobidiagonal<Precision>(a,
m,
m, tauq, taup);
8189 bidiagonal::unpackqfrombidiagonal<Precision>(a,
m,
m, tauq, ncu, u);
8190 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a,
m,
m, isupper,
w, e);
8192 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8193 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, 0);
8194 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8203 lq::lqdecomposition<Precision>(a,
m, n,
tau);
8204 lq::unpackqfromlq<Precision>(a,
m, n,
tau, nrvt, vt);
8205 for(i=1; i<=m-1; i++)
8207 for(j=i+1; j<=
m; j++)
8212 bidiagonal::tobidiagonal<Precision>(a,
m,
m, tauq, taup);
8213 bidiagonal::unpackqfrombidiagonal<Precision>(a,
m,
m, tauq, ncu, u);
8214 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a,
m,
m, isupper,
w, e);
8216 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8217 if( additionalmemory<1 )
8223 bidiagonal::multiplybypfrombidiagonal<Precision>(a,
m,
m, taup, vt,
m, n,
false,
true);
8224 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, n);
8232 bidiagonal::unpackptfrombidiagonal<Precision>(a,
m,
m, taup,
m, t2);
8233 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, t2,
m);
8234 blas::copymatrix<Precision>(vt, 1,
m, 1, n, a, 1,
m, 1, n);
8235 blas::matrixmatrixmultiply<Precision>(t2, 1,
m, 1,
m,
false, a, 1,
m, 1, n,
false,
amp::ampf<Precision>(
"1.0"), vt, 1, m, 1, n,
amp::ampf<Precision>(
"0.0"), work);
8237 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8248 bidiagonal::tobidiagonal<Precision>(a,
m, n, tauq, taup);
8249 bidiagonal::unpackqfrombidiagonal<Precision>(a,
m, n, tauq, ncu, u);
8250 bidiagonal::unpackptfrombidiagonal<Precision>(a,
m, n, taup, nrvt, vt);
8251 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a,
m, n, isupper,
w, e);
8253 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8254 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, minmn, isupper,
false, a, 0, u, nru, vt, ncvt);
8255 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8262 bidiagonal::tobidiagonal<Precision>(a,
m, n, tauq, taup);
8263 bidiagonal::unpackqfrombidiagonal<Precision>(a,
m, n, tauq, ncu, u);
8264 bidiagonal::unpackptfrombidiagonal<Precision>(a,
m, n, taup, nrvt, vt);
8265 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a,
m, n, isupper,
w, e);
8266 if( additionalmemory<2 || uneeded==0 )
8272 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, minmn, isupper,
false, u, nru, a, 0, vt, ncvt);
8281 blas::copyandtranspose<Precision>(u, 1,
m, 1, minmn, t2, 1, minmn, 1,
m);
8282 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, minmn, isupper,
false, u, 0, t2,
m, vt, ncvt);
8283 blas::copyandtranspose<Precision>(t2, 1, minmn, 1,
m, u, 1,
m, 1, minmn);
const ampf< Precision > halfpi()
void rmatrixlqunpackq(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qrows, ap::template_2d_array< amp::ampf< Precision > > &q)
const CanonicalForm int s
complex & operator+=(const double &v)
const CanonicalForm int const CFList const Variable & y
void rmatrixqr(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
void setcontent(int iLow1, int iHigh1, int iLow2, int iHigh2, const T *pContent)
static const ampf getAlgoPascalMinNumber()
void mu(int **points, int sizePoints)
int getlowbound(int iBoundNum) const
void rmatrixbdmultiplybyp(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
void tobidiagonal(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_1d_array< amp::ampf< Precision > > &taup)
template_1d_array(const template_1d_array &rhs)
void copyandtranspose(const ap::template_2d_array< amp::ampf< Precision > > &a, int is1, int is2, int js1, int js2, ap::template_2d_array< amp::ampf< Precision > > &b, int id1, int id2, int jd1, int jd2)
gmp_float exp(const gmp_float &a)
bool svddecomposition(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, int uneeded, int vtneeded, int additionalmemory, ap::template_1d_array< amp::ampf< Precision > > &w, ap::template_2d_array< amp::ampf< Precision > > &u, ap::template_2d_array< amp::ampf< Precision > > &vt)
void vmoveneg(raw_vector< T > vdst, const_raw_vector< T > vsrc)
double maxreal(double m1, double m2)
void matrixmatrixmultiply(const ap::template_2d_array< amp::ampf< Precision > > &a, int ai1, int ai2, int aj1, int aj2, bool transa, const ap::template_2d_array< amp::ampf< Precision > > &b, int bi1, int bi2, int bj1, int bj2, bool transb, amp::ampf< Precision > alpha, ap::template_2d_array< amp::ampf< Precision > > &c, int ci1, int ci2, int cj1, int cj2, amp::ampf< Precision > beta, ap::template_1d_array< amp::ampf< Precision > > &work)
void vsub(raw_vector< T > vdst, const_raw_vector< T > vsrc)
void unpackqfromqr(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
template_2d_array< int > integer_2d_array
ampf & operator*=(const T &v)
complex & operator*=(const complex &z)
static const ampf getAlgoPascalEpsilon()
const ampf< Precision > log2(const ampf< Precision > &x)
void multiplybyqfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
template_1d_array< complex > complex_1d_array
const complex operator/(const complex &lhs, const complex &rhs)
const ampf< Precision > acos(const ampf< Precision > &x)
void generatereflection(ap::template_1d_array< amp::ampf< Precision > > &x, int n, amp::ampf< Precision > &tau)
const complex operator-(const complex &lhs)
void vAdd(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
campf(const ampf< Precision > &_x)
ampf & operator/=(const T &v)
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
complex & operator+=(const complex &z)
const bool operator==(const complex &lhs, const complex &rhs)
campf(const ampf< Precision > &_x, const ampf< Precision > &_y)
gmp_float log(const gmp_float &a)
const ampf< Precision > atan2(const ampf< Precision > &y, const ampf< Precision > &x)
const ampf< Precision > minimum(const ampf< Precision > &x, const ampf< Precision > &y)
const template_1d_array & operator=(const template_1d_array &rhs)
void WerrorS(const char *s)
const double machineepsilon
template_1d_array< double > real_1d_array
int rowidxabsmax(const ap::template_2d_array< amp::ampf< Precision > > &x, int j1, int j2, int i)
void rmatrixbd(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_1d_array< amp::ampf< Precision > > &taup)
static void make_assertion(bool bClause)
campf< Precision > & operator/=(campf< Precision > &lhs, const campf< Precision > &rhs)
template_1d_array< int > integer_1d_array
void copymatrix(const ap::template_2d_array< amp::ampf< Precision > > &a, int is1, int is2, int js1, int js2, ap::template_2d_array< amp::ampf< Precision > > &b, int id1, int id2, int jd1, int jd2)
int randominteger(int maxv)
ampf & operator+=(const T &v)
const T * getcontent() const
mpfr_record * mpfr_record_ptr
Rational abs(const Rational &a)
lists getList(spectrum &spec)
void unpackptfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, int ptrows, ap::template_2d_array< amp::ampf< Precision > > &pt)
void unpackdiagonalsfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &b, int m, int n, bool &isupper, ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > &e)
template_1d_array< bool > boolean_1d_array
const ampf< Precision > sinh(const ampf< Precision > &x)
int gethighbound(int iBoundNum=0) const
int columnidxabsmax(const ap::template_2d_array< amp::ampf< Precision > > &x, int i1, int i2, int j)
bool bidiagonalsvddecompositioninternal(ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > e, int n, bool isupper, bool isfractionalaccuracyrequired, ap::template_2d_array< amp::ampf< Precision > > &u, int ustart, int nru, ap::template_2d_array< amp::ampf< Precision > > &c, int cstart, int ncc, ap::template_2d_array< amp::ampf< Precision > > &vt, int vstart, int ncvt)
complex(const double &_x)
bool bidiagonalsvddecomposition(ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > e, int n, bool isupper, bool isfractionalaccuracyrequired, ap::template_2d_array< amp::ampf< Precision > > &u, int nru, ap::template_2d_array< amp::ampf< Precision > > &c, int ncc, ap::template_2d_array< amp::ampf< Precision > > &vt, int ncvt)
const complex operator*(const complex &lhs, const complex &rhs)
template_2d_array< double > real_2d_array
bool operator>=(const ampf< Precision > &op1, const ampf< Precision > &op2)
int vectoridxabsmax(const ap::template_1d_array< amp::ampf< Precision > > &x, int i1, int i2)
amp::ampf< Precision > extsignbdsqr(amp::ampf< Precision > a, amp::ampf< Precision > b)
void generaterotation(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > &cs, amp::ampf< Precision > &sn, amp::ampf< Precision > &r)
BOOLEAN fa(leftv res, leftv args)
void vMul(ap::raw_vector< ampf< Precision > > vDst, T2 alpha)
template_2d_array< bool > boolean_2d_array
const T * getcontent() const
template_2d_array(const template_2d_array &rhs)
bool wrongIdx(int i) const
const T & operator()(int i1, int i2) const
const T * GetData() const
void vmul(raw_vector< T > vdst, T2 alpha)
void vadd(raw_vector< T > vdst, const_raw_vector< T > vsrc)
#define __AMP_BINARY_OPF(type)
gmp_float sqrt(const gmp_float &a)
void qrdecompositionunpacked(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &q, ap::template_2d_array< amp::ampf< Precision > > &r)
const template_2d_array & operator=(const template_2d_array &rhs)
campf< Precision > & operator+=(campf< Precision > &lhs, const campf< Precision > &rhs)
bool rmatrixsvd(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, int uneeded, int vtneeded, int additionalmemory, ap::template_1d_array< amp::ampf< Precision > > &w, ap::template_2d_array< amp::ampf< Precision > > &u, ap::template_2d_array< amp::ampf< Precision > > &vt)
raw_vector< T > getrow(int iRow, int iColumnStart, int iColumnEnd)
void tau(int **points, int sizePoints, int k)
raw_vector< T > getcolumn(int iColumn, int iRowStart, int iRowEnd)
void vMoveNeg(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
void vSub(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
int minint(int m1, int m2)
double minreal(double m1, double m2)
bool wrongRow(int i) const
void rmatrixlq(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
campf< Precision > & operator-=(campf< Precision > &lhs, const campf< Precision > &rhs)
void rmatrixbdmultiplybyq(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
ampf(const ampf< Precision2 > &r)
void lqdecompositionunpacked(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &l, ap::template_2d_array< amp::ampf< Precision > > &q)
campf(const campf< Prec2 > &z)
ampf(const std::string &s)
gmp_float cos(const gmp_float &a)
const double minrealnumber
const_raw_vector< T > getvector(int iStart, int iEnd) const
int maxint(int m1, int m2)
void setbounds(int iLow, int iHigh)
campf< Precision > & operator*=(campf< Precision > &lhs, const campf< Precision > &rhs)
const ampf< Precision > ldexp2(const ampf< Precision > &x, mp_exp_t exponent)
T vdotproduct(const_raw_vector< T > v1, const_raw_vector< T > v2)
const ampf< Precision > tan(const ampf< Precision > &x)
raw_vector(T *Data, int Length, int Step)
const Variable & v
< [in] a sqrfree bivariate poly
template_2d_array< complex > complex_2d_array
void rmatrixqrunpackq(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
complex & operator/=(const complex &z)
complex & operator-=(const complex &z)
void svdv2x2(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > h, amp::ampf< Precision > &ssmin, amp::ampf< Precision > &ssmax, amp::ampf< Precision > &snr, amp::ampf< Precision > &csr, amp::ampf< Precision > &snl, amp::ampf< Precision > &csl)
const ampf< Precision > twopi()
void rmatrixbdunpackpt(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, int ptrows, ap::template_2d_array< amp::ampf< Precision > > &pt)
complex(const double &_x, const double &_y)
void applyreflectionfromtheright(ap::template_2d_array< amp::ampf< Precision > > &c, amp::ampf< Precision > tau, const ap::template_1d_array< amp::ampf< Precision > > &v, int m1, int m2, int n1, int n2, ap::template_1d_array< amp::ampf< Precision > > &work)
const ampf< Precision > frac(const ampf< Precision > &x)
ampf & operator-=(const T &v)
void rmatrixbdunpackq(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
const complex csqr(const complex &z)
void applyrotationsfromtheleft(bool isforward, int m1, int m2, int n1, int n2, const ap::template_1d_array< amp::ampf< Precision > > &c, const ap::template_1d_array< amp::ampf< Precision > > &s, ap::template_2d_array< amp::ampf< Precision > > &a, ap::template_1d_array< amp::ampf< Precision > > &work)
const_raw_vector< T > getrow(int iRow, int iColumnStart, int iColumnEnd) const
bool wrongColumn(int j) const
void svd2x2(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > h, amp::ampf< Precision > &ssmin, amp::ampf< Precision > &ssmax)
void applyrotationsfromtheright(bool isforward, int m1, int m2, int n1, int n2, const ap::template_1d_array< amp::ampf< Precision > > &c, const ap::template_1d_array< amp::ampf< Precision > > &s, ap::template_2d_array< amp::ampf< Precision > > &a, ap::template_1d_array< amp::ampf< Precision > > &work)
const T & operator()(int i) const
int gethighbound(int iBoundNum) const
const ampf< Precision > asin(const ampf< Precision > &x)
const ampf< Precision > tanh(const ampf< Precision > &x)
raw_vector< T > getvector(int iStart, int iEnd)
const double maxrealnumber
void rmatrixqrunpackr(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &r)
void lqdecomposition(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
signed long floor(const ampf< Precision > &x)
bool isZero(const CFArray &A)
checks if entries of A are zero
void applyreflectionfromtheleft(ap::template_2d_array< amp::ampf< Precision > > &c, amp::ampf< Precision > tau, const ap::template_1d_array< amp::ampf< Precision > > &v, int m1, int m2, int n1, int n2, ap::template_1d_array< amp::ampf< Precision > > &work)
complex & operator/=(const double &v)
complex & operator-=(const double &v)
bool operator>(const ampf< Precision > &op1, const ampf< Precision > &op2)
void inplacetranspose(ap::template_2d_array< amp::ampf< Precision > > &a, int i1, int i2, int j1, int j2, ap::template_1d_array< amp::ampf< Precision > > &work)
#define __AMP_BINARY_OPI(type)
bool rmatrixbdsvd(ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > e, int n, bool isupper, bool isfractionalaccuracyrequired, ap::template_2d_array< amp::ampf< Precision > > &u, int nru, ap::template_2d_array< amp::ampf< Precision > > &c, int ncc, ap::template_2d_array< amp::ampf< Precision > > &vt, int ncvt)
void vMove(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
const double abscomplex(const complex &z)
std::string toString(const gfan::ZCone *const c)
int exponent(const CanonicalForm &f, int q)
int exponent ( const CanonicalForm & f, int q )
complex(const complex &z)
static gmp_randstate_t * getRandState()
const ampf< Precision > frexp2(const ampf< Precision > &x, mp_exp_t *exponent)
const ampf< Precision > maximum(const ampf< Precision > &x, const ampf< Precision > &y)
const_raw_vector(const T *Data, int Length, int Step)
complex & operator*=(const double &v)
void rmatrixbdunpackdiagonals(const ap::template_2d_array< amp::ampf< Precision > > &b, int m, int n, bool &isupper, ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > &e)
amp::ampf< Precision > upperhessenberg1norm(const ap::template_2d_array< amp::ampf< Precision > > &a, int i1, int i2, int j1, int j2, ap::template_1d_array< amp::ampf< Precision > > &work)
void qrdecomposition(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
void setcontent(int iLow, int iHigh, const T *pContent)
const complex operator+(const complex &lhs)
void multiplybypfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
void rmatrixlqunpackl(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &l)
Rational pow(const Rational &a, int e)
const_raw_vector< T > getcolumn(int iColumn, int iRowStart, int iRowEnd) const
void matrixvectormultiply(const ap::template_2d_array< amp::ampf< Precision > > &a, int i1, int i2, int j1, int j2, bool trans, const ap::template_1d_array< amp::ampf< Precision > > &x, int ix1, int ix2, amp::ampf< Precision > alpha, ap::template_1d_array< amp::ampf< Precision > > &y, int iy1, int iy2, amp::ampf< Precision > beta)
gmp_float sin(const gmp_float &a)
const ampf< Precision > log10(const ampf< Precision > &x)
T & operator()(int i1, int i2)
const ampf< Precision > atan(const ampf< Precision > &x)
const bool operator!=(const complex &lhs, const complex &rhs)
amp::ampf< Precision > vectornorm2(const ap::template_1d_array< amp::ampf< Precision > > &x, int i1, int i2)
signed long ceil(const ampf< Precision > &x)
const complex conj(const complex &z)
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
void unpackqfromlq(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qrows, ap::template_2d_array< amp::ampf< Precision > > &q)
amp::ampf< Precision > pythag2(amp::ampf< Precision > x, amp::ampf< Precision > y)
const ampf< Precision > cosh(const ampf< Precision > &x)
void unpackqfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
int getlowbound(int iBoundNum=0) const