API tools faq
paste
VNM24ix

ComplexMath

VNM24ix
Jun 9th, 2025
38
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C# 23.78 KB | Source Code | 0 0
raw download clone embed print report
  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. }
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!