1 // Copyright (c) 2014 Jamison Hope
2 // This is free software;  for terms and warranty disclaimer see ./COPYING.
3 
4 package gnu.math;
5 import java.io.*;
6 
7 /**
8  * A quaternion number using plain double values.
9  * @author Jamison Hope
10  */
11 public class DQuaternion extends Quaternion implements Externalizable {
12     double real;
13     double imag;
14     double jmag;
15     double kmag;
16 
DQuaternion()17     public DQuaternion() {
18     }
19 
DQuaternion(double real, double imag, double jmag, double kmag)20     public DQuaternion(double real, double imag, double jmag,
21                        double kmag) {
22         this.real = real;
23         this.imag = imag;
24         this.jmag = jmag;
25         this.kmag = kmag;
26     }
27 
re()28     public RealNum re () { return new DFloNum (real); }
doubleValue()29     public double doubleValue() { return real; }
im()30     public RealNum im () { return new DFloNum (imag); }
doubleImagValue()31     public double doubleImagValue () { return imag; }
jm()32     public RealNum jm () { return new DFloNum (jmag); }
doubleJmagValue()33     public double doubleJmagValue () { return jmag; }
km()34     public RealNum km () { return new DFloNum (kmag); }
doubleKmagValue()35     public double doubleKmagValue () { return kmag; }
36 
isExact()37     public boolean isExact () {
38         return false;
39     }
40 
toExact()41     public Quaternion toExact() {
42         return new CQuaternion(DFloNum.toExact(real),
43                                DFloNum.toExact(imag),
44                                DFloNum.toExact(jmag),
45                                DFloNum.toExact(kmag));
46     }
47 
equals(Object obj)48     public boolean equals(Object obj) {
49         if (obj == null || !(obj instanceof Quaternion))
50             return false;
51         Quaternion y = (Quaternion)obj;
52         return y.unit() == Unit.Empty
53             && (Double.doubleToLongBits(real)
54                 == Double.doubleToLongBits(y.reValue()))
55             && (Double.doubleToLongBits(imag)
56                 == Double.doubleToLongBits(y.imValue()))
57             && (Double.doubleToLongBits(jmag)
58                 == Double.doubleToLongBits(y.jmValue()))
59             && (Double.doubleToLongBits(kmag)
60                 == Double.doubleToLongBits(y.kmValue()));
61     }
62 
toString()63     public String toString() {
64         String reString = DFloNum.toString(real);
65 
66         if ((Double.doubleToLongBits(imag) == 0) &&
67             (Double.doubleToLongBits(jmag) == 0) &&
68             (Double.doubleToLongBits(kmag) == 0))
69             return reString;
70 
71         StringBuilder sbuf = new StringBuilder();
72         if (!reString.equals("0.0"))
73             sbuf.append(reString);
74 
75         if (Double.doubleToLongBits(imag) != 0) {
76             String imString = DFloNum.toString(imag);
77             char ch0 = imString.charAt(0);
78             if (ch0 != '-' && ch0 != '+')
79                 sbuf.append('+');
80             sbuf.append(imString);
81             sbuf.append('i');
82         }
83 
84         if (Double.doubleToLongBits(jmag) != 0) {
85             String jmString = DFloNum.toString(jmag);
86             char ch0 = jmString.charAt(0);
87             if (ch0 != '-' && ch0 != '+')
88                 sbuf.append('+');
89             sbuf.append(jmString);
90             sbuf.append('j');
91         }
92 
93         if (Double.doubleToLongBits(kmag) != 0) {
94             String kmString = DFloNum.toString(kmag);
95             char ch0 = kmString.charAt(0);
96             if (ch0 != '-' && ch0 != '+')
97                 sbuf.append('+');
98             sbuf.append(kmString);
99             sbuf.append('k');
100         }
101         return sbuf.toString();
102     }
103 
toString(int radix)104     public String toString(int radix) {
105         if (radix == 10)
106             return toString();
107         return "#d" + toString();
108     }
109 
neg()110     public final Numeric neg() {
111         return new DQuaternion(-real, -imag, -jmag, -kmag);
112     }
113 
add(Object y, int k)114     public Numeric add (Object y, int k) {
115         if (y instanceof Quaternion) {
116             Quaternion yq = (Quaternion)y;
117             if (yq.dimensions() != Dimensions.Empty)
118                 throw new ArithmeticException ("units mis-match");
119             return Quaternion.make(real + k * yq.reValue(),
120                                    imag + k * yq.imValue(),
121                                    jmag + k * yq.jmValue(),
122                                    kmag + k * yq.kmValue());
123         }
124         return ((Numeric)y).addReversed(this, k);
125     }
126 
mul(Object y)127     public Numeric mul (Object y) {
128         if (y instanceof Quaternion) {
129             Quaternion yq = (Quaternion)y;
130             if (yq.unit() == Unit.Empty) {
131                 double y_re = yq.reValue();
132                 double y_im = yq.imValue();
133                 double y_jm = yq.jmValue();
134                 double y_km = yq.kmValue();
135                 return Quaternion.make
136                     (real * y_re - imag * y_im - jmag * y_jm - kmag * y_km,
137                      real * y_im + imag * y_re + jmag * y_km - kmag * y_jm,
138                      real * y_jm - imag * y_km + jmag * y_re + kmag * y_im,
139                      real * y_km + imag * y_jm - jmag * y_im + kmag * y_re);
140             }
141             return Quaternion.times(this, yq);
142         }
143         return ((Numeric)y).mulReversed(this);
144     }
145 
div(Object y)146     public Numeric div(Object y) {
147         if (y instanceof Quaternion) {
148             Quaternion yq = (Quaternion) y;
149             return DQuaternion.div(real, imag, jmag, kmag,
150                                    yq.doubleValue(), yq.doubleImagValue(),
151                                    yq.doubleJmagValue(), yq.doubleKmagValue());
152         }
153         return ((Numeric)y).divReversed(this);
154     }
155 
hypot4(double w, double x, double y, double z)156     public static double hypot4(double w, double x, double y, double z) {
157         return Math.hypot(Math.hypot(w,x),Math.hypot(y,z));
158     }
159 
hypot3(double x, double y, double z)160     public static double hypot3(double x, double y, double z) {
161         return Math.hypot(Math.hypot(x,y),z);
162     }
163 
power(double x_re, double x_im, double x_jm, double x_km, double y_re, double y_im, double y_jm, double y_km)164     public static Quaternion power(double x_re, double x_im,
165                                    double x_jm, double x_km,
166                                    double y_re, double y_im,
167                                    double y_jm, double y_km) {
168         if (x_jm == 0.0 && x_km == 0.0 && y_jm == 0.0 && y_km == 0.0)
169             return DComplex.power(x_re, x_im, y_re, y_im);
170 
171         if (x_re == 0.0 && x_im == 0.0 && x_jm == 0.0 && x_km == 0.0)
172             if (y_re > 0.0)
173                 return DFloNum.valueOf(0.0);
174             else if (y_re == 0.0 && y_im == 0.0 &&
175                      y_jm == 0.0 && y_km == 0.0)
176                 return DFloNum.valueOf(1.0);
177 
178         // ln(x)
179         double qmag = hypot4(x_re, x_im, x_jm, x_km);
180         double vmag = hypot3(x_im, x_jm, x_km);
181 
182         double atv = Math.atan2(vmag, x_re) / vmag;
183 
184         double ln_r = Math.log(qmag);
185         double ln_i = atv * x_im;
186         double ln_j = atv * x_jm;
187         double ln_k = atv * x_km;
188 
189         // ln(x)*y
190         double p_r = ln_r * y_re - ln_i * y_im - ln_j * y_jm - ln_k * y_km;
191         double p_i = ln_r * y_im + ln_i * y_re + ln_j * y_km - ln_k * y_jm;
192         double p_j = ln_r * y_jm - ln_i * y_km + ln_j * y_re + ln_k * y_im;
193         double p_k = ln_r * y_km + ln_i * y_jm - ln_j * y_im + ln_k * y_re;
194 
195         // exp(ln(x)*y)
196         double pvmag = hypot3(p_i,p_j,p_k);
197         double sinpvmag = Math.sin(pvmag);
198         double expr = Math.exp(p_r);
199 
200         if (pvmag == 0.0 || sinpvmag == 0.0)
201             return DFloNum.valueOf(expr * Math.cos(pvmag));
202 
203         return Quaternion.make(expr * Math.cos(pvmag),
204                                expr * sinpvmag * p_i / pvmag,
205                                expr * sinpvmag * p_j / pvmag,
206                                expr * sinpvmag * p_k / pvmag);
207     }
208 
exp(double x_re, double x_im, double x_jm, double x_km)209     public static Quaternion exp(double x_re, double x_im,
210                                  double x_jm, double x_km) {
211         if (x_jm == 0.0 && x_km == 0.0)
212             return Complex.polar(Math.exp(x_re), x_im);
213 
214         double vmag = hypot3(x_im,x_jm,x_km);
215         double sinvmag = Math.sin(vmag);
216         double expr = Math.exp(x_re);
217         return Quaternion.make(expr * Math.cos(vmag),
218                                expr * sinvmag * x_im / vmag,
219                                expr * sinvmag * x_jm / vmag,
220                                expr * sinvmag * x_km / vmag);
221     }
222 
log(double x_re, double x_im, double x_jm, double x_km)223     public static Quaternion log(double x_re, double x_im,
224                                  double x_jm, double x_km) {
225         if (x_jm == 0.0 && x_km == 0.0)
226             return DComplex.log(x_re, x_im);
227 
228         double qmag = hypot4(x_re,x_im,x_jm,x_km);
229         double vmag = hypot3(x_im,x_jm,x_km);
230 
231         double atv = Math.atan2(vmag, x_re) / vmag;
232 
233         double r = Math.log(qmag);
234         double i = atv * x_im;
235         double j = atv * x_jm;
236         double k = atv * x_km;
237 
238         return Quaternion.make(r, i, j, k);
239     }
240 
div(double x_re, double x_im, double x_jm, double x_km, double y_re, double y_im, double y_jm, double y_km)241     public static Quaternion div(double x_re, double x_im,
242                                  double x_jm, double x_km,
243                                  double y_re, double y_im,
244                                  double y_jm, double y_km) {
245         if (x_jm == 0.0 && x_km == 0.0 && y_jm == 0.0 && y_km == 0.0)
246             return DComplex.div(x_re, x_im, y_re, y_im);
247 
248         double y_norm = y_re*y_re + y_im*y_im + y_jm*y_jm + y_km*y_km;
249 
250         // This computes (y^-1 * x), which is different from (x * y^-1).
251         double r = x_re*y_re + x_im*y_im + x_jm*y_jm + x_km*y_km;
252         double i = x_im*y_re - x_re*y_im + x_km*y_jm - x_jm*y_km;
253         double j = x_jm*y_re - x_re*y_jm + x_im*y_km - x_km*y_im;
254         double k = x_km*y_re - x_re*y_km + x_jm*y_im - x_im*y_jm;
255 
256         return Quaternion.make(r/y_norm, i/y_norm, j/y_norm, k/y_norm);
257     }
258 
sqrt(double x_re, double x_im, double x_jm, double x_km)259     public static Quaternion sqrt(double x_re, double x_im,
260                                   double x_jm, double x_km) {
261         if (x_jm == 0.0 && x_km == 0.0)
262             return DComplex.sqrt(x_re, x_im);
263 
264         double qmag = hypot4(x_re,x_im,x_jm,x_km);
265         double vmag = hypot3(x_im,x_jm,x_km);
266 
267         double t = Math.acos(x_re/qmag);
268 
269         double y_mag = Math.sqrt(qmag);
270         double s = Math.sin(t/2);
271 
272         return Quaternion.make(y_mag * Math.cos(t/2),
273                                y_mag * s * x_im / vmag,
274                                y_mag * s * x_jm / vmag,
275                                y_mag * s * x_km / vmag);
276     }
277 
sin(double x_re, double x_im, double x_jm, double x_km)278     public static Quaternion sin(double x_re, double x_im,
279                                  double x_jm, double x_km) {
280         if (x_jm == 0.0 && x_km == 0.0)
281             return DComplex.sin(x_re, x_im);
282 
283         double vmag = hypot3(x_im,x_jm,x_km);
284 
285         double r = Math.sin(x_re) * Math.cosh(vmag);
286         double v = Math.cos(x_re) * Math.sinh(vmag);
287 
288         return Quaternion.make(r, v * x_im/vmag, v * x_jm/vmag, v * x_km/vmag);
289     }
290 
cos(double x_re, double x_im, double x_jm, double x_km)291     public static Quaternion cos(double x_re, double x_im,
292                                  double x_jm, double x_km) {
293         if (x_jm == 0.0 && x_km == 0.0)
294             return DComplex.cos(x_re, x_im);
295 
296         double vmag = hypot3(x_im,x_jm,x_km);
297 
298         double r =  Math.cos(x_re) * Math.cosh(vmag);
299         double v = -Math.sin(x_re) * Math.sinh(vmag);
300 
301         return Quaternion.make(r, v * x_im/vmag, v * x_jm/vmag, v * x_km/vmag);
302     }
303 
tan(double x_re, double x_im, double x_jm, double x_km)304     public static Quaternion tan(double x_re, double x_im,
305                                  double x_jm, double x_km) {
306         if (x_jm == 0.0 && x_km == 0.0)
307             return DComplex.tan(x_re, x_im);
308 
309         double vmag = hypot3(x_im,x_jm,x_km);
310 
311         double sin_re = Math.sin(x_re);
312         double cos_re = Math.cos(x_re);
313         double sinh_v = Math.sinh(vmag);
314         double cosh_v = Math.cosh(vmag);
315         // tan = sin/cos
316         return DQuaternion.div(sin_re*cosh_v,
317                                cos_re*sinh_v*x_im/vmag,
318                                cos_re*sinh_v*x_jm/vmag,
319                                cos_re*sinh_v*x_km/vmag,
320                                cos_re*cosh_v,
321                                -sin_re*sinh_v*x_im/vmag,
322                                -sin_re*sinh_v*x_jm/vmag,
323                                -sin_re*sinh_v*x_km/vmag);
324     }
325 
326     /**
327      * @serialData Writes the real part, followed by the imaginary parts.
328      *   All are written as doubles (using writeDouble).
329      */
writeExternal(ObjectOutput out)330     public void writeExternal(ObjectOutput out) throws IOException {
331         out.writeDouble(real);
332         out.writeDouble(imag);
333         out.writeDouble(jmag);
334         out.writeDouble(kmag);
335     }
336 
readExternal(ObjectInput in)337     public void readExternal(ObjectInput in)
338         throws IOException, ClassNotFoundException {
339         real = in.readDouble();
340         imag = in.readDouble();
341         jmag = in.readDouble();
342         kmag = in.readDouble();
343     }
344 }
345