Square Root Function: That Scales The Kernel

The sqrt function below uses the sqrt_kernel and findExponent functions from previous posts. Note that my actual implementation is placed in a namespace so the system sqrt will not be picked up or conflict.

#include <boost/assert.hpp>  
   
 typedef double type_t;  
 const type_t epsilon = 1e-15;  
 const type_t ten   = 10;  
   
 type_t  
 pow10(type_t value){  
   return pow( ten, value );  
 }  
   
 type_t  
 sqrt(type_t n, type_t eps = epsilon) {  
   BOOST_ASSERT(n >= 0);  
   if(n == 0) return 0.0;  
   
   type_t ScalingExponent = findExponent(n);  
   n *= pow10(-ScalingExponent);  
   
   type_t result = sqrt_kernel(n,eps);  
   
   type_t Exponent = findExponent(result);  
   return result * pow10( -Exponent + (ScalingExponent/2) );  
 }  

Finding The Exponent In Scientific Notation - Harder Than It Seems

Was needing a method to find the exponent of a number (which means that value's scientific notation exponent.) I found simpler ways for several special cases but in general this seems to work for my needs.

The farther below UnitTest++ test code uses Table I.

Table I: Used to see what the exponent is,

Debug : value = 0.0001 = 1e-4 : log10 = -4
Debug : value = 0.0005 = 5e-4 : log10 = -3.30103
Debug : value = 0.001 = 1e-3 : log10 = -3
Debug : value = 0.005 = 5e-3 : log10 = -2.30103
Debug : value = 0.01 = 1e-2 : log10 = -2
Debug : value = 0.05 = 5e-2 : log10 = -1.30103
Debug : value = 0.1 = 1e-1 : log10 = -1
Debug : value = 0.5 = 5e-1 : log10 = -0.30103
Debug : value = 1 = 1e+0 : log10 = 0
Debug : value = 2 = 2e+0 : log10 = 0.30103
Debug : value = 10 = 1e+1 : log10 = 1
Debug : value = 20 = 2e+1 : log10 = 1.30103
Debug : value = 100 = 1e+2 : log10 = 2
Debug : value = 200 = 2e+2 : log10 = 2.30103
Debug : value = 1000 = 1e+3 : log10 = 3
Debug : value = 2000 = 2e+3 : log10 = 3.30103
Debug : value = 10000 = 1e+4 : log10 = 4

The function being tested is at lines 10-13. If uses the boost modf(). Note there is some C++11 in here.

1:  #include <UnitTest++.h> 
2:  #include <vector> 
3:  #include <boost/math/special_functions/modf.hpp> 
4:  namespace bm = boost::math; 
5:   
6:  SUITE(FIND_EXPONENT) { 
7:       typedef double type_t; 
8:       const type_t epsilon = 1e-15; 
9:   
10:       type_t findExponent(type_t value){ 
11:            type_t intPart; 
12:            return ( bm::modf<type_t>( log10(fabs(value)) , &intPart ) < 0 ) ? --intPart : intPart; 
13:       } 
14:       TEST(Function_findExponent_withTableI) { 
15:            std::vector<type_t> values = 
16:            { 1e-4,5e-4, 1e-3,5e-3, 1e-2,5e-2, 1e-1,5e-1, 1e+0,2e+0, 1e+1,2e+1, 1e+2,2e+2, 1e+3,2e+3, 1e+4,2e+4 }; 
17:            std::vector<type_t> exponents = 
18:            { -4,-4, -3,-3, -2,-2, -1,-1, 0,0, +1,+1, +2,+2, +3,+3, +4,+4 }; 
19:            CHECK_EQUAL(exponents.size(), values.size()); 
20:   
21:            for(unsigned i=0; i<values.size(); ++i){ 
22:                 values.at(i) = findExponent(values.at(i)); 
23:            } 
24:            CHECK_ARRAY_EQUAL(exponents,values,values.size()); 
25:       } 
26:  } 
27:   

Python For Lisp Programmers - "Interesting"

This page seems like a good reference for those interested in Python and/or Lisp. Here is a quote (the first paragraph) from the page:

'This is a brief introduction to Python for Lisp programmers. (Although it wasn't my intent, Python programmers have told me this page has helped them learn Lisp.) Basically, Python can be seen as a dialect of Lisp with "traditional" syntax (what Lisp people call "infix" or "m-lisp" syntax). One message on comp.lang.python said "I never understood why LISP was a good idea until I started playing with python." Python supports all of Lisp's essential features except macros, and you don't miss macros all that much because it does have eval, and operator overloading, and regular expression parsing, so some--but not all--of the use cases for macros are covered.'

Python For Lisp Programmers

Challenge: Random Number order of operation usage

Below is two examples. Each example should do the same thing but it turns out there is a difference in random number usage between the two. The second example differs from the first example in that the call to the linearSampling() function has been moved out of the circle.isInside() function argument list and replaced with the declaration of double X and double Y.

The two example codes do not give the same results on the same problem due to a difference in the order of random number usage. Can you tell why?

1:   //Example 1
2:   unsigned accepted = 0;
3:   for(unsigned n=0; n<nTests; ++n) {
4:     if( circle.isInside( linearSampling(Range[Shape::X].data(), rng) , linearSampling(Range[Shape::Y].data(), rng) ) )
5:       ++accepted;
6:   }
7:
8:   //Example 2
9:   unsigned accepted = 0;
10:  for(unsigned n=0; n<nTests; ++n) {
11:    double PointX=linearSampling(Range[Shape::X].data(), rng);
12:    double PointY=linearSampling(Range[Shape::Y].data(), rng);
13:    if( circle.isInside( PointX , PointY ) )
14:      ++accepted;
15:  }
16:
17:  double
18:  linearSampling(double* pCoord, rng_t& rng) {
19:    return pCoord[0] + ( (pCoord[1]-pCoord[0]) * rng() );
20:  }
21:

Using C++11 std::tuple

A snippet of C++ unit test code showing how to use std::tuple. Using enums to help index the std::get function that is used to retrieve tuple values. Note that std::tuple_size class is used to get the value (size) of the typedef Tuple_t.

 #include <UnitTest++.h>  
 #include <tuple>  
 SUITE(TUPLES){  
   TEST(Tuples) {  
     enum TupleIndex {MEAN, STDDEV, SIZE};  
     typedef std::tuple<double,double> Tuple_t;  
     CHECK_EQUAL( SIZE, std::tuple_size<Tuple_t>::value );  
   
     Tuple_t T1( 99.0, 0.5);  
     CHECK_CLOSE( 99.0, std::get<MEAN  >(T1), 1e-09);  
     CHECK_CLOSE(  0.5, std::get<STDDEV>(T1), 1e-09);  
 #if 0  
     //Out of Range ... caught at compile time.  
     //g++ error message not so helpful.  
     std::get<SIZE>(T1);  
 #endif  
   }  
 }