'■Formデザイン
'  TextBox1にテスト結果を設定するので、FormにTextBox1を配置し、
'  プロパティMultiLine=True、ScrollBars=Bothを設定しておくこと
'
Public Class Form1
    Public Structure quaternion
        Public r As Integer '通常の四元数のとき0, 球面集合のとき1
        Public a As Double                           '実数部
        Public b As Double, c As Double, d As Double '虚数部
    End Structure
    '四元数設定
    Function Quater(r, a, b, c, d) As quaternion
        Quater.r = r : Quater.a = a
        Quater.b = b : Quater.c = c : Quater.d = d
    End Function
    Function Qnorm2(x As quaternion) As Double 'ノルムの二乗の計算
        Return x.a ^ 2 + x.b ^ 2 + x.c ^ 2 + x.d ^ 2
    End Function
    Function Qnorm(x As quaternion) As Double 'ノルムの
        Return Math.Sqrt(Qnorm2(x))
    End Function
    Function Qconj(x As quaternion) As quaternion '共役四元数の計算
        Return Quater(x.r, x.a, -x.b, -x.c, -x.d)
    End Function
    Function PMstring(V As Double, S As String) As String
        If V = 1 Then Return " + " & S
        If V = -1 Then Return " - " & S
        Dim VV As Double = Math.Abs(V)
        If V > 0 Then Return " + " & VV & S
        If V < 0 Then Return " - " & VV & S
        Return ""
    End Function
    Function QtoString(x As quaternion) As String '文字列への変換
        Dim S As String
        If x.r = 1 Then Return "実部 " & x.a & ", 球面集合半径 " & x.b
        If x.a <> 0 Then S = x.a Else S = ""
        S = S & PMstring(x.b, "i")
        S = S & PMstring(x.c, "j")
        S = S & PMstring(x.d, "k")
        Return S
    End Function
    Function PMstringF(V As Double, S As String, FM As String) As String
        If V = 1 Then Return " + " & S
        If V = -1 Then Return " - " & S
        Dim VV As Double = Math.Abs(V)
        If V > 0 Then Return " + " & Format(VV, FM) & S
        If V < 0 Then Return " - " & Format(VV, FM) & S
        Return ""
    End Function
    Function QtoStringF(x As quaternion, FM As String) As String '文字列への変換
        Dim S As String
        If x.r = 1 Then
            Return "実部 " & Format(x.a, FM) & ", 虚部 球面集合半径 " & Format(x.b, FM)
        Else
            If x.a <> 0 Then S = Format(x.a, FM) Else S = ""
            S = S & PMstringF(x.b, "i", FM)
            S = S & PMstringF(x.c, "j", FM)
            S = S & PMstringF(x.d, "k", FM)
            Return S
        End If
    End Function
    Function Qsqr(x As quaternion) As quaternion '四元数の平方根
        Dim a, b As Double
        b = x.b ^ 2 + x.c ^ 2 + x.d ^ 2 '虚部のノルムの二乗
        If b < 1.0E-64 Then   '実数のとき
            Qsqr.a = 0 : Qsqr.b = 0 : Qsqr.c = 0 : Qsqr.d = 0
            If x.a < 0 Then   '負の実数のとき球面体の集合とする
                Return Quater(1, Math.Sqrt(Math.Abs(x.a)), 0, 0, 0)
            End If           '０または正の実数のとき実数の平方根
            Return Quater(0, Math.Sqrt(x.a), 0, 0, 0)
        Else     '負の実数以外のとき通常の四元数とする
            a = Math.Sqrt((x.a + Math.Sqrt(x.a ^ 2 + b)) / 2) * 2
            Return Quater(0, a, x.b / a, x.c / a, x.d / a)
        End If
    End Function
    Function Qpow2(x As quaternion) As quaternion '四元数の平方
        Qpow2.r = 0                                     '四元数とする
        Qpow2.a = x.a ^ 2 - x.b ^ 2 - x.c ^ 2 - x.d ^ 2 '実部
        Qpow2.b = 2 * x.a * x.b                         '虚部 i
        Qpow2.c = 2 * x.a * x.c                         '     j
        Qpow2.d = 2 * x.a * x.d                         '     k
    End Function
    '四元数の加算
    Function Qadd(x As quaternion, y As quaternion) As quaternion
        Qadd.r = x.r : If x.r <> y.r Then Qadd.r = 0
        Qadd.a = x.a + y.a
        Qadd.b = x.b + y.b : Qadd.c = x.c + y.c : Qadd.d = x.d + y.d
    End Function
    '四元数の減算
    Function Qsub(x As quaternion, y As quaternion) As quaternion
        Qsub.r = x.r : If x.r <> y.r Then Qsub.r = 0
        Qsub.a = x.a - y.a
        Qsub.b = x.b - y.b
        Qsub.c = x.c - y.c
        Qsub.d = x.d - y.d
    End Function
    '四元数の乗算
    Function Qmult(x As quaternion, y As quaternion) As quaternion
        If Math.Abs(x.b) < 1.0E-32 And Math.Abs(x.c) < 1.0E-32 _
                    And Math.Abs(x.d) < 1.0E-32 Then
            Qmult.r = y.r : Qmult.a = x.a * y.a
            Qmult.b = x.a * y.b
            Qmult.c = x.a * y.c
            Qmult.d = x.a * y.d
        ElseIf Math.Abs(y.b) < 1.0E-32 And Math.Abs(y.c) < 1.0E-32 _
                    And Math.Abs(y.d) < 1.0E-32 Then
            Qmult.r = x.r : Qmult.a = y.a * x.a
            Qmult.b = y.a * x.b : Qmult.c = y.a * x.c : Qmult.d = y.a * x.d
        Else
            Qmult.r = 0
            Qmult.a = x.a * y.a - x.b * y.b - x.c * y.c - x.d * y.d
            Qmult.b = x.a * y.b + x.b * y.a + x.c * y.d - x.d * y.c
            Qmult.c = x.a * y.c + x.c * y.a + x.d * y.b - x.b * y.d
            Qmult.d = x.a * y.d + x.d * y.a + x.b * y.c - x.c * y.b
        End If
    End Function
    Function Qinv(p As quaternion) As quaternion '四元数の逆数
        Dim n2 As Double : n2 = Qnorm2(p)
        On Error GoTo ERRQinv
        Qinv.r = p.r : Qinv.a = p.a / n2
        Qinv.b = -p.b / n2 : Qinv.c = -p.c / n2 : Qinv.d = -p.d / n2
QinvRet:
        Exit Function
ERRQinv:
        MsgBox("**Error** Zero Devide Or Overflow")
        Resume QinvRet
    End Function
    '四元数の除算（p・q^-1)
    Function Qdev(p As quaternion, q As quaternion) As quaternion
        Qdev = Qmult(p, Qinv(q))
    End Function
    '四元数の除算（q^-1・p)
    Function QdevBefor(q As quaternion, p As quaternion) As quaternion
        QdevBefor = Qmult(Qinv(q), p)
    End Function
    '単位四元数
    Function Qversor(q As quaternion) As quaternion
        Dim n As Double : n = Qnorm(q)
        Qversor.r = q.r : Qversor.a = q.a / n
        Qversor.b = q.b / n : Qversor.c = q.c / n : Qversor.d = q.d / n
    End Function
    'ベクトル部分のみ
    Function Qvect(q As quaternion) As quaternion
        Qvect = Quater(q.r, 0, q.b, q.c, q.d)
    End Function
    '回転ベクトルの設定
    Function QvectRot(q As quaternion, Ang As Double) As quaternion
        Dim sinA, cosA As Double
        sinA = Math.Sin(Ang / 2) : cosA = Math.Cos(Ang / 2)
        Dim V As quaternion
        V = Qmult(Qversor(Qvect(q)), Quater(0, sinA, 0, 0, 0))
        QvectRot = Quater(0, cosA, V.b, V.c, V.d)
    End Function
    '回転
    Function QRotation(A As quaternion, q As quaternion, Ang As Double) As quaternion
        Dim R As quaternion : R = QvectRot(q, Ang)
        QRotation = Qmult(Qmult(R, A), Qconj(R))
    End Function
    Function Qexp(q As quaternion) As quaternion '指数
        Dim V As quaternion
        Dim EA As Double : Dim EASV As Double
        Dim nV As Double : Dim SV As Double : Dim cosV As Double
        V = Qvect(q) : nV = Qnorm(V) : cosV = Math.Cos(nV)
        SV = 0 : If Math.Abs(nV) > 1.0E-32 Then SV = Math.Sin(nV) / nV
        EA = Math.Exp(q.a) : EASV = EA * SV
        Qexp = Quater(0, cosV * EA, V.b * EASV, V.c * EASV, V.d * EASV)
    End Function
    Function Qln(q As quaternion) As quaternion '対数
        Dim V As quaternion : Dim nV As Double
        Dim qq As Double : Dim TH As Double
        V = Qvect(q) : nV = Qnorm(V)
        If Math.Abs(nV) < 1.0E-32 Then
            Qln = Quater(0, Math.Log(q.a), 0, 0, 0)
        Else
            qq = Qnorm(q) : TH = Math.Acos(q.a / qq) / nV
            Qln = Quater(0, Math.Log(qq), V.b * TH, V.c * TH, V.d * TH)
        End If
    End Function
    Function Qpow(q As quaternion, p As quaternion) As quaternion 'べき乗
        Qpow = Qexp(Qmult(p, Qln(q)))
    End Function
    Private Sub Form1_Load(sender As System.Object, e As System.EventArgs) _
                                                        Handles MyBase.Load
        Dim A As quaternion, B As quaternion
        Dim D As quaternion, C As quaternion
        Dim X As quaternion
        Dim S As String = "", Pi6 As Double = Math.PI / 6
        A = Quater(0, 5, 0, 0, 0)
        S = vbCrLf & " 実数の表示 " & QtoString(A)
        A = Quater(0, 1, 2, 3, 4)
        X = Quater(0, 1, 0.4, 0.5, 0.3)
        B = Quater(0, 2, 3, 4, 5)
        S = S & vbCrLf & " X = " & QtoString(X)
        S = S & vbCrLf & " A = " & QtoString(A)
        S = S & vbCrLf & " B = " & QtoString(B)
        D = Qsqr(A)
        S = S & vbCrLf & " C=sqr(A)  " & QtoStringF(D, "#0.0000")
        C = Qpow2(D)
        S = S & vbCrLf & " C^2       " & QtoString(C)
        C = Qadd(A, B)
        S = S & vbCrLf & "  A+B      " & QtoString(C)
        C = Qsub(A, B)
        S = S & vbCrLf & "  A-B      " & QtoString(C)
        C = Qmult(A, B)
        S = S & vbCrLf & " C= A*B    " & QtoString(C)
        D = Qdev(C, B)
        S = S & vbCrLf & "  C *1/B   " & QtoString(D)
        D = QdevBefor(A, C)
        S = S & vbCrLf & "  1/A *C   " & QtoString(D)
        C = QRotation(A, B, Pi6)
        S = S & vbCrLf & " rot A     " & QtoStringF(C, "#0.0000")
        D = QRotation(C, B, -Pi6)
        S = S & vbCrLf & " rot rev   " & QtoString(D)
        C = Qexp(X)
        S = S & vbCrLf & " C=exp(X)  " & QtoStringF(C, "#0.0000")
        D = Qln(C)
        S = S & vbCrLf & " log (C)   " & QtoString(D)
        C = Qpow(A, X)
        S = S & vbCrLf & " C= A^X    " & QtoStringF(C, "#0.0000")
        D = Qpow(C, Qinv(X))
        S = S & vbCrLf & " C^(1/X)   " & QtoString(D)
        TextBox1.Text = S
    End Sub
End Class