Advertisement
VNM24ix

Untitled

Jun 9th, 2025
11
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 23.78 KB | Source Code | 0 0
  1. using System;
  2. using System.Linq;
  3. using ComplexMath;
  4. using CSharpShellCore;
  5. using System.Numerics;
  6. using System.Collections.Generic;
  7. using System.IO;
  8.  
  9. namespace ComplexMath {
  10. delegate double func(double t);
  11. delegate Cplx cplxFn(Cplx z);
  12. class Polyn {
  13. Cplx[] a;
  14. public int L => a.Length;
  15. public int N => a.Length - 1;
  16. public Cplx this[int i] { get { return a[i]; } set { a[i] = value; } }
  17. public Polyn(params Cplx[] a) {
  18. this.a = CopyNoZeros(Cplx.ArrayCopy(a));
  19.  
  20. }
  21. public Polyn(Polyn p) : this(p.a) { }
  22. public Polyn() : this(new Cplx[] { }) { }
  23. public static Polyn Inv(params Cplx[] a) {
  24. Cplx[] aa = new Cplx[a.Length];
  25. for (int i = 0; i < a.Length; i++)
  26. aa[i] = a[a.Length - 1 - i];
  27. return new Polyn(aa);
  28. }
  29. public static Polyn Zero => new Polyn() { };
  30.  
  31. static Cplx[] CopyNoZeros(Cplx[] a) {
  32. int s = 0;
  33. for (int i = a.Length - 1; i >= 0; i--)
  34. if (a[i].IsZero) s++;
  35. else break;
  36. Cplx[] aa = new Cplx[a.Length - s];
  37. for (int i = 0; i < aa.Length; i++) aa[i] = a[i].Copy;
  38. return aa;
  39. }
  40. public void Show(string name, int total = 16, int frac = 12) {
  41. Console.WriteLine(name + ":");
  42. for (int i = 0; i < L; i++) {
  43. a[i].Show(total, frac, true, false);
  44. Console.Write(" * X" + i + "\n");
  45. }
  46. Console.Write("\n");
  47. }
  48. public void ShowInv(string name, int total = 16, int frac = 12) {
  49. Polyn.Inv(a).Show(name, total, frac);
  50. }
  51. public Polyn Deriv {
  52. get {
  53. if (a.Length < 2) return new Polyn(new Cplx[0]);
  54. Cplx[] aa = new Cplx[a.Length - 1];
  55. for (int i = 0; i < aa.Length; i++) {
  56. aa[i] = a[i + 1] * (i + 1);
  57. }
  58. return new Polyn(aa);
  59. }
  60. }
  61. public Cplx Val(Cplx z) {
  62. Cplx w = 0;
  63. for (int i = a.Length - 1; i >= 0; i--) {
  64. w *= z; w += a[i];
  65. }
  66. return w;
  67. }
  68. public static Polyn operator -(Polyn p) {
  69. Polyn q = new Polyn(p.a);
  70. for (int i = 0; i < q.L; i++) q[i] = -q[i];
  71. return q;
  72. }
  73. public static Polyn operator +(Polyn p, Polyn q) {
  74. if (p.L == 0) return new Polyn(q.a);
  75. else if (q.L == 0) return new Polyn(p.a);
  76. int l = Math.Max(p.L, q.L);
  77. //Console.WriteLine(l);
  78. Cplx[] a = new Cplx[l];
  79. for (int i = 0; i < l; i++) {
  80. a[i] = (i < p.L ? p[i] : 0) + (i < q.L ? q[i] : 0);
  81. }
  82. return new Polyn(a);
  83. }
  84. public static Polyn operator -(Polyn p, Polyn q) => p + -q;
  85. public static Polyn operator *(Polyn p, Cplx λ) {
  86. if (λ.IsZero) return new Polyn(new Cplx[0]);
  87. Polyn q = new Polyn(p.a);
  88. for (int i = 0; i < q.L; i++) q[i] *= λ;
  89. return q;
  90. }
  91. public static Polyn operator *(Cplx λ, Polyn p) => p * λ;
  92. public static Polyn operator *(Polyn p, Polyn q) {
  93. if (p.L == 0 || q.L == 0) return Zero;
  94. //Console.WriteLine(p.L+" "+q.L+" "+p.N+" "+ q.N);
  95. Cplx[] a = new Cplx[p.L + q.L - 1];
  96. for (int i = 0; i < a.Length; i++) a[i] = 0;
  97. for (int i = 0; i < p.L; i++)
  98. for (int j = 0; j < q.L; j++) {
  99. //Console.WriteLine(i+j-1);
  100. a[i + j] += p[i] * q[j];
  101. }
  102. return new Polyn(a);
  103. }
  104. public Polyn Pow(int n) {
  105. if (n == 0) return new Polyn(new Cplx[] { 1 });
  106. Polyn pow = new Polyn(a);
  107. for (int i = 2; i <= n; i++) pow *= this;
  108. return pow;
  109. }
  110. //суперпозиция полиномов
  111. public Polyn SP(Polyn q) {
  112. Polyn res = Polyn.Zero;
  113. for (int i = N; i >= 0; i--) {
  114. res *= q;
  115. res += new Polyn(new Cplx[] { a[i] });
  116. }
  117. return res;
  118. }
  119. public static Polyn[] Div(Polyn p, Polyn q) {
  120. //Console.WriteLine("In Polyn.Div"); p.Show("p");
  121. //Console.WriteLine("In Polyn.Div"); q.Show("q");
  122. if (q.L == 0) return null;
  123. else if (p.L < q.L) return new Polyn[] { Polyn.Zero, new Polyn(p) };
  124. Cplx[] pa = Cplx.ArrayCopy(p.a);
  125. Cplx[] qa = Cplx.ArrayCopy(q.a);
  126. Cplx[] u = new Cplx[p.L - q.L + 1];
  127. //int k = 0;
  128. for (int n = p.N, k = 0; n >= q.N; n--, k++) {
  129. Cplx g = pa[n] / qa[q.N];
  130. u[u.Length - 1 - k] = g.Copy;
  131. for (int i = 0; i < q.L; i++) {
  132. //Console.WriteLine((n - k - i) + " " + (q.N - i));
  133. pa[n - i] -= qa[q.N - i] * g;
  134. }
  135. }
  136. //Console.WriteLine("In Polyn.Div"); new Polyn(u).Show("u");
  137. //Console.WriteLine("In Polyn.Div"); new Polyn(pa).Show("pa");
  138. return new Polyn[] { new Polyn(u), new Polyn(pa) };
  139. }
  140. public static Polyn GCD(Polyn p, Polyn q) {
  141. //Console.WriteLine("In GCD:");
  142. Polyn pp = new Polyn(p);
  143. Polyn qq = new Polyn(q);
  144. if (pp.L == 1 || qq.L == 1)
  145. return new Polyn(new Cplx[] { 1 });
  146. else {
  147. Polyn r = Div(pp, qq)[1];
  148. //Console.WriteLine("Div(pp,qq)[1]:");
  149. //r.Show();
  150. if (r.L == 0) return qq;
  151. else {
  152. pp = new Polyn(qq);
  153. qq = r;
  154. return GCD(pp, qq);
  155. }
  156. }
  157. }
  158. public Cplx[] Roots0() {
  159. //Console.WriteLine("In Polyn Roots0: Start");
  160. Cplx[] r = new Cplx[N];
  161. Polyn q = new Polyn(this);
  162. for (int i = 0; i < N; i++) {
  163. r[i] = Cplx.Newton(q, new Cplx(2, 1));
  164. //Console.WriteLine($"In Polyn.Roots0: Found root{i}"); r[i].Show();
  165. Polyn[] div = Div(q, new Polyn(new Cplx[] { -r[i], 1 }));
  166. //Console.WriteLine("In Polyn.Div"); div[0].Show("div[0]");
  167. q = div[0];
  168. //Console.WriteLine("In Polyn Roots0"); q.Show("Polyn q");
  169. }
  170. //Console.WriteLine("In Polyn Roots0 OK");
  171. return r;
  172. }
  173. public Cplx[] Roots() {
  174. Polyn q = GCD(this, Deriv);
  175. //Console.WriteLine("In Polyn.Roots");
  176. //q.Show("GCD");
  177. if (q.L == 1) return Roots0();
  178. else {
  179. Polyn p = Div(this, q)[0];
  180. //p.Show(15, 9);
  181. Cplx[] r0 = p.Roots0();
  182. Cplx[] r1 = q.Roots();
  183. Cplx[] r = new Cplx[r0.Length + r1.Length];
  184. for (int i = 0; i < r.Length; i++) {
  185. r[i] = (i < r0.Length ? r0[i] : r1[i - r0.Length]).Copy;
  186. }
  187. return r;
  188. }
  189. }
  190. }
  191.  
  192. class Cplx {
  193. //delegate double DD(double t);
  194. static Random rnd = new Random();
  195. public const double π = Math.PI;
  196. const double ε = 1e-6;
  197. func sqr = t => t * t;
  198. func sqrt = Math.Sqrt;
  199. func sin = Math.Sin;
  200. func cos = Math.Cos;
  201. func cbrt = t => Math.Pow(t, (double)1 / 3);
  202. func exp = t => Math.Exp(t);
  203. func asin = Math.Asin;
  204. func acos = Math.Acos;
  205. func ln = Math.Log;
  206. static double pow(double a, double b) => Math.Pow(a, b);
  207. public double X { get; private set; }
  208. public double Y { get; private set; }
  209. public Cplx(double x, double y) {
  210. X = x;
  211. Y = y;
  212. }
  213. public Cplx Copy => new Cplx(X, Y);
  214. public bool IsZero => R < ε;
  215. public bool isReal => Math.Abs(Y) < ε;
  216. public static implicit operator Cplx(double x) => new Cplx(x, 0);
  217. public static Cplx operator +(Cplx z, Cplx w) {
  218. double x = z.X + w.X;
  219. double y = z.Y + w.Y;
  220. return new Cplx(x, y);
  221. }
  222. public static Cplx I = new Cplx(0, 1);
  223. public Cplx Iz => new Cplx(-Y, X);
  224. public double R2 => sqr(X) + sqr(Y);
  225. public double R => sqrt(R2);
  226. public double Φ {
  227. get {
  228. double φ = acos(X / R);
  229. return Y < 0 ? -φ : φ;
  230. }
  231. }
  232. public static Cplx operator -(Cplx z) => new Cplx(-z.X, -z.Y);
  233. public static Cplx operator -(Cplx z, Cplx w) => z + -w;
  234. public static Cplx operator *(Cplx z, double q) => new Cplx(z.X * q, z.Y * q);
  235. public static Cplx operator *(double q, Cplx z) => z * q;
  236. public static Cplx operator *(Cplx z, Cplx w) {
  237. double x = z.X * w.X - z.Y * w.Y;
  238. double y = z.X * w.Y + z.Y * w.X;
  239. return new Cplx(x, y);
  240. }
  241. public static Cplx operator ~(Cplx z) => new Cplx(z.X, -z.Y);
  242. public static Cplx operator /(Cplx z, double q) => new Cplx(z.X / q, z.Y / q);
  243. public static Cplx operator /(Cplx z, Cplx w) => z * ~w / w.R2;
  244. public Cplx PowN(int n) {
  245. if (n > 0) {
  246. if (n == 1) return this.Copy;
  247. else {
  248. int nn = n >> 1;
  249. Cplx z = (PowN(nn)), zz = z * z;
  250. return (n & 1) == 1 ? zz * this : zz;
  251. }
  252. } else return n == 0 ? 1 : 1 / PowN(-n);
  253. }
  254. public Cplx Root(int n) {
  255. double r = pow(R, (double)1 / n);
  256. double φ = Φ / n;
  257. return new Cplx(r * cos(φ), r * sin(φ));
  258. }
  259. public Cplx[] Cbrt {
  260. get {
  261. double eps = 1e-9;
  262. Cplx w = new Cplx(5, 4);
  263. bool f = false;
  264. while (true) {
  265. Cplx ww = (this / (w * w) + w + w) / 3;
  266. if (f) { w = ww; break; } else {
  267. f = ((w - ww) / ww).R < eps;
  268. w = ww.Copy;
  269. }
  270. }
  271. Cplx[] res = new Cplx[3];
  272. res[0] = w;
  273. Cplx q = new Cplx(-1, sqrt(3)) / 2;
  274. res[1] = w * q;
  275. res[2] = w * ~q;
  276. return res;
  277. }
  278. }
  279. public Cplx Exp {
  280. get {
  281. double r = exp(X);
  282. return r * new Cplx(cos(Y), sin(Y));
  283. }
  284. }
  285. public Cplx Cos => (Iz.Exp + (-Iz).Exp) / 2;
  286. public Cplx Sin => ((-Iz).Exp - Iz.Exp).Iz / 2;
  287. public Cplx Tg => Sin / Cos;
  288. public Cplx Ln => new Cplx(ln(R), Φ);
  289. public Cplx Pow(Cplx z) => (z * Ln).Exp;
  290. public static Cplx Newton(Polyn p, Cplx w0) {
  291. //Console.WriteLine("In Newton"); p.Show("p");
  292. w0 = new Cplx(0.5 + rnd.NextDouble(), 0.5 + rnd.NextDouble());
  293. Cplx z = w0.Copy;
  294. bool flag = false;
  295. double ε1 = 1e-2, ε2 = 1e-4;
  296. while (true) {
  297. Cplx δ = p.Val(z) / p.Deriv.Val(z);
  298. double λ = δ.R > ε1 ? 0.25 : (δ.R > ε2 ? 0.5 : 1);
  299. Cplx w = z - δ * λ;
  300. if (flag) return w;
  301. flag = (w - z).R < 1e-9;
  302. z = w.Copy;
  303. //Console.WriteLine("In Cplx.Newton"); z.Show(15, 9);
  304. }
  305. }
  306. public static cplxFn D1(cplxFn f) {
  307. double ε = 1e-5;
  308. double φ = π * (2 * rnd.NextDouble() - 1) / 2;
  309. Cplx δ = ε * (Math.Cos(φ) + Cplx.I * Math.Sin(φ));
  310. return z => (f(z + δ / 2) - f(z - δ / 2)) / δ;
  311. }
  312.  
  313. public static Cplx Newton(cplxFn f, Cplx z0) {
  314. cplxFn df = D1(f);
  315. Cplx z = z0.Copy;
  316. bool flag1 = false, flag2 = false;
  317. int cnt = 0;
  318. double k = 0.5;
  319. while (true) {
  320. Cplx zz = z - k * f(z) / df(z);
  321. if (flag2) return zz;
  322. else {
  323. double δ = (zz - z).R;
  324. flag1 = δ < 1e-3;
  325. if (flag1) k = 1;
  326. flag2 = δ < 1e-9;
  327. z = zz.Copy;
  328. cnt++;
  329. if (cnt > 50) return null;
  330. }
  331. }
  332. }
  333.  
  334. //public string Str(int total, int frac, )
  335. public void Show1(int total = 15, int frac = 9) {
  336. string s = "";
  337. for (int i = 0; i < 2; i++) {
  338. s += "{" + i + "," + total + ":f" + frac + "}";
  339. }
  340. s = "(" + s + ")\n";
  341. Console.ForegroundColor = isReal ? ConsoleColor.Red : ConsoleColor.Cyan;
  342. Console.Write(s, X, Y);
  343. Console.ForegroundColor = ConsoleColor.White;
  344. }
  345. public void Show(int total = 15, int frac = 9, bool color = true, bool newLine = true) {
  346. Console.ForegroundColor = ConsoleColor.White;
  347. Console.Write("(");
  348. Console.ForegroundColor = ConsoleColor.Red;
  349. string frm = "{0," + total + ":f" + frac + "}";
  350. Console.Write(frm, X);
  351. Console.ForegroundColor = ConsoleColor.White;
  352. Console.Write(" : ");
  353. Console.ForegroundColor = isReal ? ConsoleColor.DarkGray : ConsoleColor.Cyan;
  354. Console.Write(frm, Y);
  355. Console.ForegroundColor = ConsoleColor.White;
  356. Console.Write(")");
  357. if (newLine) Console.Write("\n");
  358. }
  359. public static void Show(Cplx[] z, int total, int frac) {
  360. for (int i = 0; i < z.Length; i++) {
  361. z[i].Show(total, frac);
  362. }
  363. Console.Write("\n");
  364. }
  365. public static void Show(Cplx[,] a, int total, int frac, bool real = false) {
  366. string form = "{0," + total + ":f" + frac + "}";
  367. for (int i = 0; i < a.GetLength(0); i++) {
  368. Console.ForegroundColor = ConsoleColor.Red;
  369. for (int j = 0; j < a.GetLength(1); j++) Console.Write(form, a[i, j].X);
  370. Console.WriteLine();
  371. if (!real) {
  372. for (int j = 0; j < a.GetLength(1); j++) {
  373. Console.ForegroundColor = a[i, j].isReal ? ConsoleColor.DarkGray : ConsoleColor.Cyan;
  374. Console.Write(form, a[i, j].Y);
  375. }
  376. Console.WriteLine("\n");
  377. }
  378. Console.ForegroundColor = ConsoleColor.White;
  379. }
  380. if (real) Console.WriteLine();
  381. }
  382. public static Cplx[] ArrayCopy(Cplx[] z) {
  383. Cplx[] w = new Cplx[z.Length];
  384. for (int i = 0; i < z.Length; i++) w[i] = z[i].Copy;
  385. return w;
  386. }
  387. public static Cplx[,] ArrayCopy(Cplx[,] z) {
  388. Cplx[,] w = new Cplx[z.GetLength(0), z.GetLength(1)];
  389. for (int i = 0; i < z.GetLength(0); i++)
  390. for (int j = 0; j < z.GetLength(1); j++) w[i, j] = z[i, j].Copy;
  391. return w;
  392. }
  393. public static Cplx[] GetRow(Cplx[,] a, int n) {
  394. Cplx[] b = new Cplx[a.GetLength(1)];
  395. for (int i = 0; i < b.Length; i++) b[i] = a[n, i].Copy;
  396. return b;
  397. }
  398. public static Cplx[] GetCol(Cplx[,] a, int n) {
  399. Cplx[] b = new Cplx[a.GetLength(0)];
  400. for (int i = 0; i < b.Length; i++) b[i] = a[i, n].Copy;
  401. return b;
  402. }
  403. public static void SetRow(Cplx[,] a, int n, Cplx[] b) {
  404. for (int i = 0; i < b.Length; i++) a[n, i] = b[i].Copy;
  405. }
  406.  
  407. public static void SetCol(Cplx[,] a, int n, Cplx[] b) {
  408. for (int i = 0; i < b.Length; i++) a[i, n] = b[i].Copy;
  409. }
  410. /*
  411. public static (Cplx[,] arr, int rnk) RankArr(Cplx[,] a0) {
  412. Cplx[,] a = Cplx.ArrayCopy(a0);
  413. int n0 = a.GetLength(0), n1 = a.GetLength(1);
  414. for (int i = 0; i < n0; i++)
  415. for (int j = 0; j < n1; j++)
  416. a[i, j] = a0[i, j].Copy;
  417. int cnt = 0;
  418. while (cnt < n0) {
  419. int j0 = -1;
  420. for (int j = 0; j < n1; j++) if (!a[cnt, j].IsZero) { j0 = j; break; }
  421. if (j0 == -1) {
  422. Cplx[,] tmp = new Cplx[n0 - 1, n1];
  423. for (int i = 0; i < n0 - 1; i++)
  424. for (int j = 0; j < n1; j++)
  425. tmp[i, j] = i < cnt ? a[i, j].Copy : a[i + 1, j].Copy;
  426. n0--;
  427. a = tmp;
  428. continue;
  429. }
  430. Cplx w = a[cnt, j0].Copy;
  431. for (int j = 0; j < n1; j++) a[cnt, j] /= w;
  432. for (int i = 0; i < n0; i++)
  433. if (i != cnt) {
  434. w = a[i, j0].Copy;
  435. for (int j = 0; j < n1; j++) a[i, j] -= a[cnt, j] * w;
  436. }
  437. cnt++;
  438. }
  439. return (a, Math.Min(a.GetLength(0), a.GetLength(1)));
  440. }
  441. */
  442. }
  443. class CplxMatrix {
  444. static Random rnd = new Random();
  445. public readonly Cplx[,] a;
  446. public int N { get; private set; }
  447. public Cplx this[int i, int j] { get { return a[i, j].Copy; } set { a[i, j] = value.Copy; } }
  448. public CplxMatrix(Cplx[,] a) {
  449. N = a.GetLength(0);
  450. this.a = new Cplx[N, N];
  451. for (int i = 0; i < N; i++)
  452. for (int j = 0; j < N; j++)
  453. this.a[i, j] = a[i, j].Copy;
  454. }
  455. public CplxMatrix Copy() => new CplxMatrix(a);
  456. public static CplxMatrix Zero(int n) {
  457. Cplx[,] a = new Cplx[n, n];
  458. for (int i = 0; i < n; i++)
  459. for (int j = 0; j < n; j++) a[i, j] = 0;
  460. return new CplxMatrix(a);
  461. }
  462. public bool IsZero() {
  463. foreach (Cplx z in a) if (!z.IsZero) return false;
  464. return true;
  465. }
  466. public static CplxMatrix I(int n) {
  467. Cplx[,] a = new Cplx[n, n];
  468. for (int i = 0; i < n; i++)
  469. for (int j = 0; j < n; j++) a[i, j] = (i == j ? 1 : 0);
  470. return new CplxMatrix(a);
  471. }
  472. public static CplxMatrix RandInt(int n, int m) {
  473. Cplx[,] a = new Cplx[n, n];
  474. for (int i = 0; i < n; i++)
  475. for (int j = 0; j < n; j++)
  476. a[i, j] = rnd.Next(-m, m + 1);
  477. return new CplxMatrix(a);
  478. }
  479. public static CplxMatrix RandErmit(int n, int m, bool real = false) {
  480. Cplx[,] a = new Cplx[n, n];
  481. for (int i = 0; i < n; i++) {
  482. int x = rnd.Next(-m, m + 1);
  483. a[i, i] = x;
  484. for (int j = i + 1; j < n; j++) {
  485. x = rnd.Next(-m, m + 1);
  486. int y = real ? 0 : rnd.Next(-m, m + 1);
  487. a[i, j] = new Cplx(x, y);
  488. a[j, i] = ~a[i, j];
  489. }
  490. }
  491. return new CplxMatrix(a);
  492. }
  493. public double Max() {
  494. double max = 0;
  495. foreach (Cplx z in a) if (z.R > max) max = z.R;
  496. return max;
  497. }
  498. public CplxMatrix Tr() {
  499. Cplx[,] a = new Cplx[N, N];
  500. for (int i = 0; i < N; i++)
  501. for (int j = 0; j < N; j++)
  502. a[i, j] = this.a[j, i].Copy;
  503. return new CplxMatrix(a);
  504. }
  505. public CplxMatrix Inv() {
  506. int n = N, nn = N + N;
  507. Cplx[,] a = new Cplx[n, nn];
  508. for (int i = 0; i < n; i++) {
  509. for (int j = 0; j < nn; j++)
  510. a[i, j] = j < n ? this[i, j] : j == N + i ? 1 : 0;
  511. }
  512. for (int i = 0; i < n; i++) {
  513. if (a[i, i].IsZero) {
  514. int ii = i;
  515. for (int j = i + 1; j < n; j++) if (!a[j, i].IsZero) {
  516. ii = j;
  517. break;
  518. }
  519. if (ii == i) return null;
  520. for (int k = i; k < nn; k++) { Cplx tmp = a[i, k].Copy; a[i, k] = a[ii, k].Copy; a[ii, k] = tmp; }
  521. }
  522. Cplx aii = a[i, i].Copy;
  523. for (int j = 0; j < nn; j++) a[i, j] /= aii;
  524. for (int j = 0; j < n; j++)
  525. if (j != i) {
  526. Cplx aji = a[j, i].Copy;
  527. for (int k = i; k < nn; k++) a[j, k] -= aji * a[i, k];
  528. }
  529. }
  530. Cplx[,] b = new Cplx[n, n];
  531. for (int i = 0; i < n; i++)
  532. for (int j = 0; j < n; j++)
  533. b[i, j] = a[i, n + j].Copy;
  534. return new CplxMatrix(b);
  535. }
  536. public static CplxMatrix operator +(CplxMatrix a, CplxMatrix b) {
  537. CplxMatrix res = new CplxMatrix(a.a);
  538. for (int i = 0; i < res.N; i++)
  539. for (int j = 0; j < res.N; j++)
  540. res[i, j] += b[i, j];
  541. return res;
  542. }
  543. public static CplxMatrix operator -(CplxMatrix a) {
  544. CplxMatrix res = new CplxMatrix(a.a);
  545. for (int i = 0; i < res.N; i++)
  546. for (int j = 0; j < res.N; j++)
  547. res[i, j] = -res[i, j];
  548. return res;
  549. }
  550. public static CplxMatrix operator -(CplxMatrix a, CplxMatrix b) => a + -b;
  551. public static CplxMatrix operator *(CplxMatrix a, Cplx μ) {
  552. CplxMatrix b = new CplxMatrix(a.a);
  553. for (int i = 0; i < b.N; i++)
  554. for (int j = 0; j < b.N; j++) b[i, j] *= μ;
  555. return b;
  556. }
  557. public static CplxMatrix operator *(CplxMatrix a, CplxMatrix b) {
  558. CplxMatrix res = Zero(a.N);
  559. for (int i = 0; i < a.N; i++)
  560. for (int j = 0; j < a.N; j++)
  561. for (int k = 0; k < a.N; k++) res[i, j] += a[i, k] * b[k, j];
  562. return res;
  563. }
  564. public CplxMatrix PowN(int n) {
  565. if (n > 1) {
  566. CplxMatrix r = PowN(n / 2), rr = r * r;
  567. return (n & 1) == 1 ? rr * this : rr;
  568. } else if (n == 0) return I(N);
  569. else if (n == 1) return Copy();
  570. else return Inv().PowN(-n);
  571. }
  572. public static Cplx[] operator *(CplxMatrix a, Cplx[] b) {
  573. if (a.N != b.Length) return null;
  574. Cplx[] res = new Cplx[a.N];
  575. for (int i = 0; i < a.N; i++) {
  576. res[i] = 0;
  577. for (int j = 0; j < a.N; j++)
  578. res[i] += a[i, j] * b[j];
  579. }
  580. return res;
  581. }
  582. public static Cplx[,] operator *(CplxMatrix a, Cplx[,] b) {
  583. if (a.N != b.GetLength(0)) return null;
  584. Cplx[,] res = new Cplx[b.GetLength(0), b.GetLength(1)];
  585. for (int i = 0; i < res.GetLength(0); i++)
  586. for (int j = 0; j < res.GetLength(1); j++) {
  587. res[i, j] = 0;
  588. for (int k = 0; k < a.N; k++) res[i, j] += a[i, k] * b[k, j];
  589. }
  590. return res;
  591. }
  592. public static Cplx[] LinEq(Cplx[,] a0) {
  593. int n = a0.GetLength(0);
  594. Cplx[,] a = new Cplx[n, n];
  595. for (int i = 0; i < n; i++)
  596. for (int j = 0; j < n; j++)
  597. a[i, j] = a0[i, j].Copy;
  598. CplxMatrix mx = new CplxMatrix(a).Inv();
  599. Cplx[] z = new Cplx[n];
  600. for (int i = 0; i < n; i++) {
  601. z[i] = 0;
  602. for (int j = 0; j < n; j++)
  603. z[i] += mx[i, j] * a0[j, n];
  604. }
  605. return z;
  606. }
  607. public void Show(string name, int total, int frac, bool real = false) {
  608. Console.WriteLine(name + ":");
  609. Cplx.Show(a, total, frac, real);
  610. }
  611. Polyn[,] M() {
  612. Polyn[,] p = new Polyn[N, N];
  613. for (int i = 0; i < N; i++)
  614. for (int j = 0; j < N; j++) {
  615. p[i, j] = (i == j) ? new Polyn(a[i, j], -1) :
  616. new Polyn(a[i, j]);
  617. }
  618. return p;
  619. }
  620. static void Swap(int[] q, int i, int j) {
  621. if (i != j) { q[i] ^= q[j]; q[j] ^= q[i]; q[i] ^= q[j]; }
  622. }
  623. public static int[] Permut(int n, int k) {
  624. int[] p = new int[n];
  625. for (int i = 0; i < n; i++) p[i] = i;
  626. for (int i = 0; i < n; i++) {
  627. int j = k % (n - i);
  628. Swap(p, i, i + j);
  629. k /= (n - i);
  630. }
  631. return p;
  632. }
  633. static int Σ(int[] p) {
  634. int σ = 1;
  635. for (int i = 0; i < p.Length - 1; i++)
  636. for (int j = i + 1; j < p.Length; j++)
  637. if (p[i] > p[j]) σ = -σ;
  638. return σ;
  639. }
  640. static int Fact(int n) => n == 0 ? 1 : n * Fact(n - 1);
  641. }
  642. class Program {
  643. delegate bool db(Cplx z);
  644. delegate int di(Polyn p);
  645. delegate Cplx func(double p);
  646. static Random rnd = new Random();
  647. public static void Main() {
  648. void Show(Polyn p) {
  649. string s = "";
  650. string Pow(int n) {
  651. string pw = "⁰¹²³⁴⁵⁶⁷⁸⁹";
  652. string res = "";
  653. for (int i = 0; i < n.ToString().Length; i++) {
  654. res += "⁰¹²³⁴⁵⁶⁷⁸⁹"[n.ToString()[i] - '0'];
  655. }
  656. return res;
  657. }
  658. for (int i = p.N; i >= 0; i--) {
  659.  
  660. //int a0 = (int)Math.Round(p[i].X);
  661. //int a = Math.Abs((int)Math.Round(p[i].X));
  662. double a0 = Math.Round(p[i].X, 3);
  663. double a = Math.Abs(a0);
  664. int σ = Math.Sign(a0);
  665. if (a == 0) continue;
  666. s += s != "" || s == "" && σ == -1 ? (σ == 1 ? " + " : " – ") : "";
  667. if (i == 0 || i > 0 && a != 1)
  668. s += a;
  669. if (i > 0) {
  670. s += "x";
  671. if (i > 1) {
  672. s += Pow(i);
  673. //s += "^";
  674. //s += i;
  675. }
  676. }
  677.  
  678. }
  679. Console.WriteLine(s);
  680. }
  681. //var w = new Cplx(2, 1);
  682. //var ww = w.Pow(w);
  683. //ww.Show();
  684. cplxFn f = w => w * w.Ln;
  685. var w0 = Cplx.Newton(w => f(w + 1) - f(w - 1) + new Cplx(0, 4 * Math.PI), new Cplx(rnd.NextDouble() * 10 - 5, rnd.NextDouble() * 10 - 5));
  686. //var w0 = Cplx.Newton(w => w.Pow(w) - w, new Cplx(2.86, 3.22));
  687. w0.Show(16, 12);
  688. f(w0 - 1).Show(18, 14);
  689. f(w0 + 1).Show(18, 14);
  690. //((w0 - 1) * w0.Ln).Show();
  691. //w0.Ln.Show();
  692. return;
  693. Cplx sum = 0;
  694. int n = 100000;
  695. double π = Math.PI;
  696. double r = 0.51;
  697. int σ = 1;
  698. for (int i = 0; i < n; i++) {
  699. double φ0 = 2 * π * i / n, φ1 = 2 * π * (i + 1) / n;
  700. Cplx z0 = r * (Cplx.I * φ0).Exp;
  701. Cplx z1 = r * (Cplx.I * φ1).Exp;
  702. Cplx dz = σ * (z1 - z0);
  703. double φ = 2 * π * (i + 0.5) / n;
  704. Cplx z = r * (Cplx.I * φ).Exp;
  705. sum += (-1d / 2) / (z * z - 5 * Cplx.I * z / 2 - 1) * dz;
  706. //sum += 1 / (z - Cplx.I / 2) * dz;
  707. }
  708. sum.Show();
  709. }
  710. }
  711. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement