--- ./boost/math/octonion.hpp.org	2006-01-27 06:14:27.000000000 +0900
+++ ./boost/math/octonion.hpp	2010-02-25 21:30:26.000000000 +0900
@@ -4720,7 +4720,159 @@
                 return(pow(octonion<T>(1)/o,-n));
             }
         }
+         
         
+        template<typename T>
+        inline octonion  <T>                    pow(octonion  <T> const & q,
+                                                    T const & p)
+        {
+            // q^p = e^(p log q)
+            return(((abs(q) != static_cast<T>(0)) ? exp(p*log(q)) : octonion  <T>(0)));
+        }
+        
+        
+        template<typename T>
+        inline octonion  <T>                    pow(octonion  <T> const & q,
+                                                    octonion  <T> const & p)
+        {
+            // q^p = e^(p log q)
+            return(((abs(q) != static_cast<T>(0)) ? exp(p*log(q)) : octonion  <T>(0)));
+        }
+        
+        
+        template<typename T>
+        inline octonion  <T>                    sqrt(octonion  <T> const & q)
+        {
+            return(pow(q, static_cast<T>(0.5)));
+        }
+
+
+        // inverse transcendentals
+        // (please see the documentation)
+        
+        
+        template<typename T>
+        inline octonion  <T>                    log(octonion  <T> const & q)
+        {
+            using ::std::log;
+            using ::std::atan2;
+            T    z = abs(unreal(q));
+
+            // log(q) = log(|q|) + Vq/|Vq| atan(|Vq|/Sq)
+            return(log(abs(q)) + ((z != static_cast<T>(0)) ? unreal(q)/z*(atan2(z, real(q))) : octonion  <T>(0, (real(q)<0)?M_PI:0, 0, 0, 0, 0, 0, 0)));
+        }
+        
+        
+        template<typename T>
+        inline octonion  <T>                    acos(octonion  <T> const & q)
+        {
+            octonion  <T> u = (abs(unreal(q)) != static_cast<T>(0)) ? unreal(q)/abs(unreal(q)) : octonion  <T>(0, 1, 0, 0, 0, 0, 0, 0);
+#if 1
+            // acos(q) = 1/UVq log(q + sqrt(q^2 - 1))
+            //  = -UVq log(q + sqrt(q^2 - 1))
+            //return(-u*log(q + sqrt(pow(q, 2) - static_cast<T>(1))));
+            //  = -UVq log(q + sqrt(UVq^2(1 - q^2)))
+            //  = -UVq log(q + UVq sqrt((1 - q^2)))
+            //return(-u*log(q + (u)*sqrt(static_cast<T>(1) - pow(q, 2)))); // the signs of imagpart are different with complex ones when |Vq| == 0 && Sq < 0
+            octonion  <T> r = -u*log(q + (u)*sqrt(static_cast<T>(1) - pow(q, 2)));
+            //return(r.real() + r.unreal()*static_cast<T>((abs(unreal(q)) != static_cast<T>(0) || real(q) > 0) ? 1 : -1)); // the values of realpart are different |Vq| == 0 && when Sq > 1
+            return(((abs(unreal(q)) != static_cast<T>(0) || real(q) < 1) ? r.real() : 0) + r.unreal()*static_cast<T>((abs(unreal(q)) != static_cast<T>(0) || real(q) > 0) ? 1 : -1));
+#else
+            // acos(q) = 1/UVq acosh(q)
+            //return(static_cast<T>(1)/u*acosh(q)); // the signs of imagpart are different with complex ones when |Vq| == 0
+            octonion  <T> r = static_cast<T>(1)/u*acosh(q);
+            return(r.real() + r.unreal()*static_cast<T>((abs(unreal(q)) != static_cast<T>(0)) ? 1 : -1));
+#endif
+        }
+        
+        
+        template<typename T>
+        inline octonion  <T>                    asin(octonion  <T> const & q)
+        {
+            octonion  <T> u = (abs(unreal(q)) != static_cast<T>(0)) ? unreal(q)/abs(unreal(q)) : octonion  <T>(0, 1, 0, 0, 0, 0, 0, 0);
+#if 1
+            // asin(q) = 1/UVq log(UVq q + sqrt(1 - q^2))
+            //  = -UVq log(UVq q + sqrt(1 - q^2))
+            //return(-u*log((u)*q + sqrt(static_cast<T>(1) - pow(q, 2)))); // the signs of imagpart are different with complex ones when |Vq| == 0 && Sq > 0
+            octonion  <T> r = -u*log((u)*q + sqrt(static_cast<T>(1) - pow(q, 2)));
+            return(r.real() + r.unreal()*static_cast<T>((abs(unreal(q)) != static_cast<T>(0) || real(q) < 0) ? 1 : -1));
+#else
+            // asin(z) = 1/UVq asinh(UVq z)
+            //return(static_cast<T>(1)/u*asinh(u*q)); // the signs of imagpart are different with complex ones when |Vq| == 0 && Sq > 0
+            octonion  <T> r = static_cast<T>(1)/u*asinh(u*q);
+            return(r.real() + r.unreal()*static_cast<T>((abs(unreal(q)) != static_cast<T>(0) || real(q) < 0) ? 1 : -1));
+#endif
+        }
+        
+        
+        template<typename T>
+        inline octonion  <T>                    atan(octonion  <T> const & q)
+        {
+            static const T inf = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : static_cast<T>(HUGE_VAL);
+            if (q == unreal(q)/abs(unreal(q)) || q == abs(unreal(q))/unreal(q))
+                return(octonion  <T>((q.R_component_1()==static_cast<T>(0)) ? 0 : q.R_component_1()*inf,
+                                     (q.R_component_2()==static_cast<T>(0)) ? 0 : q.R_component_2()*inf,
+                                     (q.R_component_3()==static_cast<T>(0)) ? 0 : q.R_component_3()*inf,
+                                     (q.R_component_4()==static_cast<T>(0)) ? 0 : q.R_component_4()*inf,
+                                     (q.R_component_5()==static_cast<T>(0)) ? 0 : q.R_component_5()*inf,
+                                     (q.R_component_6()==static_cast<T>(0)) ? 0 : q.R_component_6()*inf,
+                                     (q.R_component_7()==static_cast<T>(0)) ? 0 : q.R_component_7()*inf,
+                                     (q.R_component_8()==static_cast<T>(0)) ? 0 : q.R_component_8()*inf));
+            octonion  <T> u = (abs(unreal(q)) != static_cast<T>(0)) ? unreal(q)/abs(unreal(q)) : octonion  <T>(0, 1, 0, 0, 0, 0, 0, 0);
+#if 1
+            // atan(q) = 1/(2UVq) log( (1+UVq.q)/(1-UVq.q) )
+            //  = -UVq/2 log( (1+UVq.q)/(1-UVq.q) )
+            //return(-u/static_cast<T>(2)*log((static_cast<T>(1)+(u)*q)/(static_cast<T>(1)-(u)*q)));
+            octonion  <T> d = static_cast<T>(1)-(u)*q, c = (static_cast<T>(1)+(u)*q)/d, r = -u/static_cast<T>(2)*log(c);
+            //return((abs(d) != static_cast<T>(0) && abs(c) != static_cast<T>(0)) ? r : octonion  <T>(0));
+            return((abs(d) != static_cast<T>(0) && abs(c) != static_cast<T>(0)) ? ((abs(unreal(q)) != static_cast<T>(0)) ? r : octonion  <T>(real(r))) : octonion  <T>(0));
+            // the signs of realpart are different with complex ones when Sq == 0 && Vq < -1
+#else
+            // atan(q) = -UVq atanh(UVq.q)
+            return(static_cast<T>(1)/u*atanh(u*q)); // the signs of realpart are different with complex ones when Sq == 0 && Vq < -1
+#endif
+        }
+        
+        
+        template<typename T>
+        inline octonion  <T>                    acosh(octonion  <T> const & q)
+        {
+            // acosh(q) = log(q + sqrt(q^2 - 1))
+            //return(log(q + sqrt(pow(q, 2) - static_cast<T>(1))));
+            octonion  <T> s = sqrt(pow(q, 2) - static_cast<T>(1)), r = log(q + s);
+            return((real(r) < 0) ? log(q - s) : r);
+        }
+        
+        
+        template<typename T>
+        inline octonion  <T>                    asinh(octonion  <T> const & q)
+        {
+            // asinh(q) = log(q + sqrt(q^2 + 1))
+            return(log(q + sqrt(pow(q, 2) + static_cast<T>(1)))); // the signs of realpart are different with complex ones when Sq == 0 && Vq > 1
+        }
+        
+        
+        template<typename T>
+        inline octonion  <T>                    atanh(octonion  <T> const & q)
+        {
+            static const T inf = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : static_cast<T>(HUGE_VAL);
+            if (q == -static_cast<T>(1) || q == static_cast<T>(1))
+                return(octonion  <T>((q.R_component_1()==static_cast<T>(0)) ? 0 : q.R_component_1()*inf,
+                                     (q.R_component_2()==static_cast<T>(0)) ? 0 : q.R_component_2()*inf,
+                                     (q.R_component_3()==static_cast<T>(0)) ? 0 : q.R_component_3()*inf,
+                                     (q.R_component_4()==static_cast<T>(0)) ? 0 : q.R_component_4()*inf,
+                                     (q.R_component_5()==static_cast<T>(0)) ? 0 : q.R_component_5()*inf,
+                                     (q.R_component_6()==static_cast<T>(0)) ? 0 : q.R_component_6()*inf,
+                                     (q.R_component_7()==static_cast<T>(0)) ? 0 : q.R_component_7()*inf,
+                                     (q.R_component_8()==static_cast<T>(0)) ? 0 : q.R_component_8()*inf));
+            // atanh(q) = 1/2 log( (1+q)/(1-q) )
+            //return(static_cast<T>(1)/static_cast<T>(2)*log((static_cast<T>(1)+q)/(static_cast<T>(1)-q)));
+            octonion  <T> d = (static_cast<T>(1)-q), c = (static_cast<T>(1)+q)/d, r = static_cast<T>(1)/static_cast<T>(2)*log(c);
+            //return((abs(d) != static_cast<T>(0) && abs(c) != static_cast<T>(0)) ? r : octonion  <T>(0));
+            return((abs(d) != static_cast<T>(0) && abs(c) != static_cast<T>(0)) ? ((real(q) != static_cast<T>(0)) ? r : unreal(r)) : octonion  <T>(0));
+        }
+
+       
         
         // helper templates for converting copy constructors (definition)
         
--- ./boost/math/quaternion.hpp.org	2005-05-04 22:40:34.000000000 +0900
+++ ./boost/math/quaternion.hpp	2010-02-25 21:30:26.000000000 +0900
@@ -1790,6 +1790,8 @@
             
             T    w = sinc_pi(z);
             
+            // e^q = e^(Sq+Vq) = e^(Sq) (cos(|Vq|) + UVq sin(|Vq|))
+            //  = e^(Sq) (cos(|Vq|) + Vq sin(|Vq|)/|Vq|)
             return(u*quaternion<T>(cos(z),
                 w*q.R_component_2(), w*q.R_component_3(),
                 w*q.R_component_4()));
@@ -1809,6 +1811,10 @@
             
             T    w = -sin(q.real())*sinhc_pi(z);
             
+            // cos(q) = cos(Sq+Vq)
+            //  = cos(Sq)cos(Vq) - sin(Sq)sin(Vq)
+            //  = cos(Sq)cosh(|Vq|) - UVq sin(Sq)sinh(|Vq|)
+            //  = cos(Sq)cosh(|Vq|) - Vq sin(Sq)sinh(|Vq|)/|Vq|
             return(quaternion<T>(cos(q.real())*cosh(z),
                 w*q.R_component_2(), w*q.R_component_3(),
                 w*q.R_component_4()));
@@ -1828,6 +1834,10 @@
             
             T    w = +cos(q.real())*sinhc_pi(z);
             
+            // sin(q) = sin(Sq+Vq)
+            //  = sin(Sq)cos(Vq) + cos(Sq)sin(Vq)
+            //  = sin(Sq)cosh(|Vq|) + UVq cos(Sq)sinh(|Vq|)
+            //  = sin(Sq)cosh(|Vq|) + Vq cos(Sq)sinh(|Vq|)/|Vq|
             return(quaternion<T>(sin(q.real())*cosh(z),
                 w*q.R_component_2(), w*q.R_component_3(),
                 w*q.R_component_4()));
@@ -1844,14 +1854,24 @@
         template<typename T>
         inline quaternion<T>                    cosh(quaternion<T> const & q)
         {
+            // sinh(q) = (e^q + e^(-q))/2
             return((exp(+q)+exp(-q))/static_cast<T>(2));
+            // cosh(q) = cosh(Sq+Vq)
+            //  = cosh(Sq)cosh(Vq) + sinh(Sq)sinh(Vq)
+            //  = cosh(Sq)cos(-|Vq|) - UVq sinh(Sq)sin(-|Vq|)
+            //  = cosh(Sq)cos(|Vq|) + Vq sinh(Sq)sin(|Vq|)/|Vq|
         }
         
         
         template<typename T>
         inline quaternion<T>                    sinh(quaternion<T> const & q)
         {
+            // sinh(q) = (e^q - e^(-q))/2
             return((exp(+q)-exp(-q))/static_cast<T>(2));
+            // sinh(q) = sinh(Sq+Vq)
+            //  = sinh(Sq)cosh(Vq) + cosh(Sq)sinh(Vq)
+            //  = sinh(Sq)cos(-|Vq|) - UVq cosh(Sq)sin(-|Vq|)
+            //  = sinh(Sq)cos(|Vq|) + Vq cosh(Sq)sin(|Vq|)/|Vq|
         }
         
         
@@ -1896,6 +1916,150 @@
         }
         
         
+        template<typename T>
+        inline quaternion<T>                    pow(quaternion<T> const & q,
+                                                    T const & p)
+        {
+            // q^p = e^(p log q)
+            return(((abs(q) != static_cast<T>(0)) ? exp(p*log(q)) : quaternion<T>(0)));
+        }
+        
+        
+        template<typename T>
+        inline quaternion<T>                    pow(quaternion<T> const & q,
+                                                    quaternion<T> const & p)
+        {
+            // q^p = e^(p log q)
+            return(((abs(q) != static_cast<T>(0)) ? exp(p*log(q)) : quaternion<T>(0)));
+        }
+        
+        
+        template<typename T>
+        inline quaternion<T>                    sqrt(quaternion<T> const & q)
+        {
+            return(pow(q, static_cast<T>(0.5)));
+        }
+
+
+        // inverse transcendentals
+        // (please see the documentation)
+        
+        
+        template<typename T>
+        inline quaternion<T>                    log(quaternion<T> const & q)
+        {
+            using ::std::log;
+            using ::std::atan2;
+            T    z = abs(unreal(q));
+
+            // log(q) = log(|q|) + Vq/|Vq| atan(|Vq|/Sq)
+            return(log(abs(q)) + ((z != static_cast<T>(0)) ? unreal(q)/z*(atan2(z, real(q))) : quaternion<T>(0, (real(q)<0)?M_PI:0, 0, 0)));
+        }
+        
+        
+        template<typename T>
+        inline quaternion<T>                    acos(quaternion<T> const & q)
+        {
+            quaternion<T> u = (abs(unreal(q)) != static_cast<T>(0)) ? unreal(q)/abs(unreal(q)) : quaternion<T>(0, 1, 0, 0);
+#if 1
+            // acos(q) = 1/UVq log(q + sqrt(q^2 - 1))
+            //  = -UVq log(q + sqrt(q^2 - 1))
+            //return(-u*log(q + sqrt(pow(q, 2) - static_cast<T>(1))));
+            //  = -UVq log(q + sqrt(UVq^2(1 - q^2)))
+            //  = -UVq log(q + UVq sqrt((1 - q^2)))
+            //return(-u*log(q + (u)*sqrt(static_cast<T>(1) - pow(q, 2)))); // the signs of imagpart are different with complex ones when |Vq| == 0 && Sq < 0
+            quaternion<T> r = -u*log(q + (u)*sqrt(static_cast<T>(1) - pow(q, 2)));
+            //return(r.real() + r.unreal()*static_cast<T>((abs(unreal(q)) != static_cast<T>(0) || real(q) > 0) ? 1 : -1)); // the values of realpart are different |Vq| == 0 && when Sq > 1
+            return(((abs(unreal(q)) != static_cast<T>(0) || real(q) < 1) ? r.real() : 0) + r.unreal()*static_cast<T>((abs(unreal(q)) != static_cast<T>(0) || real(q) > 0) ? 1 : -1));
+#else
+            // acos(q) = 1/UVq acosh(q)
+            //return(static_cast<T>(1)/u*acosh(q)); // the signs of imagpart are different with complex ones when |Vq| == 0
+            quaternion<T> r = static_cast<T>(1)/u*acosh(q);
+            return(r.real() + r.unreal()*static_cast<T>((abs(unreal(q)) != static_cast<T>(0)) ? 1 : -1));
+#endif
+        }
+        
+        
+        template<typename T>
+        inline quaternion<T>                    asin(quaternion<T> const & q)
+        {
+            quaternion<T> u = (abs(unreal(q)) != static_cast<T>(0)) ? unreal(q)/abs(unreal(q)) : quaternion<T>(0, 1, 0, 0);
+#if 1
+            // asin(q) = 1/UVq log(UVq q + sqrt(1 - q^2))
+            //  = -UVq log(UVq q + sqrt(1 - q^2))
+            //return(-u*log((u)*q + sqrt(static_cast<T>(1) - pow(q, 2)))); // the signs of imagpart are different with complex ones when |Vq| == 0 && Sq > 0
+            quaternion<T> r = -u*log((u)*q + sqrt(static_cast<T>(1) - pow(q, 2)));
+            return(r.real() + r.unreal()*static_cast<T>((abs(unreal(q)) != static_cast<T>(0) || real(q) < 0) ? 1 : -1));
+#else
+            // asin(z) = 1/UVq asinh(UVq z)
+            //return(static_cast<T>(1)/u*asinh(u*q)); // the signs of imagpart are different with complex ones when |Vq| == 0 && Sq > 0
+            quaternion<T> r = static_cast<T>(1)/u*asinh(u*q);
+            return(r.real() + r.unreal()*static_cast<T>((abs(unreal(q)) != static_cast<T>(0) || real(q) < 0) ? 1 : -1));
+#endif
+        }
+        
+        
+        template<typename T>
+        inline quaternion<T>                    atan(quaternion<T> const & q)
+        {
+            static const T inf = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : static_cast<T>(HUGE_VAL);
+            if (q == unreal(q)/abs(unreal(q)) || q == abs(unreal(q))/unreal(q))
+                return(quaternion<T>((q.R_component_1()==static_cast<T>(0)) ? 0 : q.R_component_1()*inf,
+                                     (q.R_component_2()==static_cast<T>(0)) ? 0 : q.R_component_2()*inf,
+                                     (q.R_component_3()==static_cast<T>(0)) ? 0 : q.R_component_3()*inf,
+                                     (q.R_component_4()==static_cast<T>(0)) ? 0 : q.R_component_4()*inf));
+            quaternion<T> u = (abs(unreal(q)) != static_cast<T>(0)) ? unreal(q)/abs(unreal(q)) : quaternion<T>(0, 1, 0, 0);
+#if 1
+            // atan(q) = 1/(2UVq) log( (1+UVq.q)/(1-UVq.q) )
+            //  = -UVq/2 log( (1+UVq.q)/(1-UVq.q) )
+            //return(-u/static_cast<T>(2)*log((static_cast<T>(1)+(u)*q)/(static_cast<T>(1)-(u)*q)));
+            quaternion<T> d = static_cast<T>(1)-(u)*q, c = (static_cast<T>(1)+(u)*q)/d, r = -u/static_cast<T>(2)*log(c);
+            //return((abs(d) != static_cast<T>(0) && abs(c) != static_cast<T>(0)) ? r : quaternion<T>(0));
+            return((abs(d) != static_cast<T>(0) && abs(c) != static_cast<T>(0)) ? ((abs(unreal(q)) != static_cast<T>(0)) ? r : quaternion<T>(real(r))) : quaternion<T>(0));
+            // the signs of realpart are different with complex ones when Sq == 0 && Vq < -1
+#else
+            // atan(q) = -UVq atanh(UVq.q)
+            return(static_cast<T>(1)/u*atanh(u*q)); // the signs of realpart are different with complex ones when Sq == 0 && Vq < -1
+#endif
+        }
+        
+        
+        template<typename T>
+        inline quaternion<T>                    acosh(quaternion<T> const & q)
+        {
+            // acosh(q) = log(q + sqrt(q^2 - 1))
+            //return(log(q + sqrt(pow(q, 2) - static_cast<T>(1))));
+            quaternion<T> s = sqrt(pow(q, 2) - static_cast<T>(1)), r = log(q + s);
+            return((real(r) < 0) ? log(q - s) : r);
+        }
+        
+        
+        template<typename T>
+        inline quaternion<T>                    asinh(quaternion<T> const & q)
+        {
+            // asinh(q) = log(q + sqrt(q^2 + 1))
+            return(log(q + sqrt(pow(q, 2) + static_cast<T>(1)))); // the signs of realpart are different with complex ones when Sq == 0 && Vq > 1
+        }
+        
+        
+        template<typename T>
+        inline quaternion<T>                    atanh(quaternion<T> const & q)
+        {
+            static const T inf = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : static_cast<T>(HUGE_VAL);
+            if (q == -static_cast<T>(1) || q == static_cast<T>(1))
+                return(quaternion<T>((q.R_component_1()==static_cast<T>(0)) ? 0 : q.R_component_1()*inf,
+                                     (q.R_component_2()==static_cast<T>(0)) ? 0 : q.R_component_2()*inf,
+                                     (q.R_component_3()==static_cast<T>(0)) ? 0 : q.R_component_3()*inf,
+                                     (q.R_component_4()==static_cast<T>(0)) ? 0 : q.R_component_4()*inf));
+            // atanh(q) = 1/2 log( (1+q)/(1-q) )
+            //return(static_cast<T>(1)/static_cast<T>(2)*log((static_cast<T>(1)+q)/(static_cast<T>(1)-q)));
+            quaternion<T> d = (static_cast<T>(1)-q), c = (static_cast<T>(1)+q)/d, r = static_cast<T>(1)/static_cast<T>(2)*log(c);
+            //return((abs(d) != static_cast<T>(0) && abs(c) != static_cast<T>(0)) ? r : quaternion<T>(0));
+            return((abs(d) != static_cast<T>(0) && abs(c) != static_cast<T>(0)) ? ((real(q) != static_cast<T>(0)) ? r : unreal(r)) : quaternion<T>(0));
+        }
+
+
+        
         // helper templates for converting copy constructors (definition)
         
         namespace detail