项目中需要快速求解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;
} }
}
}
}