Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

mod for floats

3 views
Skip to first unread message

Maarten van Reeuwijk

unread,
Feb 18, 2005, 8:38:09 AM2/18/05
to
Hi,

I am trying to get Maple to calculate things like 10.3 mod 1 = 0.3 but this
function only seems to work for integers. Is there a corresponding function
to easily do this?

Thx Maarten
--
===================================================================
Maarten van Reeuwijk Thermal and Fluids Sciences
Phd student dept. of Multiscale Physics
www.ws.tn.tudelft.nl Delft University of Technology

dev...@math.udel.edu

unread,
Feb 18, 2005, 9:07:42 AM2/18/05
to
frac(10.3) = 0.3

Maarten van Reeuwijk

unread,
Feb 18, 2005, 11:23:26 AM2/18/05
to
dev...@math.udel.edu wrote:

> frac(10.3) = 0.3

Thanks! But that only gives correct results for values > 0 as

> frac(-1.7) = -0.7

What I am trying to do is to define an iterative map such as

x_n+1 = f(x_n) and e.g. f(x) = a * x mod 1

so that 0 < x < 1.

frac would work for a > 0 but not for a < 0.

Any further suggestions?

dev...@math.udel.edu

unread,
Feb 18, 2005, 11:34:01 AM2/18/05
to
mod1:= x-> x-floor(x):
mod1(-1.7);
0.3

Joe Riel

unread,
Feb 18, 2005, 9:27:43 PM2/18/05
to
dev...@math.udel.edu writes:

The floor is a bit expensive. This will speed it up by about 10x for
x > 0, 5x for x < 0:

fastmod1 := proc(x)
local y;
if x >= 0 then x - trunc(x)
else
y := x - trunc(x) + 1;
y - trunc(y)
fi
end:


Joe

Alec Mihailovs

unread,
Feb 19, 2005, 10:48:57 PM2/19/05
to
"Joe Riel" <jo...@k-online.com> wrote in message
news:87oeehs...@k-online.com...

It can be also slightly improved,

fastmod2:=proc(x) local y;
y:=x-trunc(x);
if y<0 then y+1 else y fi end;

Alec

Alec Mihailovs

unread,
Feb 19, 2005, 11:32:25 PM2/19/05
to

"Alec Mihailovs" <al...@mihailovs.com> wrote in message
news:JMTRd.8258$A23....@news02.roc.ny...

>
> It can be also slightly improved,
>
> fastmod2:=proc(x) local y;
> y:=x-trunc(x);
> if y<0 then y+1 else y fi end;

It seems as if it is faster with 0. and 1. instead of 0 and 1,

fastmod2:=proc(x) local y;
y:=x-trunc(x);

if y<0. then y+1. else y fi end;

Alec


dev...@math.udel.edu

unread,
Feb 20, 2005, 1:05:08 AM2/20/05
to
Even faster:

mod1:= proc(x)
option inline;
x-trunc(x)+`if`(x>0., 0, 1.)
end proc;

Joe Riel

unread,
Feb 20, 2005, 12:32:26 PM2/20/05
to
Joe Riel <jo...@k-online.com> writes:

> Except that that doesn't work with negative integers,
> that is, it returns 1. instead of 0. I couldn't figure
> out a faster technique then the one I proposed, though
> I suspect that there is one.

Note that checking for integer doesn't work, because
the input might be, say, -2.0. The following handles
the special case, and is marginally faster.

fastmod3 := proc(x)
local y;
y := x-trunc(x);
if y<0. then
y := y+1.;
if y = 1. then 0. else y
fi
else y fi
end:

Making my original inline (per Carl's suggestion), I get the following,
which is faster for positives and essentially the same for negatives:

fastmod4 := proc(x)
option inline;
`if`(x>=0.
,x-trunc(x)
,x-trunc(x)+1-trunc(x-trunc(x)+1))
end:

Seems as though it could be improved.


Joe

Joe Riel

unread,
Feb 20, 2005, 12:14:18 PM2/20/05
to
"Alec Mihailovs" <al...@mihailovs.com> writes:

Except that that doesn't work with negative integers,
that is, it returns 1. instead of 0. I couldn't figure
out a faster technique then the one I proposed, though
I suspect that there is one.

Joe

Alec Mihailovs

unread,
Feb 20, 2005, 1:16:24 PM2/20/05
to
"Joe Riel" <jo...@k-online.com> wrote in message
news:87k6p3r...@k-online.com...

>> fastmod2:=proc(x) local y;
>> y:=x-trunc(x);
>> if y<0. then y+1. else y fi end;
>
>
> Except that that doesn't work with negative integers,
> that is, it returns 1. instead of 0. I couldn't figure
> out a faster technique then the one I proposed, though
> I suspect that there is one.

It works OK on my Maple,

> fastmod2(-3);
0

> fastmod2(-3.);

0.
Alec


Alec Mihailovs

unread,
Feb 20, 2005, 1:25:31 PM2/20/05
to
<dev...@math.udel.edu> wrote in message
news:1108879508.1...@c13g2000cwb.googlegroups.com...

I tried that, too. In my calculations, fastmod2 is faster. For example,

> tt:=time():
> for n to 1000000 do mod1(n-0.7) od:
> time()-tt;

5.781

> tt:=time():
> for n to 1000000 do fastmod2(n-0.7) od:
> time()-tt;

5.359
Same with positive numbers,

> tt:=time():
> for n to 1000000 do mod1(n+0.7) od:
> time()-tt;

5.750

> tt:=time():
> for n to 1000000 do fastmod2(n+0.4) od:
> time()-tt;

5.360

Generally speaking, `if` is slower than if ... fi.

Alec


Joe Riel

unread,
Feb 20, 2005, 1:38:35 PM2/20/05
to
"Alec Mihailovs" <al...@mihailovs.com> writes:

I stand corrected, thanks. I have a bunch of these in one file,
probably I tested one of my old ones with the same name. From
the code it is obvious that yours has to work.

Joe

Alec Mihailovs

unread,
Feb 20, 2005, 1:45:15 PM2/20/05
to

"Joe Riel" <jo...@k-online.com> wrote in message
news:87fyzrr...@k-online.com...

>> Except that that doesn't work with negative integers,
>> that is, it returns 1. instead of 0. I couldn't figure
>> out a faster technique then the one I proposed, though
>> I suspect that there is one.
>
> Note that checking for integer doesn't work, because
> the input might be, say, -2.0.

I don't see why you have problem with that. As I said in another message, it
works OK without checking for integer,

> fastmod2(-2.0);

0.

> The following handles
> the special case, and is marginally faster.
>
> fastmod3 := proc(x)
> local y;
> y := x-trunc(x);
> if y<0. then
> y := y+1.;
> if y = 1. then 0. else y
> fi
> else y fi
> end:
>
> Making my original inline (per Carl's suggestion), I get the following,
> which is faster for positives and essentially the same for negatives:
>
> fastmod4 := proc(x)
> option inline;
> `if`(x>=0.
> ,x-trunc(x)
> ,x-trunc(x)+1-trunc(x-trunc(x)+1))
> end:
>
> Seems as though it could be improved.

In my time calculations, my fastmod2 was always faster than yours fastmod1,
fastmod3, fastmod4, and Carl's modp1. I just tried it again,

> tt:=time():
> for n to 1000000 do mod1(n+0.4) od:
> time()-tt;

5.938

> tt:=time():
> for n to 1000000 do fastmod1(n+0.4) od:
> time()-tt;

6.188

> tt:=time():
> for n to 1000000 do fastmod2(n+0.4) od:
> time()-tt;

5.375

> tt:=time():
> for n to 1000000 do fastmod3(n+0.4) od:
> time()-tt;

5.469

> tt:=time():
> for n to 1000000 do fastmod4(n+0.4) od:
> time()-tt;

5.531

> tt:=time():
> for n to 1000000 do mod1(n-0.3) od:
> time()-tt;

5.907

> tt:=time():
> for n to 1000000 do fastmod1(n-0.3) od:
> time()-tt;

6.125

> tt:=time():
> for n to 1000000 do fastmod2(n-0.3) od:
> time()-tt;

5.390

> tt:=time():
> for n to 1000000 do fastmod3(n-0.3) od:
> time()-tt;

5.453

> tt:=time():
> for n to 1000000 do fastmod4(n-0.3) od:
> time()-tt;

5.578

Evidently, fastmod2 is faster.

Alec


Alec Mihailovs

unread,
Feb 20, 2005, 1:57:44 PM2/20/05
to
"Alec Mihailovs" <al...@mihailovs.com> wrote in message
news:%U4Sd.8996$US3....@news02.roc.ny...

> In my time calculations, my fastmod2 was always faster than yours
> fastmod1, fastmod3, fastmod4, and Carl's modp1. I just tried it again,

I just did it again for negative numbers,

> tt:=time():
> for n to 1000000 do mod1(0.3-n) od:
> time()-tt;

6.156

> tt:=time():
> for n to 1000000 do fastmod1(0.3-n) od:
> time()-tt;

8.298

> tt:=time():
> for n to 1000000 do fastmod2(0.3-n) od:
> time()-tt;

6.109

> tt:=time():
> for n to 1000000 do fastmod3(0.3-n) od:
> time()-tt;

7.250

> tt:=time():
> for n to 1000000 do fastmod4(0.3-n) od:
> time()-tt;

9.312

Again, fastmod2 was fastest.

Alec

Joe Riel

unread,
Feb 20, 2005, 1:57:27 PM2/20/05
to
"Alec Mihailovs" <al...@mihailovs.com> writes:

> I tried that, too. In my calculations, fastmod2 is faster. For example,
>
>> tt:=time():
>> for n to 1000000 do mod1(n-0.7) od:
>> time()-tt;
>
> 5.781
>
>> tt:=time():
>> for n to 1000000 do fastmod2(n-0.7) od:
>> time()-tt;
>
> 5.359


> Same with positive numbers,

Your prevous ones are also positive. You need to change
to mod1(-n-0.7). But I didn't see much difference
in the speed.

My timing has Carl's as being slightly faster
(5 - 10%) but I consider that too small to matter.

Joe

Alec Mihailovs

unread,
Feb 20, 2005, 2:34:33 PM2/20/05
to
"Joe Riel" <jo...@k-online.com> wrote in message
news:871xbbr...@k-online.com...

> Your prevous ones are also positive. You need to change
> to mod1(-n-0.7). But I didn't see much difference
> in the speed.
>
> My timing has Carl's as being slightly faster
> (5 - 10%) but I consider that too small to matter.

My timing still shows that fastmod2 is faster,

> tt:=time():
> for n to 1000000 do mod1(0.3-n) od:
> time()-tt;

6.234
> tt:=time():
> for n to 1000000 do fastmod2(0.3-n) od:
> time()-tt;

6.110

Perhaps, mod1 is faster on Linux?

This way of doing time calculations is not the most precise though. I tried
it with setting gc high and including gc() after the calculations as Carl
suggested in another thread. This time modp1 was slightly faster for
negative numbers but slower for positive,

> kernelopts(gcfreq= 2^22):
> gc();tt:=time():
> for n to 1000000 do mod1(0.3-n) od:gc();
> time()-tt;

5.702
> gc();tt:=time():
> for n to 1000000 do fastmod2(0.3-n) od:gc();
> time()-tt;

5.750
> gc();tt:=time():
> for n to 1000000 do mod1(n+0.4) od:gc();
> time()-tt;

5.625
> gc();tt:=time():
> for n to 1000000 do fastmod2(n+0.4) od:gc();
> time()-tt;

4.938
Alec

Alec Mihailovs

unread,
Feb 20, 2005, 2:43:45 PM2/20/05
to
<dev...@math.udel.edu> wrote in message
news:1108879508.1...@c13g2000cwb.googlegroups.com...

However, it gives wrong results for negative integers (similarly to Joe
Riel's note),

> mod1(0);
1.

> mod1(-3);

1.

> mod1(-3.0);

1.

Alec


dev...@math.udel.edu

unread,
Feb 20, 2005, 4:52:47 PM2/20/05
to
Alec Mihailovs wrote:
> "Joe Riel" <jo...@k-online.com> wrote in message
> > My timing has Carl's as being slightly faster
> > (5 - 10%) but I consider that too small to matter.

Joe, I'd like to see your test method.

> Perhaps, mod1 is faster on Linux?

I am using WindowsXP, Pentium IV.

> This way of doing time calculations is not the most precise though. I
tried
> it with setting gc high and including gc() after the calculations as
Carl
> suggested in another thread. This time modp1 was slightly faster for
> negative numbers but slower for positive,
>
> > kernelopts(gcfreq= 2^22):
> > gc();tt:=time():
> > for n to 1000000 do mod1(0.3-n) od:gc();
> > time()-tt;
> 5.702
> > gc();tt:=time():
> > for n to 1000000 do fastmod2(0.3-n) od:gc();
> > time()-tt;
> 5.750

I can confirm that fastmod2 is faster on my machine when the test is
done as above. I don't know why that is. Perhaps it has something to
do with repetition of data. Here are three tests which seem to me to
be a bit more realistic that show that mod1 is faster. Please confirm
these on your machine.

.> restart;
.> kernelopts(version);
Maple 9.52, IBM INTEL NT, Dec 17 2004 Build ID 175529
.> mod1:= proc(x)
.> option
.> inline
.> ,`Carl's entry`
.> ;
.> x-trunc(x)+`if`(x>0., 0, 1.)
.> end proc:
.> fastmod2:= proc(x)
.> option `Alec's entry`;
.> local y;
.> y:=x-trunc(x);
.> if y<0. then y+1. else y fi
.> end proc:

Generate test data:
.> Num:= rand(-2^29..2^29): Den:= rand(1..2^29):
.> Data:= [seq](evalf[14](Num()/Den()), i= 1..2^18):

Test applying 'seq' to the data (tests done twice):
> kernelopts(gcfreq= 2^22):
> gc(): st:= time(): L:= [seq](mod1(x), x= Data): eval(L): gc():
time()-st;
5.516
> gc(): st:= time(): L:= [seq](fastmod2(x), x= Data): eval(L): gc():
time()-st;
7.702
> gc(): st:= time(): L:= [seq](mod1(x), x= Data): gc(): eval(L):
time()-st;
5.828
> gc(): st:= time(): L:= [seq](fastmod2(x), x= Data): eval(L): gc():
time()-st;
7.796

Test them in a for loop:
> gc(): st:= time(): for n to 2^18 do fastmod2(Data[n]) od: gc():
time()-st;
7.890
> gc(): st:= time(): for n to 2^18 do mod1(Data[n]) od: gc():
time()-st;
5.593

Test applying map to the data.
> gc(): st:= time(): L:= map(mod1, Data): eval(L): gc(): time()-st;
6.702
> gc(): st:= time(): L:= map(fastmod2, Data): eval(L): gc(): time()-st;
7.672

Alec Mihailovs

unread,
Feb 20, 2005, 5:42:10 PM2/20/05
to
<dev...@math.udel.edu> wrote in message
news:1108936367.5...@c13g2000cwb.googlegroups.com...

>
> I can confirm that fastmod2 is faster on my machine when the test is
> done as above. I don't know why that is. Perhaps it has something to
> do with repetition of data. Here are three tests which seem to me to
> be a bit more realistic that show that mod1 is faster. Please confirm
> these on your machine.

Yes, in these tests mod1 is clearly faster than fastmod2. Note, however,
that mod1 should be fixed so that it give 0 for negative integers.

Alec


Joe Riel

unread,
Feb 20, 2005, 5:58:56 PM2/20/05
to
dev...@math.udel.edu writes:

> Alec Mihailovs wrote:
>> "Joe Riel" <jo...@k-online.com> wrote in message
>> > My timing has Carl's as being slightly faster
>> > (5 - 10%) but I consider that too small to matter.
>
> Joe, I'd like to see your test method.

I'm now confused about what I tested, that is whether I
tested your procedure or a modified version (with inline option) that
works properly for nonpositive integers (i.e. returns 0, not 1).
Probably the latter. My first tests made no attempt to vary the
input, just keep looping on the same function call. Using random,
but the same, data is probably better.

Here are the results running your test:

Carl's Time Joe's Time Diff
----------- ---------- ----
1. 5.516 6.409 +16%
2. 7.702 7.720 0.3
3. 5.828 5.990 2.8
4. 7.796 7.850 0.7
5. 7.890 8.201 3.9
6. 5.593 6.510 16
7. 6.702 7.660 14
8. 7.672 8.220 7

This is Maple 9.51, Linux, Pentium IV, 1.6GHz.

Joe

Joe Riel

unread,
Feb 20, 2005, 6:00:05 PM2/20/05
to
Joe Riel <jo...@k-online.com> writes:

> This is Maple 9.51, Linux, Pentium IV, 1.6GHz.

I meant 9.52.

> kernelopts(version);
Maple 9.52, IBM INTEL LINUX, Dec 17 2004 Build ID 175529

Joe

Joe Riel

unread,
Feb 20, 2005, 6:08:08 PM2/20/05
to
"Alec Mihailovs" <al...@mihailovs.com> writes:

Tthere is a bug in Alec's routine:

fastmod2(0.99999999999999999999999);
1.000000000

If you look at Maple's floor routine, you can see that they
use the maximum of Digits or the precision in a call to evalf
to correctly handle this case.

Joe

Alec Mihailovs

unread,
Feb 20, 2005, 6:22:08 PM2/20/05
to

"Alec Mihailovs" <al...@mihailovs.com> wrote in message
news:6n8Sd.9158$6p7....@news01.roc.ny...

It seems as if the method of accessing the data is important. Accessing data
in lists works faster with option inline. For the data calculated in place,
opposite. For example, with Carl's setting,

> gc(): st:= time(): for n to 2^16 do fastmod2(sin(2.*n)) od: gc():
> time()-st;

22.703

> gc(): st:= time(): for n to 2^16 do mod1(sin(2.*n)) od: gc():
> time()-st;
>

24.141

This is with non-modified version of mod1 and in Maple 9.52 on Windows XP.
Joe mentioned the modified version of mod1, but I don't know what it is.

Alec


Alec Mihailovs

unread,
Feb 20, 2005, 6:52:24 PM2/20/05
to

"Joe Riel" <jo...@k-online.com> wrote in message
news:87ll9ir...@k-online.com...

> Tthere is a bug in Alec's routine:
>
> fastmod2(0.99999999999999999999999);
> 1.000000000

Well, mod1 also gives the same answer,

> mod1(0.99999999999999999999999);

1.000000000

That explains by the following Maple calculation,

> 0.99999999999999999999-0;

1.000000000

I didn't check that, but Joe's routines also should have that bug.

There is more serious bug though,in both facmod2 and mod1, and, probably, in
Joe's routines,

> facmod2(123456789012.9);
0.

> mod1(123456789012.9);

0.

It happens because of the following calculation (bug?) in Maple,

> 123456789012.9-123456789012;

0.

> If you look at Maple's floor routine, you can see that they
> use the maximum of Digits or the precision in a call to evalf
> to correctly handle this case.

By the way, Carl's initial routine also has both of these bugs in spite of
using floor,

> fmod1:=x->x-floor(x):

> fmod1(0.99999999999999999999);

1.000000000

> fmod1(123456789012.9);

0.
Alec


dev...@math.udel.edu

unread,
Feb 20, 2005, 7:03:09 PM2/20/05
to
Alec Mihailovs wrote:
> It seems as if the method of accessing the data is important.
Accessing data
> in lists works faster with option inline. For the data calculated in
place,
> opposite. For example, with Carl's setting,
>
> > gc(): st:= time(): for n to 2^16 do fastmod2(sin(2.*n)) od: gc():
> > time()-st;
> 22.703
> > gc(): st:= time(): for n to 2^16 do mod1(sin(2.*n)) od: gc():
> > time()-st;
> 24.141

I believe that the reason for the difference here is the repetition of
the input in inline code, rather than a list versus non-list thing.
When the expression mod1(sin(2.*n)) is expanded, it becomes
sin(2.*n) - trunc(sin(2.*n)) + `if`(sin(2*.n)>0., 0, 1.)
so that sin(2.*n) is computed three times. Most of that is done with
remember tables, but still there is the time for the table lookups.

.> forget(evalf): gc(): st:= time():
.> for n to 2^15 do fastmod2(sin(2.*n)) od:
.> gc(): time()-st;
16.250
.> forget(evalf): gc(): st:= time():
.> for n to 2^15 do x:= sin(2.*n); mod1(x) od:
.> gc(): time()-st;
16.282

So, the moral is, when you apply an inline procedure to a compound
expression, you need to be wary of how many time the parameter is
replicated in the expression.

Oh, I forgot to add something to my last post, an important and
little-known fact. Note that seq is faster than map. This seems
generally true, and is very useful to know because seq can be used
instead of map in most cases. In a case where the user constructs an
arrow expression for map, seq is MUCH faster. For example,

L:= [$1..2^18]:
[seq](x^2, x in L)
is a little more than twice as fast as
map(x-> x^2, L)

In a case where the user constructs an `@` expression, the difference
is extreme. For example.
[seq](f(g(x)), x in L)
is about 12 times as fast as
map(f@g, L).
This is because `@` is a high level procedure, and the above map
generates 2^18 calls to it. Check out showstat(`@`).

dev...@math.udel.edu

unread,
Feb 20, 2005, 7:09:34 PM2/20/05
to
All these examples with 0.999.... and subtracting two numbers with many
digits are not bugs. They are artifacts of the floating-point number
system. They will not occur in practice because the input numbers
will then be at the same precision as the computation. Note that
0.999999999999999999 IS 1.0
.> evalb(0.999999999999999999 = 1.0);
true
and 1234567899012.9 IS 123456789012.
.> evalb(123456789012.9 = 123456789012.);
true

Alec Mihailovs

unread,
Feb 20, 2005, 7:47:18 PM2/20/05
to

<dev...@math.udel.edu> wrote in message
news:1108944188....@l41g2000cwc.googlegroups.com...

> Alec Mihailovs wrote:
>> It seems as if the method of accessing the data is important.
> Accessing data
>> in lists works faster with option inline. For the data calculated in
> place,
>> opposite. For example, with Carl's setting,
>>
>> > gc(): st:= time(): for n to 2^16 do fastmod2(sin(2.*n)) od: gc():
>> > time()-st;
>> 22.703
>> > gc(): st:= time(): for n to 2^16 do mod1(sin(2.*n)) od: gc():
>> > time()-st;
>> 24.141
>
> I believe that the reason for the difference here is the repetition of
> the input in inline code, rather than a list versus non-list thing.
> When the expression mod1(sin(2.*n)) is expanded, it becomes
> sin(2.*n) - trunc(sin(2.*n)) + `if`(sin(2*.n)>0., 0, 1.)
> so that sin(2.*n) is computed three times. Most of that is done with
> remember tables, but still there is the time for the table lookups.

> So, the moral is, when you apply an inline procedure to a compound


> expression, you need to be wary of how many time the parameter is
> replicated in the expression.

You are right, moving calculation of sin(2*n) out makes mod1 faster while
fastmod2 has about the same timing,

> gc(): st:= time(): for n to 2^16 do x:=sin(2.*n);fastmod2(x) od: gc():
> time()-st;

22.766

> gc(): st:= time(): for n to 2^16 do x:=sin(2.*n);mod1(x) od: gc():
> time()-st;
>

22.564

Still, in calculations with lists timing is significantly different,

> L:=[seq](sin(2.*n),n=1..2^16):

> gc(): st:= time(): for n to 2^16 do fastmod2(L[n]) od: gc():
> time()-st;
>

0.828

> gc(): st:= time(): for n to 2^16 do mod1(L[n]) od: gc():
> time()-st;
>

0.687

Alec


Joe Riel

unread,
Feb 20, 2005, 8:33:39 PM2/20/05
to
dev...@math.udel.edu writes:

Well, I'd consider them potential bugs. If the output is supposed to
be in the real range [0,1) and we get 1.0, that's not a good thing.
However, we never rigorously specified the requirements. Note that in
your test data you use evalf[14], so the list of data has fourteen
digits, which means the bug could potentially arise (actually, because
the denominator is limited to 2^29 it isn't going to, but I wonder
whether you considered that when generating the data). A simple fix,
which only hurts the performance about 10%, is the following:

slowmod3 := proc(x)
local y;
description "Alec's method, modified to handle arbirary precision";
y := evalf(x-trunc(x));
if y<0. then y := y+1. fi;
if y=1. then 0. else y fi
end:

Note that the output is now always a float, unlike the originals,
which return rationals given rationals.

Joe

No one

unread,
Feb 20, 2005, 9:55:58 PM2/20/05
to
Joe Riel wrote:

[...]


> Well, I'd consider them potential bugs. If the output is supposed to
> be in the real range [0,1) and we get 1.0, that's not a good thing.

you're aware of "ANSI/IEEE Standard 754-1985", right? as Mr Devore
wrote, this is a consequence of attempting precision arithmetic on a
finite-precision machine. as long as maplesoft did a diligent job, as
long as you are aware of this and can program around it, you're safe.

dev...@math.udel.edu

unread,
Feb 20, 2005, 11:29:45 PM2/20/05
to

Indeed, Maple does provide an environment variable called Rounding
which can be used to control this behavior. Rounding can be set to
'nearest', '0', 'infinity', '-infinity'. The default is 'nearest',
which is appropriate for most computations. 'nearest' means that
numbers are rounded to the nearest representable number, with ties
broken by going to the one that ends in an even digit. This is also
called "round to even", and is specified in the ANSI standard. '0'
means round towards 0, '-infinity' means round down, and 'infinity'
means round up.

In some computations, making an error in one direction is more serious
than making an error in the other direction. Say, for example,
measuring the dose of a danderous medicine: You can always give a
little more if you need to, but if you gave too much the it is too
late. Such is the case with the computation under discussion in this
thread. I think setting Rounding to 0 is the best choice in this case.
Then 0.999999999999 becomes 0.9999999999 (assuming 10 digits
precision, the default) rather than 1.0, and mod1(-0.9999999999) is
0.0000000001. But numbers larger than 10^10 (assuming 10 digits
precision) in absolute value are going to look like integers no matter
what the Rounding is set at (though it may change which integer that
is), and hence they will appear to have 0 fractional part.

Joe Riel

unread,
Feb 21, 2005, 12:59:53 AM2/21/05
to
dev...@math.udel.edu writes:

> Indeed, Maple does provide an environment variable called Rounding
> which can be used to control this behavior. Rounding can be set to
> 'nearest', '0', 'infinity', '-infinity'. The default is 'nearest',
> which is appropriate for most computations. 'nearest' means that
> numbers are rounded to the nearest representable number, with ties
> broken by going to the one that ends in an even digit. This is also
> called "round to even", and is specified in the ANSI standard. '0'
> means round towards 0, '-infinity' means round down, and 'infinity'
> means round up.

I'd forgotten about that. We can use it to ensure that Alec's original
procedure stays in [0,1):

fastmod2 := proc(x)
local y;
description "Alec's entry";
Rounding := 0;
y := x-trunc(x);


if y<0. then
y+1.
else y fi

end:

Digits := 2: # to make it easy to see...
fastmod2(0.999);
0.99

fastmod2(-0.999); # not ideal, but...
0.01

fastmod2(-1.001);
0.


Joe

Joe Riel

unread,
Feb 21, 2005, 12:46:17 AM2/21/05
to
No one <no...@thisplace.thanks> writes:

Yes. But that isn't the cause of the problem---Maple is not a fixed
precision machine. That is how, for example, Maple's built-in trunc
procedure gives the proper result for an input like trunc(0.999...99)
(0) and why its floor and ceil procedures, which are not built-in, go
to some lengths to give the proper result. I considered the operation
of Carl's procedure to be buggy not because it wouldn't work in the
application (which isn't defined), but because it doesn't work with
the usual range of Maple inputs.

It isn't that hard to get the "right" result, what makes it interesting
is doing it without too much a speed penalty. Something like the
following should work (I'm not claiming that this is bug free):

slowmod2 := proc(x)


local y;
description "Alec's method, modified to handle arbirary precision";

if type(y,integer) then return 0 fi;
Digits := max(Digits,length(op(1,x)));
y := x-trunc(x);


if y<0. then y := y+1. fi;
if y=1. then 0. else y fi;
end:

With that, if x is specified to a greater precision than Digits,
it will use that as the precision for the computation. This
is basically what Maple does with the ceil and floor precision
with numerical input, though those routines are consideraly more
general and a bit cleverer.

slowmod2(0.99999999999999999999999999999999999999999999999);

0.99999999999999999999999999999999999999999999999

slowmod2

slowmod2(1-1/10^20);
99999999999999999999
---------------------
100000000000000000000

That's not to say that this is the right approach for the problem
at hand. Consider

evalf(slowmod2(1-1/10^20));

1.000000000

What happens here is that, while slowmod2 does the conversion
properly, the call to evalf rounds the result to one; the extra effort
to get the "right" result is wasted.


Joe

Alec Mihailovs

unread,
Feb 21, 2005, 2:27:51 AM2/21/05
to
"Vladimir Bondarenko" <v...@cybertester.com> wrote in message
news:1108966961.2...@g14g2000cwa.googlegroups.com...
>
> In ANY Maple version since 1994
>
>> limit(sin(z)^2+cos(z)^2, z=infinity);
>
> 0..2

As usual, a proper approach is giving the correct answer,

> limit(simplify(sin(z)^2+cos(z)^2),z=infinity);

1

The ?limit,return help page says, "If limit returns a numeric range it means
that the value of the limiting expression is known to lie in that range for
arguments restricted to some neighborhood of the limit point. It does not
necessarily imply that the limiting expression is known to achieve every
value infinitely often in this range."

Evidently, 1 is in the range 0..2. What's wrong with that?

Alec


Vladimir Bondarenko

unread,
Feb 21, 2005, 12:27:32 AM2/21/05
to
No one <n...@thisplace.thanks> writes on Sun, Feb 20 2005 6:55 pm

NO> as long as maplesoft did a diligent job


ROTFLMAO :)


In ANY Maple version since 1994

> limit(sin(z)^2+cos(z)^2, z=infinity);

0..2


Whilst Derive, MACSYMA, Mathematica, Maxima, and MuPAD return the
correct answer, known already to Pythagoras, 1.

They does hate Gauss...

2 years ago, at the entry page of www.cybertester.com, I pointed
out PUBLICLY to Maplesoft how to spell Gauss... and what we see?

Open in Maple 9.5.2 ?examples,LA_NAG and read about 'fraction-
free Guassian elimination'.


http://groups-beta.google.com/group/comp.soft-sys.math.maple/msg/76ab0802cb1450ab

PL> Did you watch Maplesoft over the years? The price went up,
PL> their hype went up, their arrogance went up -
PL> and the quality went down.


Several minutes ago the VM machine showed me yet another piece
of Maple bunk... or quality degradation... name it as you wish.


int(int(arctan((2*x*y)/(x^2 - y^2)), x=0..1), y=0..1);


-------------------- (2004) Maple 9.5.1 ----------------------

0

-------------------- (2004) Maple 9.5 ------------------------

0

-------------------- (2003) Maple 9 --------------------------

0

-------------------- (2002) Maple 8 --------------------------

Error, (in assuming) when calling `limit`. Received: 'invalid
limiting point'

-------------------- (2001) Maple 7 --------------------------

int(int(arctan(2*x*y/(x^2-y^2)),x = 0 .. 1),y = 0 .. 1)

-------------------- (2000) Maple 6 --------------------------

int(-1/2/y/(-1/y^2)^(1/2)*Sum(2^(2*_k1)*(-1)^_k1*y^(2*_k1)*(-y^
2)^(-2*_k1)*(GAMMA(1+_k1)*Pi*csc(Pi*_k1)/GAMMA(1-_k1)*(-1/y^2)^
(-1/2-_k1)-1/_k1*(-1/y^2)^(-1/2-2*_k1)*hypergeom([1+2*_k1, _k1]
,[1+_k1],y^2)*GAMMA(1+2*_k1))/(1/2+_k1)/GAMMA(1+2*_k1),_k1 = 0
.. infinity),y = 0 .. 1)

-------------------- (1997) Maple V Rel 5 --------------------

1/2*Pi

-------------------- (1995) Maple V Rel 4 --------------------

1/2*Pi

-------------------- (1994) Maple V Rel 3 --------------------

1/2*Pi

---------------------------------------------------------------


Only old Maple of 1994-1997 can calculate the integral correctly.


To my (and my purse) sadness, new-fangled Maplesoft of the last
years and 'diligent job' are incompatible essences.

Maple is bug-ridden, it is a colossus with clay feet.

Read our joint Review on the verge of being published.


The VM machine + Vladimir Bondarenko,
The # 1 world's connoisseurs of Maple bugs

http://www.cybertester.com/
http://maple.bug-list.org/
http://www.CAS-testing.org/

Message has been deleted

Vladimir Bondarenko

unread,
Feb 21, 2005, 2:15:56 AM2/21/05
to
No one <n...@thisplace.thanks> writes on Sun, Feb 20 2005 6:55 pm

http://groups-beta.google.com/group/comp.soft-sys.math.maple/msg/d49b9d8aba9f5ea0

NO> as long as maplesoft did a diligent job

In ANY Maple version since 1994

> limit(sin(z)^2+cos(z)^2, z=infinity);

0..2


Whilst Derive, MACSYMA, Mathematica, Maxima, and MuPAD return the
correct answer, known already to Pythagoras, 1.

They does hate Gauss...

2 years ago, at the entry page of www.cybertester.com, I pointed
out PUBLICLY to Maplesoft how to spell Gauss... and what we see?

Open in Maple 9.5.2 ?examples,LA_NAG and read about 'fraction-
free Guassian elimination'.


http://groups-beta.google.com/group/comp.soft-sys.math.maple/msg/76ab0802cb1450ab

PL> Did you watch Maplesoft over the years? The price went up,
PL> their hype went up, their arrogance went up -
PL> and the quality went down.


Several minutes ago the VM machine showed me yet another piece
of Maple bunk... or quality degradation... name it as you wish.


int(abs(z)^2-z^2, z);

-------------------- (2004) Maple 9.5.1 ----------------------

0

-------------------- (2004) Maple 9.5 ------------------------

0

-------------------- (2003) Maple 9 --------------------------

0

-------------------- (2002) Maple 8 --------------------------

0

-------------------- (2001) Maple 7 --------------------------

0

-------------------- (2000) Maple 6 --------------------------

0

-------------------- (1997) Maple V Rel 5 --------------------

1/3*abs(z)^2*z-1/3*z^3

-------------------- (1995) Maple V Rel 4 --------------------

1/3*abs(z)^2*z-1/3*z^3

-------------------- (1994) Maple V Rel 3 --------------------

int(abs(z)^2-z^2,z)

---------------------------------------------------------------


Only old Maple of 1995-1997 can calculate the integral correctly.


To my (and my purse) sadness, new-fangled Maplesoft of the last
years and 'diligent job' are incompatible essences.

Maple is bug-ridden, it is a colossus with clay feet.

Read our joint Review on the verge of being published.


The VM machine + Vladimir Bondarenko,

The # 1 connoisseurs of Maple bugs

Maarten van Reeuwijk

unread,
Feb 21, 2005, 5:34:09 AM2/21/05
to
Maarten van Reeuwijk wrote:

Thank you for all the responses! I see that there are plenty ways of doing
it. Yet it is strange that such an elementary operation is not implemented
in Maple's core.

Maarten

> Hi,
>
> I am trying to get Maple to calculate things like 10.3 mod 1 = 0.3 but
> this function only seems to work for integers. Is there a corresponding
> function to easily do this?
>
> Thx Maarten

--
===================================================================
Maarten van Reeuwijk Thermal and Fluids Sciences
Phd student dept. of Multiscale Physics
www.ws.tn.tudelft.nl Delft University of Technology

dev...@math.udel.edu

unread,
Feb 21, 2005, 2:37:29 AM2/21/05
to
Alec Mihailovs wrote:
> Still, in calculations with lists timing is significantly different,

Yes, I did other tests (shown below) that show that the difference
between the inline and non-inline procedure is more pronounced when the
data is from a list. Maybe it has something to do with memory
cacheing.

When you test the inline proc, do
x:= List[n]; mod1(x);
instead of
mod1(List[n]);
because in the latter, the list lookup is repeated three time per call
to mod1.

Here are my tests:

.> Num:= rand(-2^29..2^29): Den:= rand(1..2^29):

.> N:= 2^18:
.> ListData:= [seq](evalf[14](Num()/Den()), i= 1..N):
.> for n to N do TableData[n]:= ListData[n] od:
.> ProcData:= ()-> evalf[14](Num()/Den()):
.> VecData:= Vector(N, ListData, datatype= float[8]):
.> kernelopts(gcfreq= 2^22):
.> gc(): st:= time():
.> for n to N do x:= ListData[n]; mod1(x) od:
.> gc(): time()-st;
8.767

.> gc(): st:= time():
.> for n to N do fastmod2(ListData[n]) od:
.> gc(): time()-st;
10.077

.> gc(): st:= time():
.> for n to N do x:= TableData[n]; mod1(x) od:
.> gc(): time()-st;
9.953

.> gc(): st:= time():
.> for n to N do fastmod2(TableData[n]) od:
.> gc(): time()-st;
10.452

.> gc(): st:= time():
.> for n to N do x:= ProcData(): mod1(x) od:
.> gc(): time()-st;
21.797

.> gc(): st:= time():
.> for n to N do fastmod2(ProcData()) od:
.> gc(): time()-st;
22.124

The difference fastmod2 minus mod1 for list, table, and procedure,
respectively:

> 10.077-8.767, 10.452-9.953, 22.124-21.797;

1.310, 0.499, 0.327

Joe Riel

unread,
Feb 21, 2005, 11:56:35 AM2/21/05
to
Maarten van Reeuwijk <maarten@remove_this_ws.tn.tudelft.nl> writes:

> Thank you for all the responses! I see that there are plenty ways of doing
> it. Yet it is strange that such an elementary operation is not implemented
> in Maple's core.

I don't know that it would be useful enough to be a standard Maple
procedure. A more interesting, and potentially useful, procedure
would be a two argument version of modp that handles floats. The
following simplistic procedure, based on scaling the input and output
and using the previous mod1, is poor, it loses digits of precision:

modf1 := proc(x,m)
local y;
y := x/m;
y := y - trunc(y);
if y < 0. then
m*(y+1)
else
m*y
fi
end proc:

> Digits := 3;
> modf1(3,2.9);
0.087
> modf1(3,2.99);
0.

Better is something along the lines of

modf2 := proc(x,m)
local y;
y := frem(x,m);
if y < 0. then y := y+m fi;
if y=m then 0. else y fi
end:

> modf2(3,2.9);
0.1
> modf2(3,2.99);
0.01

Joe

Alec Mihailovs

unread,
Feb 21, 2005, 5:07:16 PM2/21/05
to

"Joe Riel" <jo...@k-online.com> wrote in message
news:87zmxxq...@k-online.com...

> I don't know that it would be useful enough to be a standard Maple
> procedure.

That would be very useful for many of my calculations. I use a procedure
written in assembler and compiled to dll though instead. Nothing can be
faster. For not very intensive calculations, the following simple procedure
is pretty good and avoids problems with precision,

fmod1:=proc(x) local y;
Digits:=length(x);
y:=x-trunc(x);
if y<0 then y+1 else y fi end:

For example,

fmod1(357892567.333333333333334444444444455555555);

0.333333333333334444444444455555555

fmod1(-35789256733333333333333.4444444444455555555);

0.5555555555544444445

fmod1(-3);

0

fmod1(-1234567890123456789012345678/8888876487);

383252572
---------
987652943
Alec


Alec Mihailovs

unread,
Feb 21, 2005, 5:16:08 PM2/21/05
to
"Alec Mihailovs" <al...@mihailovs.com> wrote in message
news:oYsSd.9812$6p7....@news01.roc.ny...

> For not very intensive calculations, the following simple procedure is
> pretty good and avoids problems with precision,
>
> fmod1:=proc(x) local y;
> Digits:=length(x);
> y:=x-trunc(x);
> if y<0 then y+1 else y fi end:

Well, just found a bug in it as well, fmod1(0) gives an error :-)

Alec


Alec Mihailovs

unread,
Feb 21, 2005, 5:29:41 PM2/21/05
to

"Alec Mihailovs" <al...@mihailovs.com> wrote in message
news:I4tSd.9814$6p7....@news01.roc.ny...

>>
>> fmod1:=proc(x) local y;
>> Digits:=length(x);
>> y:=x-trunc(x);
>> if y<0 then y+1 else y fi end:
>
> Well, just found a bug in it as well, fmod1(0) gives an error :-)

The following seems working without much of an overhead,

fmod1:=proc(x) local y;
try Digits:=length(x);


y:=x-trunc(x);
if y<0 then y+1 else y fi

catch: return 0 end end:

For example,

fmod1(0);
0

fmod1(infinity);

undefined

fmod1(z);

0

Alec


0 new messages