Loren on the Art of MATLAB

May 24th, 2006

Programming for Multiple Datatypes

How do you write M-files in MATLAB that you would like to have work correctly whether the inputs are double precision or single? By correctly, I mean that the program gives an answer appropriate to the precision of the input.

Contents

Special Constants and Functions in MATLAB to Help with Precision

There are some fairly new functions in MATLAB to help with precision, and some updates to existing functions. Let me first list some "constants" that might be interesting for working with different precision data:

  • eps
  • zeros, ones, Inf

There are several places in the documentation that are also particularly relevant:

And Cleve wrote this article for News and Notes in 1996 on floating point.

Problem Description

Let's write a program that calculates a result either in single or double precision, depending on the input. Fibonacci numbers are a sequence of numbers with some interesting characteristics. As you calculate more numbers in the sequence, the ratio between successive Fibonacci numbers approaches the golden mean or golden ratio. In other words,

My current web search found 42,000,000 pages for golden mean, 9,490,000 for golden ratio - wow!

Golden Mean

The golden mean can expressed as

What is its value, in double precision?

format long
GoldenMeanDouble = (1+sqrt(5))/2
GoldenMeanDouble =

   1.61803398874989

What about single precision?

GoldenMeanSingle = (1+sqrt(single(5)))/2
GoldenMeanSingle =

   1.6180340

How Many Fibonacci Terms Should We Compute?

How many terms in the Fibonacci sequence should we compute so that the ratio between successive terms doesn't change our estimate of the golden mean to within the appropriate precision? First, stop and take a guess for both double and single precision. Let's start with double, using the filter code from the previous blog to calculate the Fibonacci numbers. 

format
nf = 100;
x = [1 zeros(1,nf-1)];
a = [1 -1 -1];
b = 1;
FibDouble = filter(b, a, x);

Calculate a few of the ratios and compare to the golden mean.

RatioDouble = FibDouble(2:11)./FibDouble(1:10)
ResidualDouble = FibDouble(2:11)./FibDouble(1:10) - GoldenMeanDouble
RatioDouble =

  Columns 1 through 6 

    1.0000    2.0000    1.5000    1.6667    1.6000    1.6250

  Columns 7 through 10 

    1.6154    1.6190    1.6176    1.6182


ResidualDouble =

  Columns 1 through 6 

   -0.6180    0.3820   -0.1180    0.0486   -0.0180    0.0070

  Columns 7 through 10 

   -0.0026    0.0010   -0.0004    0.0001

Calculate Number of Terms

Here are the number of terms needed for double precision. Most people expect the number of required terms to be considerably higher. 

numdouble = goldFibDouble
ResidDouble = FibDouble(numdouble)/FibDouble(numdouble-1) - GoldenMeanDouble
numdouble =

    41


ResidDouble =

     0

Double Algorithm

Let's look at the algorithm. In this M-file, I do not worry about efficiency for calculating the Fibonacci sequence, but

type goldFibDouble
function nterms = goldFibDouble
% goldFib Number of Fibonacci terms to reach golden mean ratio.

fcurrent = 1;
fnext = 1;
goldenMean = (1+sqrt(5))/2;
tol = eps(goldenMean);  % Calculate increment for next closest number.
nterms = 2;
while abs(fnext/fcurrent - goldenMean) >= tol
    nterms = nterms + 1;
    temp  = fnext;
    fnext = fnext + fcurrent;
    fcurrent = temp;
end

Floating Point Algorithm

What do I need to change in goldFibDouble to allow it to also calculate the number of terms for single precision? There are really only a few things I need to change. First, I probably want to calculate the Fibonacci numbers themselves in single, and the golden mean and tolerance. Finally, I need to change the function so I can specify what precision I want the output to be. Let's now look at goldFib where I made the necessary changes. 

dbtype goldFib
1     function nterms = goldFib(dtype)
2     % goldFib Number of Fibonacci terms to reach golden mean ratio.
3     
4     fcurrent = ones(dtype);
5     fnext = fcurrent;
6     goldenMean = (1+sqrt(cast(5,dtype)))/2;
7     tol = eps(goldenMean);
8     nterms = 2;
9     while abs(fnext/fcurrent - goldenMean) >= tol
10        nterms = nterms + 1;
11        temp  = fnext;
12        fnext = fnext + fcurrent;
13        fcurrent = temp;
14    end

Number of Terms for Single and Double

Now let's see how many terms we need for both single and double.

numSingle = goldFib('single')
numDouble = goldFib('double')
numSingle =

    19


numDouble =

    41

You'll notice that on line 7 of goldFib, we calculate the relative accuracy for the golden mean in the correct precision. You can see the values here. 

TolSingle = eps(GoldenMeanSingle)
TolDouble = eps(GoldenMeanDouble)
TolSingle =

  1.1921e-007


TolDouble =

  2.2204e-016

Tools in MATLAB for Writing "Generic" Floating Point Programs

MATLAB has tools so you can create programs that work correctly on both single and double precision. You can see with this example that we only needed to modify a small number of lines and terms in the M-file in order to convert the file from handling double precision only to being able to handle both single and double precision calculations appropriately. Do you need to do something similar? Is it because you have large datasets and use single precision to help manage the memory? Let me know.


Published with MATLAB® 7.2

By Loren Shure

08:08 UTC | Posted in Best Practice, New Feature, Robustness | Permalink |

You can follow any responses to this entry through the RSS 2.0 feed. You can skip to the end and leave a response. Pinging is currently not allowed.

9 Responses to “Programming for Multiple Datatypes”

  1. David Jones replied on July 6th, 2006 at 15:07 UTC :

    Hi Loren,

    I found your article about single / double data types useful.
    Since you invited comments, … I am running up against the current
    memory limitations in MATLAB and I might need to switch from double
    to single to try to reduce memory consumption. Do you know when Mathworks
    plans to upgrade MATLAB to have an address space that is larger than 32 bits?
    I am running Mac OS X 10.4 on a G5 with 4 GB of RAM, but poor MATLAB can only address 2 GB at a time.

    thanks,
    David Jones

  2. Leslie Mehrez replied on July 7th, 2006 at 09:51 UTC :

    Hi David,
    We don’t have any plans to add 64-bit support for PowerMacs. We are focusing our Mac resources on Intel Mac support. Once Apple formally announces their 64-bit roadmap for Intel Macs we’ll be able to start planning for 64-bit Mac support.

    Leslie Mehrez
    Product Manager, Platforms, Installation & Licensing
    The MathWorks

  3. Andrew Catellier replied on July 25th, 2006 at 10:05 UTC :

    I have an interesting problem. I’m trying to make MATLAB behave like a C program using 32 bit floats. I’ve written the code I need and I’m in the process of verifying it. I realized that the accuracy of my code is near perfect when I cast all my variables to singles. However, I have many, many variables, and was wondering if there was a way to force MATLAB to create and use single precision variables. This, of course, would allow me to avoid the sheer pain of converting the variables in nearly 60 .m files. Any insight or pointers would be of great interest to me.

    Thanks for your time,

    Andrew

  4. Loren replied on July 25th, 2006 at 10:51 UTC :

    Andrew-

    There is no way to force MATLAB to make all variables single.

    In MATLAB itself, (i.e., I’m NOT covering toolboxes with this statement), we have tried to make sure the functions that are most likely to be used with singles work correctly with them. If you write code “generically,” such as we have done in many of our math functions (especially look in the ones in $MATLABROOT/toolbox/matlab/funfun), then singles should propogate through correctly.

    It would be interesting in your application to see what happens if you start off with singles and then see where the first double creeps in. If it’s in a MATLAB function (one we supply), it would be a candidate for updating to be more generic.

  5. Andrew Catellier replied on July 25th, 2006 at 16:18 UTC :

    So far the MATLAB functions have been very well behaved. Of course, hindsight is 20/20, and I should have imported all of my arrays and constants as singles in the first place. However, in the interest of time, I was wondering if there was a way I could salvage some of my work. I now know that my first idea does not exist, and I will find another way to finish the work.

    Thanks for your insight,
    Andrew

  6. Andrew Catellier replied on July 28th, 2006 at 17:01 UTC :

    As I made my way further into my program, it turns out that the function “bitshift” is not defined for use with singles. For now I can recast them, but I thought this might be useful for future type expansion.

    Regards,

    Andrew

  7. StephenLL replied on January 12th, 2007 at 12:50 UTC :

    I think this is where MATLAB shines vs other languages. When we were on the topic of scalar expansion I wrote code in c++ to do the expansion. My problem is that I am not a strong c/c++ programmer and my code only handled double precision numbers. Yes, it is possible to write c++ to handle all the possibilities, but I found MATLAB to naturally handle them.

    One of the main reasons I think twice about recoding bottlenecks in c or c++ is the issue with dealing with multiple precisions and data types. This advise I pass on to anyone looking to use MEX.

    Thankfully the MATLAB language evolves over time and those bottlenecks go away.

    Stephen

  8. Chris Mejia replied on October 19th, 2008 at 09:14 UTC :

    I find myself in a situation where I wrote a MEX-file for efficiency and now need a single-precision version. My current solution is to have the gateway routine select between two computational routines, and the two computational routines are identical except I’ve replaced all of the “double” keywords with “float”. Is there a more elegant solution?

    Thanks in advance,
    –Chris

  9. Loren replied on October 19th, 2008 at 12:38 UTC :

    Chris-

    I think you can use function pointers in Fortran somehow but my Fortran is rusty. I don’t know if the newest Fortran language versions have a concept of templates. If so, they may help. Otherwise, I think you’re stuck doing something similar to what you are currently doing.

    Maybe someone else who is more current on Fortran can add suggestions.

    –Loren

Leave a Reply

Wrap code fragments inside <pre> tags, like this:

<pre class="code">
a = magic(3);
sum(a)
</pre>

If you have a "<" character in your code, either follow it with a space or replace it with "&lt;" (including the semicolon).

About

Loren Shure works on design of the MATLAB language at The MathWorks. She writes here about once a week on MATLAB programming and related topics.

Subscribe


Recent Comments

  • Jun: I totally can not believe it, Loren. You are really helpful. Thank you so much, MATLAB master!
  • Loren: Wow folks- Always lots of interest when there’s a quickie to try out! I will only make 2 general...
  • Loren: Jun- ismember is your friend here: >> [aa,ind] = ismember(Array2,Arra y1) aa = 1 1 1 1 1 1 1 ind = 1 2 1 4 4 3...
  • Dan: I like the first way better than the second way. Combining the arrays into one and running any is nice, although...
  • James Myatt: How about I = (a == 0 | b == 0); a(I) = []; b(I) = [];
  • Tunc: Hello Loren, love your blog because of such inspiring and challenging comments to such ’small’...
  • Pekka Kumpulainen: Here is my tradeoff. I usually want to keep the original variables as they are most probably...
  • Iain: Followup: Of course, to allow NaNs (counting them as non-zero): mask = (a~=0) & (b~=0); The mask says “a...
  • Matt Fig: I would usually go with something like this: y = a&b; x = a(y); y = b(y); But I was surprised to find...
  • kk: c=all([a;b]) a(c) a(b)

Blog Admin

These postings are the author's and don't necessarily represent the opinions of The MathWorks.