求歐拉范數是一個比較簡單的算法，似乎沒有什么可說的，一般代碼如下：
/** \fn double enorm1(long n, const double* x)
* \brief 求歐拉范數，簡單算法。
* \param [in] n 向量長度
* \param [in] x 向量值
* \return 歐拉范數
*/
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時，發現它求歐拉范數的函數enorm則要復雜得多。仔細比較發現上面算法有不足，它沒有考慮溢出的情況，當x[i]為很小或者很大的數時，x[i]*x[i]則會下溢或者上溢，最后結果可能不準確。但x[i]*x[i]只是中間結果，最后的范數數量級應該和x[i]相同。它溢出的可能性要小得多。改進中間結果則可以改進算法。
/** \fn double enorm2(long n, const double* x)
* \brief 求歐拉范數，考慮溢出的算法。
* \param [in] n 向量長度
* \param [in] x 向量值
* \return 歐拉范數
*/
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)，降低了溢出的可能。當然該算法沒有區分可能溢出和不溢出的數，計算量較大。下面的算法是仿照minpack的enorm函數編寫：
/** \fn double enorm3(long n, const double* x)
* \brief 求歐拉范數，仿照minpack的enorm函數
* \param [in] n 向量長度
* \param [in] x 向量值
* \return 歐拉范數
*/
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
導航
|
|
| 27 | 28 | 29 | 30 | 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 | 6 | 7 |
常用鏈接
留言簿(4)
隨筆檔案
文章檔案
搜索
最新隨筆
最新評論
閱讀排行榜
評論排行榜
久久久久九九精品影院|
国产精品久久影院|
久久成人小视频|
久久婷婷国产综合精品|
99久久久精品|
香港aa三级久久三级老师2021国产三级精品三级在
|
伊人久久大香线蕉精品|
久久se这里只有精品|
亚洲精品无码久久久影院相关影片|
久久国产精品99精品国产|
久久精品国产亚洲av瑜伽|
熟妇人妻久久中文字幕|
久久国产综合精品五月天|
无码国产69精品久久久久网站|
国内精品久久久久影院网站
|
久久久精品久久久久影院|
精品久久久久久中文字幕|
久久笫一福利免费导航
|
色狠狠久久综合网|
99久久精品无码一区二区毛片
|
欧美久久精品一级c片片|
久久久久久午夜精品|
国产三级精品久久|
97久久超碰国产精品旧版|
综合久久国产九一剧情麻豆
|
国产精自产拍久久久久久蜜|
2021少妇久久久久久久久久|
久久综合亚洲色HEZYO社区
|
亚洲中文字幕无码久久综合网|
国产—久久香蕉国产线看观看|
久久精品国产亚洲麻豆|
人妻精品久久无码区|
人妻丰满AV无码久久不卡|
亚洲第一极品精品无码久久|
久久久无码精品亚洲日韩京东传媒
|
热久久视久久精品18|
欧洲性大片xxxxx久久久|
亚洲午夜福利精品久久|
亚洲国产天堂久久综合|
色狠狠久久综合网|
亚洲va中文字幕无码久久不卡|