Hello,
I was modifying "for" loop implementation in Freemat and realized that neither Freemat 3.5 nor Octave 3.0 match the behavior of "for" loops in matlab. If you run the script below you will see that the first loop performs 20 iterations in matlab, but only 19 in Octave. "Loop 2" below performs 20 iterations in both systems. For the second test you should change "if 1" to "if 0". Plot "for loop vs. my implementation" In matlab is constant zero. However, in Octave the plot shows the difference of more than 20*eps. A side note, in matlab "for loop vs. my implementation" is constant zero, but "vectorized vs. my implementation" shows difference of up to 4*eps. This means that "k=[first:step:last]" and "for k=first:step:last" use different algorithms and produce different results!!!! Compatible "for" loop implementation will go in the next release of Freemat, and Freemat already has almost correct implementation of "k=[first:step:last]" construct (a fairly complicated algorithm). Hope this helps Octave too. Cheers, Eugene ================================================================================================================== clear if 1 first1 = 1.e16; last1 = (1.e16+2*pi); step1 = (1./pi); else first1 = 0; last1 = 5.8; step1 = .1; end % loop 1 % matlab for loop output n1=0; for k=first1:step1:last1 mat_k(n1+1)=k; n1=n1+1; end % loop 2 % matching internal implementation nst = round((last1-first1)/step1)+ 1; %precompute number of steps for n = 0:(nst-1) my_k(n+1)=first1+n*step1; end [ length(my_k) length(mat_k) ] dif = my_k-mat_k; k = [first1:step1:last1]; dif1 = my_k - k; subplot(1,2,1);plot(dif./eps); title('for loop vs. my implementation');xlabel('step');ylabel('difference / eps') subplot(1,2,2);plot(dif1./eps); title('vectorized vs. implementation');xlabel('step');ylabel('difference / eps') |
On 11-Feb-2008, EI wrote:
| Hello, | | I was modifying "for" loop implementation in Freemat and realized that | neither Freemat 3.5 nor Octave 3.0 match the behavior of "for" loops in | matlab. | | If you run the script below you will see that the first loop performs 20 | iterations in matlab, but only 19 in Octave. "Loop 2" below performs 20 | iterations in both systems. | | For the second test you should change "if 1" to "if 0". Plot "for loop vs. | my implementation" In matlab is constant zero. However, in Octave the plot | shows the difference of more than 20*eps. | | A side note, in matlab "for loop vs. my implementation" is constant zero, | but "vectorized vs. my implementation" shows difference of up to 4*eps. This | means that "k=[first:step:last]" and "for k=first:step:last" use different | algorithms and produce different results!!!! In Octave, we were using base + i * increment to compute the values when doing something like this: [ base:increment:limit ] which converts a range (stored as base, increment, and limit values only) to an array with all values expanded, but the for loop implmentation was using tmp = base; for (i = 0; i < nelem; i++, tmp += increment) ... which obviously fails for cases like this when base + increment == base. I've fixed this problem so both calculations are performed using multiplication. I think this change fixes the precision problem you see. | Compatible "for" loop implementation will go in the next release of Freemat, | and Freemat already has almost correct implementation of | "k=[first:step:last]" construct (a fairly complicated algorithm). Would you care to share the code for this? Octave already goes to some lengths to get the number of elements right. Thanks, jwe |
John W. Eaton wrote:
> On 11-Feb-2008, EI wrote: > > | Hello, > | > | I was modifying "for" loop implementation in Freemat and realized that > | neither Freemat 3.5 nor Octave 3.0 match the behavior of "for" loops in > | matlab. > | > | If you run the script below you will see that the first loop performs 20 > | iterations in matlab, but only 19 in Octave. "Loop 2" below performs 20 > | iterations in both systems. > | > | For the second test you should change "if 1" to "if 0". Plot "for loop vs. > | my implementation" In matlab is constant zero. However, in Octave the plot > | shows the difference of more than 20*eps. > | > | A side note, in matlab "for loop vs. my implementation" is constant zero, > | but "vectorized vs. my implementation" shows difference of up to 4*eps. This > | means that "k=[first:step:last]" and "for k=first:step:last" use different > | algorithms and produce different results!!!! > > In Octave, we were using > > base + i * increment > > to compute the values when doing something like this: > > [ base:increment:limit ] > > which converts a range (stored as base, increment, and limit values > only) to an array with all values expanded, but the for loop > implmentation was using > > tmp = base; > for (i = 0; i < nelem; i++, tmp += increment) > ... > > which obviously fails for cases like this when base + increment == base. > I've fixed this problem so both calculations are performed using > multiplication. I think this change fixes the precision problem you > see. > > | Compatible "for" loop implementation will go in the next release of Freemat, > | and Freemat already has almost correct implementation of > | "k=[first:step:last]" construct (a fairly complicated algorithm). > > Would you care to share the code for this? > > Octave already goes to some lengths to get the number of elements > right. > > Thanks, > > jwe > > Octave already has good matlab compatibility here, and does quite complex for loops like looping over matrices and cell arrays. For example a = eye (3); for i = a, disp(i); printf("\n"); end 1 0 0 0 1 0 0 0 1 Precision issues are never going to go away. Change the underlying FPU and the result will probably change. All we can hope for is to be inconsistent in the same manner as Matlab itself :-) D. -- David Bateman [hidden email] Motorola Labs - Paris +33 1 69 35 48 04 (Ph) Parc Les Algorithmes, Commune de St Aubin +33 6 72 01 06 33 (Mob) 91193 Gif-Sur-Yvette FRANCE +33 1 69 35 77 01 (Fax) The information contained in this communication has been classified as: [x] General Business Information [ ] Motorola Internal Use Only [ ] Motorola Confidential Proprietary |
Free forum by Nabble | Edit this page |