目的
SwiftのDecimalで数学の関数演算
SwiftのDecimal型で数学の関数演算（その１）
こちらの記事で数学の関数演算を行うプログラムを書いていた。そのとき、
$$\tan x= \frac{\sin x}{\cos x}$$
としていたら、精度が直接マクローリン展開で求めている$\sin$や$\cos$よりも精度が少し落ちた。
tanのマクローリン展開の公式
tanは以下の公式で求められる
$$\tan x = \sum_{n=1}^{\infty} \frac{B_{2n}(-4)^n(1-4^n)}{(2n)!}x^{2n-1}$$
この計算式をプログラムで利用するときに難しくしているのは$B_{2n}$のベルヌーイ数である。
このベルヌーイ数、一般項にすると$\sum$が2重になっており、計算量が多くなる。マクローリン展開の公式にある$\sum$と1つまとめられるが面倒。なので以前tanのプログラムを書いたときに楽してsinとcosから計算した。
今回、ベルヌーイ数をなんとかする目処が立ったので、この公式を利用して直接tanを計算してみた。
計算方法
ベルヌーイ数
とはいえ、ベルヌーイ数をプログラムで逐一求めるのはやはり簡単ではない。
しかし、精度的にtanのマクローリン展開では偶数番目のベルヌーイ数しか使わないし、$B_{200}$も使わないだろうと思い、そのくらいなら最初からベルヌーイ数をプログラムに打ち込んでおいてそれを使おうということになった。
ちなみに、ベルヌーイ数の求め方は以下を見てほしい。
ベルヌーイ数を分数で求める
（分数で求める計算にしたのは、こうすればとりあえず有効桁数を考えなくてもいいと思ったから）
各項の係数
ベルヌーイ数を最初から打ち込むなら、もういっそマクローリン展開の各項の係数を打ち込んでしまえと思った。係数は、ベルヌーイ数を求めるPythonプログラムを書き換えて求めた。そして、それを元につくった係数を返すSwiftのメソッドが以下のものである。
private func tan_coefficients(_ n: Int) -> Decimal {
switch n {
case 0: return Decimal(string: "0")!
case 1: return Decimal(string: "1")!
case 2: return Decimal(string: "3.33333333333333333333333333333333333333e-1")!
case 3: return Decimal(string: "1.33333333333333333333333333333333333333e-1")!
case 4: return Decimal(string: "5.39682539682539682539682539682539682540e-2")!
case 5: return Decimal(string: "2.18694885361552028218694885361552028219e-2")!
case 6: return Decimal(string: "8.86323552990219656886323552990219656886e-3")!
case 7: return Decimal(string: "3.59212803657248101692546136990581435026e-3")!
case 8: return Decimal(string: "1.45583438705131826824948518070211191904e-3")!
case 9: return Decimal(string: "5.90027440945585981378075993700021088754e-4")!
case 10: return Decimal(string: "2.39129114243552481485731458885112809321e-4")!
case 11: return Decimal(string: "9.69153795692945032559587500038880937786e-5")!
case 12: return Decimal(string: "3.92783238833168340533708080931233350952e-5")!
case 13: return Decimal(string: "1.59189050693289647407442798165700430553e-5")!
case 14: return Decimal(string: "6.45168921565543076319084231530255804463e-6")!
case 15: return Decimal(string: "2.61477115129075455426359425640974109737e-6")!
case 16: return Decimal(string: "1.05972683201046543509135539412477761900e-6")!
case 17: return Decimal(string: "4.29491107827380585482035128039727224262e-7")!
case 18: return Decimal(string: "1.74066189635716477798622923133027292194e-7")!
case 19: return Decimal(string: "7.05463694640096832521451854448055656984e-8")!
case 20: return Decimal(string: "2.85913666230525390833109453414164694357e-8")!
case 21: return Decimal(string: "1.15876444327988522001979320290917903963e-8")!
case 22: return Decimal(string: "4.69629539823090162881904347246767238237e-9")!
case 23: return Decimal(string: "1.90333683393127592131970659546238403474e-9")!
case 24: return Decimal(string: "7.71393363535906229047893307216272854899e-10")!
case 25: return Decimal(string: "3.12633954589208704577812001344089995544e-10")!
case 26: return Decimal(string: "1.26705769303054010600962880812234386249e-10")!
case 27: return Decimal(string: "5.13519140803936774007604951305429517699e-11")!
case 28: return Decimal(string: "2.08121468677004741989722719506097108605e-11")!
case 29: return Decimal(string: "8.43484541909433829384914818855870645342e-12")!
case 30: return Decimal(string: "3.41851408681115581403868029248568322278e-12")!
case 31: return Decimal(string: "1.38547157429484689939005446535238823491e-12")!
case 32: return Decimal(string: "5.61510479241468008325910118777644799993e-13")!
case 33: return Decimal(string: "2.27571625537287484554337456537684152713e-13")!
case 34: return Decimal(string: "9.22313058513953173504879603107797769815e-14")!
case 35: return Decimal(string: "3.73799403109673888502238947405794090155e-14")!
case 36: return Decimal(string: "1.51495191871486050751108883110796643109e-14")!
case 37: return Decimal(string: "6.13986886261681382451746533696988528729e-15")!
case 38: return Decimal(string: "2.48839512227627893824726042472915416004e-15")!
case 39: return Decimal(string: "1.00850855663540964870384529296132611377e-15")!
case 40: return Decimal(string: "4.08733122686901376190242207655547998999e-16")!
case 41: return Decimal(string: "1.65653295137862836727701494037447856591e-16")!
case 42: return Decimal(string: "6.71367517504871165607064488754700646572e-17")!
case 43: return Decimal(string: "2.72095006130445879103141466410559076155e-17")!
case 44: return Decimal(string: "1.10275952337223714588178395502199383518e-17")!
case 45: return Decimal(string: "4.46931600723741246634893781304951464441e-18")!
case 46: return Decimal(string: "1.81134555169976670766703620664393894875e-18")!
case 47: return Decimal(string: "7.34110701134014754553807111838677641765e-19")!
case 48: return Decimal(string: "2.97523860653619392921492989609583791722e-19")!
case 49: return Decimal(string: "1.20581878893050734904690845296116419207e-19")!
case 50: return Decimal(string: "4.88699947810437052169625792765250532302e-20")!
case 51: return Decimal(string: "1.98062628632408709016622156485819526210e-20")!
case 52: return Decimal(string: "8.02717598733937200765001813965609787179e-21")!
case 53: return Decimal(string: "3.25329188937030611595031956866135424751e-21")!
case 54: return Decimal(string: "1.31850953986006721125916811820106428292e-21")!
case 55: return Decimal(string: "5.34371788889344575062001458049629755636e-22")!
case 56: return Decimal(string: "2.16572728621367329750300435520822931587e-22")!
case 57: return Decimal(string: "8.77736208342709631440828592909340662994e-23")!
case 58: return Decimal(string: "3.55733086220083717062142133170547112850e-23")!
case 59: return Decimal(string: "1.44173189426242713612942108221539763739e-23")!
case 60: return Decimal(string: "5.84311928086315527456101200576160877141e-24")!
case 61: return Decimal(string: "2.36812704680168148272645722819731848943e-24")!
case 62: return Decimal(string: "9.59765741589520735830425747092434974659e-25")!
case 63: return Decimal(string: "3.88978403828865315672657090738009855536e-25")!
case 64: return Decimal(string: "1.57647009148991537547411178428131301522e-25")!
case 65: return Decimal(string: "6.38919262586010974913948804380105312352e-26")!
case 66: return Decimal(string: "2.58944223748416985965318963120096859901e-26")!
case 67: return Decimal(string: "1.04946140990143836024495875591479161918e-26")!
case 68: return Decimal(string: "4.25330688952681405256164016942044644830e-27")!
case 69: return Decimal(string: "1.72380035376386711484848849851727760915e-27")!
case 70: return Decimal(string: "6.98629968825742385528396636132349440140e-28")!
case 71: return Decimal(string: "2.83144061477735026377654296605129707184e-28")!
case 72: return Decimal(string: "1.14753965800891865921189020355489745354e-28")!
case 73: return Decimal(string: "4.65080305703948518684171840710561728678e-29")!
case 74: return Decimal(string: "1.88489948250657435169961472446836671826e-29")!
case 75: return Decimal(string: "7.63920986457584206322630627772664723545e-30")!
case 76: return Decimal(string: "3.09605514228418860428398709661919214923e-30")!
case 77: return Decimal(string: "1.25478388655272188777344516732156695068e-30")!
case 78: return Decimal(string: "5.08544754403418635484654998884796004877e-31")!
case 79: return Decimal(string: "2.06105425805025358458071826134632994067e-31")!
case 80: return Decimal(string: "8.35313827906927862322509137772487564807e-32")!
case 81: return Decimal(string: "3.38539943025367870320680181809775111926e-32")!
case 82: return Decimal(string: "1.37205070950159456714938511390305308291e-32")!
case 83: return Decimal(string: "5.56071207616043587916786180407146774513e-33")!
case 84: return Decimal(string: "2.25367171780326745485958130404835763519e-33")!
case 85: return Decimal(string: "9.13378743956350754286489209726672013556e-34")!
case 86: return Decimal(string: "3.70178461805636902657526156093255700422e-34")!
case 87: return Decimal(string: "1.50027679636188220953889816781117707439e-34")!
case 88: return Decimal(string: "6.08039283194081822004757735016215631004e-35")!
case 89: return Decimal(string: "2.46429039497051947686151655652059735907e-35")!
case 90: return Decimal(string: "9.98739278627428336064283895489871204858e-36")!
case 91: return Decimal(string: "4.04773783442502470979059680289066571219e-36")!
case 92: return Decimal(string: "1.64048635383126611341544034483999711043e-36")!
case 93: return Decimal(string: "6.64864076477147254412542294982917302839e-37")!
case 94: return Decimal(string: "2.69459260759737365918818258309735624519e-37")!
case 95: return Decimal(string: "1.09207724974263718942093784346277317742e-37")!
case 96: return Decimal(string: "4.42602238291171649664684266331137827319e-38")!
case 97: return Decimal(string: "1.79379930665638181520445814233093501434e-38")!
case 98: return Decimal(string: "7.26999475868917716285758069211738577864e-39")!
case 99: return Decimal(string: "2.94641789609591676456961619152728021094e-39")!
default: return .nan
}
}
プログラム
extension Decimal {
func tan() -> Decimal {
if isNaN { return .nan }
var angle = nomalize(unit: .pi)
if angle == .pi_two { return .nan }
var sign: Decimal = 1
if angle > .pi_two { angle = .pi - angle; sign = -1 }
let square = angle.mul(angle)
var x: Decimal = angle
var result: Decimal = tan_coefficients(1) * x
var next: Decimal = result
var i = 2
while true {
x = x.mul(square)
next += tan_coefficients(i).mul(x)
if next == result || next.isNaN {
break
}
result = next
i += 1
}
return result * sign
}
}
結果
let v = (Decimal.pi / 3)
print(v.sin() / v.cos())
> 1.73205080756887729352744634150587236692
print(v.tan())
> 1.7320508075688772935274463415058705858
あれ？
ちなみに
$$\tan \frac{\pi}{3} = \sqrt{3} = 1.73205080756887729352744634150587236694$$
なので精度が全然よろしくない。
結論
手間の割に精度が全然良くない。なので、tanを求めるプログラムを書くときはsin / cosを使うことをお勧めする。
ちなみに、なぜここまで精度が悪いのかはよくわかっていない。なぜ精度が悪いのかや、精度をよくする方法を知ってる人は教えてほしい。