MicroPythonで有効数字6桁以上の数値計算をする方法

　先日、Raspberry Pi Picoに太陽誘電のGPSモジュールGYSFDMAXBを接続して2地点間の距離を計算しようとしたところ、Raspberry Pi PicoのMicroPythonでは、float型は有効数字6桁以上では誤差が生じることがわかりました。まぁ、500円ちょいのマイコンですから、そんなものかもしれません。とはいえ、GPS関係では精密な計算が必要なので、有効数字6桁以上で計算する方法を調べてみました。

MicroPython向けのDecimalモジュール
DecimalNumberでヒュベニの式を計算してみる

MicroPython向けのDecimalモジュール

　PythonにはDecimalというモジュールがありますが、MicroPythonでは標準では使えません。探してみたところ、GitHubでMicroPython向けのDecimalモジュールが作成・公開されていました。
　使用方法は小数を文字列型(str)で指定し、DecimalNumberという型にします。DecimalNumberでは四則演算、べき乗、平方根、比較、指数関数、三角関数などが使え、円周率やネイピア数なども準備されています。上記のサイトで公開されているmpy_decimal.pyというコードを、Raspberry Pi Picoにコピーし、自分のソースコードから呼び出します。

DecimalNumberでヒュベニの式を計算してみる

　東京都庁の目の前の横断歩道にgoogle mapのカーソルを合わせてみると、緯度経度は都庁側が（35.68958, 139.69225）、都民広場側が（35.68962, 139.69243）でした。この2点の距離を国土地理院の測量計算サイトで計算すると、16.886（m）でした。同じ座標を先日公開したfloat型で計算するヒュベニの式に入れると、17.21543（m）となり、約2％の誤差が生じていました（小さい数を扱うときの計算誤差なので、もう少し長い距離の測量では誤差は小さい）。

　MicroPython向けのDecimalモジュールを利用して、以下のようなプログラムを作成しました。

from mpy_decimal import *

#ヒュベニの式で2点間の距離を求める
def hubenyDist(lat1,lon1,lat2,lon2):      
    #赤道半径（WGS84）
    Rx=DecimalNumber("6378137")
    #極半径（WGS84）
    Ry=DecimalNumber("6356752.314")
    #緯度の差（ラジアン）
    Dy = lat1-lat2
    if Dy < 0:
        Dy*=-1
    Dy=Dy/180*DecimalNumber.pi()
    #経度の差（ラジアン）
    Dx = lon1-lon2
    if Dx < 0:
        Dx*=-1
    Dx=Dx/180*DecimalNumber.pi()
    #平均緯度（ラジアン）
    P=(lat1+lat2)/2/180*DecimalNumber.pi()
    #離心率
    E=(Rx*Rx-Ry*Ry)/Rx/Rx
    E=E.square_root()
    W=1-E*E*P.sin()*P.sin()
    W=W.square_root()
    #子午線曲率半径
    M=Rx*(1-E*E)/W/W/W
    #卯酉線曲率半径
    N=Rx/W
    #2点間の距離
    D=Dy*Dy*M*M+Dx*Dx*N*N*P.cos()*P.cos()
    return D.square_root()

print(hubenyDist(DecimalNumber("35.68958"),DecimalNumber("139.69225"),DecimalNumber("35.68962"),DecimalNumber("139.69243")))

　計算結果は、16.8864883775618749となり、無事、ラズパイピコで国土地理院の測量計算サイトと同等程度の精度の計算ができました。良かった良かった。