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