1+ #ifndef CORE_H
2+ #define CORE_H
3+ #include < iostream>
4+ #include < iomanip>
5+ #include < cmath>
6+ #include < functional>
7+ #include < vector>
8+ using std::cerr;
9+ using std::cout;
10+ using std::function;
11+ using std::vector;
12+
13+ namespace cor {
14+
15+ // -----常量部分-----
16+ const double Machine_Epsilon = pow(2 , -52 );
17+
18+ // -----声明部分-----
19+
20+ // 杂例I
21+ template <class DM >
22+ double Norm (const vector<DM> &, double );
23+ template <class DJ > DJ Abs (DJ);
24+ template <class DJ > DJ Sign (DJ);
25+ template <class DJ > vector<DJ> Range (DJ, DJ,DJ);
26+ template <class DJ > vector<DJ> Range (DJ, DJ);
27+ template <class DJ > vector<DJ> Rand (int );
28+ template <class DJ > DJ Limit (std::function<DJ(DJ)>, DJ);
29+ template <class DJ > DJ N (DJ,int );
30+ inline int Sig (){return rand ()%2 ?1 :-1 ;}// random Sign(+/-)
31+
32+ template <class DD > DD Iteration (std::function<DD(DD)>, DD);
33+ template <class DD > DD Sqrt (DD);
34+ template <> double Sqrt (double );
35+
36+
37+ // -----定义部分-----
38+
39+
40+ // 杂例I
41+ template <class DM > double Norm (const vector<DM> &a, double p)
42+ {
43+ if (p==INFINITY)
44+ {
45+ double max = 0 ;
46+ for (long long unsigned int i = 0 ; i < a.size ();++i)
47+ {
48+ if (Abs (a[i])>max)
49+ max = Abs (a[i]);
50+ }
51+ return max;
52+ }
53+ double s = 0 ;
54+ for (long long unsigned int i = 0 ; i < a.size ();++i)
55+ {
56+ s += pow (Abs (a[i]), p);
57+ }
58+ s = pow (s, 1 / p);
59+ return s;
60+ }
61+ template <class DJ > DJ Abs (DJ a)
62+ {
63+ if (a>0 )
64+ return a;
65+ return -a;
66+ }
67+ template <class DJ > DJ Sign (DJ a)
68+ {
69+ if (a>0 )
70+ return 1 ;
71+ if (a<0 )
72+ return -1 ;
73+ return 0 ;
74+ }
75+ template <class DJ > vector<DJ> Range (DJ a, DJ b,DJ c)
76+ {
77+ vector<DJ> r;
78+ DJ s = a;
79+ while (s<=b)
80+ {
81+ r.emplace_back (s);
82+ s = s + DJ (c);
83+ }
84+ return r;
85+ }
86+ template <class DJ > vector<DJ> Range (DJ a, DJ b)
87+ {
88+ return Range (a, b, DJ (1 ));
89+ }
90+ template <class DJ > vector<DJ> Rand (int a)
91+ {
92+ vector<DJ> r (a);
93+ for (int i = 0 ;i<a;++i)
94+ r[i]=std::rand ()*Sig ();
95+ return r;
96+ }
97+ template <class DJ > DJ Limit (std::function<DJ(DJ)> f, DJ x0)
98+ {
99+ const DJ h = 1e-7 ;
100+ return 0.5 * (f (x0+h)+f (x0-h));
101+ }
102+ template <class DJ > DJ N (DJ a,int n)
103+ {
104+ cout << std::setprecision (n) << a << ' \n ' ;
105+ return a;
106+ }
107+
108+ template <class DD > DD Iteration (std::function<DD(DD)> phi, DD a)
109+ {
110+ int times = 100 ;// 迭代次数
111+ int k = 0 ;
112+ DD pre =a;
113+ DD x = phi (a);
114+ while (Abs (x-pre )>1e-14 && k < times)
115+ {
116+ pre = x;
117+ x = phi (x);
118+ ++k;
119+ }
120+ if (k>=times)
121+ {
122+ cerr<<" 警告:未在指定次数内达到预期精度。" <<' \n ' ;
123+ cerr << " 可能是初值的选择使迭代法不收敛,或者迭代函数的选择使收敛速度很慢。" << ' \n ' ;
124+ }
125+ return x;
126+ }
127+ template <class DD > DD Sqrt (DD c)
128+ {
129+ std::function<DD (DD)> phi = [&](DD x) { return (x + c / x) / 2 ; };
130+ return Iteration (phi,c);
131+ }
132+ template <> double Sqrt<double >(double c)
133+ {
134+ if (c<0 )
135+ {
136+ cerr << " 错误:根号下不能是负数。" << ' \n ' ;
137+ return NAN;
138+ }
139+ double pre =0 ;
140+ double x = 0.5 *c;
141+ while (fabs (pre -x)>1e-14 )
142+ {
143+ pre = x;
144+ x = 0.5 * (x + c / x);
145+ }
146+ return x;
147+ }
148+
149+ }
150+ #endif
0 commit comments