using System; using System.Linq; using System.Collections.Generic; using System.IO; //using Float.Flt; //using Float.Cplx; //using Float.Equations; //using Float.Program; using FuncSqrt; //using FuncSqrt.Flt; namespace FuncSqrt; delegate Flt func(Flt x); delegate Flt funcN(Flt[] x); delegate Flt func2(Flt x, Flt y); static class Arr { static Random rnd = new Random(); public static int NZ0(uint[] a) { int res = a.Length; for (int i = a.Length - 1; i >= 0; i--) if (a[i] == 0) res--; else break; return res; } public static int NZ(uint[] a) { int n0 = NZ0(a); if (n0 == 0) return 0; int n1 = 4; for (int i = 3; i >= 0; i--) if ((a[n0 - 1] & (0xff << (i << 3))) == 0) n1--; else break; uint b = (a[n0 - 1] >> ((n1 - 1) << 3)) & 0xff; int n2 = 8; for (int i = 7; i >= 0; i--) if ((b & (1 << i)) == 0) n2--; else break; return ((n0 - 1) << 5) + ((n1 - 1) << 3) + n2; } public static uint[] Copy(uint[] a) { uint[] aa = new uint[a.Length]; a.CopyTo(aa, 0); return aa; } public static uint[] Rand(int size) { uint[] a = new uint[size]; for (int i = 0; i < size; i++) for (int j = 0; j < 4; j++) a[i] |= (uint)rnd.Next(0x100) << (j << 3); return a; } static string StrB(uint x) { string s = ""; for (int i = 31; i >= 0; i--) s += ((x & (1u << i)) > 0 ? '1' : '0'); return s; } public static string StrB(uint[] a, bool nz = false) { if (NZ0(a) == 0 && nz) return "0"; string s = ""; for (int i = a.Length - 1; i >= 0; i--) s += StrB(a[i]); if (nz) { int k = 0; for (int i = 0; i < s.Length; i++) if (s[i] == '0') k++; else break; s = s.Substring(k, s.Length - k); } return s; } public static bool GetBit(uint[] a, int n) { if (n > (a.Length << 5) - 1) return false; int n0 = n >> 5, n1 = n & 0x1f; return (a[n0] & (1 << n1)) > 0; } public static void SetBit(ref uint[] a, int n, bool f) { int n0 = n >> 5, n1 = n & 0x1f; if (f && n0 > a.Length - 1) { uint[] aa = new uint[n0 + 1]; a.CopyTo(aa, 0); a = aa; } a[n0] = f ? (a[n0] | (1u << n1)) : (a[n0] & ~(1u << n1)); } public static void SetLen(ref uint[] a, int l) { uint[] aa = new uint[l]; if (l >= a.Length) a.CopyTo(aa, 0); else for (int i = 0; i < l; i++) aa[i] = a[i]; a = aa; } public static int Cmp(uint[] a, uint[] b) { int m = Math.Max(a.Length, b.Length); for (int i = m - 1; i >= 0; i--) { uint aa = i < a.Length ? a[i] : 0; uint bb = i < b.Length ? b[i] : 0; if (aa < bb) return -1; else if (aa > bb) return 1; } return 0; } public static uint[] ShiftL(uint[] a, int n) { //Console.WriteLine("In ShiftL"); if (n == 0 || NZ0(a) == 0) return Copy(a); int p = NZ(a) - 1 + n; //Console.WriteLine($"p={p} NZ(a)-1={NZ(a) - 1}"); int l = Math.Max((p >> 5) + 1, a.Length); //Console.WriteLine($"l={l} a.Length={a.Length}"); uint[] aa = new uint[l]; int n0 = n >> 5, n1 = n & 0x1f; for (int i = Arr.NZ0(a) - 1; i >= 0; i--) aa[i + n0] = a[i]; if (n1 > 0) { //Console.WriteLine("In Sh"); for (int i = aa.Length - 1; i >= 0; i--) { ulong b = ((ulong)aa[i] << 32) | (i > 0 ? aa[i - 1] : 0); b <<= n1; aa[i] = (uint)(b >> 32); } //Console.WriteLine("aa after n1"); //Console.WriteLine(Arr.StrB(aa)); } return aa; } public static uint[] ShiftR(uint[] a, int n) { if (n == 0) return Copy(a); bool f = GetBit(a, n - 1); int n0 = n >> 5, n1 = n & 0x1f; uint[] aa = new uint[a.Length]; for (int i = n0; i < a.Length; i++) aa[i - n0] = a[i]; if (n1 > 0) for (int i = 0; i < NZ0(aa); i++) { ulong b = aa[i] | ((ulong)(i < aa.Length - 1 ? aa[i + 1] : 0) << 32); aa[i] = (uint)(b >> n1); } if (f) aa = Add(aa, new uint[] { 1 }); return aa; } public static uint[] Add(uint[] a, uint[] b) { uint[] res = new uint[Math.Max(a.Length, b.Length) + 1]; uint c = 0; for (int i = 0; i < res.Length; i++) { ulong aa = i < a.Length ? a[i] : 0; ulong s = aa + (i < b.Length ? b[i] : 0) + c; res[i] = (uint)s; c = (s > uint.MaxValue ? 1u : 0); } return res; } public static void DelZeros(ref uint[] a) => SetLen(ref a, NZ0(a)); public static void DelFirstZeros(ref uint[] a) { int k = 0; for (int i = 0; i < a.Length; i++) if (a[i] == 0) k++; else break; uint[] aa = new uint[a.Length - k]; for (int i = 0; i < aa.Length; i++) aa[i] = a[k + i]; a = aa; } public static uint[] Sub(uint[] a, uint[] b) { uint[] res = new uint[Math.Max(a.Length, b.Length)]; uint c = 0; for (int i = 0; i < res.Length; i++) { ulong s = 0x100000000ul + (i < a.Length ? a[i] : 0) - (i < b.Length ? b[i] : 0) - c; res[i] = (uint)s; c = s > uint.MaxValue ? 0 : 1u; } return res; } public static uint[] Mul(uint[] a, uint[] b) { uint[] res = new uint[a.Length + b.Length]; for (int i = 0; i < NZ0(a); i++) for (int j = 0; j < NZ0(b); j++) { ulong w = (ulong)a[i] * b[j]; int k0 = i + j; uint[] u = new uint[] { (uint)w, (uint)(w >> 32) }; for (int t = 0; t < 2; t++) { int k = k0 + t; do { ulong v = (ulong)res[k] + u[t]; res[k++] = (uint)v; u[t] = (uint)(v >> 32); } while (u[t] > 0); } } return res; } } class Flt { public static readonly int N = Program.N; uint[] a; public uint[] A => Arr.Copy(a); public int exp { get; private set; } public int sign { get; private set; } public Flt(uint[] a, int exp, int sign) { if (sign != 0) { var v = Norm(a, exp); this.a = v.a; this.exp = v.exp; this.sign = Arr.NZ0(this.a) == 0 ? 0 : Math.Sign(sign); } else { this.a = new uint[N]; this.exp = 0; this.sign = 0; } } public Flt(Flt x) : this(x.a, x.exp, x.sign) { } public static implicit operator Flt(long n) { if (n == 0) return Zero; int sign = Math.Sign(n); ulong u = (ulong)Math.Abs(n); uint[] a = new uint[N]; a[^2] = (uint)u; a[^1] = (uint)(u >> 32); int d = (N << 5) - Arr.NZ(a); int exp = 64 - d; a = Arr.ShiftL(a, d); return new Flt(a, exp, sign); } public static implicit operator Flt(double x) { string str = x.ToString(); bool neg = str[0] == '-'; if (neg) str = str.Substring(1, str.Length - 1); int p = str.IndexOf(','); if (p == -1) { long n = long.Parse(str); return (Flt)(neg ? -n : n); } else { string str0 = str.Substring(0, p); string str1 = str.Substring(p + 1, str.Length - p - 1); Flt res = (Flt)long.Parse(str0) + (Flt)long.Parse(str1) / ((Flt)10).Pow(str1.Length); return neg ? -res : res; } } public Flt Copy => new Flt(a, exp, sign); public static Flt[] CopyArr(Flt[] x) { Flt[] xx = new Flt[x.Length]; for (int i = 0; i < x.Length; i++) xx[i] = x[i].Copy; return xx; } public Flt Abs => sign >= 0 ? Copy : new Flt(a, exp, 1); public static Flt Zero => new Flt(new uint[N], 0, 0); public bool IsZero => sign == 0; public static (uint[] a, int exp) Norm(uint[] a0, int exp0) { //Console.WriteLine("In norm"); int n = Arr.NZ(a0); if (n == 0) return (new uint[N], 0); if (n == (N << 5)) { //Console.WriteLine("Here"); uint[] a = Arr.Copy(a0); Arr.SetLen(ref a, N); return (a, exp0); } else if (n < (N << 5)) return (Arr.ShiftL(a0, (N << 5) - n), exp0 + n - (N << 5)); else { uint[] a = Arr.ShiftR(a0, n - (N << 5)); Arr.SetLen(ref a, N); int exp = exp0 + n - (N << 5); return (a, exp); } } public static Flt operator *(Flt x, Flt y) { if (x.IsZero || y.IsZero) return Zero; uint[] z = Arr.Mul(x.A, y.A); z = Arr.ShiftR(z, (N << 5)); //Console.WriteLine("in * before norm\n" + Arr.StrB(z)); //var v = Norm(z, x.exp + y.exp - (N << 5)); var v = Norm(z, x.exp + y.exp); //Console.WriteLine("in * after norm\n" + Arr.StrB(v.a)); Arr.SetLen(ref v.a, N); return new Flt(v.a, v.exp, x.sign * y.sign); } public bool Eq(Flt x) => sign == x.sign && exp == x.exp && Arr.Cmp(a, x.a) == 0; public bool NEq(Flt x) => !Eq(x); public Flt Pow(int n) { if (n > 1) { Flt x = Pow(n / 2); Flt x2 = x * x; return (n & 1) == 0 ? x2 : x2 * this; } else if (n == 1) return Copy; else if (n == 0) return 1; else return 1 / Pow(-n); } public static int CmpAbs(Flt x, Flt y) { if (x.exp != y.exp) return Math.Sign(x.exp - y.exp); else return Arr.Cmp(x.a, y.a); } static int Cmp(Flt x, Flt y) { if (x.sign != y.sign) return Math.Sign(x.sign - y.sign); else return CmpAbs(x, y) * x.sign; } public static bool operator <(Flt x, Flt y) { if (x.sign != y.sign) return x.sign < y.sign; else if (x.sign == 0) return false; else { int s = CmpAbs(x, y); return x.sign == -s; } } public static bool operator >(Flt x, Flt y) => y < x; public Flt Root(int n) { Flt x = 1; bool f = false; while (true) { Flt xx = (this / x.Pow(n - 1) + x * (n - 1)) / n; if (f) return xx; f = x.exp == xx.exp && x.sign == xx.sign; for (int i = x.A.Length - 1; i > 0; i--) f &= x.A[i] == xx.A[i]; x = xx; } } public static Flt operator /(Flt x, Flt y) { if (y.IsZero) return null; else if (x.IsZero) return Zero; uint[] xa = x.A; uint[] ya = y.A; Arr.SetLen(ref xa, N + 1); Arr.SetLen(ref ya, N + 1); xa = Arr.ShiftL(xa, 31); ya = Arr.ShiftL(ya, 16); uint[] z = new uint[1]; while (Arr.NZ(z) <= (N << 5)) { ulong u = ((ulong)xa[^1] << 16) | (xa[^2] >> 16); uint v = (ya[^1] << 16) | (ya[^2] >> 16); uint w = (uint)(u / v); uint[] yw = Arr.Mul(ya, new uint[] { w }); Arr.SetLen(ref yw, N + 1); if (Arr.Cmp(xa, yw) < 0) { w--; yw = Arr.Sub(yw, ya); } xa = Arr.Sub(xa, yw); xa = Arr.ShiftL(xa, 16); z = Arr.ShiftL(z, 16); z[0] |= w; } int n = (Arr.NZ0(z) << 5) - Arr.NZ(z); //Console.WriteLine($"n={n}"); z = Arr.ShiftL(z, n); //Console.WriteLine($"zlen={z.Length}"); uint[] zz = new uint[N]; for (int i = 0; i < N; i++) zz[^(1 + i)] = z[^(1 + i)]; //Console.WriteLine(Arr.StrB(zz)); return new Flt(zz, x.exp - y.exp + 1 - (n & 1), x.sign * y.sign); } public static Flt operator +(Flt x, Flt y) { if (x.sign == 0) return y.Copy; else if (y.sign == 0) return x.Copy; else { uint[] xa = x.A; uint[] ya = y.A; int expX = x.exp, expY = y.exp; bool f = expX > expY || expX == expY && Arr.Cmp(xa, ya) >= 0; int d = Math.Max(expX, expY) - Math.Min(expX, expY); uint[] r; if (f) ya = Arr.ShiftR(ya, d); else xa = Arr.ShiftR(xa, d); //Console.WriteLine("In + xa\n" + Arr.StrB(xa)); //Console.WriteLine("In + ya\n" + Arr.StrB(ya)); if (x.sign == y.sign) r = Arr.Add(xa, ya); else r = f ? Arr.Sub(xa, ya) : Arr.Sub(ya, xa); //Console.WriteLine("In + xa+ya\n" + Arr.StrB(r)); var v = Norm(r, Math.Max(x.exp, y.exp)); int sign = f ? x.sign : y.sign; return new Flt(v.a, v.exp, sign); } } public static Flt operator -(Flt x) => new Flt(x.a, x.exp, -x.sign); public static Flt operator -(Flt x, Flt y) => x + -y; public static explicit operator Flt(string str) { int indP(string str) { int i = str.IndexOf('.'); return i > -1 ? i : str.IndexOf(','); } Flt IntToFlt(string str) { Flt x = 0; for (int i = 0; i < str.Length; i++) { x *= 10; x += (int)str[i] - (int)'0'; } return x; } Flt FracToFlt(string str) { int p = indP(str); if (p == -1) return IntToFlt(str); string s0 = str.Substring(0, p); string s1 = str.Substring(p + 1, str.Length - (p + 1)); return IntToFlt(s0) + IntToFlt(s1) / ((Flt)10).Pow(s1.Length); } if (str[0] == '-') return -(Flt)str.Substring(1, str.Length - 1); else if (str[0] == '+') return (Flt)str.Substring(1, str.Length - 1); int indE = str.IndexOf('E'); if (indE == -1) indE = str.IndexOf('e'); if (indE == -1) return FracToFlt(str); else { string s0 = str.Substring(0, indE); string s1 = str.Substring(indE + 1, str.Length - (indE + 1)); return FracToFlt(s0) * ((Flt)10).Pow(int.Parse(s1)); } } public void ShowB(string name = "") { if (name != "") Console.WriteLine(name + " =\n"); string s = Arr.StrB(a); for (int i = 0; i < N; i++) { string s1 = s.Substring((i << 5), 16); Console.ForegroundColor = ConsoleColor.Yellow; Console.Write(s1); string s2 = s.Substring((i << 5) + 16, 16); Console.ForegroundColor = ConsoleColor.Cyan; Console.Write(s2); } Console.Write("\n"); Console.ForegroundColor = ConsoleColor.White; Console.Write($"exp={exp} sign={sign}\n\n"); } public string Str10 => new Dcm(this).Str(); public string Str10Round(int n) => new Dcm(this).Str(n); public void Show10(string name = "") { if (name != "") Console.WriteLine(name + " ="); Console.WriteLine(Str10); } public void Show10Round(int n, string name = "") { if (name != "") Console.WriteLine(name + " ="); string str = Str10Round(n); if (str[0] != '-') str = " " + str; Console.WriteLine(str); } public static void ShowArray(Flt[] a, int prec, string name, ConsoleColor clr = ConsoleColor.White) { Console.ForegroundColor = clr; for (int i = 0; i < a.Length; i++) { if (name != "") Console.WriteLine($"{name}[{i}] = "); a[i].Show10Round(prec); } Console.ForegroundColor = ConsoleColor.White; } public class Dcm { Flt x; public Dcm(Flt x) { this.x = x.Copy; } public uint[] Int() { if (x.exp <= 0) return new uint[0]; uint[] a0 = Arr.ShiftL(x.a, x.exp); uint[] a = new uint[a0.Length - N]; for (int i = 0; i < a.Length; i++) a[i] = a0[N + i]; return a; } public uint[] Frac() { //Console.WriteLine("In Frac"); if (x.exp == 0) return Arr.Copy(x.a); else if (x.exp > 0) { //Console.WriteLine("Aaaa"); uint[] a = Arr.ShiftL(x.a, x.exp); Arr.SetLen(ref a, N); return a; } else { //Console.WriteLine("Bbbb"); int n = ((-x.exp - 1) >> 5) + 1; //Console.WriteLine("n =" + n); uint[] a = x.A; Arr.SetLen(ref a, N + n); a = Arr.ShiftL(a, (n << 5) + x.exp); return a; } } static (uint[] q, uint r) Div(uint[] a, uint b) { int n0 = Arr.NZ0(a); uint[] q; uint r; if (n0 < 3) { ulong u = n0 == 0 ? 0 : (n0 == 1 ? a[0] : ((ulong)a[1] << 32) | a[0]); ulong q0 = u / b; r = (uint)(u % b); q = (q0 > uint.MaxValue ? new uint[] { (uint)q0, (uint)(q0 >> 32) } : new uint[] { (uint)q0 }); } else { uint[] aa = Arr.Copy(a); Arr.DelZeros(ref aa); uint[] bb0 = new uint[] { b }; uint[] bb = Arr.ShiftL(bb0, Arr.NZ(aa) - Arr.NZ(bb0) + 1); q = new uint[] { 0 }; while (Arr.Cmp(bb, bb0) == 1) { bb = Arr.ShiftR(bb, 1); q = Arr.ShiftL(q, 1); if (Arr.Cmp(aa, bb) >= 0) { aa = Arr.Sub(aa, bb); q[0]++; } } r = aa[0]; } return (q, r); } static string DcmInt(uint[] a) { var l = new List(); uint[] aa = Arr.Copy(a); while (Arr.NZ0(aa) > 0) { var v = Div(aa, 1000000000u); l.Add(v.r); aa = v.q; } if (l.Count == 0) return "0"; string res = ""; for (int i = l.Count - 1; i >= 0; i--) { string s = l[i].ToString(); if (i < l.Count - 1) s = new String('0', 9 - s.Length) + s; res += s; } return res; } public string DcmInt() => DcmInt(Int()); static string DcmFrac(uint[] a) { //Console.WriteLine("In DcmFrac"); uint[] aa = Arr.Copy(a); var l = new List(); while (Arr.NZ0(aa) > 0) { Arr.DelFirstZeros(ref aa); Arr.SetLen(ref aa, aa.Length + 1); uint c = 0; for (int i = 0; i < aa.Length; i++) { ulong s = (ulong)aa[i] * 1000000000u + c; aa[i] = (uint)s; c = (uint)(s >> 32); } l.Add(aa[aa.Length - 1]); Arr.SetLen(ref aa, aa.Length - 1); } string res = ""; for (int i = 0; i < l.Count; i++) { string s = l[i].ToString(); s = new String('0', 9 - s.Length) + s; res += s; } return res == "" ? "0" : res; } public string DcmFrac() => DcmFrac(Frac()); public string Str() { if (x.IsZero) return "0"; int l = (int)Math.Round(N * 9.6) - 1; string res; string sInt = DcmInt(); string sFrac = DcmFrac(); //if (sInt.Length >= l) return sInt; if (sInt.Length >= l) res = sInt; else if (sInt != "0") { if (sInt.Length + sFrac.Length <= l) res = sInt + "." + sFrac; else res = sInt + "." + sFrac.Substring(0, l - sInt.Length); } else { int k = 0; for (int i = 0; i < sFrac.Length; i++) if (sFrac[i] != '0') break; else k++; //return "0." + (sFrac.Length < k + l ? sFrac : sFrac.Substring(0, k + l)); res = "0." + (sFrac.Length < k + l ? sFrac : sFrac.Substring(0, k + l)); } return (x.sign < 0 ? "-" : "") + res; } public string Str(int n) { string str = Str(); //Console.WriteLine("str=" + str); bool neg = str[0] == '-'; //Console.WriteLine("neg=" + neg); //Console.WriteLine(neg + " " + str[0]); if (neg) str = str.Substring(1, str.Length - 1); //Console.WriteLine("str0=" + str); int p = str.IndexOf('.'); if (str.Length <= p + 1 + n) return (neg ? "-" : "") + str;//str; bool f = str[p + 1 + n] > '4'; str = str.Substring(0, p + 1 + n); if (!f) return (neg ? "-" : "") + str; int k = -1; for (int i = str.Length - 1; i >= 0; i--) if (Char.IsDigit(str[i]) && str[i] < '9') { k = i; break; } if (k == -1) { str = "0" + str; k = 0; } //Console.WriteLine("str=" + str); string str1 = str.Substring(0, k); str1 += (char)((int)str[k] + 1); for (int i = k + 1; i < str.Length; i++) str1 += (str[i] == '.' ? '.' : '0'); //Console.WriteLine("check neg"); return neg ? "-" + str1 : str1; } } public class Polyn { Flt[] a; public int L => a.Length; public int N => L - 1; public Flt[] A { get { Flt[] aa = new Flt[L]; for (int i = 0; i < L; i++) aa[i] = a[i].Copy; return aa; } } public Flt this[int i] { get { return a[i].Copy; } set { a[i] = value.Copy; } } public Polyn(params Flt[] a) { this.a = new Flt[a.Length]; for (int i = 0; i < a.Length; i++) this.a[i] = a[i].Copy; } //public static Polyn Zero => new Polyn(new Flt[0]); public static Polyn operator -(Polyn p) { Flt[] a = new Flt[p.L]; for (int i = 0; i < p.L; i++) a[i] = -p[i]; return new Polyn(a); } public static Polyn operator +(Polyn p1, Polyn p2) { Flt[] a = new Flt[Math.Max(p1.L, p2.L)]; for (int i = 0; i < a.Length; i++) a[i] = (i < p1.L ? p1[i] : 0) + (i < p2.L ? p2[i] : 0); return new Polyn(a); } public static Polyn operator -(Polyn p1, Polyn p2) => p1 + -p2; public static Polyn operator *(Polyn p, Flt μ) { Polyn pp = new Polyn(p.a); for (int i = 0; i < pp.L; i++) pp[i] *= μ; return pp; } public static Polyn operator *(Flt μ, Polyn p) => p * μ; public static Polyn operator *(Polyn p1, Polyn p2) { Flt[] a = new Flt[p1.L + p2.L - 1]; for (int i = 0; i < p1.L; i++) for (int j = 0; j < p2.L; j++) { int k = i + j; if (a[k] == null) a[k] = 0; a[k] += p1[i] * p2[j]; } return new Polyn(a); } public Flt Value(Flt x) { Flt res = 0; for (int i = N; i >= 0; i--) { res *= x; res += a[i]; } return res; } public Polyn D1 { get { if (L == 0) return new Polyn(new Flt[0]); Flt[] aa = new Flt[N]; for (int i = 0; i < N; i++) aa[i] = a[i + 1] * (i + 1); return new Polyn(aa); } } public Polyn D2 => D1.D1; public Polyn Intg { get { Flt[] aa = new Flt[L + 1]; aa[0] = 0; for (int i = 1; i <= L; i++) aa[i] = a[i - 1] / i; return new Polyn(aa); } } } static class M { public static readonly Flt π = GetPi(); static Flt GetPi() => (ArcTg0((Flt)1 / 5) * 4 - ArcTg0((Flt)1 / 239)) * 4; //ArcSin0((Flt)1 / 2) * 6; public static readonly Flt ε0 = GetEps(-Flt.N * 16); public static readonly Flt ε1 = GetEps(-Flt.N * 11); public static readonly Flt ε2 = GetEps(-Flt.N * 22); static Flt GetEps(int n) { uint[] a = new uint[Flt.N]; Arr.SetBit(ref a, (Flt.N << 5) - 1, true); return new Flt(a, n, 1); } public static readonly Flt ln2 = GetLn2(); static Flt GetLn2() => Ln0((Flt)4 / 3) - Ln0((Flt)2 / 3); static Flt Ln0(Flt x) { Flt y = x - 1; Flt res = y.Copy; Flt z = y.Copy; int n = 1; while (true) { Flt t = res + (z = -z * y) / (++n); if (t.Eq(res)) return res; res = t; } } public static Flt Ln(Flt x) { if (x.sign != 1) return null; Flt y = new Flt(x.A, 1, 1); Flt z = (y - 1) / (y + 1); return Ln0(1 + z) - Ln0(1 - z) + (x.exp - 1) * ln2; } public static Flt Sin(Flt x) { Flt t = x.Copy; while (t.Abs > M.π) if (t.sign == 1) t -= M.π * 2; else t += M.π * 2; Flt xx = t * t; Flt y = t.Copy; Flt s = t.Copy; int n = 0; //bool f = false; while (true) { n += 2; y = -y * xx / (n * (n + 1)); Flt st = s + y; if (st.Eq(s)) return st; s = st; } } public static Flt Cos(Flt x) { Flt t = x.Copy; while (t.Abs > M.π) if (t.sign == 1) t -= M.π * 2; else t += M.π * 2; Flt xx = t * t; Flt y = 1; Flt s = 1; int n = 0; while (true) { n += 2; y = -y * xx / (n * (n - 1)); Flt st = s + y; if (st.Eq(s)) break; s = st; } return s; } public static Flt Tg(Flt x) => Sin(x) / Cos(x); public static Flt Exp(Flt x) { if (x.IsZero) return 1; else if (x.sign == -1) return 1 / Exp(-x); else { Flt y = x.exp > 0 ? new Flt(x.A, 0, 1) : x.Copy; Flt res = 1, z = 1; int n = 0; while (true) { Flt t = res + (z *= y / ++n); if (t.Eq(res)) break; res = t; } for (int i = 0; i < x.exp; i++) res = res.Pow(2); return res; } } public static Flt ArcTg0(Flt x) { Flt y = x.Copy; Flt xx = x * x; Flt res = x.Copy; int n = 1; while (true) { n += 2; Flt t = res + (y = -y * xx) / n; if (t.Eq(res)) break; res = t; } return res; } public static Flt ArcTg(Flt x) { if (x < 0) return -ArcTg(-x); else if (x > 1) return π / 2 - ArcTg(1 / x); else return Newton(t => Tg(t) - x, (Flt)11 / 28); } public static Flt ArcSin0(Flt x) => Newton(t => Sin(t) - x, (Flt)11 / 28); public static Flt Pow(Flt x, Flt y) => Exp(y * Ln(x)); public static func D1(func F) => x => (F(x + ε1 / 2) - F(x - ε1 / 2)) / ε1; public static Flt Newton(func F, Flt x0) { func Fx = D1(F); int cnt = 0; Flt x = x0.Copy; bool f = false; while (true) { Flt xx = x - F(x) / Fx(x); if (cnt == 5) return (xx + x) / 2; f = (x - xx).Abs < ε0; if (f) cnt++; else cnt = 0; x = xx; } } public static Flt Intg(func F, Flt a, Flt b, int n = 1) { if (n == 1) { Flt h = b - a, hh = h / 6; Flt x0 = a; Flt x6 = b; Flt x3 = (a + b) / 2; Flt x1 = a + hh; Flt x2 = x3 - hh; Flt x4 = x3 + hh; Flt x5 = x6 - hh; Flt y3 = F(x3); Flt y24 = (F(x2) + F(x4)) / 2; Flt y15 = (F(x1) + F(x5)) / 2; Flt y06 = (F(x0) + F(x6)) / 2; return (y3 * 136 + y24 * 27 + y15 * 216 + y06 * 41) * h / 420; } else { Flt res = 0; Flt h = (b - a) / n; for (int i = 0; i < n; i++) { Flt aa = a + h * i; Flt bb = aa + h; res += Intg(F, aa, bb); } return res; } } public static Flt Dihotom(func F, Flt a, Flt b, Flt ε) { Flt aa = a.Copy, bb = b.Copy; Flt fa = F(aa), fb = F(bb); while ((aa - bb).Abs > ε) { Flt c = (aa + bb) / 2, fc = F(c); if (fc.Eq(0)) return c; if (fc.sign == fa.sign) aa = c.Copy; else bb = c.Copy; } return (fb * aa - fa * bb) / (fb - fa); } } /* public static Flt[] Grad(funcN f, Flt[] x) { var gr = new Flt[x.Length]; for (int i = 0; i < gr.Length; i++) gr[i] = D1(f, i, x); return gr; } */ class Equations { public static Flt D1(funcN f, Flt[] x, int n) { Flt[] x1 = Flt.CopyArr(x); x1[n] -= M.ε1 / 2; Flt[] x2 = Flt.CopyArr(x); x2[n] += M.ε1 / 2; return (f(x2) - f(x1)) / M.ε1; } public static Flt[] Grad(funcN f, Flt[] x) { var gr = new Flt[x.Length]; for (int i = 0; i < gr.Length; i++) gr[i] = D1(f, x, i); return gr; } static (Flt[,] a, Flt[] b) Mx(funcN[] f, Flt[] x) { var a = new Flt[f.Length, f.Length]; var b = new Flt[f.Length]; for (int i = 0; i < f.Length; i++) { b[i] = -f[i](x); for (int j = 0; j < x.Length; j++) { a[i, j] = D1(f[i], x, j); } } return (a, b); } public static (Flt[] x, Flt max) Xnew(funcN[] f, Flt[] x0, double μ = 1) { var mx = Mx(f, x0); Flt[] dx = Matr.SolveLinEq(mx.a, mx.b); Flt[] x = new Flt[x0.Length]; for (int i = 0; i < x.Length; i++) x[i] = x0[i] + dx[i] * μ; Flt max = 0; for (int i = 0; i < dx.Length; i++) if (dx[i].Abs > max) max = dx[i].Abs; return (x, max); } public static Flt[] Solve(funcN[] f, Flt[] x0) { double μ = 0.5; bool s = false; Flt[] x = Flt.CopyArr(x0); while (true) { var v = Xnew(f, x, μ); if (s) return v.x; else if (v.max < M.ε2) s = true; else if (v.max < 0.001) μ = 1; x = v.x; //////////////////////////// Console.Clear(); for (int i = 0; i < x.Length; i++) x[i].Show10(); //Console.WriteLine(); //Console.ReadLine(); //////////////////////////// } } public static class Matr { static bool IsZero(Flt x) => x.Abs < M.ε2; static Flt[,] Inv(Flt[,] a0) { int n = a0.GetLength(0), nn = n * 2; Flt[,] a = new Flt[n, nn]; for (int i = 0; i < n; i++) for (int j = 0; j < nn; j++) a[i, j] = j < n ? a0[i, j].Copy : (j == i + n ? 1 : 0); for (int i = 0; i < n; i++) { if (IsZero(a[i, i])) { int ii = i; for (int j = i + 1; j < n; j++) if (!IsZero(a[j, i])) { ii = j; break; } if (ii == i) return null; for (int j = i; j < nn; j++) { Flt tmp = a[i, j].Copy; a[i, j] = a[ii, j].Copy; a[ii, j] = tmp; } } Flt aii = a[i, i].Copy; for (int j = i; j < nn; j++) a[i, j] /= aii; for (int j = 0; j < n; j++) if (j != i) { Flt aji = a[j, i].Copy; for (int k = i; k < nn; k++) a[j, k] -= aji * a[i, k]; } } Flt[,] aa = new Flt[n, n]; for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) aa[i, j] = a[i, j + n].Copy; return aa; } public static Flt[] Mul(Flt[,] a, Flt[] b) { int n = b.Length; Flt[] c = new Flt[n]; for (int i = 0; i < n; i++) { c[i] = 0; for (int j = 0; j < n; j++) c[i] += a[i, j] * b[j]; } return c; } public static Flt[] SolveLinEq(Flt[,] a, Flt[] b) { Flt[,] aa = Inv(a); if (aa == null) return null; return Mul(aa, b); } } } class Cplx { Flt x, y; public Flt X { get { return x.Copy; } private set { x = value.Copy; } } public Flt Y { get { return y.Copy; } private set { y = value.Copy; } } public Cplx(Flt x, Flt y) { this.x = x; this.y = y; } public Cplx(Cplx z) : this(z.x, z.y) { } public Cplx Copy => new Cplx(this); public Cplx(Flt x) : this(x, 0) { } public static implicit operator Cplx(Flt x) => new Cplx(x); public Flt R2 => x.Pow(2) + y.Pow(2); public Flt R => R2.Root(2); public static Cplx Zero => (Flt)0; public bool IsZero => R2.exp < -(Flt.N << 32); public static Cplx operator -(Cplx z) => new Cplx(-z.x, -z.y); public static Cplx operator ~(Cplx z) => new Cplx(z.x, -z.y); public static Cplx operator +(Cplx z1, Cplx z2) => new Cplx(z1.x + z2.x, z1.y + z2.y); public static Cplx operator -(Cplx z1, Cplx z2) => z1 + -z2; public static Cplx operator *(Cplx z1, Cplx z2) => new Cplx(z1.x * z2.x - z1.y * z2.y, z1.x * z2.y + z1.y * z2.x); public static Cplx operator /(Cplx z, Flt a) => new Cplx(z.x / a, z.y / a); public static Cplx operator /(Cplx z1, Cplx z2) => z1 * ~z2 / z2.R2; public Cplx Pow(int n) { if (n < 0) return (Flt)1 / Pow(-n); else if (n == 0) return (Flt)1; else if (n == 1) return Copy; else { int m = n / 2; Cplx z = Pow(m); Cplx res = z * z; if ((n & 1) == 1) res *= this; return res; } } public static Cplx[] RootsOf1(int n) { Cplx[] z = new Cplx[n]; Flt φ = M.π * 2 / n; Cplx ζ = new Cplx(M.Cos(φ), M.Sin(φ)); for (int i = 0; i < n; i++) z[i] = ζ.Pow(i); return z; } public void Show(string name = "") { if (name != "") Console.WriteLine(name + " ="); Console.ForegroundColor = ConsoleColor.Red; x.Show10(); if (!y.IsZero) { Console.ForegroundColor = ConsoleColor.Cyan; y.Show10(); } Console.ForegroundColor = ConsoleColor.White; } public class Polyn { Cplx[] a; public Cplx this[int i] { get { return a[i].Copy; } set { a[i] = value.Copy; } } public int L => a.Length; public int N => L - 1; public Polyn(params Cplx[] a) { int k = a.Length; for (int i = a.Length - 1; i >= 0; i--) if (a[i].IsZero) k--; else break; this.a = new Cplx[k]; for (int i = 0; i < k; i++) this.a[i] = a[i].IsZero ? (Flt)0 : a[i].Copy; } public static Polyn operator -(Polyn p) { Cplx[] a = new Cplx[p.L]; for (int i = 0; i < p.L; i++) a[i] = -p[i]; return new Polyn(a); } public static Polyn operator +(Polyn p1, Polyn p2) { Cplx[] a = new Cplx[Math.Max(p1.L, p2.L)]; for (int i = 0; i < a.Length; i++) a[i] = (i < p1.L ? p1[i] : (Flt)0) + (i < p2.L ? p2[i] : (Flt)0); return new Polyn(a); } public static Polyn operator -(Polyn p1, Polyn p2) => p1 + -p2; public static Polyn operator *(Polyn p, Cplx b) { Cplx[] a = new Cplx[p.L]; for (int i = 0; i < p.L; i++) a[i] *= b; return new Polyn(a); } public static Polyn operator *(Cplx b, Polyn p) => p * b; public static Polyn operator *(Polyn p1, Polyn p2) { Cplx[] a = new Cplx[p1.L + p2.L - 1]; for (int i = 0; i < a.Length; i++) a[i] = Zero; for (int i = 0; i < p1.L; i++) for (int j = 0; j < p2.L; j++) { a[i + j] += p1[i] * p2[j]; } return new Polyn(a); } } } public static class Program { public const int N = 32; static Random rnd = new Random(); class Pt { Flt x, y; public Flt X { get { return x; } private set { x = value.Copy; } } public Flt Y { get { return y; } private set { y = value.Copy; } } public Pt(Flt x, Flt y) { this.x = x.Copy; this.y = y.Copy; } public Pt Copy => new Pt(x, y); public static Pt operator -(Pt p) => new Pt(-p.x, -p.y); public static Pt operator +(Pt p1, Pt p2) => new Pt(p1.x + p2.x, p1.y + p2.y); public static Pt operator -(Pt p1, Pt p2) => p1 + -p2; public static Pt operator *(Pt p, Flt μ) => new Pt(p.x * μ, p.y * μ); public static Pt operator *(Flt μ, Pt p) => p * μ; public static Pt operator /(Pt p, Flt μ) => new Pt(p.x / μ, p.y / μ); public static Flt operator *(Pt p1, Pt p2) => p1.x * p2.x + p1.y * p2.y; public Flt R2 => this * this; public Flt R => R2.Root(2); public Pt E => this / R; } static Flt.Polyn Intpol(Pt[] p) { Flt.Polyn res = new Polyn(new Flt[0]); for (int i = 0; i < p.Length; i++) { Flt.Polyn q = new Flt.Polyn(new Flt[] { 1 }); Flt μ = 1; for (int j = 0; j < p.Length; j++) if (i != j) { q *= new Flt.Polyn(-p[j].X, 1); μ *= p[i].X - p[j].X; } //q *= p[i].Y / μ; res += q * (p[i].Y / μ); } return res; } static Flt.Polyn Pln(Flt[] prm) { Flt[] q = new Flt[prm.Length + 1]; q[0] = 0; q[1] = 0; for (int i = 2; i < q.Length; i++) q[i] = prm[i - 1]; return new Flt.Polyn(q); } static Flt I(Pt[] pt, int n) { Flt res = 0; int m = (pt.Length - 1) / n; for (int i = 0; i < m; i++) { Pt[] pt0 = new Pt[n + 1]; // Flt.Polyn p = Parab(pt0).Intg; for (int j = 0; j <= n; j++) pt0[j] = pt[i * n + j]; Flt.Polyn p = Intpol(pt0).Intg; res += p.Value(pt0[^1].X) - p.Value(pt0[0].X); } return res; } static Flt.Polyn Pln1(Pt[] pt) { var f = new funcN[pt.Length - 1]; for (int i = 0; i < f.Length; i++) { int j = i; f[j] = t => { var a = new Flt[pt.Length + 1]; a[0] = pt[0].Y.Copy; a[1] = 0; for (int k = 2; k < a.Length; k++) a[k] = t[k - 2].Copy; var pln = new Flt.Polyn(a); return pln.Value(pt[j + 1].X) - pt[j + 1].Y; }; } var a0 = new Flt[f.Length]; for (int i = 0; i < a0.Length; i++) a0[i] = 0; var a1 = Equations.Solve(f, a0); //Console.WriiteLine("OK"); //Console.ReadLine(); var aa = new Flt[a1.Length + 2]; aa[0] = pt[0].Y.Copy; aa[1] = 0; for (int i = 2; i < aa.Length; i++) aa[i] = a1[i - 2].Copy; return new Flt.Polyn(aa); } public static void Main() { int prec = 50; void Show(Flt x, string name = "") => x.Show10Round(prec, name); { Flt g(Flt[] a) { Flt res = 1; for (int i = a.Length - 1; i >= 0; i--) res = M.Pow(a[i], res); return res; } int nt = 9; Flt r0 = 2; var ft0 = new funcN[nt]; var ft = new funcN[nt]; for (int i = 0; i < ft.Length; i++) { int j = i; ft0[j] = t => { var x = new Flt[ft.Length]; for (int k = 0; k < x.Length; k++) { x[k] = t[(j + k) % x.Length]; } return g(x); }; ft[j] = t => ft0[j](t) - (r0 + j); } /* Flt R(Flt[] t) { Flt res = 0; for (int i = 0; i < ft.Length; i++) res += ft[i](t).Pow(2); return res; }*/ Flt R(Flt[] t) { Flt res = 0; for (int i = 0; i < ft.Length; i++) { var s = ft0[i](t) / (r0 + i); res += (s - 1).Pow(2) + (1 / s - 1).Pow(2); } return res; } /* var q0 = new Flt[] { 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2 }; */ var q0 = new Flt[nt]; /* for (int i = 0; i < nt; i++) q0[i] = 1.3; */ string path = "q2.txt"; if (File.Exists(path)) { var str = File.ReadAllLines(path); for (int i = 0; i < q0.Length; i++) q0[i] = (Flt)str[i]; } else { q0[0] = (Flt)"1.306026805891242807652545950399"; q0[1] = (Flt)"1.422496123930371511508945347115"; q0[2] = (Flt)"1.481968384211597922114406899939"; q0[3] = (Flt)"1.517880134636269721047143119619"; q0[4] = (Flt)"1.537379894011106094871255325070"; q0[5] = (Flt)"1.428813343342501327902427160487"; q0[6] = (Flt)"1.298260843634120034650137509060"; q0[7] = (Flt)"1.251605981176236997085156615874"; q0[8] = (Flt)"3.351161933057958748208486588071"; } var t = Equations.Solve(ft, q0); //Flt.ShowArray(t, 30, "solution"); var res1 = new Flt[ft.Length]; for (int i = 0; i < ft.Length; i++) res1[i] = ft[i](t) + r0 + i; Flt.ShowArray(t, 60, "q0", ConsoleColor.Yellow); Flt.ShowArray(res1, 60, "f", ConsoleColor.Cyan); return; //Console.WriteLine("Start"); Flt δ0 = R(q0); //Console.WriteLine("δ0 OK"); Flt hmax = 0.0005; Flt h = 0.00025; int cnt = 0; while (δ0 > ((Flt)10).Pow(-50)) { //Show(h); var grd = Equations.Grad(R, q0); var qt = new Flt[grd.Length]; for (int i = 0; i < grd.Length; i++) qt[i] = q0[i] - grd[i] * h; var δt = R(qt); if (δt < δ0) { q0 = Flt.CopyArr(qt); δ0 = δt.Copy; cnt++; Console.Clear(); Show(h); var res = new Flt[ft.Length]; for (int i = 0; i < ft.Length; i++) res[i] = ft[i](q0) + r0 + i; Flt.ShowArray(q0, prec, "q0", ConsoleColor.Yellow); Flt.ShowArray(res, prec, "f", ConsoleColor.Cyan); Show(δ0); if (cnt == 5) { var str = new string[q0.Length]; for (int i = 0; i < q0.Length; i++) str[i] = q0[i].Str10; File.WriteAllLines(path, str); Console.WriteLine("Saved"); if (h < hmax) h *= 1.25; cnt = 0; } } else { h *= 0.7; Console.ForegroundColor = ConsoleColor.Red; Show(h); Console.ForegroundColor = ConsoleColor.White; cnt = 0; } } } /* //Polyn p = new Polyn(1, 2, 3); func F0 = x => x * (x - 1) + 1; //func F0 = M.Exp; Flt FPln(Flt x, Flt[] q, Flt a) { Flt xx = (x - a).Pow(2); Flt res = 0; for (int i = q.Length - 1; i >= 0; i--) { res *= xx; res += q[i]; } return res; } Flt x1 = -0.1, x2 = 1.1; Flt xm = (x1 + x2) / 2; int n = 13; Flt Δ(int i, Flt[] q, Flt x1, Flt x2) { //Flt x = x1 + (x2 - x1) * i / (q.Length - 1); Flt x = xm + (x2 - xm) * i / (q.Length - 1); //var pln = new Polyn(q); //return pln.Value(pln.Value(x)) - F0(x); return FPln(FPln(x, q, 0.5), q, 0.5) - F0(x); } var f = new funcN[n]; for (int i = 0; i < n; i++) { int j = i; f[j] = q => Δ(j, q, x1, x2); } funcN R = q => { Flt r = 0; for (int i = 0; i < q.Length; i++) { Flt δ = Δ(i, q, x1, x2); r += δ * δ; } return r; }; var q0 = new Flt[n]; for (int i = 0; i < n; i++) q0[i] = 0; q0[0] = 0.69; q0[1] = 1.78; q0[2] = -5.9; var q = Equations.Solve(f, q0); for (int i = 0; i < q.Length; i++) { Show(q[i], i.ToString()); } Console.WriteLine(); //var pln = new Polyn(q); Flt dx = (x2 - xm) / (2 * (n - 1)); for (int i = 0; i <= 2 * (n - 1); i++) { var x = xm + dx * i; Show(x); //Show(pln.Value(x)); Show(FPln(x, q, 0.5)); Show(FPln(FPln(x, q, 0.5), q, 0.5)); //var t = pln.Value(pln.Value(x)) / F0(x); var t = FPln(FPln(x, q, 0.5), q, 0.5) / F0(x); Show(t); Console.WriteLine(); } return; */ } } }