求歐拉范數(shù)是一個比較簡單的算法，似乎沒有什么可說的，一般代碼如下：
/** \fn double enorm1(long n, const double* x)
* \brief 求歐拉范數(shù)，簡單算法。
* \param [in] n 向量長度
* \param [in] x 向量值
* \return 歐拉范數(shù)
*/
double enorm1(long n, const double* x)
{
double ret = 0.0;
long i = 0;
for (i=0; i<n; ++i)
ret += x[i]*x[i];
return sqrt(ret);
}
然而近日在學習minpack時，發(fā)現(xiàn)它求歐拉范數(shù)的函數(shù)enorm則要復雜得多。仔細比較發(fā)現(xiàn)上面算法有不足，它沒有考慮溢出的情況，當x[i]為很小或者很大的數(shù)時，x[i]*x[i]則會下溢或者上溢，最后結果可能不準確。但x[i]*x[i]只是中間結果，最后的范數(shù)數(shù)量級應該和x[i]相同。它溢出的可能性要小得多。改進中間結果則可以改進算法。
/** \fn double enorm2(long n, const double* x)
* \brief 求歐拉范數(shù)，考慮溢出的算法。
* \param [in] n 向量長度
* \param [in] x 向量值
* \return 歐拉范數(shù)
*/
double enorm2(long n, const double *x)
{
double ret = 0.0;
long i= 0;
double xmax = 0.0;
for (i=0; i<n; ++i)
{
double xabs = fabs(x[i]);
/* 這個比較方式不需要考慮xabs和xmax為0的情況 */
if (xabs < xmax)
{
ret += (xabs/xmax)*(xabs/xmax);
}
else if (xabs == xmax)
{
ret += 1;
}
else
{
ret = 1 + ret*(xmax/xabs)*(xmax/xabs);
xmax = xabs;
}
}
return sqrt(ret) * xmax;
}
將中間結果改成(x[i]/xmax)*(x[i]/xmax)，降低了溢出的可能。當然該算法沒有區(qū)分可能溢出和不溢出的數(shù)，計算量較大。下面的算法是仿照minpack的enorm函數(shù)編寫：
/** \fn double enorm3(long n, const double* x)
* \brief 求歐拉范數(shù)，仿照minpack的enorm函數(shù)
* \param [in] n 向量長度
* \param [in] x 向量值
* \return 歐拉范數(shù)
*/
double enorm3(long n, const double *x)
{
double s1 = 0.0;
double s2 = 0.0;
double s3 = 0.0;
/* 上溢和下溢的邊界，不一定要十分精確 */
const double dwarf = 1.483e-154; /* 下溢的邊界 */
const double giant = 1.341e154 / n; /* 上溢的邊界 */
double x1max = 0.0;
double x3max = 0.0;
long i = 0;
for (i=0; i<n; ++i)
{
double xabs = fabs(x[i]);
if (xabs > dwarf && xabs < giant)
{
s2 += xabs*xabs;
}
else if (xabs <= dwarf)
{
/* 這個比較方式不需要考慮xabs和xmax為0的情況 */
if (xabs < x3max)
{
s3 += (xabs/x3max)*(xabs/x3max);
}
else if (xabs == x3max)
{
s3 += 1;
}
else
{
s3 =1 + s3*(x3max/xabs)*(x3max/xabs);
x3max = xabs;
}
}
else /* if (xabs >= giant) */
{
/* 不需要考慮xabs和xmax為0的情況 */
if (xabs <= x1max)
{
s1 += (xabs/x1max)*(xabs/x1max);
}
else if (xabs == x1max)
{
s1 += 1;
}
else
{
s1 =1 + s1*(x1max/xabs)*(x1max/xabs);
x1max = xabs;
}
}
}
if (s1 != 0.0)
{
return x1max*sqrt(s1 + s2/x1max/x1max);
}
else if (s2 != 0.0)
{
return sqrt(s2 + x3max*x3max*s3);
/* 下面為minpack中enorm的代碼，好像沒有必要
if (s2 >= x3max)
return sqrt(s2 * (1 + x3max/s2*x3max*s3));
else
return sqrt(x3max * (s2/x3max + x3max*s2));
*/
}
else
{
return x3max * sqrt(s3);
}
}
下面是測試結果：
|
向量
|
enorm1
|
enorm2
|
enorm3
|
|
{1, 2, 3}
|
3.74166
|
3.74166
|
3.74166
|
|
{1e200, 2e200, 3e200}
|
1.#INF
|
<td style="BORDER-RIGHT: rgb(0,0,0) 0.5pt solid; PADDING-RIGHT: 5.4pt; BORDER-TOP: medium none; PADDING-LEFT: 5.4pt; PADDING-BOTTOM: 0pt; BO
導航
|
|
| 26 | 27 | 28 | 29 | 30 | 31 | 1 |
| 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| 9 | 10 | 11 | 12 | 13 | 14 | 15 |
| 16 | 17 | 18 | 19 | 20 | 21 | 22 |
| 23 | 24 | 25 | 26 | 27 | 28 | 29 |
| 30 | 31 | 1 | 2 | 3 | 4 | 5 |
常用鏈接
留言簿(4)
隨筆檔案
文章檔案
搜索
最新隨筆
最新評論
閱讀排行榜
評論排行榜
久久久久久国产精品免费无码|
狠狠色婷婷综合天天久久丁香|
欧美大战日韩91综合一区婷婷久久青草|
四虎国产永久免费久久|
欧美激情精品久久久久久久九九九|
久久这里只精品99re66|
99久久免费国产精品热|
看全色黄大色大片免费久久久
|
久久精品国产一区二区|
久久青青草视频|
免费国产99久久久香蕉|
久久SE精品一区二区|
岛国搬运www久久|
色欲久久久天天天综合网|
国产99久久久久久免费看|
亚洲人成伊人成综合网久久久
|
亚洲综合久久久|
国产精品久久久久久福利69堂|
无码任你躁久久久久久老妇|
久久se精品一区精品二区|
伊人久久综合无码成人网|
久久播电影网|
97超级碰碰碰碰久久久久|
国产精品美女久久久久|
亚洲Av无码国产情品久久|
99久久亚洲综合精品成人|
99国产精品久久久久久久成人热|
久久热这里只有精品在线观看|
久久精品无码一区二区日韩AV|
色综合久久最新中文字幕|
久久精品aⅴ无码中文字字幕重口|
国产精品一区二区久久精品涩爱|
久久综合伊人77777麻豆|
欧美大战日韩91综合一区婷婷久久青草
|
久久99精品国产麻豆蜜芽|
91精品国产高清久久久久久91
|
丰满少妇人妻久久久久久|
99久久国语露脸精品国产|
欧美一区二区三区久久综|
久久国产精品无码HDAV|
国内精品久久久久久99蜜桃|