summaryrefslogtreecommitdiffstats log msg author committer range
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180  \section{The Pythagorean sum} In this section, we analyse the computation of the Pythagorean sum, which is used in two different computations, that is the norm of a complex number and the 2-norm of a vector of real values. In the first part, we briefly present the mathematical formulas for these two computations. We then present the naïve algorithm based on these mathematical formulas. In the second part, we make some experiments in Scilab and compare our naïve algorithm with the \emph{abs} and \emph{norm} Scilab primitives. In the third part, we analyse why and how floating point numbers must be taken into account when the Pythagorean sum is to compute. \subsection{Theory} \subsection{Experiments} % TODO : compare both abs and norm. \lstset{language=Scilab} \lstset{numbers=left} \lstset{basicstyle=\footnotesize} \lstset{keywordstyle=\bfseries} \begin{lstlisting} // Straitforward implementation function mn2 = mynorm2(a,b) mn2 = sqrt(a^2+b^2) endfunction // With scaling function mn2 = mypythag1(a,b) if (a==0.0) then mn2 = abs(b); elseif (b==0.0) then mn2 = abs(a); else if (abs(b)>abs(a)) then r = a/b; t = abs(b); else r = b/a; t = abs(a); end mn2 = t * sqrt(1 + r^2); end endfunction // With Moler & Morrison's // At most 7 iterations are required. function mn2 = mypythag2(a,b) p = max(abs(a),abs(b)) q = min(abs(a),abs(b)) //index = 0 while (q<>0.0) //index = index + 1 //mprintf("index = %d, p = %e, q = %e\n",index,p,q) r = (q/p)^2 s = r/(4+r) p = p + 2*s*p q = s * q end mn2 = p endfunction function compare(x) mprintf("Re(x)=%e, Im(x)=%e\n",real(x),imag(x)); p = abs(x); mprintf("%20s : %e\n","Scilab",p); p = mynorm2(real(x),imag(x)); mprintf("%20s : %e\n","Naive",p); p = mypythag1(real(x),imag(x)); mprintf("%20s : %e\n","Scaling",p); p = mypythag2(real(x),imag(x)); mprintf("%20s : %e\n","Moler & Morrison",p); endfunction // Test #1 : all is fine x = 1 + 1 * %i; compare(x); // Test #2 : more difficult when x is large x = 1.e200 + 1 * %i; compare(x); // Test #3 : more difficult when x is small x = 1.e-200 + 1.e-200 * %i; compare(x); \end{lstlisting} \begin{verbatim} *************************************** Example #1 : simple computation with Scilab 5.1 x(1)=1.000000e+000, x(2)=1.000000e+000 Scilab : 1.414214e+000 Naive : 1.414214e+000 Scaling : 1.414214e+000 Moler & Morrison : 1.414214e+000 *************************************** Example #2 : with large numbers ? Scilab : Inf Naive : Inf Scaling : 1.000000e+200 Moler & Morrison : 1.000000e+200 *************************************** Example #3 : with small numbers ? x(1)=1.000000e-200, x(2)=1.000000e-200 Scilab : 0.000000e+000 Naive : 0.000000e+000 Scaling : 1.414214e-200 Moler & Morrison : 1.414214e-200 *************************************** > Conclusion : Scilab is naive ! Octave 3.0.3 *************************************** octave-3.0.3.exe:29> compare(x); *************************************** x(1)=1.000000e+000, x(2)=1.000000e+000 Octave : 1.414214e+000 Naive : 1.414214e+000 Scaling : 1.414214e+000 Moler & Morrison : 1.414214e+000 *************************************** x(1)=1.000000e+200, x(2)=1.000000e+000 Octave : 1.000000e+200 Naive : Inf Scaling : 1.000000e+200 Moler & Morrison : 1.000000e+200 *************************************** octave-3.0.3.exe:33> compare(x) x(1)=1.000000e-200, x(2)=1.000000e-200 Octave : 1.414214e-200 Naive : 0.000000e+000 Scaling : 1.414214e-200 Moler & Morrison : 1.414214e-200 *************************************** > Conclusion : Octave is not naive ! With complex numbers. *************************************** Re(x)=1.000000e+000, Im(x)=1.000000e+000 Scilab : 1.414214e+000 Naive : 1.414214e+000 Scaling : 1.414214e+000 Moler & Morrison : 1.414214e+000 *************************************** Re(x)=1.000000e+200, Im(x)=1.000000e+000 Scilab : 1.000000e+200 Naive : Inf Scaling : 1.000000e+200 Moler & Morrison : 1.000000e+200 *************************************** Re(x)=1.000000e-200, Im(x)=1.000000e-200 Scilab : 1.414214e-200 Naive : 0.000000e+000 Scaling : 1.414214e-200 Moler & Morrison : 1.414214e-200 *************************************** > Conclusion : Scilab is not naive ! \end{verbatim} \subsection{Explanations} \subsection{References} The paper by Moler and Morrisson 1983 \cite{journals/ibmrd/MolerM83} gives an algorithm to compute the Pythagorean sum $a\oplus b = \sqrt{a^2 + b^2}$ without computing their squares or their square roots. Their algorithm is based on a cubically convergent sequence. The BLAS linear algebra suite of routines \cite{900236} includes the SNRM2, DNRM2 and SCNRM2 routines which conpute the euclidian norm of a vector. These routines are based on Blue \cite{355771} and Cody \cite{Cody:1971:SEF}. In his 1978 paper \cite{355771}, James Blue gives an algorithm to compute the Euclidian norm of a n-vector $\|x\| = \sqrt{\sum_{i=1,n}x_i^2}$. The exceptionnal values of the \emph{hypot} operator are defined as the Pythagorean sum in the IEEE 754 standard \cite{P754:2008:ISF,ieee754-1985}. The \emph{\_\_ieee754\_hypot(x,y)} C function is implemented in the Fdlibm software library \cite{fdlibm} developed by Sun Microsystems and available at netlib. This library is used by Matlab \cite{matlab-hypot} and its \emph{hypot} command.