实现 Math.Asin 迈克劳林(泰勒)展开式,结果比Math.Asin 慢一倍
项目中需要快速求解Asin(x) 的近似值,原以为用泰勒展开式会快一些,结果比原生的慢一倍。
Math.ASin
Time Elapsed: 9ms
Gen 0: 0
Gen 1: 0
Gen 2: 0
Maclaurin.ASin
Time Elapsed: 17ms
Gen 0: 4
Gen 1: 0
Gen 2: 0
各位,谁有能力改进?
附:
http://www.mathportal.org/formulas/pdf/taylor-series-formulas.pdf
http://pages.pacificcoast.net/~cazelais/260/maclaurin.pdf
using System; using System.Collections.Generic; using System.Linq; using System.Runtime.Remoting.Messaging; using System.Text; using Diagnostics; namespace Asin { class Program { static void Main(string[] args) { int count = 100000; List<double> values = new List<double>(count); Random r = new Random(); for (var i = 0; i <= count; ++i) { values .Add(r.NextDouble() * 2 - 1); } CodeTime.Init(); int? iter = 0; CodeTime.Timer("Math.ASin", count, () => { var i = iter.Value + 1; iter = i; Math.Asin(values[i]); }); iter = 0; CodeTime.Timer("Maclaurin.ASin", count, () => { var i = iter.Value + 1; iter = i; Maclaurin.Asin(values[i],3); }); while (true) { iter = 0; CodeTime.Timer("Math.ASin", count, () => { var i = iter.Value + 1; iter = i; Math.Asin(values[i]); }); iter = 0; CodeTime.Timer("Maclaurin.ASin", count, () => { var i = iter.Value + 1; iter = i; Maclaurin.Asin(values[i], 3); }); } //var ret = Maclaurin.Asin(0.5, 3); //var ret2 = Math.Asin(0.5); //Console.WriteLine(ret); //Console.WriteLine(ret2); Console.ReadLine(); } } class Maclaurin { class ASinImpl { private List<double> quotieties = new List<double>(); private IEnumerator<double> computeQuotieties = null; public ASinImpl() { this.computeQuotieties = ComputeQuotiety(); } public double Calc(double v, int precision = 2) { if (quotieties.Count < precision) { for (var i = quotieties.Count; i < precision; ++i) { computeQuotieties.MoveNext(); quotieties.Add(computeQuotieties.Current); } } double ret = 0; var values = ComputeValues(v); for (int i = 0; i < precision; ++i) { values.MoveNext(); ret += quotieties[i]*values.Current; } return ret; } private IEnumerator<double> ComputeValues(double v) { double ret = 1; double q = v*v; for(int i = 0;;++i) { if (i == 0) { ret = v; yield return ret; } else { ret *= q; yield return ret; } } throw new NotImplementedException(); } private IEnumerator<double> ComputeQuotiety() { for (int i = 0;; i++) { double up = Factorial(2*i); double down = Math.Pow(Math.Pow(2, i)*Factorial(i), 2)*(2*i + 1); double quotiety = up/down; yield return quotiety; } throw new NotImplementedException(); } private long Factorial(long v ) { if( v < 0) throw new ArgumentOutOfRangeException("v"); if (v == 0) return 1; if (v == 1) return 1; long ret = 1; for (int i = 2; i <= v; ++i) { ret *= i; } return ret; } } private static ASinImpl asinImpl = new ASinImpl(); /// <summary> /// /// </summary> /// <param name="v"></param> /// <param name="precision"></param> /// <returns></returns> public static double Asin(double v, int precision) { if (v < -1 || v > 1) { throw new ArgumentOutOfRangeException("v"); } return asinImpl.Calc(v, precision); } } }
通过一下优化:基本持平
class ASinImpl { private readonly int _precision; private double[] _quotieties = null; private long[] _factorials =null; public ASinImpl(int precision = 3) { _precision = precision; _quotieties = new double[precision + 1]; _factorials = new long[(precision + 2)*2 + 1]; Factorial(); ComputeQuotiety(); } public double Calc(double v) { unchecked { double retVal = 0; double vVal = 1; double q = v * v; unsafe { fixed (double* pq = _quotieties) { for (int i = 0; i < _precision; ++i) { if (i == 0) { vVal = v; //yield return ret; retVal += pq[i] * vVal; } else { vVal *= q; //yield return ret; retVal += pq[i] * vVal; } } return retVal; } } } } private void ComputeQuotiety() { unchecked { int precision = _quotieties.Length; for (int i = 0; i < precision; i++) { double up = _factorials[2*i]; double down = Math.Pow(Math.Pow(2, i)*_factorials[i], 2)*(2*i + 1); double quotiety = up/down; _quotieties[i] = quotiety; } } } private void Factorial() { unchecked { int precision = _factorials.Length ; long ret = 1; for (long v = 0; v < precision; ++v) { if (v == 0) { this._factorials[v] = 1; } else if (v == 1){ this._factorials[v] = 1; } else { ret *= v; this._factorials[v] = ret; } } } } }
QQ:273352165
evlon#126.com
转载请注明出处。