//■四元数プログラム例 // //@演算子%は、通常は剰余の意味を持つが、ここでは、前方向からの除算の // 意味に用いている。すなわち、四元数の場合、剰余という概念はないの // で演算子を借用している。 //A確認テスト結果はtextBox1に設定されるので、フォームにtextBox1を配置し、 // MultiLineプロパティをTrue, ScrollbarsプロパティをBothに設定して実行する。 // using System; using System.Collections.Generic; using System.ComponentModel; using System.Data; using System.Drawing; using System.Linq; using System.Text; using System.Windows.Forms; namespace quaternion { public partial class Form1 : Form { public Form1() { InitializeComponent(); } private void Form1_Load(object sender, EventArgs e) { Quater q, q2, qr, qr2, qe; string s = ""; double Pi6 = Math.PI / 6; q= new Quater(5); s += "\r\n" + q.ToString("0"); q = new Quater(1, 2, 3, 4); qe = new Quater(1, 0.4, 0.5, 0.3); q2 = new Quater(2, 3, 4, 5) ; s += "\r\n " + q.ToString("0"); qr = Quater.Sqrt(q) ; s += "\r\n sqr " + qr.ToString("0.0000"); qr2 = Quater.Pow2(qr) ; s += "\r\n pow2 " + qr2.ToString("0.0000"); qr2 = q + q2 ; s += "\r\n add " + qr2.ToString("0.0000"); qr2 = q - q2 ; s += "\r\n sub " + qr2.ToString("0.0000"); qr2 = q * q2 ; s += "\r\n mult " + qr2.ToString("0.0000"); qr = qr2 / q2 ; s += "\r\n div " + qr.ToString("0.0000"); qr = q % qr2 ; s += "\r\n div befor " + qr.ToString("0.0000"); qr2 = Quater.Rot(q, q2, Pi6) ; s += "\r\n rot " + qr2.ToString("0.0000"); qr = Quater.Rot(qr2, q2, -Pi6) ; s += "\r\n rot rev " + qr.ToString("0.0000"); qr2 = Quater.Exp(qe) ; s += "\r\n exp " + qr2.ToString("0.0000"); qr = Quater.Log(qr2) ; s += "\r\n log " + qr.ToString("0.0000"); qr2 = Quater.Pow(q, qe) ; s += "\r\n pow " + qr2.ToString("0.0000"); qr = Quater.Pow(qr2, Quater.Inv(qe)); s += "\r\n pow rev " + qr.ToString("0.0000"); textBox1.Text=s; } } public class Quater { int r; double a, b, c, d;// 球面上の集合のとき r = 1, 一般の四元数のとき r = 0 public Quater(int R, double A, double B, double C,double D){//四元数 r=R; a=A; b=B; c=C; d=D; } public Quater(double A, double B, double C,double D){//球面上の集合以外 r=0; a=A; b=B; c=C; d=D; } public Quater(double R){r = 1; a = 0; b = R; c = 0; d = 0;}//球面上の集合 public static string I2S(double d, string si, string Fm) { if(d==0) return ""; string S= d<0? " - ": " + "; if(Math.Abs(d)==1) return S + si; string S2 = (Fm == "" ? Math.Abs(d).ToString() : Math.Abs(d).ToString(Fm)); return S + S2 + " " + si +" "; } public override string ToString() {//四元数を文字列に変換 if (r == 1) return ("実部 " + a.ToString() + ", 虚部球面集合半径 " + b.ToString()); return a.ToString() + " " + I2S(b, "i", "") + I2S(c, "j", "") + I2S(b, "k", ""); } public string ToString(string F) {//四元数を文字列に変換 if (r == 1) return ("実部 " + a.ToString(F) + ", 虚部球面集合半径 " + b.ToString(F)); return a.ToString(F) + " " + I2S(b, "i", F) + I2S(c, "j", F) + I2S(b, "k", F); } public static string ToString(Quater p) { return p.ToString(); } public static string ToString(Quater p, string F){return p.ToString(F); } public static double Norm2(Quater x)//ノルムの２乗 { return x.a * x.a + x.b * x.b + x.c * x.c + x.d * x.d; } public static double Norm(Quater x)//ノルム { return Math.Sqrt(x.a * x.a + x.b * x.b + x.c * x.c + x.d * x.d); } public static Quater Conj(Quater x) //共役四元数 { return new Quater(x.r, x.a, -x.b, -x.c, -x.d); } public static Quater Sqrt(Quater x) { try { double bb = x.b * x.b + x.c * x.c + x.d * x.d; //虚部のノルムの２乗 if (bb < 1E-64) {//実数のとき（虚部の大きさ0） if (x.a < 0) return new Quater(Math.Abs(x.a));//負の実数のとき return D2Q(Math.Sqrt(x.a)); //正の実数のとき } //実数以外のとき、通常の四元数とする。 double aa = Math.Sqrt((x.a + Math.Sqrt(x.a * x.a + bb)) / 2); double a2 = aa * 2; return new Quater(aa, x.b / a2, x.c / a2, x.d / a2); } catch (Exception e) { MessageBox.Show(e.Message, "Error", MessageBoxButtons.OK, MessageBoxIcon.Error); return D2Q(0); } } public static int globe(int r1, int r2) { return (r1 != r2) ? 0 : r1; } public static Quater Pow2(Quater x) { return new Quater(x.a * x.a - x.b * x.b - x.c * x.c - x.d * x.d, 2 * x.a * x.b, 2 * x.a * x.c, 2 * x.a * x.d); } public static Quater Add(Quater x, Quater y) { return new Quater(globe(y.r, x.r), x.a + y.a, x.b + y.b, x.c + y.c, x.d + y.d); } public static Quater Sub(Quater x, Quater y) { return new Quater(globe(y.r, x.r), x.a - y.a, x.b - y.b, x.c - y.c, x.d - y.d); } public static bool IsReal(Quater x) { return (Math.Abs(x.b) < 1E-32 && Math.Abs(x.c) < 1E-32 && Math.Abs(x.d) < 1E-32); } public static Quater MltR(Quater x, double R) { return new Quater(x.r, R * x.a, R * x.b, R * x.c, R * x.d); } public static Quater Mlt(Quater x, Quater y) { if (IsReal(x)) return MltR(y, x.a);// xが実数のときyを実数倍 if (IsReal(y)) return MltR(x, y.a); // yが実数のときxを実数倍 return new Quater(x.a * y.a - x.b * y.b - x.c * y.c - x.d * y.d,//一般の四元数 x.a * y.b + x.b * y.a + x.c * y.d - x.d * y.c, x.a * y.c + x.c * y.a + x.d * y.b - x.b * y.d, x.a * y.d + x.d * y.a + x.b * y.c - x.c * y.b); } public static Quater Inv(Quater p) { try { double n2 = Norm2(p); return new Quater(p.r, p.a / n2, -p.b / n2, -p.c / n2, -p.d / n2);} catch (Exception e) { MessageBox.Show(e.Message, "Error", MessageBoxButtons.OK, MessageBoxIcon.Error); return D2Q(1); } } public static Quater Div(Quater p, Quater q){return p * Inv(q);}// p * q^(-1)　除算 public static Quater DivB(Quater p, Quater q){return Inv(p)*q;} // p^(-1) * q public static Quater Versor(Quater p) { try { double n = Norm(p); return new Quater(p.r, p.a / n, p.b / n, p.c / n, p.d / n); } catch (Exception e) { MessageBox.Show(e.Message, "Error", MessageBoxButtons.OK, MessageBoxIcon.Error); return D2Q(1); } } public Quater Vector(){ return new Quater(r, 0, b, c, d);}// ベクトル部分のみ public static Quater Vector(Quater p) {return p.Vector(); } public static Quater RotVector(Quater q, double angle)//回転ベクトル設定 { Quater V = Versor(Vector(q)) * D2Q(Math.Sin(angle/2)); return new Quater(Math.Cos(angle/2), V.b, V.c, V.d); } public static Quater Rot(Quater A, Quater axis, double angle)//軸（axis）周りの回転 { Quater q=RotVector(axis,angle); return q * A * Conj(q); } public static Quater Exp(Quater x) { Quater V=Vector(x); double nV=Norm(V), cosV=Math.Cos(nV), SV=0; if(Math.Abs(nV)>1E-32) SV = Math.Sin(nV) / nV; double EA = Math.Exp(x.a), ES = EA * SV; return new Quater(cosV * EA, x.b * ES, x.c * ES, x.d * ES); } public static Quater Log(Quater x) { Quater V=Vector(x); double nV=Norm(V); if(Math.Abs(nV)<1E-32) return D2Q(Math.Log(x.a)); double qq=Norm(x), TH=Math.Acos(x.a/qq)/nV; return new Quater(Math.Log(qq), x.b * TH, x.c * TH, x.d * TH); } public static bool Eq(Quater x,Quater y) { return (x.r == y.r && x.a == y.a && x.b == y.b && x.c == y.c && x.d == y.d); } public static Quater Pow(Quater q, Quater p){return Exp(p * Log(q));}// べき乗 public static Quater D2Q(double x) { return new Quater(x, 0, 0, 0); } public static Quater operator +(Quater x, Quater y) { return Add(x, y); } public static Quater operator +(double x, Quater y) { return Add(D2Q(x), y);} public static Quater operator +(Quater x, double y) { return Add(x, D2Q(y)); } public static Quater operator +(Quater x) { return x; } public static Quater operator -(Quater x, Quater y) { return Sub(x, y); } public static Quater operator -(double x, Quater y) { return Sub(D2Q(x), y); } public static Quater operator -(Quater x, double y) { return Sub(x, D2Q(y)); } public static Quater operator -(Quater x) { return new Quater(-x.a, -x.b, -x.c, -x.d); } public static Quater operator *(Quater x, Quater y) { return Mlt(x, y); } public static Quater operator *(double x, Quater y) { return Mlt(D2Q(x), y); } public static Quater operator *(Quater x, double y) { return Mlt(x, D2Q(y)); } public static Quater operator /(Quater x, Quater y) { return Div(x, y); } public static Quater operator /(double x, Quater y) { return Div(D2Q(x), y); } public static Quater operator /(Quater x, double y) { return Div(x, D2Q(y)); } public static Quater operator %(Quater x, Quater y) { return DivB(x, y); } public static Quater operator %(double x, Quater y) { return DivB(D2Q(x), y); } public static Quater operator %(Quater x, double y) { return DivB(x, D2Q(y)); } public static bool operator ==(Quater x, Quater y) { return Eq(x, y); } public static bool operator ==(double x, Quater y) { return Eq(D2Q(x), y); } public static bool operator ==(Quater x, double y) { return Eq(x, D2Q(y)); } public static bool operator !=(Quater x, Quater y) { return !Eq(x, y); } public static bool operator !=(double x, Quater y) { return !Eq(D2Q(x), y); } public static bool operator !=(Quater x, double y) { return !Eq(x, D2Q(y)); } public override bool Equals(object obj){ return base.Equals(obj);} public override int GetHashCode(){ return base.GetHashCode();} } }