Re: Calculation issue with octave (data unexpectedly go to zero at some point)

[Juan Pablo Carbajal]

On Fri, Feb 28, 2014 at 3:11 AM, Dmitri A. Sergatskov
<address@hidden> wrote:
>
> Back to the original problem.
> I think if we pre-scale argument such that X is about MAXINT, we get the max
> possible accuracy in this situation. In this particular case
>
>
>
> Gamma = 1.62e7;
> duration = 10/Gamma;
> dt = 0.0025/Gamma;
> t   = 0:dt:duration;
> sc = 2^31 / duration ;
> y = mod (sc*t, sc*0.2/Gamma)/sc;
>
> [a,b] = find(abs(y) < eps);
>
> sum(a)
> ans =  51
>
> b =
>
>  Columns 1 through 11:
>
>       1     81    161    241    321    401    481    561    641    721
> 801
>
>  Columns 12 through 22:
>
>     881    961   1041   1121   1201   1281   1361   1441   1521   1601
> 1681
>
>  Columns 23 through 33:
>
>    1761   1841   1921   2001   2081   2161   2241   2321   2401   2481
> 2561
>
>  Columns 34 through 44:
>
>    2641   2721   2801   2881   2961   3041   3121   3201   3281   3361
> 3441
>
>  Columns 45 through 51:
>
>    3521   3601   3681   3761   3841   3921   4001
>
>
> Perhaps we should include such pre-scaling in the code of rem and mod
> if we try for matlab compatibility...
>
> Dmitri.
> --
>
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://mailman.cae.wisc.edu/listinfo/help-octave
>

Hi I am opening a bug for mod, flagged as incompatible result. I get
the following comaring betwene octave 3.8.1 and matlab2008b
Gamma = 1.62e7;
duration = 10/Gamma;
dt = 0.0025/Gamma;
t   = 0:dt:duration;
y = mod (t, 0.2/Gamma);
find (y==0,3,'first')

octave
1   241   401

r2008b
1    81   161

Reading the help of mod in matlab it says
MOD(x,y) is x - n.*y where n = floor(x./y) if y ~= 0.  If y is not an
    integer and the quotient x./y is within roundoff error of an integer,
    then n is that integer.

So indeed matlab is giving a result considerig roundoff error, I
assume they do something like
function m = mod_ml(x,y)
  if fix(y) != y
   err      = abs (x./y - round(x./y)) < sqrt (eps);
   m       = mod (x,y);
   m(err) = 0;
  endif
endfunction

Shall I fix our mod function?

@Renaud: you can use mod_ml for the moment and tell us if it fixes our issues