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);
1021 #include <stdexcept> 1035 class incorrectPrecision :
public exception {};
1036 class overflow :
public exception {};
1037 class divisionByZero :
public exception {};
1038 class sqrtOfNegativeNumber :
public exception {};
1039 class invalidConversion :
public exception {};
1040 class invalidString :
public exception {};
1041 class internalError :
public exception {};
1042 class domainError :
public exception {};
1049 unsigned int refCount;
1050 unsigned int Precision;
1063 static mpfr_record* newMpfr(
unsigned int Precision);
1064 static void deleteMpfr(mpfr_record* ref);
1066 static gmp_randstate_t* getRandState();
1068 static mpfr_record_ptr&
getList(
unsigned int Precision);
1074 class mpfr_reference
1078 mpfr_reference(
const mpfr_reference& r);
1079 mpfr_reference& operator= (
const mpfr_reference &r);
1085 mpfr_srcptr getReadPtr()
const;
1086 mpfr_ptr getWritePtr();
1094 template<
unsigned int Precision>
1104 if( rval->refCount==0 )
1105 mpfr_storage::deleteMpfr(rval);
1114 ampf (
long double v) { InitializeAsDouble(v); }
1115 ampf (
double v) { InitializeAsDouble(v); }
1116 ampf (
float v) { InitializeAsDouble(v); }
1117 ampf (
signed long v) { InitializeAsSLong(v); }
1118 ampf (
unsigned long v) { InitializeAsULong(v); }
1119 ampf (
signed int v) { InitializeAsSLong(v); }
1120 ampf (
unsigned int v) { InitializeAsULong(v); }
1121 ampf (
signed short v) { InitializeAsSLong(v); }
1122 ampf (
unsigned short v) { InitializeAsULong(v); }
1123 ampf (
signed char v) { InitializeAsSLong(v); }
1124 ampf (
unsigned char v) { InitializeAsULong(v); }
1131 ampf (
const char *
s) { InitializeAsString(s); }
1141 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS 1142 template<
unsigned int Precision2>
1146 rval = mpfr_storage::newMpfr(Precision);
1147 mpfr_set(getWritePtr(), r.getReadPtr(), GMP_RNDN);
1154 ampf& operator= (
long double v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
return *
this; }
1155 ampf& operator= (
double v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
return *
this; }
1156 ampf& operator= (
float v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
return *
this; }
1157 ampf& operator= (
signed long v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN);
return *
this; }
1158 ampf& operator= (
unsigned long v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
return *
this; }
1159 ampf& operator= (
signed int v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN);
return *
this; }
1160 ampf& operator= (
unsigned int v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
return *
this; }
1161 ampf& operator= (
signed short v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN);
return *
this; }
1162 ampf& operator= (
unsigned short v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
return *
this; }
1163 ampf& operator= (
signed char v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN);
return *
this; }
1164 ampf& operator= (
unsigned char v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
return *
this; }
1165 ampf& operator= (
const char *
s) { mpfr_strtofr(getWritePtr(), s,
NULL, 0, GMP_RNDN);
return *
this; }
1166 ampf& operator= (
const std::string &
s) { mpfr_strtofr(getWritePtr(), s.c_str(),
NULL, 0, GMP_RNDN);
return *
this; }
1175 if( rval->refCount==0 )
1176 mpfr_storage::deleteMpfr(rval);
1182 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS 1183 template<
unsigned int Precision2>
1186 if( (
void*)
this==(
void*)(&r) )
1188 mpfr_set(getWritePtr(), r.getReadPtr(), GMP_RNDN);
1205 mpfr_srcptr getReadPtr()
const;
1206 mpfr_ptr getWritePtr();
1211 bool isFiniteNumber()
const;
1212 bool isPositiveNumber()
const;
1214 bool isNegativeNumber()
const;
1215 const ampf getUlpOf();
1220 double toDouble()
const;
1229 static const ampf getUlpOf(
const ampf &
x);
1230 static const ampf getUlp();
1231 static const ampf getUlp256();
1232 static const ampf getUlp512();
1233 static const ampf getMaxNumber();
1234 static const ampf getMinNumber();
1235 static const ampf getAlgoPascalEpsilon();
1236 static const ampf getAlgoPascalMaxNumber();
1237 static const ampf getAlgoPascalMinNumber();
1238 static const ampf getRandom();
1240 void CheckPrecision();
1241 void InitializeAsZero();
1242 void InitializeAsSLong(
signed long v);
1243 void InitializeAsULong(
unsigned long v);
1244 void InitializeAsDouble(
long double v);
1245 void InitializeAsString(
const char *
s);
1257 template<
unsigned int Precision>
1262 WerrorS(
"incorrectPrecision");
1265 template<
unsigned int Precision>
1269 rval = mpfr_storage::newMpfr(Precision);
1270 mpfr_set_ui(getWritePtr(), 0, GMP_RNDN);
1273 template<
unsigned int Precision>
1277 rval = mpfr_storage::newMpfr(Precision);
1278 mpfr_set_si(getWritePtr(), sv, GMP_RNDN);
1281 template<
unsigned int Precision>
1285 rval = mpfr_storage::newMpfr(Precision);
1286 mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
1289 template<
unsigned int Precision>
1293 rval = mpfr_storage::newMpfr(Precision);
1294 mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
1297 template<
unsigned int Precision>
1301 rval = mpfr_storage::newMpfr(Precision);
1302 mpfr_strtofr(getWritePtr(), s,
NULL, 0, GMP_RNDN);
1305 template<
unsigned int Precision>
1311 template<
unsigned int Precision>
1314 if( rval->refCount==1 )
1316 mpfr_record *newrval = mpfr_storage::newMpfr(Precision);
1317 mpfr_set(newrval->value, rval->value, GMP_RNDN);
1323 template<
unsigned int Precision>
1326 return mpfr_number_p(getReadPtr())!=0;
1329 template<
unsigned int Precision>
1332 if( !isFiniteNumber() )
1334 return mpfr_sgn(getReadPtr())>0;
1337 template<
unsigned int Precision>
1340 return mpfr_zero_p(getReadPtr())!=0;
1343 template<
unsigned int Precision>
1346 if( !isFiniteNumber() )
1348 return mpfr_sgn(getReadPtr())<0;
1351 template<
unsigned int Precision>
1354 return getUlpOf(*
this);
1357 template<
unsigned int Precision>
1360 return mpfr_get_d(getReadPtr(), GMP_RNDN);
1363 template<
unsigned int Precision>
1369 if( !isFiniteNumber() )
1374 ptr = mpfr_get_str(
NULL, &_e, 16, 0, getReadPtr(), GMP_RNDN);
1385 signed long iexpval;
1389 ptr = mpfr_get_str(
NULL, &expval, 16, 0, getReadPtr(), GMP_RNDN);
1392 if( iexpval!=expval )
1395 sprintf(buf_e,
"%ld",
long(iexpval));
1405 mpfr_free_str(ptr2);
1409 template<
unsigned int Precision>
1417 if( !isFiniteNumber() )
1422 ptr = mpfr_get_str(
NULL, &_e, 10, 0, getReadPtr(), GMP_RNDN);
1433 signed long iexpval;
1437 ptr = mpfr_get_str(
NULL, &expval, 10, 0, getReadPtr(), GMP_RNDN);
1440 if( iexpval!=expval )
1443 sprintf(buf_e,
"%ld",
long(iexpval));
1453 mpfr_free_str(ptr2);
1456 template<
unsigned int Precision>
1459 char *toString_Block=(
char *)
omAlloc(256);
1463 if( !isFiniteNumber() )
1467 ptr = mpfr_get_str(
NULL, &_e, 10, 0, getReadPtr(), GMP_RNDN);
1468 strcpy(toString_Block, ptr);
1470 return toString_Block;
1478 signed long iexpval;
1482 ptr = mpfr_get_str(
NULL, &expval, 10, 0, getReadPtr(), GMP_RNDN);
1485 if( iexpval!=expval )
1488 sprintf(buf_e,
"%ld",
long(iexpval));
1492 sprintf(toString_Block,
"-0.%sE%s",ptr,buf_e);
1495 sprintf(toString_Block,
"0.%sE%s",ptr,buf_e);
1496 mpfr_free_str(ptr2);
1497 return toString_Block;
1500 template<
unsigned int Precision>
1503 if( !x.isFiniteNumber() )
1508 mpfr_nextabove(r.getWritePtr());
1509 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1513 mpfr_get_exp(x.getReadPtr()),
1523 template<
unsigned int Precision>
1527 mpfr_nextabove(r.getWritePtr());
1528 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1532 template<
unsigned int Precision>
1536 mpfr_nextabove(r.getWritePtr());
1537 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1546 template<
unsigned int Precision>
1550 mpfr_nextabove(r.getWritePtr());
1551 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1560 template<
unsigned int Precision>
1564 mpfr_nextbelow(r.getWritePtr());
1565 mpfr_set_exp(r.getWritePtr(),mpfr_get_emax());
1569 template<
unsigned int Precision>
1573 mpfr_set_exp(r.getWritePtr(),mpfr_get_emin());
1577 template<
unsigned int Precision>
1583 template<
unsigned int Precision>
1587 mp_exp_t e1 = mpfr_get_emax();
1588 mp_exp_t e2 = -mpfr_get_emin();
1589 mp_exp_t e = e1>e2 ? e1 : e2;
1590 mpfr_set_exp(r.getWritePtr(), e-5);
1594 template<
unsigned int Precision>
1598 mp_exp_t e1 = mpfr_get_emax();
1599 mp_exp_t e2 = -mpfr_get_emin();
1600 mp_exp_t e = e1>e2 ? e1 : e2;
1601 mpfr_set_exp(r.getWritePtr(), 2-(e-5));
1605 template<
unsigned int Precision>
1616 template<
unsigned int Precision>
1619 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())==0;
1622 template<
unsigned int Precision>
1625 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())!=0;
1628 template<
unsigned int Precision>
1631 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())<0;
1634 template<
unsigned int Precision>
1637 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())>0;
1640 template<
unsigned int Precision>
1643 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())<=0;
1646 template<
unsigned int Precision>
1649 return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())>=0;
1655 template<
unsigned int Precision>
1661 template<
unsigned int Precision>
1664 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1665 mpfr_neg(v->value, op1.getReadPtr(), GMP_RNDN);
1669 template<
unsigned int Precision>
1672 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1673 mpfr_add(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1677 template<
unsigned int Precision>
1680 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1681 mpfr_sub(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1686 template<
unsigned int Precision>
1689 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1690 mpfr_mul(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1694 template<
unsigned int Precision>
1697 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1698 mpfr_div(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1705 template<
unsigned int Precision>
1710 mpfr_sqr(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
1714 template<
unsigned int Precision>
1717 int s = mpfr_sgn(x.getReadPtr());
1725 template<
unsigned int Precision>
1730 mpfr_abs(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
1734 template<
unsigned int Precision>
1739 mpfr_max(res.getWritePtr(), x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
1743 template<
unsigned int Precision>
1748 mpfr_min(res.getWritePtr(), x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
1752 template<
unsigned int Precision>
1757 mpfr_sqrt(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
1761 template<
unsigned int Precision>
1766 mpfr_trunc(tmp.getWritePtr(), x.getReadPtr());
1767 if( mpfr_integer_p(tmp.getReadPtr())==0 )
1770 mpfr_clear_erangeflag();
1771 r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
1772 if( mpfr_erangeflag_p()!=0 )
1778 template<
unsigned int Precision>
1783 mpfr_frac(r.getWritePtr(), x.getReadPtr(), GMP_RNDN);
1787 template<
unsigned int Precision>
1792 mpfr_floor(tmp.getWritePtr(), x.getReadPtr());
1793 if( mpfr_integer_p(tmp.getReadPtr())==0 )
1796 mpfr_clear_erangeflag();
1797 r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
1798 if( mpfr_erangeflag_p()!=0 )
1804 template<
unsigned int Precision>
1809 mpfr_ceil(tmp.getWritePtr(), x.getReadPtr());
1810 if( mpfr_integer_p(tmp.getReadPtr())==0 )
1813 mpfr_clear_erangeflag();
1814 r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
1815 if( mpfr_erangeflag_p()!=0 )
1821 template<
unsigned int Precision>
1826 mpfr_round(tmp.getWritePtr(), x.getReadPtr());
1827 if( mpfr_integer_p(tmp.getReadPtr())==0 )
1830 mpfr_clear_erangeflag();
1831 r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
1832 if( mpfr_erangeflag_p()!=0 )
1838 template<
unsigned int Precision>
1843 if( !x.isFiniteNumber() )
1853 *exponent = mpfr_get_exp(r.getReadPtr());
1854 mpfr_set_exp(r.getWritePtr(),0);
1858 template<
unsigned int Precision>
1863 mpfr_mul_2si(r.getWritePtr(), x.getReadPtr(),
exponent, GMP_RNDN);
1870 #define __AMP_BINARY_OPI(type) \ 1871 template<unsigned int Precision> const ampf<Precision> operator+(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \ 1872 template<unsigned int Precision> const ampf<Precision> operator+(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \ 1873 template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const signed type& op2) { return op1+ampf<Precision>(op2); } \ 1874 template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const unsigned type& op2) { return op1+ampf<Precision>(op2); } \ 1875 template<unsigned int Precision> const ampf<Precision> operator-(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \ 1876 template<unsigned int Precision> const ampf<Precision> operator-(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \ 1877 template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const signed type& op2) { return op1-ampf<Precision>(op2); } \ 1878 template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const unsigned type& op2) { return op1-ampf<Precision>(op2); } \ 1879 template<unsigned int Precision> const ampf<Precision> operator*(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \ 1880 template<unsigned int Precision> const ampf<Precision> operator*(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \ 1881 template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const signed type& op2) { return op1*ampf<Precision>(op2); } \ 1882 template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const unsigned type& op2) { return op1*ampf<Precision>(op2); } \ 1883 template<unsigned int Precision> const ampf<Precision> operator/(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \ 1884 template<unsigned int Precision> const ampf<Precision> operator/(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \ 1885 template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const signed type& op2) { return op1/ampf<Precision>(op2); } \ 1886 template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const unsigned type& op2) { return op1/ampf<Precision>(op2); } \ 1887 template<unsigned int Precision> bool operator==(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \ 1888 template<unsigned int Precision> bool operator==(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \ 1889 template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const signed type& op2) { return op1==ampf<Precision>(op2); } \ 1890 template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const unsigned type& op2) { return op1==ampf<Precision>(op2); } \ 1891 template<unsigned int Precision> bool operator!=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \ 1892 template<unsigned int Precision> bool operator!=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \ 1893 template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const signed type& op2) { return op1!=ampf<Precision>(op2); } \ 1894 template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const unsigned type& op2) { return op1!=ampf<Precision>(op2); } \ 1895 template<unsigned int Precision> bool operator<=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \ 1896 template<unsigned int Precision> bool operator<=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \ 1897 template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const signed type& op2) { return op1<=ampf<Precision>(op2); } \ 1898 template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const unsigned type& op2) { return op1<=ampf<Precision>(op2); } \ 1899 template<unsigned int Precision> bool operator>=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \ 1900 template<unsigned int Precision> bool operator>=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \ 1901 template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const signed type& op2) { return op1>=ampf<Precision>(op2); } \ 1902 template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const unsigned type& op2) { return op1>=ampf<Precision>(op2); } \ 1903 template<unsigned int Precision> bool operator<(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \ 1904 template<unsigned int Precision> bool operator<(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \ 1905 template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const signed type& op2) { return op1<ampf<Precision>(op2); } \ 1906 template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const unsigned type& op2) { return op1<ampf<Precision>(op2); } \ 1907 template<unsigned int Precision> bool operator>(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \ 1908 template<unsigned int Precision> bool operator>(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \ 1909 template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const signed type& op2) { return op1>ampf<Precision>(op2); } \ 1910 template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const unsigned type& op2) { return op1>ampf<Precision>(op2); } 1915 #undef __AMP_BINARY_OPI 1916 #define __AMP_BINARY_OPF(type) \ 1917 template<unsigned int Precision> const ampf<Precision> operator+(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \ 1918 template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const type& op2) { return op1+ampf<Precision>(op2); } \ 1919 template<unsigned int Precision> const ampf<Precision> operator-(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \ 1920 template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const type& op2) { return op1-ampf<Precision>(op2); } \ 1921 template<unsigned int Precision> const ampf<Precision> operator*(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \ 1922 template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const type& op2) { return op1*ampf<Precision>(op2); } \ 1923 template<unsigned int Precision> const ampf<Precision> operator/(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \ 1924 template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const type& op2) { return op1/ampf<Precision>(op2); } \ 1925 template<unsigned int Precision> bool operator==(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \ 1926 template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const type& op2) { return op1==ampf<Precision>(op2); } \ 1927 template<unsigned int Precision> bool operator!=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \ 1928 template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const type& op2) { return op1!=ampf<Precision>(op2); } \ 1929 template<unsigned int Precision> bool operator<=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \ 1930 template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const type& op2) { return op1<=ampf<Precision>(op2); } \ 1931 template<unsigned int Precision> bool operator>=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \ 1932 template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const type& op2) { return op1>=ampf<Precision>(op2); } \ 1933 template<unsigned int Precision> bool operator<(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \ 1934 template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const type& op2) { return op1<ampf<Precision>(op2); } \ 1935 template<unsigned int Precision> bool operator>(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \ 1936 template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const type& op2) { return op1>ampf<Precision>(op2); } 1940 #undef __AMP_BINARY_OPF 1945 template<
unsigned int Precision>
1948 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1949 mpfr_const_pi(v->value, GMP_RNDN);
1953 template<
unsigned int Precision>
1956 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1957 mpfr_const_pi(v->value, GMP_RNDN);
1958 mpfr_mul_2si(v->value, v->value, -1, GMP_RNDN);
1962 template<
unsigned int Precision>
1965 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1966 mpfr_const_pi(v->value, GMP_RNDN);
1967 mpfr_mul_2si(v->value, v->value, +1, GMP_RNDN);
1971 template<
unsigned int Precision>
1974 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1975 mpfr_sin(v->value, x.getReadPtr(), GMP_RNDN);
1979 template<
unsigned int Precision>
1982 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1983 mpfr_cos(v->value, x.getReadPtr(), GMP_RNDN);
1987 template<
unsigned int Precision>
1990 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1991 mpfr_tan(v->value, x.getReadPtr(), GMP_RNDN);
1995 template<
unsigned int Precision>
1998 mpfr_record *v = mpfr_storage::newMpfr(Precision);
1999 mpfr_asin(v->value, x.getReadPtr(), GMP_RNDN);
2003 template<
unsigned int Precision>
2006 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2007 mpfr_acos(v->value, x.getReadPtr(), GMP_RNDN);
2011 template<
unsigned int Precision>
2014 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2015 mpfr_atan(v->value, x.getReadPtr(), GMP_RNDN);
2019 template<
unsigned int Precision>
2022 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2023 mpfr_atan2(v->value, y.getReadPtr(), x.getReadPtr(), GMP_RNDN);
2027 template<
unsigned int Precision>
2030 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2031 mpfr_log(v->value, x.getReadPtr(), GMP_RNDN);
2035 template<
unsigned int Precision>
2038 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2039 mpfr_log2(v->value, x.getReadPtr(), GMP_RNDN);
2043 template<
unsigned int Precision>
2046 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2047 mpfr_log10(v->value, x.getReadPtr(), GMP_RNDN);
2051 template<
unsigned int Precision>
2054 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2055 mpfr_exp(v->value, x.getReadPtr(), GMP_RNDN);
2059 template<
unsigned int Precision>
2062 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2063 mpfr_sinh(v->value, x.getReadPtr(), GMP_RNDN);
2067 template<
unsigned int Precision>
2070 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2071 mpfr_cosh(v->value, x.getReadPtr(), GMP_RNDN);
2075 template<
unsigned int Precision>
2078 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2079 mpfr_tanh(v->value, x.getReadPtr(), GMP_RNDN);
2083 template<
unsigned int Precision>
2086 mpfr_record *v = mpfr_storage::newMpfr(Precision);
2087 mpfr_pow(v->value, x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
2094 template<
unsigned int Precision>
2113 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS 2114 template<
unsigned int Prec2>
2118 campf& operator= (
long double v) { x=
v; y=0;
return *
this; }
2119 campf& operator= (
double v) { x=
v; y=0;
return *
this; }
2120 campf& operator= (
float v) { x=
v; y=0;
return *
this; }
2121 campf& operator= (
signed long v) { x=
v; y=0;
return *
this; }
2122 campf& operator= (
unsigned long v) { x=
v; y=0;
return *
this; }
2123 campf& operator= (
signed int v) { x=
v; y=0;
return *
this; }
2124 campf& operator= (
unsigned int v) { x=
v; y=0;
return *
this; }
2125 campf& operator= (
signed short v) { x=
v; y=0;
return *
this; }
2126 campf& operator= (
unsigned short v) { x=
v; y=0;
return *
this; }
2127 campf& operator= (
signed char v) { x=
v; y=0;
return *
this; }
2128 campf& operator= (
unsigned char v) { x=
v; y=0;
return *
this; }
2129 campf& operator= (
const char *s) { x=
s; y=0;
return *
this; }
2137 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS 2138 template<
unsigned int Precision2>
2153 template<
unsigned int Precision>
2155 {
return lhs.
x==rhs.
x && lhs.
y==rhs.
y; }
2157 template<
unsigned int Precision>
2159 {
return lhs.
x!=rhs.
x || lhs.
y!=rhs.
y; }
2161 template<
unsigned int Precision>
2165 template<
unsigned int Precision>
2167 { lhs.
x += rhs.
x; lhs.
y += rhs.
y;
return lhs; }
2169 template<
unsigned int Precision>
2173 template<
unsigned int Precision>
2177 template<
unsigned int Precision>
2179 { lhs.
x -= rhs.
x; lhs.
y -= rhs.
y;
return lhs; }
2181 template<
unsigned int Precision>
2185 template<
unsigned int Precision>
2194 template<
unsigned int Precision>
2198 template<
unsigned int Precision>
2208 result.
x = (lhs.
x+lhs.
y*e)/f;
2209 result.
y = (lhs.
y-lhs.
x*e)/f;
2215 result.
x = (lhs.
y+lhs.
x*e)/f;
2216 result.
y = (-lhs.
x+lhs.
y*e)/f;
2221 template<
unsigned int Precision>
2228 template<
unsigned int Precision>
2235 w = xabs>yabs ? xabs : yabs;
2236 v = xabs<yabs ? xabs : yabs;
2246 template<
unsigned int Precision>
2252 template<
unsigned int Precision>
2262 #define __AMP_BINARY_OPI(type) \ 2263 template<unsigned int Precision> const campf<Precision> operator+ (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \ 2264 template<unsigned int Precision> const campf<Precision> operator+ (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \ 2265 template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \ 2266 template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \ 2267 template<unsigned int Precision> const campf<Precision> operator- (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \ 2268 template<unsigned int Precision> const campf<Precision> operator- (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \ 2269 template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \ 2270 template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \ 2271 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); } \ 2272 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); } \ 2273 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); } \ 2274 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); } \ 2275 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; } \ 2276 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; } \ 2277 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); } \ 2278 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); } \ 2279 template<unsigned int Precision> bool operator==(const signed type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \ 2280 template<unsigned int Precision> bool operator==(const unsigned type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \ 2281 template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const signed type& op2) { return op1.x==op2 && op1.y==0; } \ 2282 template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const unsigned type& op2) { return op1.x==op2 && op1.y==0; } \ 2283 template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const signed type& op2) { return op1.x!=op2 || op1.y!=0; } \ 2284 template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const unsigned type& op2) { return op1.x!=op2 || op1.y!=0; } \ 2285 template<unsigned int Precision> bool operator!=(const signed type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } \ 2286 template<unsigned int Precision> bool operator!=(const unsigned type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } 2291 #undef __AMP_BINARY_OPI 2292 #define __AMP_BINARY_OPF(type) \ 2293 template<unsigned int Precision> const campf<Precision> operator+ (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \ 2294 template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \ 2295 template<unsigned int Precision> const campf<Precision> operator- (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \ 2296 template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \ 2297 template<unsigned int Precision> const campf<Precision> operator* (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \ 2298 template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \ 2299 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; } \ 2300 template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \ 2301 template<unsigned int Precision> bool operator==(const type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \ 2302 template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const type& op2) { return op1.x==op2 && op1.y==0; } \ 2303 template<unsigned int Precision> bool operator!=(const type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } \ 2304 template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const type& op2) { return op1.x!=op2 || op1.y!=0; } 2309 #undef __AMP_BINARY_OPF 2314 template<
unsigned int Precision>
2318 int i, cnt = v1.GetLength();
2319 const ampf<Precision> *p1 = v1.GetData();
2320 const ampf<Precision> *p2 = v2.GetData();
2321 mpfr_record *r =
NULL;
2322 mpfr_record *t =
NULL;
2325 r = mpfr_storage::newMpfr(Precision);
2326 t = mpfr_storage::newMpfr(Precision);
2327 mpfr_set_ui(r->value, 0, GMP_RNDN);
2328 for(i=0; i<cnt; i++)
2330 mpfr_mul(t->value, p1->getReadPtr(), p2->getReadPtr(), GMP_RNDN);
2331 mpfr_add(r->value, r->value, t->value, GMP_RNDN);
2335 mpfr_storage::deleteMpfr(t);
2348 template<
unsigned int Precision>
2352 int i, cnt = vDst.GetLength();
2353 ampf<Precision> *pDst = vDst.GetData();
2354 const ampf<Precision> *pSrc = vSrc.GetData();
2357 for(i=0; i<cnt; i++)
2360 pDst += vDst.GetStep();
2361 pSrc += vSrc.GetStep();
2365 template<
unsigned int Precision>
2369 int i, cnt = vDst.GetLength();
2370 ampf<Precision> *pDst = vDst.GetData();
2371 const ampf<Precision> *pSrc = vSrc.GetData();
2372 for(i=0; i<cnt; i++)
2375 mpfr_ptr v = pDst->getWritePtr();
2376 mpfr_neg(v, v, GMP_RNDN);
2377 pDst += vDst.GetStep();
2378 pSrc += vSrc.GetStep();
2382 template<
unsigned int Precision,
class T2>
2386 int i, cnt = vDst.GetLength();
2387 ampf<Precision> *pDst = vDst.GetData();
2388 const ampf<Precision> *pSrc = vSrc.GetData();
2389 ampf<Precision> a(alpha);
2390 for(i=0; i<cnt; i++)
2393 mpfr_ptr v = pDst->getWritePtr();
2394 mpfr_mul(v, v, a.getReadPtr(), GMP_RNDN);
2395 pDst += vDst.GetStep();
2396 pSrc += vSrc.GetStep();
2400 template<
unsigned int Precision>
2404 int i, cnt = vDst.GetLength();
2405 ampf<Precision> *pDst = vDst.GetData();
2406 const ampf<Precision> *pSrc = vSrc.GetData();
2407 for(i=0; i<cnt; i++)
2409 mpfr_ptr v = pDst->getWritePtr();
2410 mpfr_srcptr vs = pSrc->getReadPtr();
2411 mpfr_add(v, v, vs, GMP_RNDN);
2412 pDst += vDst.GetStep();
2413 pSrc += vSrc.GetStep();
2417 template<
unsigned int Precision,
class T2>
2421 int i, cnt = vDst.GetLength();
2422 ampf<Precision> *pDst = vDst.GetData();
2423 const ampf<Precision> *pSrc = vSrc.GetData();
2424 ampf<Precision> a(alpha), tmp;
2425 for(i=0; i<cnt; i++)
2427 mpfr_ptr v = pDst->getWritePtr();
2428 mpfr_srcptr vs = pSrc->getReadPtr();
2429 mpfr_mul(tmp.getWritePtr(), a.getReadPtr(), vs, GMP_RNDN);
2430 mpfr_add(v, v, tmp.getWritePtr(), GMP_RNDN);
2431 pDst += vDst.GetStep();
2432 pSrc += vSrc.GetStep();
2436 template<
unsigned int Precision>
2440 int i, cnt = vDst.GetLength();
2441 ampf<Precision> *pDst = vDst.GetData();
2442 const ampf<Precision> *pSrc = vSrc.GetData();
2443 for(i=0; i<cnt; i++)
2445 mpfr_ptr v = pDst->getWritePtr();
2446 mpfr_srcptr vs = pSrc->getReadPtr();
2447 mpfr_sub(v, v, vs, GMP_RNDN);
2448 pDst += vDst.GetStep();
2449 pSrc += vSrc.GetStep();
2453 template<
unsigned int Precision,
class T2>
2456 vAdd(vDst, vSrc, -alpha);
2459 template<
unsigned int Precision,
class T2>
2462 int i, cnt = vDst.GetLength();
2463 ampf<Precision> *pDst = vDst.GetData();
2464 ampf<Precision> a(alpha);
2465 for(i=0; i<cnt; i++)
2467 mpfr_ptr v = pDst->getWritePtr();
2468 mpfr_mul(v, a.getReadPtr(),
v, GMP_RNDN);
2469 pDst += vDst.GetStep();
2516 template<
unsigned int Precision>
2520 template<
unsigned int Precision>
2529 template<
unsigned int Precision>
2580 template<
unsigned int Precision>
2610 mx = amp::maximum<Precision>(amp::abs<Precision>(
x(j)), mx);
2617 xnorm = xnorm+amp::sqr<Precision>(
x(j)/mx);
2619 xnorm = amp::sqrt<Precision>(xnorm)*mx;
2634 mx = amp::maximum<Precision>(amp::abs<Precision>(
alpha), amp::abs<Precision>(xnorm));
2635 beta = -mx*amp::sqrt<Precision>(amp::sqr<Precision>(alpha/mx)+amp::sqr<Precision>(xnorm/mx));
2640 tau = (beta-
alpha)/beta;
2675 template<
unsigned int Precision>
2688 if( tau==0 || n1>n2 || m1>m2 )
2696 for(i=n1; i<=n2; i++)
2700 for(i=m1; i<=m2; i++)
2703 ap::vadd(work.getvector(n1, n2), c.getrow(i, n1, n2), t);
2709 for(i=m1; i<=m2; i++)
2712 ap::vsub(c.getrow(i, n1, n2), work.getvector(n1, n2), t);
2745 template<
unsigned int Precision>
2760 if( tau==0 || n1>n2 || m1>m2 )
2769 for(i=m1; i<=m2; i++)
2778 for(i=m1; i<=m2; i++)
2781 ap::vsub(c.getrow(i, n1, n2), v.getvector(1, vm), t);
2828 template<
unsigned int Precision>
2834 template<
unsigned int Precision>
2841 template<
unsigned int Precision>
2851 template<
unsigned int Precision>
2858 template<
unsigned int Precision>
2868 template<
unsigned int Precision>
2875 template<
unsigned int Precision>
2881 template<
unsigned int Precision>
2888 template<
unsigned int Precision>
2898 template<
unsigned int Precision>
2905 template<
unsigned int Precision>
2915 template<
unsigned int Precision>
2974 template<
unsigned int Precision>
3004 tauq.setbounds(0, n-1);
3005 taup.setbounds(0, n-1);
3009 tauq.setbounds(0, m-1);
3010 taup.setbounds(0, m-1);
3018 for(i=0; i<=n-1; i++)
3025 reflections::generatereflection<Precision>(t, m-
i, ltau);
3027 ap::vmove(a.getcolumn(i, i, m-1), t.getvector(1, m-i));
3033 reflections::applyreflectionfromtheleft<Precision>(a, ltau, t,
i, m-1, i+1, n-1, work);
3041 ap::vmove(t.getvector(1, n-i-1), a.getrow(i, i+1, n-1));
3042 reflections::generatereflection<Precision>(t, n-1-
i, ltau);
3044 ap::vmove(a.getrow(i, i+1, n-1), t.getvector(1, n-1-i));
3050 reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1, m-1, i+1, n-1, work);
3064 for(i=0; i<=m-1; i++)
3071 reflections::generatereflection<Precision>(t, n-
i, ltau);
3073 ap::vmove(a.getrow(i, i, n-1), t.getvector(1, n-i));
3079 reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1, m-1,
i, n-1, work);
3087 ap::vmove(t.getvector(1, m-1-i), a.getcolumn(i, i+1, m-1));
3088 reflections::generatereflection<Precision>(t, m-1-
i, ltau);
3090 ap::vmove(a.getcolumn(i, i+1, m-1), t.getvector(1, m-1-i));
3096 reflections::applyreflectionfromtheleft<Precision>(a, ltau, t, i+1, m-1, i+1, n-1, work);
3128 template<
unsigned int Precision>
3142 if( m==0 || n==0 || qcolumns==0 )
3150 q.setbounds(0, m-1, 0, qcolumns-1);
3151 for(i=0; i<=m-1; i++)
3153 for(j=0; j<=qcolumns-1; j++)
3169 rmatrixbdmultiplybyq<Precision>(qp,
m, n, tauq, q,
m, qcolumns,
false,
false);
3202 template<
unsigned int Precision>
3222 if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
3272 reflections::applyreflectionfromtheright<Precision>(z, tauq(i),
v, 0, zrows-1,
i, m-1, work);
3276 reflections::applyreflectionfromtheleft<Precision>(z, tauq(i),
v,
i, m-1, 0, zcolumns-1, work);
3280 while( i!=i2+istep );
3320 reflections::applyreflectionfromtheright<Precision>(z, tauq(i),
v, 0, zrows-1, i+1, m-1, work);
3324 reflections::applyreflectionfromtheleft<Precision>(z, tauq(i),
v, i+1, m-1, 0, zcolumns-1, work);
3328 while( i!=i2+istep );
3355 template<
unsigned int Precision>
3369 if( m==0 || n==0 || ptrows==0 )
3377 pt.setbounds(0, ptrows-1, 0, n-1);
3378 for(i=0; i<=ptrows-1; i++)
3380 for(j=0; j<=n-1; j++)
3396 rmatrixbdmultiplybyp<Precision>(qp,
m, n, taup, pt, ptrows, n,
true,
true);
3429 template<
unsigned int Precision>
3449 if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
3503 reflections::applyreflectionfromtheright<Precision>(z, taup(i),
v, 0, zrows-1, i+1, n-1, work);
3507 reflections::applyreflectionfromtheleft<Precision>(z, taup(i),
v, i+1, n-1, 0, zcolumns-1, work);
3511 while( i!=i2+istep );
3550 reflections::applyreflectionfromtheright<Precision>(z, taup(i),
v, 0, zrows-1,
i, n-1, work);
3554 reflections::applyreflectionfromtheleft<Precision>(z, taup(i),
v,
i, n-1, 0, zcolumns-1, work);
3558 while( i!=i2+istep );
3585 template<
unsigned int Precision>
3603 d.setbounds(0, n-1);
3604 e.setbounds(0, n-1);
3605 for(i=0; i<=n-2; i++)
3610 d(n-1) =
b(n-1,n-1);
3614 d.setbounds(0, m-1);
3615 e.setbounds(0, m-1);
3616 for(i=0; i<=m-2; i++)
3621 d(m-1) =
b(m-1,m-1);
3630 template<
unsigned int Precision>
3654 taup.setbounds(1, minmn);
3655 tauq.setbounds(1, minmn);
3670 reflections::generatereflection<Precision>(t, mmip1, ltau);
3672 ap::vmove(a.getcolumn(i, i, m), t.getvector(1, mmip1));
3678 reflections::applyreflectionfromtheleft<Precision>(a, ltau, t,
i,
m, i+1, n, work);
3688 ap::vmove(t.getvector(1, nmi), a.getrow(i, ip1, n));
3689 reflections::generatereflection<Precision>(t, nmi, ltau);
3691 ap::vmove(a.getrow(i, ip1, n), t.getvector(1, nmi));
3697 reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1,
m, i+1, n, work);
3719 reflections::generatereflection<Precision>(t, nmip1, ltau);
3721 ap::vmove(a.getrow(i, i, n), t.getvector(1, nmip1));
3727 reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1,
m,
i, n, work);
3737 ap::vmove(t.getvector(1, mmi), a.getcolumn(i, ip1, m));
3738 reflections::generatereflection<Precision>(t, mmi, ltau);
3740 ap::vmove(a.getcolumn(i, ip1, m), t.getvector(1, mmi));
3746 reflections::applyreflectionfromtheleft<Precision>(a, ltau, t, i+1,
m, i+1, n, work);
3761 template<
unsigned int Precision>
3778 if( m==0 || n==0 || qcolumns==0 )
3786 q.setbounds(1, m, 1, qcolumns);
3795 for(j=1; j<=qcolumns; j++)
3814 reflections::applyreflectionfromtheleft<Precision>(q, tauq(i),
v,
i,
m, 1, qcolumns, work);
3819 for(i=
ap::minint(m-1, qcolumns-1); i>=1; i--)
3825 reflections::applyreflectionfromtheleft<Precision>(q, tauq(i),
v, i+1,
m, 1, qcolumns, work);
3835 template<
unsigned int Precision>
3857 if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
3908 reflections::applyreflectionfromtheright<Precision>(z, tauq(i),
v, 1, zrows,
i,
m, work);
3912 reflections::applyreflectionfromtheleft<Precision>(z, tauq(i),
v,
i,
m, 1, zcolumns, work);
3916 while( i!=i2+istep );
3958 reflections::applyreflectionfromtheright<Precision>(z, tauq(i),
v, 1, zrows, i+1,
m, work);
3962 reflections::applyreflectionfromtheleft<Precision>(z, tauq(i),
v, i+1,
m, 1, zcolumns, work);
3966 while( i!=i2+istep );
3976 template<
unsigned int Precision>
3993 if( m==0 || n==0 || ptrows==0 )
4001 pt.setbounds(1, ptrows, 1, n);
4008 for(i=1; i<=ptrows; i++)
4030 reflections::applyreflectionfromtheright<Precision>(pt, taup(i),
v, 1, ptrows, i+1, n, work);
4040 reflections::applyreflectionfromtheright<Precision>(pt, taup(i),
v, 1, ptrows,
i, n, work);
4050 template<
unsigned int Precision>
4072 if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
4128 reflections::applyreflectionfromtheright<Precision>(z, taup(i),
v, 1, zrows, i+1, n, work);
4132 reflections::applyreflectionfromtheleft<Precision>(z, taup(i),
v, i+1, n, 1, zcolumns, work);
4136 while( i!=i2+istep );
4176 reflections::applyreflectionfromtheright<Precision>(z, taup(i),
v, 1, zrows,
i, n, work);
4180 reflections::applyreflectionfromtheleft<Precision>(z, taup(i),
v,
i, n, 1, zcolumns, work);
4184 while( i!=i2+istep );
4193 template<
unsigned int Precision>
4213 for(i=1; i<=n-1; i++)
4224 for(i=1; i<=m-1; i++)
4276 template<
unsigned int Precision>
4281 template<
unsigned int Precision>
4288 template<
unsigned int Precision>
4293 template<
unsigned int Precision>
4298 template<
unsigned int Precision>
4305 template<
unsigned int Precision>
4351 template<
unsigned int Precision>
4372 tau.setbounds(0, minmn-1);
4378 for(i=0; i<=k-1; i++)
4385 reflections::generatereflection<Precision>(t, m-
i, tmp);
4387 ap::vmove(a.getcolumn(i, i, m-1), t.getvector(1, m-i));
4395 reflections::applyreflectionfromtheleft<Precision>(a,
tau(i), t,
i, m-1, i+1, n-1, work);
4421 template<
unsigned int Precision>
4438 if( m<=0 || n<=0 || qcolumns<=0 )
4448 q.setbounds(0, m-1, 0, qcolumns-1);
4451 for(i=0; i<=m-1; i++)
4453 for(j=0; j<=qcolumns-1; j++)
4469 for(i=k-1; i>=0; i--)
4477 reflections::applyreflectionfromtheleft<Precision>(q,
tau(i),
v,
i, m-1, 0, qcolumns-1, work);
4497 template<
unsigned int Precision>
4512 r.setbounds(0, m-1, 0, n-1);
4513 for(i=0; i<=n-1; i++)
4517 for(i=1; i<=m-1; i++)
4519 ap::vmove(r.getrow(i, 0, n-1), r.getrow(0, 0, n-1));
4521 for(i=0; i<=k-1; i++)
4523 ap::vmove(r.getrow(i, i, n-1), a.getrow(i, i, n-1));
4531 template<
unsigned int Precision>
4549 tau.setbounds(1, minmn);
4563 reflections::generatereflection<Precision>(t, mmip1, tmp);
4565 ap::vmove(a.getcolumn(i, i, m), t.getvector(1, mmip1));
4573 reflections::applyreflectionfromtheleft<Precision>(a,
tau(i), t,
i,
m, i+1, n, work);
4582 template<
unsigned int Precision>
4600 if( m==0 || n==0 || qcolumns==0 )
4610 q.setbounds(1, m, 1, qcolumns);
4615 for(j=1; j<=qcolumns; j++)
4640 reflections::applyreflectionfromtheleft<Precision>(q,
tau(i),
v,
i,
m, 1, qcolumns, work);
4648 template<
unsigned int Precision>
4669 q.setbounds(1, m, 1, m);
4670 r.setbounds(1, m, 1, n);
4675 qrdecomposition<Precision>(a,
m, n,
tau);
4686 ap::vmove(r.getrow(i, 1, n), r.getrow(1, 1, n));
4690 ap::vmove(r.getrow(i, i, n), a.getrow(i, i, n));
4696 unpackqfromqr<Precision>(a,
m, n,
tau,
m, q);
4736 template<
unsigned int Precision>
4741 template<
unsigned int Precision>
4748 template<
unsigned int Precision>
4753 template<
unsigned int Precision>
4758 template<
unsigned int Precision>
4765 template<
unsigned int Precision>
4807 template<
unsigned int Precision>
4826 tau.setbounds(0, minmn-1);
4828 for(i=0; i<=k-1; i++)
4835 reflections::generatereflection<Precision>(t, n-
i, tmp);
4837 ap::vmove(a.getrow(i, i, n-1), t.getvector(1, n-i));
4845 reflections::applyreflectionfromtheright<Precision>(a,
tau(i), t, i+1, m-1,
i, n-1, work);
4871 template<
unsigned int Precision>
4888 if( m<=0 || n<=0 || qrows<=0 )
4898 q.setbounds(0, qrows-1, 0, n-1);
4901 for(i=0; i<=qrows-1; i++)
4903 for(j=0; j<=n-1; j++)
4919 for(i=k-1; i>=0; i--)
4927 reflections::applyreflectionfromtheright<Precision>(q,
tau(i),
v, 0, qrows-1,
i, n-1, work);
4947 template<
unsigned int Precision>
4961 l.setbounds(0, m-1, 0, n-1);
4962 for(i=0; i<=n-1; i++)
4966 for(i=1; i<=m-1; i++)
4968 ap::vmove(
l.getrow(i, 0, n-1),
l.getrow(0, 0, n-1));
4970 for(i=0; i<=m-1; i++)
4973 ap::vmove(
l.getrow(i, 0, k), a.getrow(i, 0, k));
4982 template<
unsigned int Precision>
5002 tau.setbounds(1, minmn);
5016 reflections::generatereflection<Precision>(t, nmip1, tmp);
5018 ap::vmove(a.getrow(i, i, n), t.getvector(1, nmip1));
5026 reflections::applyreflectionfromtheright<Precision>(a,
tau(i), t, i+1,
m,
i, n, work);
5036 template<
unsigned int Precision>
5054 if( m==0 || n==0 || qrows==0 )
5064 q.setbounds(1, qrows, 1, n);
5067 for(i=1; i<=qrows; i++)
5094 reflections::applyreflectionfromtheright<Precision>(q,
tau(i),
v, 1, qrows,
i, n, work);
5102 template<
unsigned int Precision>
5118 q.setbounds(1, n, 1, n);
5119 l.setbounds(1, m, 1, n);
5124 lqdecomposition<Precision>(a,
m, n,
tau);
5147 unpackqfromlq<Precision>(a,
m, n,
tau, n, q);
5187 template<
unsigned int Precision>
5191 template<
unsigned int Precision>
5195 template<
unsigned int Precision>
5200 template<
unsigned int Precision>
5205 template<
unsigned int Precision>
5212 template<
unsigned int Precision>
5223 template<
unsigned int Precision>
5230 template<
unsigned int Precision>
5241 template<
unsigned int Precision>
5256 template<
unsigned int Precision>
5259 template<
unsigned int Precision>
5282 template<
unsigned int Precision>
5303 result = amp::abs<Precision>(
x(i1));
5308 for(ix=i1; ix<=i2; ix++)
5312 absxi = amp::abs<Precision>(
x(ix));
5315 ssq = 1+ssq*amp::sqr<Precision>(scl/absxi);
5320 ssq = ssq+amp::sqr<Precision>(absxi/scl);
5324 result = scl*amp::sqrt<Precision>(ssq);
5329 template<
unsigned int Precision>
5340 a = amp::abs<Precision>(
x(result));
5341 for(i=i1+1; i<=i2; i++)
5343 if( amp::abs<Precision>(
x(i))>amp::abs<Precision>(
x(result)) )
5352 template<
unsigned int Precision>
5364 a = amp::abs<Precision>(
x(result,j));
5365 for(i=i1+1; i<=i2; i++)
5367 if( amp::abs<Precision>(
x(i,j))>amp::abs<Precision>(
x(result,j)) )
5376 template<
unsigned int Precision>
5388 a = amp::abs<Precision>(
x(i,result));
5389 for(j=j1+1; j<=j2; j++)
5391 if( amp::abs<Precision>(
x(i,j))>amp::abs<Precision>(
x(i,result)) )
5400 template<
unsigned int Precision>
5414 for(j=j1; j<=j2; j++)
5418 for(i=i1; i<=i2; i++)
5422 work(j) = work(j)+amp::abs<Precision>(a(i,j));
5426 for(j=j1; j<=j2; j++)
5428 result = amp::maximum<Precision>(
result, work(j));
5434 template<
unsigned int Precision>
5450 if( is1>is2 || js1>js2 )
5456 for(isrc=is1; isrc<=is2; isrc++)
5458 idst = isrc-is1+id1;
5459 ap::vmove(
b.getrow(idst, jd1, jd2), a.getrow(isrc, js1, js2));
5464 template<
unsigned int Precision>
5479 if( i1>i2 || j1>j2 )
5484 for(i=i1; i<=i2-1; i++)
5490 ap::vmove(work.getvector(1, l), a.getcolumn(j, ips, i2));
5491 ap::vmove(a.getcolumn(j, ips, i2), a.getrow(i, jps, j2));
5492 ap::vmove(a.getrow(i, jps, j2), work.getvector(1, l));
5497 template<
unsigned int Precision>
5513 if( is1>is2 || js1>js2 )
5519 for(isrc=is1; isrc<=is2; isrc++)
5521 jdst = isrc-is1+jd1;
5522 ap::vmove(
b.getcolumn(jdst, id1, id2), a.getrow(isrc, js1, js2));
5527 template<
unsigned int Precision>
5553 if( i1>i2 || j1>j2 )
5565 for(i=iy1; i<=iy2; i++)
5578 for(i=i1; i<=i2; i++)
5581 y(iy1+i-i1) =
y(iy1+i-i1)+alpha*
v;
5590 if( i1>i2 || j1>j2 )
5602 for(i=iy1; i<=iy2; i++)
5615 for(i=i1; i<=i2; i++)
5617 v = alpha*
x(ix1+i-i1);
5618 ap::vadd(y.getvector(iy1, iy2), a.getrow(i, j1, j2),
v);
5624 template<
unsigned int Precision>
5635 xabs = amp::abs<Precision>(
x);
5636 yabs = amp::abs<Precision>(
y);
5637 w = amp::maximum<Precision>(xabs, yabs);
5638 z = amp::minimum<Precision>(xabs, yabs);
5645 result = w*amp::sqrt<Precision>(1+amp::sqr<Precision>(z/
w));
5651 template<
unsigned int Precision>
5712 if( arows<=0 || acols<=0 || brows<=0 || bcols<=0 )
5733 for(i=ci1; i<=ci2; i++)
5735 for(j=cj1; j<=cj2; j++)
5743 for(i=ci1; i<=ci2; i++)
5752 if( !transa && !transb )
5754 for(l=ai1; l<=ai2; l++)
5756 for(r=bi1; r<=bi2; r++)
5758 v = alpha*a(l,aj1+r-bi1);
5760 ap::vadd(c.getrow(k, cj1, cj2),
b.getrow(r, bj1, bj2),
v);
5769 if( !transa && transb )
5771 if( arows*acols<brows*bcols )
5773 for(r=bi1; r<=bi2; r++)
5775 for(l=ai1; l<=ai2; l++)
5778 c(ci1+l-ai1,cj1+r-bi1) = c(ci1+l-ai1,cj1+r-bi1)+alpha*
v;
5785 for(l=ai1; l<=ai2; l++)
5787 for(r=bi1; r<=bi2; r++)
5790 c(ci1+l-ai1,cj1+r-bi1) = c(ci1+l-ai1,cj1+r-bi1)+alpha*
v;
5800 if( transa && !transb )
5802 for(l=aj1; l<=aj2; l++)
5804 for(r=bi1; r<=bi2; r++)
5806 v = alpha*a(ai1+r-bi1,l);
5808 ap::vadd(c.getrow(k, cj1, cj2),
b.getrow(r, bj1, bj2),
v);
5817 if( transa && transb )
5819 if( arows*acols<brows*bcols )
5821 for(r=bi1; r<=bi2; r++)
5823 for(i=1; i<=crows; i++)
5827 for(l=ai1; l<=ai2; l++)
5829 v = alpha*
b(r,bj1+l-ai1);
5831 ap::vadd(work.getvector(1, crows), a.getrow(l, aj1, aj2),
v);
5833 ap::vadd(c.getcolumn(k, ci1, ci2), work.getvector(1, crows));
5839 for(l=aj1; l<=aj2; l++)
5842 ap::vmove(work.getvector(1, k), a.getcolumn(l, ai1, ai2));
5843 for(r=bi1; r<=bi2; r++)
5846 c(ci1+l-aj1,cj1+r-bi1) = c(ci1+l-aj1,cj1+r-bi1)+alpha*
v;
5897 template<
unsigned int Precision>
5907 template<
unsigned int Precision>
5917 template<
unsigned int Precision>
5950 template<
unsigned int Precision>
5968 if( m1>m2 || n1>n2 )
5984 for(j=m1; j<=m2-1; j++)
5988 if( ctemp!=1 || stemp!=0 )
5991 ap::vmove(work.getvector(n1, n2), a.getrow(jp1, n1, n2), ctemp);
5992 ap::vsub(work.getvector(n1, n2), a.getrow(j, n1, n2), stemp);
5993 ap::vmul(a.getrow(j, n1, n2), ctemp);
5994 ap::vadd(a.getrow(j, n1, n2), a.getrow(jp1, n1, n2), stemp);
5995 ap::vmove(a.getrow(jp1, n1, n2), work.getvector(n1, n2));
6005 for(j=m1; j<=m2-1; j++)
6009 if( ctemp!=1 || stemp!=0 )
6012 a(j+1,n1) = ctemp*temp-stemp*a(j,n1);
6013 a(j,n1) = stemp*temp+ctemp*a(j,n1);
6026 for(j=m2-1; j>=m1; j--)
6030 if( ctemp!=1 || stemp!=0 )
6033 ap::vmove(work.getvector(n1, n2), a.getrow(jp1, n1, n2), ctemp);
6034 ap::vsub(work.getvector(n1, n2), a.getrow(j, n1, n2), stemp);
6035 ap::vmul(a.getrow(j, n1, n2), ctemp);
6036 ap::vadd(a.getrow(j, n1, n2), a.getrow(jp1, n1, n2), stemp);
6037 ap::vmove(a.getrow(jp1, n1, n2), work.getvector(n1, n2));
6047 for(j=m2-1; j>=m1; j--)
6051 if( ctemp!=1 || stemp!=0 )
6054 a(j+1,n1) = ctemp*temp-stemp*a(j,n1);
6055 a(j,n1) = stemp*temp+ctemp*a(j,n1);
6088 template<
unsigned int Precision>
6118 for(j=n1; j<=n2-1; j++)
6122 if( ctemp!=1 || stemp!=0 )
6125 ap::vmove(work.getvector(m1, m2), a.getcolumn(jp1, m1, m2), ctemp);
6126 ap::vsub(work.getvector(m1, m2), a.getcolumn(j, m1, m2), stemp);
6127 ap::vmul(a.getcolumn(j, m1, m2), ctemp);
6128 ap::vadd(a.getcolumn(j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
6129 ap::vmove(a.getcolumn(jp1, m1, m2), work.getvector(m1, m2));
6139 for(j=n1; j<=n2-1; j++)
6143 if( ctemp!=1 || stemp!=0 )
6146 a(m1,j+1) = ctemp*temp-stemp*a(m1,j);
6147 a(m1,j) = stemp*temp+ctemp*a(m1,j);
6160 for(j=n2-1; j>=n1; j--)
6164 if( ctemp!=1 || stemp!=0 )
6167 ap::vmove(work.getvector(m1, m2), a.getcolumn(jp1, m1, m2), ctemp);
6168 ap::vsub(work.getvector(m1, m2), a.getcolumn(j, m1, m2), stemp);
6169 ap::vmul(a.getcolumn(j, m1, m2), ctemp);
6170 ap::vadd(a.getcolumn(j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
6171 ap::vmove(a.getcolumn(jp1, m1, m2), work.getvector(m1, m2));
6181 for(j=n2-1; j>=n1; j--)
6185 if( ctemp!=1 || stemp!=0 )
6188 a(m1,j+1) = ctemp*temp-stemp*a(m1,j);
6189 a(m1,j) = stemp*temp+ctemp*a(m1,j);
6205 template<
unsigned int Precision>
6234 r = amp::sqrt<Precision>(amp::sqr<Precision>(f1)+amp::sqr<Precision>(g1));
6237 if( amp::abs<Precision>(f)>amp::abs<Precision>(
g) && cs<0 )
6290 template<
unsigned int Precision>
6295 bool isfractionalaccuracyrequired,
6302 template<
unsigned int Precision>
6307 bool isfractionalaccuracyrequired,
6314 template<
unsigned int Precision>
6319 bool isfractionalaccuracyrequired,
6329 template<
unsigned int Precision>
6332 template<
unsigned int Precision>
6338 template<
unsigned int Precision>
6428 template<
unsigned int Precision>
6433 bool isfractionalaccuracyrequired,
6453 result = bidiagonalsvddecompositioninternal<Precision>(d1, e1, n, isupper, isfractionalaccuracyrequired, u, 0, nru, c, 0, ncc, vt, 0, ncvt);
6454 ap::vmove(d.getvector(0, n-1), d1.getvector(1, n));
6471 template<
unsigned int Precision>
6476 bool isfractionalaccuracyrequired,
6487 result = bidiagonalsvddecompositioninternal<Precision>(d, e, n, isupper, isfractionalaccuracyrequired, u, 1, nru, c, 1, ncc, vt, 1, ncvt);
6495 template<
unsigned int Precision>
6500 bool isfractionalaccuracyrequired,
6557 bool matrixsplitflag;
6585 ap::vmul(vt.getrow(vstart, vstart, vstart+ncvt-1), -1);
6611 for(i=1; i<=n-1; i++)
6616 for(i=1; i<=n-1; i++)
6635 for(i=1; i<=n-1; i++)
6637 rotations::generaterotation<Precision>(d(i), e(i), cs, sn, r);
6650 rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, 1+ustart-1, n+ustart-1, work0, work1, u, utemp);
6654 rotations::applyrotationsfromtheleft<Precision>(fwddir, 1+cstart-1, n+cstart-1, cstart, cend, work0, work1, c, ctemp);
6663 tolmul = amp::maximum<Precision>(10, amp::minimum<Precision>(100, amp::pow<Precision>(eps, -
amp::ampf<Precision>(
"0.125"))));
6665 if( !isfractionalaccuracyrequired )
6676 smax = amp::maximum<Precision>(smax, amp::abs<Precision>(d(i)));
6678 for(i=1; i<=n-1; i++)
6680 smax = amp::maximum<Precision>(smax, amp::abs<Precision>(e(i)));
6689 sminoa = amp::abs<Precision>(d(1));
6695 mu = amp::abs<Precision>(d(i))*(mu/(mu+amp::abs<Precision>(e(i-1))));
6696 sminoa = amp::minimum<Precision>(sminoa,
mu);
6703 sminoa = sminoa/amp::sqrt<Precision>(n);
6704 thresh = amp::maximum<Precision>(tol*sminoa, maxitr*n*n*unfl);
6712 thresh = amp::maximum<Precision>(amp::abs<Precision>(tol)*smax, maxitr*n*n*unfl);
6752 if( tol<0 && amp::abs<Precision>(d(m))<=thresh )
6756 smax = amp::abs<Precision>(d(m));
6758 matrixsplitflag =
false;
6759 for(lll=1; lll<=m-1; lll++)
6762 abss = amp::abs<Precision>(d(ll));
6763 abse = amp::abs<Precision>(e(ll));
6764 if( tol<0 && abss<=thresh )
6770 matrixsplitflag =
true;
6773 smin = amp::minimum<Precision>(smin, abss);
6774 smax = amp::maximum<Precision>(smax, amp::maximum<Precision>(abss, abse));
6776 if( !matrixsplitflag )
6808 svdv2x2<Precision>(d(m-1), e(m-1), d(m), sigmn, sigmx, sinr, cosr, sinl, cosl);
6819 mm1 = m-1+(vstart-1);
6822 ap::vmul(vt.getrow(mm0, vstart, vend), cosr);
6823 ap::vsub(vt.getrow(mm0, vstart, vend), vt.getrow(mm1, vstart, vend), sinr);
6832 ap::vmul(u.getcolumn(mm0, ustart, uend), cosl);
6833 ap::vsub(u.getcolumn(mm0, ustart, uend), u.getcolumn(mm1, ustart, uend), sinl);
6842 ap::vmul(c.getrow(mm0, cstart, cend), cosl);
6843 ap::vsub(c.getrow(mm0, cstart, cend), c.getrow(mm1, cstart, cend), sinl);
6860 if( idir==1 && amp::abs<Precision>(d(ll))<
amp::ampf<Precision>(
"1.0E-3")*amp::abs<Precision>(d(m)) )
6864 if( idir==2 && amp::abs<Precision>(d(m))<
amp::ampf<Precision>(
"1.0E-3")*amp::abs<Precision>(d(ll)) )
6868 if( ll!=oldll || m!=oldm || bchangedir )
6870 if( amp::abs<Precision>(d(ll))>=amp::abs<Precision>(d(m)) )
6898 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 )
6910 mu = amp::abs<Precision>(d(ll));
6913 for(lll=ll; lll<=m-1; lll++)
6915 if( amp::abs<Precision>(e(lll))<=tol*
mu )
6922 mu = amp::abs<Precision>(d(lll+1))*(mu/(mu+amp::abs<Precision>(e(lll))));
6923 sminl = amp::minimum<Precision>(sminl,
mu);
6938 if( amp::abs<Precision>(e(ll))<=amp::abs<Precision>(tol)*amp::abs<Precision>(d(ll)) || tol<0 && amp::abs<Precision>(e(ll))<=thresh )
6950 mu = amp::abs<Precision>(d(m));
6953 for(lll=m-1; lll>=ll; lll--)
6955 if( amp::abs<Precision>(e(lll))<=tol*
mu )
6962 mu = amp::abs<Precision>(d(lll))*(mu/(mu+amp::abs<Precision>(e(lll))));
6963 sminl = amp::minimum<Precision>(sminl,
mu);
6978 if( tol>=0 && n*tol*(sminl/smax)<=amp::maximum<Precision>(eps,
amp::ampf<Precision>(
"0.01")*tol) )
6994 sll = amp::abs<Precision>(d(ll));
6995 svd2x2<Precision>(d(m-1), e(m-1), d(m), shift, r);
6999 sll = amp::abs<Precision>(d(m));
7000 svd2x2<Precision>(d(ll), e(ll), d(ll+1), shift, r);
7008 if( amp::sqr<Precision>(shift/sll)<eps )
7034 for(i=ll; i<=m-1; i++)
7036 rotations::generaterotation<Precision>(d(i)*cs, e(i), cs, sn, r);
7041 rotations::generaterotation<Precision>(oldcs*r, d(i+1)*sn, oldcs, oldsn, tmp);
7045 work2(i-ll+1) = oldcs;
7046 work3(i-ll+1) = oldsn;
7057 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work0, work1, vt, vttemp);
7061 rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work2, work3, u, utemp);
7065 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work2, work3, c, ctemp);
7071 if( amp::abs<Precision>(e(m-1))<=thresh )
7085 for(i=m; i>=ll+1; i--)
7087 rotations::generaterotation<Precision>(d(i)*cs, e(i-1), cs, sn, r);
7092 rotations::generaterotation<Precision>(oldcs*r, d(i-1)*sn, oldcs, oldsn, tmp);
7096 work2(i-ll) = oldcs;
7097 work3(i-ll) = -oldsn;
7108 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work2, work3, vt, vttemp);
7112 rotations::applyrotationsfromtheright<Precision>(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work0, work1, u, utemp);
7116 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work0, work1, c, ctemp);
7122 if( amp::abs<Precision>(e(ll))<=thresh )
7141 f = (amp::abs<Precision>(d(ll))-shift)*(extsignbdsqr<Precision>(1, d(ll))+shift/d(ll));
7143 for(i=ll; i<=m-1; i++)
7145 rotations::generaterotation<Precision>(
f,
g, cosr, sinr, r);
7150 f = cosr*d(i)+sinr*e(i);
7151 e(i) = cosr*e(i)-sinr*d(i);
7153 d(i+1) = cosr*d(i+1);
7154 rotations::generaterotation<Precision>(
f,
g, cosl, sinl, r);
7156 f = cosl*e(i)+sinl*d(i+1);
7157 d(i+1) = cosl*d(i+1)-sinl*e(i);
7161 e(i+1) = cosl*e(i+1);
7163 work0(i-ll+1) = cosr;
7164 work1(i-ll+1) = sinr;
7165 work2(i-ll+1) = cosl;
7166 work3(i-ll+1) = sinl;
7175 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work0, work1, vt, vttemp);
7179 rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work2, work3, u, utemp);
7183 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work2, work3, c, ctemp);
7189 if( amp::abs<Precision>(e(m-1))<=thresh )
7201 f = (amp::abs<Precision>(d(m))-shift)*(extsignbdsqr<Precision>(1, d(m))+shift/d(m));
7203 for(i=m; i>=ll+1; i--)
7205 rotations::generaterotation<Precision>(
f,
g, cosr, sinr, r);
7210 f = cosr*d(i)+sinr*e(i-1);
7211 e(i-1) = cosr*e(i-1)-sinr*d(i);
7213 d(i-1) = cosr*d(i-1);
7214 rotations::generaterotation<Precision>(
f,
g, cosl, sinl, r);
7216 f = cosl*e(i-1)+sinl*d(i-1);
7217 d(i-1) = cosl*d(i-1)-sinl*e(i-1);
7221 e(i-2) = cosl*e(i-2);
7224 work1(i-ll) = -sinr;
7226 work3(i-ll) = -sinl;
7233 if( amp::abs<Precision>(e(ll))<=thresh )
7243 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work2, work3, vt, vttemp);
7247 rotations::applyrotationsfromtheright<Precision>(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work0, work1, u, utemp);
7251 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work0, work1, c, ctemp);
7276 ap::vmul(vt.getrow(i+vstart-1, vstart, vend), -1);
7285 for(i=1; i<=n-1; i++)
7293 for(j=2; j<=n+1-
i; j++)
7313 ap::vmove(vt.getrow(isub+vstart-1, vstart, vend), vt.getrow(j+vstart-1, vstart, vend));
7320 ap::vmove(u.getcolumn(isub+ustart-1, ustart, uend), u.getcolumn(j+ustart-1, ustart, uend));
7327 ap::vmove(c.getrow(isub+cstart-1, cstart, cend), c.getrow(j+cstart-1, cstart, cend));
7336 template<
unsigned int Precision>
7345 result = amp::abs<Precision>(a);
7349 result = -amp::abs<Precision>(a);
7355 template<
unsigned int Precision>
7373 fa = amp::abs<Precision>(
f);
7374 ga = amp::abs<Precision>(
g);
7375 ha = amp::abs<Precision>(
h);
7376 fhmn = amp::minimum<Precision>(
fa, ha);
7377 fhmx = amp::maximum<Precision>(
fa, ha);
7387 ssmax = amp::maximum<Precision>(fhmx, ga)*amp::sqrt<Precision>(1+amp::sqr<Precision>(amp::minimum<Precision>(fhmx, ga)/amp::maximum<Precision>(fhmx, ga)));
7395 at = (fhmx-fhmn)/fhmx;
7396 au = amp::sqr<Precision>(ga/fhmx);
7397 c = 2/(amp::sqrt<Precision>(aas*aas+au)+amp::sqrt<Precision>(at*at+au));
7412 ssmin = fhmn*fhmx/ga;
7418 at = (fhmx-fhmn)/fhmx;
7419 c = 1/(amp::sqrt<Precision>(1+amp::sqr<Precision>(aas*au))+amp::sqrt<Precision>(1+amp::sqr<Precision>(at*au)));
7421 ssmin = ssmin+ssmin;
7429 template<
unsigned int Precision>
7468 fa = amp::abs<Precision>(ft);
7470 ha = amp::abs<Precision>(
h);
7495 ga = amp::abs<Precision>(gt);
7558 s = amp::sqrt<Precision>(tt+mm);
7561 r = amp::abs<Precision>(
m);
7565 r = amp::sqrt<Precision>(l*l+mm);
7578 t = extsignbdsqr<Precision>(2, ft)*extsignbdsqr<Precision>(1, gt);
7582 t = gt/extsignbdsqr<Precision>(d, ft)+m/t;
7587 t = (m/(s+t)+m/(r+l))*(1+a);
7589 l = amp::sqrt<Precision>(t*t+4);
7592 clt = (crt+srt*
m)/a;
7617 tsign = extsignbdsqr<Precision>(1, csr)*extsignbdsqr<Precision>(1, csl)*extsignbdsqr<Precision>(1,
f);
7621 tsign = extsignbdsqr<Precision>(1, snr)*extsignbdsqr<Precision>(1, csl)*extsignbdsqr<Precision>(1,
g);
7625 tsign = extsignbdsqr<Precision>(1, snr)*extsignbdsqr<Precision>(1, snl)*extsignbdsqr<Precision>(1,
h);
7627 ssmax = extsignbdsqr<Precision>(ssmax, tsign);
7628 ssmin = extsignbdsqr<Precision>(ssmin, tsign*extsignbdsqr<Precision>(1,
f)*extsignbdsqr<Precision>(1, h));
7686 template<
unsigned int Precision>
7692 int additionalmemory,
7696 template<
unsigned int Precision>
7702 int additionalmemory,
7758 template<
unsigned int Precision>
7764 int additionalmemory,
7801 w.setbounds(1, minmn);
7808 u.setbounds(0, nru-1, 0, ncu-1);
7814 u.setbounds(0, nru-1, 0, ncu-1);
7822 vt.setbounds(0, nrvt-1, 0, ncvt-1);
7828 vt.setbounds(0, nrvt-1, 0, ncvt-1);
7843 qr::rmatrixqr<Precision>(a,
m, n,
tau);
7844 for(i=0; i<=n-1; i++)
7846 for(j=0; j<=i-1; j++)
7851 bidiagonal::rmatrixbd<Precision>(a, n, n, tauq, taup);
7852 bidiagonal::rmatrixbdunpackpt<Precision>(a, n, n, taup, nrvt, vt);
7853 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, n, n, isupper,
w, e);
7854 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, n, isupper,
false, u, 0, a, 0, vt, ncvt);
7863 qr::rmatrixqr<Precision>(a,
m, n,
tau);
7864 qr::rmatrixqrunpackq<Precision>(a,
m, n,
tau, ncu, u);
7865 for(i=0; i<=n-1; i++)
7867 for(j=0; j<=i-1; j++)
7872 bidiagonal::rmatrixbd<Precision>(a, n, n, tauq, taup);
7873 bidiagonal::rmatrixbdunpackpt<Precision>(a, n, n, taup, nrvt, vt);
7874 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, n, n, isupper,
w, e);
7875 if( additionalmemory<1 )
7881 bidiagonal::rmatrixbdmultiplybyq<Precision>(a, n, n, tauq, u,
m, n,
true,
false);
7882 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, n, isupper,
false, u,
m, a, 0, vt, ncvt);
7891 bidiagonal::rmatrixbdunpackq<Precision>(a, n, n, tauq, n, t2);
7892 blas::copymatrix<Precision>(u, 0, m-1, 0, n-1, a, 0, m-1, 0, n-1);
7893 blas::inplacetranspose<Precision>(t2, 0, n-1, 0, n-1, work);
7894 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, n, isupper,
false, u, 0, t2, n, vt, ncvt);
7895 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);
7913 lq::rmatrixlq<Precision>(a,
m, n,
tau);
7914 for(i=0; i<=m-1; i++)
7916 for(j=i+1; j<=m-1; j++)
7921 bidiagonal::rmatrixbd<Precision>(a,
m,
m, tauq, taup);
7922 bidiagonal::rmatrixbdunpackq<Precision>(a,
m,
m, tauq, ncu, u);
7923 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a,
m,
m, isupper,
w, e);
7925 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7926 result = bdsvd::rmatrixbdsvd<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, 0);
7927 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7936 lq::rmatrixlq<Precision>(a,
m, n,
tau);
7937 lq::rmatrixlqunpackq<Precision>(a,
m, n,
tau, nrvt, vt);
7938 for(i=0; i<=m-1; i++)
7940 for(j=i+1; j<=m-1; j++)
7945 bidiagonal::rmatrixbd<Precision>(a,
m,
m, tauq, taup);
7946 bidiagonal::rmatrixbdunpackq<Precision>(a,
m,
m, tauq, ncu, u);
7947 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a,
m,
m, isupper,
w, e);
7949 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7950 if( additionalmemory<1 )
7956 bidiagonal::rmatrixbdmultiplybyp<Precision>(a,
m,
m, taup, vt,
m, n,
false,
true);
7957 result = bdsvd::rmatrixbdsvd<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, n);
7965 bidiagonal::rmatrixbdunpackpt<Precision>(a,
m,
m, taup,
m, t2);
7966 result = bdsvd::rmatrixbdsvd<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, t2,
m);
7967 blas::copymatrix<Precision>(vt, 0, m-1, 0, n-1, a, 0, m-1, 0, n-1);
7968 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);
7970 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7981 bidiagonal::rmatrixbd<Precision>(a,
m, n, tauq, taup);
7982 bidiagonal::rmatrixbdunpackq<Precision>(a,
m, n, tauq, ncu, u);
7983 bidiagonal::rmatrixbdunpackpt<Precision>(a,
m, n, taup, nrvt, vt);
7984 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a,
m, n, isupper,
w, e);
7986 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7987 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, minmn, isupper,
false, a, 0, u, nru, vt, ncvt);
7988 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
7995 bidiagonal::rmatrixbd<Precision>(a,
m, n, tauq, taup);
7996 bidiagonal::rmatrixbdunpackq<Precision>(a,
m, n, tauq, ncu, u);
7997 bidiagonal::rmatrixbdunpackpt<Precision>(a,
m, n, taup, nrvt, vt);
7998 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a,
m, n, isupper,
w, e);
7999 if( additionalmemory<2 || uneeded==0 )
8005 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, minmn, isupper,
false, u, nru, a, 0, vt, ncvt);
8014 blas::copyandtranspose<Precision>(u, 0, m-1, 0, minmn-1, t2, 0, minmn-1, 0, m-1);
8015 result = bdsvd::rmatrixbdsvd<Precision>(
w, e, minmn, isupper,
false, u, 0, t2,
m, vt, ncvt);
8016 blas::copyandtranspose<Precision>(t2, 0, minmn-1, 0, m-1, u, 0, m-1, 0, minmn-1);
8026 template<
unsigned int Precision>
8032 int additionalmemory,
8069 w.setbounds(1, minmn);
8076 u.setbounds(1, nru, 1, ncu);
8082 u.setbounds(1, nru, 1, ncu);
8090 vt.setbounds(1, nrvt, 1, ncvt);
8096 vt.setbounds(1, nrvt, 1, ncvt);
8111 qr::qrdecomposition<Precision>(a,
m, n,
tau);
8114 for(j=1; j<=i-1; j++)
8119 bidiagonal::tobidiagonal<Precision>(a, n, n, tauq, taup);
8120 bidiagonal::unpackptfrombidiagonal<Precision>(a, n, n, taup, nrvt, vt);
8121 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, n, n, isupper,
w, e);
8122 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, n, isupper,
false, u, 0, a, 0, vt, ncvt);
8131 qr::qrdecomposition<Precision>(a,
m, n,
tau);
8132 qr::unpackqfromqr<Precision>(a,
m, n,
tau, ncu, u);
8135 for(j=1; j<=i-1; j++)
8140 bidiagonal::tobidiagonal<Precision>(a, n, n, tauq, taup);
8141 bidiagonal::unpackptfrombidiagonal<Precision>(a, n, n, taup, nrvt, vt);
8142 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, n, n, isupper,
w, e);
8143 if( additionalmemory<1 )
8149 bidiagonal::multiplybyqfrombidiagonal<Precision>(a, n, n, tauq, u,
m, n,
true,
false);
8150 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, n, isupper,
false, u,
m, a, 0, vt, ncvt);
8159 bidiagonal::unpackqfrombidiagonal<Precision>(a, n, n, tauq, n, t2);
8160 blas::copymatrix<Precision>(u, 1,
m, 1, n, a, 1,
m, 1, n);
8161 blas::inplacetranspose<Precision>(t2, 1, n, 1, n, work);
8162 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, n, isupper,
false, u, 0, t2, n, vt, ncvt);
8163 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);
8181 lq::lqdecomposition<Precision>(a,
m, n,
tau);
8182 for(i=1; i<=m-1; i++)
8184 for(j=i+1; j<=
m; j++)
8189 bidiagonal::tobidiagonal<Precision>(a,
m,
m, tauq, taup);
8190 bidiagonal::unpackqfrombidiagonal<Precision>(a,
m,
m, tauq, ncu, u);
8191 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a,
m,
m, isupper,
w, e);
8193 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8194 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, 0);
8195 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8204 lq::lqdecomposition<Precision>(a,
m, n,
tau);
8205 lq::unpackqfromlq<Precision>(a,
m, n,
tau, nrvt, vt);
8206 for(i=1; i<=m-1; i++)
8208 for(j=i+1; j<=
m; j++)
8213 bidiagonal::tobidiagonal<Precision>(a,
m,
m, tauq, taup);
8214 bidiagonal::unpackqfrombidiagonal<Precision>(a,
m,
m, tauq, ncu, u);
8215 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a,
m,
m, isupper,
w, e);
8217 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8218 if( additionalmemory<1 )
8224 bidiagonal::multiplybypfrombidiagonal<Precision>(a,
m,
m, taup, vt,
m, n,
false,
true);
8225 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, n);
8233 bidiagonal::unpackptfrombidiagonal<Precision>(a,
m,
m, taup,
m, t2);
8234 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, t2,
m);
8235 blas::copymatrix<Precision>(vt, 1,
m, 1, n, a, 1,
m, 1, n);
8236 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);
8238 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8249 bidiagonal::tobidiagonal<Precision>(a,
m, n, tauq, taup);
8250 bidiagonal::unpackqfrombidiagonal<Precision>(a,
m, n, tauq, ncu, u);
8251 bidiagonal::unpackptfrombidiagonal<Precision>(a,
m, n, taup, nrvt, vt);
8252 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a,
m, n, isupper,
w, e);
8254 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8255 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, minmn, isupper,
false, a, 0, u, nru, vt, ncvt);
8256 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
8263 bidiagonal::tobidiagonal<Precision>(a,
m, n, tauq, taup);
8264 bidiagonal::unpackqfrombidiagonal<Precision>(a,
m, n, tauq, ncu, u);
8265 bidiagonal::unpackptfrombidiagonal<Precision>(a,
m, n, taup, nrvt, vt);
8266 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a,
m, n, isupper,
w, e);
8267 if( additionalmemory<2 || uneeded==0 )
8273 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, minmn, isupper,
false, u, nru, a, 0, vt, ncvt);
8282 blas::copyandtranspose<Precision>(u, 1,
m, 1, minmn, t2, 1, minmn, 1,
m);
8283 result = bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, minmn, isupper,
false, u, 0, t2,
m, vt, ncvt);
8284 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