The mysterious J function and the magic number 3

[Niels Henrik Bruun]

Hi
I'm trying to develop a sort of an iterator - generate_numbers.
When it gets a #.# value the integer value is repeated the decimal number of times.
And it can handle vectors as well.

So generate_numbers((1, 2.2, 3.3)) should return 1,2,2,3,3,3

But something strange happens when I want to repeat a number 3 times (I only get it repeated 2 times!!)
I've tried a lot of other things without luck eg assigning rep[r] before it's use.

:         mata clear

:         mata set matalnum on

:         //mata set strict on
: 
:         function generate_numbers(real rowvector base)
>         {
>                 real rowvector values, rep
>                 real scalar r
>                 
>                 values = floor(base)
>                 rep = (d=10 * (base - values)) :+ (d :== 0)
>                 out = J(1,0,0)
>                 for(r=1; r<=cols(values); r++) {
>                         out = out, J(1, rep[r], values[r])
>                 }
>                 return(out)
>         }

:         
:         x = generate_numbers((3.2, 6.3, 5, 7.4))

:         // Should have 6 3 times, but only 2 occur
:         1..cols(x) \ x
       1   2   3   4   5   6   7   8   9
    +-------------------------------------+
  1 |  1   2   3   4   5   6   7   8   9  |
  2 |  3   3   6   6   5   7   7   7   7  |
    +-------------------------------------+

:         // 6 4 times is no problem
:         x = generate_numbers((3.2, 6.4, 5, 7.4))

:         1..cols(x) \ x
        1    2    3    4    5    6    7    8    9   10   11
    +--------------------------------------------------------+
  1 |   1    2    3    4    5    6    7    8    9   10   11  |
  2 |   3    3    6    6    6    6    5    7    7    7    7  |
    +--------------------------------------------------------+

:         // Should have 5 3 times, but only 2 occur
:         x = generate_numbers((3.2, 6.4, 5.3, 7.4))

:         1..cols(x) \ x
        1    2    3    4    5    6    7    8    9   10   11   12
    +-------------------------------------------------------------+
  1 |   1    2    3    4    5    6    7    8    9   10   11   12  |
  2 |   3    3    6    6    6    6    5    5    7    7    7    7  |
    +-------------------------------------------------------------+

The only thing I can find is that it is the J function that is the problem.


The J function works in almost similar cases like:

: values = 1..5

: repeat = 5..1

: for(r=1; r<=cols(repeat); r++) J(1, repeat[r], values[r])
       1   2   3   4   5
    +---------------------+
  1 |  1   1   1   1   1  |
    +---------------------+
       1   2   3   4
    +-----------------+
  1 |  2   2   2   2  |
    +-----------------+
       1   2   3
    +-------------+
  1 |  3   3   3  |
    +-------------+
       1   2
    +---------+
  1 |  4   4  |
    +---------+
  5

:

Is there a workaround for this?

[Robert Picard]

It's a precision issue

:         base = (3.2, 6.3, 5, 7.4)

:         values = floor(base)

:         rep = (d=10 * (base - values)) :+ (d :== 0)

:         out = J(1,0,0)

:         for(r=1; r<=cols(values); r++) {
>                 printf("The reps is %21.0g\n", rep[r])
>                 out = out, J(1, rep[r], values[r])
>         }
The reps is  2.000000000000001776
The reps is  2.999999999999998224
The reps is                     1
The reps is  4.000000000000003553

:         out
       1   2   3   4   5   6   7   8   9
    +-------------------------------------+
  1 |  3   3   6   6   5   7   7   7   7  |
    +-------------------------------------+

[Niels Henrik Bruun]

Ah, of course
Why hadn't I thought of that?
Thank you very much