maxima-discuss

[Maxima-discuss] Problem with bfloat

[couraud <l.c...@gm...>]

Hi everyone,

 

 

It seems there is a problem when I try to solve this:

 

eq :
(400099987013675381332258793266100212757877348865572395547868673444799488095
5893635730834603909005391274797631516180749919175073156505741219342752612513
9943062483732353369059902143126883533765447106960495457745029916237144786993
152

 *r6^4

 
-635169949538565681576982326160549436240883584439338432673525994990913080807
9502608605649463006744118010506365410372078424922565514273721530399556185069
4536691922119800765043952479301200067938845077085523967714959372759043100513
4028800

  *r6^3

 
+377358294448778745583055149280325938821134900815134080413155213818809381538
5623045926247772506066538790464420557034914421623994643207226503772289573557
4746770514637506165721077202702808030741639574712238081225104648285021754007
48646400000

  *r6^2

 
-994827962362789386985928920959408041604180841307656030461241481458711871699
1627738173731600116720301970551077968613578942034278006729503574412645544866
2299286973331692985051776142386278947698580925838624655188543051738657008016
94924800000000

  *r6

 
+982292276311002270931173985237493254470108298967760469537118227582552304036
5218991617759102209164976234430910614344753842949200199390345412495844846726
1547191777272062929190479200745823489671558337870426293270315863437665782915
39968000000000000)

 /582076609134674072265625$

 

trace(bfloat)$

 

solve(eq, r6);

1 Enter bfloat [0.0 + 1..#INF00e+000E+2678440]

 2 Enter bfloat [0.0]

 2 Exit  bfloat 0.0b0

 2 Enter bfloat [1..#INF00e+000E+69]

bfloat: attempted conversion of floating-point infinity.

 -- an error. To debug this try: debugmode(true);

 

 

[Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[Robert D. <rob...@gm...>]

On 2014-08-30, couraud <l.c...@gm...> wrote:

> solve(eq, r6);
>
> 1 Enter bfloat [0.0 + 1..#INF00e+000E+2678440]
>
>  2 Enter bfloat [0.0]
>
>  2 Exit  bfloat 0.0b0
>
>  2 Enter bfloat [1..#INF00e+000E+69]
>
> bfloat: attempted conversion of floating-point infinity.

I think what's going on here is that Maxima is attempting to determine
the sign of some constant expression (involving ratios of enormous
integers), and tries to compute a floating point approximation.
I think there is some provision in the code for failure to compute a
float, and then it goes to a bigfloat approximation. That scheme
works OK for Lisp implementations (Clisp, SBCL, Clozure CL) which can't
compute floating point infinity by default, but it fails for those
which do (GCL, ECL) -- after computing floating point infinity, Maxima
tries to convert it to bigfloat (I don't know what's the purpose of
that) and fails.

The code of interest is in SIGN1 (src/compar.lisp) at lines 1243--1259.
I guess one way to fix it is to account for non-finite floats.
Can someone fix it?

best

Robert Dodier

Re: [Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[David B. <dbm...@gm...>]

On 31/08/2014 7:28 AM, Robert Dodier wrote:
> On 2014-08-30, couraud <l.c...@gm...> wrote:
>
>> solve(eq, r6);
>>
>> 1 Enter bfloat [0.0 + 1..#INF00e+000E+2678440]
>>
>>   2 Enter bfloat [0.0]
>>
>>   2 Exit  bfloat 0.0b0
>>
>>   2 Enter bfloat [1..#INF00e+000E+69]
>>
>> bfloat: attempted conversion of floating-point infinity.
> I think what's going on here is that Maxima is attempting to determine
> the sign of some constant expression (involving ratios of enormous
> integers), and tries to compute a floating point approximation.
> I think there is some provision in the code for failure to compute a
> float, and then it goes to a bigfloat approximation. That scheme
> works OK for Lisp implementations (Clisp, SBCL, Clozure CL) which can't
> compute floating point infinity by default, but it fails for those
> which do (GCL, ECL) -- after computing floating point infinity, Maxima
> tries to convert it to bigfloat (I don't know what's the purpose of
> that) and fails.
>
> The code of interest is in SIGN1 (src/compar.lisp) at lines 1243--1259.
> I guess one way to fix it is to account for non-finite floats.
> Can someone fix it?
>
> best
>
> Robert Dodier
Here is a reduced test case.  Nothing magic about the constants, they 
just need to be large.

    I:87^611$
    M:91^211$
    sign(sqrt(I)-M);

Maxima 5.34post http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.10 (a.k.a. GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) display2d:false;

(%o1) false
(%i2) trace(?sign1,bfloat);

(%o2) [sign1,bfloat]
(%i3) I:87^611$
(%i4) M:91^211$
(%i5) sign(sqrt(I)-M);

1 Enter sign1 [87]
1 Exit  sign1 pos
1 Enter sign1
[87^(611/2)-227894466256574889420211635345748034449970012014651359835345211477247577813962279418947874919384990943608901614822595179322674507119231045537829696304840485591301824628190054337112380920376337397871582364269303971018654072412430217008671593811253258418371291241902259707966572789245722710797385164506264434277218952320905773492038425282472059016207718540274567812052960443332201900642178741918587353798310622059491]

  1 Enter bfloat [9.327379053088815^611]
   2 Enter bfloat [9.327379053088815]
   2 Exit  bfloat 9.327379053088816b0
  1 Exit  bfloat 3.335275205992107b592
  1 Enter bfloat [9.327379053088815^611]
   2 Enter bfloat [9.327379053088815]
   2 Exit  bfloat 9.327379053088816b0
  1 Exit  bfloat 3.335275205992107b592
  1 Enter bfloat [-i.nfE+4445842+3.335275205992107b592]
   2 Enter bfloat [-i.nfE+4445842]
bfloat: attempted conversion of floating-point infinity.
  -- an error. To debug this try: debugmode(true);

Re: [Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[Richard F. <fa...@be...>]

On 9/1/2014 7:52 AM, David Billinghurst wrote:
> ... snip...



> , Maxima tries to convert it to bigfloat (I don't know what's the 
> purpose of that) and fails.
..

If you need to know the sign of  v  : sqrt(10^500+1)-10^250, how do you 
propose to
do this?

Here's one way...

fpprec:500$
sign(bfloat(v));

note that fpprec:499   is not enough, and gives "zero"  for the sign.


In general, the question of whether an expression is different from zero 
is recursively
undecidable.  This doesn't prevent us from solving the problem , in 
practice, most of the
time.  Just that you can expect some hairy problems to remain unsolved.

Re: [Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[Robert D. <rob...@gm...>]

On 2014-09-01, David Billinghurst <dbm...@gm...> wrote:

>     I:87^611$
>     M:91^211$
>     sign(sqrt(I)-M);

Thanks for the example. Here's a proposed fix. This causes $FLOAT to
always fail (even for Lisps which otherwise allow float infinity)
and therefore SIGN1 takes a different route and sign returns 'pos.

The patch assumes that (float inf) minus (float inf) != 0.
I'm pretty sure that's required by IEEE 754 (in fact I think the
result is supposed to be NaN) but I can't quote chapter & verse
at the moment. 

best

Robert Dodier

PS.
$ git diff src/comm.lisp
diff --git a/src/comm.lisp b/src/comm.lisp
index 42bdc12..1d6d178 100644
--- a/src/comm.lisp
+++ b/src/comm.lisp
@@ -1178,7 +1178,7 @@
            (car result))))))
            
 (defmfun $float (e)
-  (cond ((numberp e) (float e))
+  (cond ((numberp e) (let ((e1 (float e))) (if (= (- e1 e1) 0) e1 (signal 'floating-point-overflow))))
        ((and (symbolp e) (mget e '$numer)))
        ((or (atom e) (member 'array (cdar e) :test #'eq)) e)
        ((eq (caar e) 'rat) (fpcofrat e))

Re: [Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[David B. <dbm...@gm...>]

On 3/09/2014 8:21 AM, Robert Dodier wrote:
> On 2014-09-01, David Billinghurst <dbm...@gm...> wrote:
>
>>      I:87^611$
>>      M:91^211$
>>      sign(sqrt(I)-M);
> Thanks for the example. Here's a proposed fix. This causes $FLOAT to
> always fail (even for Lisps which otherwise allow float infinity)
> and therefore SIGN1 takes a different route and sign returns 'pos.
>
> The patch assumes that (float inf) minus (float inf) != 0.
> I'm pretty sure that's required by IEEE 754 (in fact I think the
> result is supposed to be NaN) but I can't quote chapter & verse
> at the moment.
>
> best
>
> Robert Dodier
>
> PS.
> $ git diff src/comm.lisp
> diff --git a/src/comm.lisp b/src/comm.lisp
> index 42bdc12..1d6d178 100644
> --- a/src/comm.lisp
> +++ b/src/comm.lisp
> @@ -1178,7 +1178,7 @@
>              (car result))))))
>              
>   (defmfun $float (e)
> -  (cond ((numberp e) (float e))
> +  (cond ((numberp e) (let ((e1 (float e))) (if (= (- e1 e1) 0) e1 (signal 'floating-point-overflow))))
>          ((and (symbolp e) (mget e '$numer)))
>          ((or (atom e) (member 'array (cdar e) :test #'eq)) e)
>          ((eq (caar e) 'rat) (fpcofrat e))
This works for me with gcl-2.6.10 on Windows XP.  Thanks.  All 
regression tests pass.

Re: [Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[Raymond T. <toy...@gm...>]

>>>>> "Robert" == Robert Dodier <rob...@gm...> writes:

    Robert> On 2014-09-01, David Billinghurst <dbm...@gm...> wrote:
    >> I:87^611$
    >> M:91^211$
    >> sign(sqrt(I)-M);

    Robert> Thanks for the example. Here's a proposed fix. This causes $FLOAT to
    Robert> always fail (even for Lisps which otherwise allow float infinity)
    Robert> and therefore SIGN1 takes a different route and sign returns 'pos.

    Robert> The patch assumes that (float inf) minus (float inf) != 0.
    Robert> I'm pretty sure that's required by IEEE 754 (in fact I think the
    Robert> result is supposed to be NaN) but I can't quote chapter & verse
    Robert> at the moment. 

Yes, inf - inf is NaN, in IEEE754 arithmetic.  As is inf*0. I would
add a comment to the code saying so for the next user who isn't so
familiar with IEEE754.

On the other hand, didn't you also write some kind of float-infinity-p
function already? If not, we should probably add a float-infinity-p
and float-nan-p functions just for cases like this.

--
Ray

Re: [Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[Richard F. <fa...@be...>]

A solution to this problem of inf-inf (and a number of others)
is to tag each inf uniquely.  So then inf-inf is zero  if those infs are
identically tagged.

Thus

k :  limit(1/x,x,0);

k-k

returns 0

  On 9/3/2014 9:08 AM, Raymond Toy wrote:
>>>>>> "Robert" == Robert Dodier <rob...@gm...> writes:
>      Robert> On 2014-09-01, David Billinghurst <dbm...@gm...> wrote:
>      >> I:87^611$
>      >> M:91^211$
>      >> sign(sqrt(I)-M);
>
>      Robert> Thanks for the example. Here's a proposed fix. This causes $FLOAT to
>      Robert> always fail (even for Lisps which otherwise allow float infinity)
>      Robert> and therefore SIGN1 takes a different route and sign returns 'pos.
>
>      Robert> The patch assumes that (float inf) minus (float inf) != 0.
>      Robert> I'm pretty sure that's required by IEEE 754 (in fact I think the
>      Robert> result is supposed to be NaN) but I can't quote chapter & verse
>      Robert> at the moment.
>
> Yes, inf - inf is NaN, in IEEE754 arithmetic.  As is inf*0. I would
> add a comment to the code saying so for the next user who isn't so
> familiar with IEEE754.
>
> On the other hand, didn't you also write some kind of float-infinity-p
> function already? If not, we should probably add a float-infinity-p
> and float-nan-p functions just for cases like this.
>
> --
> Ray
>
>
> ------------------------------------------------------------------------------
> Slashdot TV.
> Video for Nerds.  Stuff that matters.
> http://tv.slashdot.org/
> _______________________________________________
> Maxima-discuss mailing list
> Max...@li...
> https://lists.sourceforge.net/lists/listinfo/maxima-discuss

Re: [Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[Raymond T. <toy...@gm...>]

>>>>> "Richard" == Richard Fateman <fa...@be...> writes:

    Richard> A solution to this problem of inf-inf (and a number of others)
    Richard> is to tag each inf uniquely.  So then inf-inf is zero  if those infs are
    Richard> identically tagged.

I think this "inf" is the IEEE754 infinity value, not maxima's inf.

--
Ray

Re: [Maxima-discuss] Problem with bfloat

[Raymond T. <toy...@gm...>]

I didn't have any bfloat issues when trying to solve this with cmucl and a
very recent git version of maxima. I don't see why solve should use bfloats
since it's supposed to be doing a symbolic solution.

However, if you just want the numerical roots, you can try allroots or
bfallroots.



On Sat, Aug 30, 2014 at 6:47 AM, couraud <l.c...@gm...> wrote:

>  Hi everyone,
>
>
>
>
>
> It seems there is a problem when I try to solve this:
>
>
>
> eq :
> (40009998701367538133225879326610021275787734886557239554786867344479948809558936357308346039090053912747976315161807499191750731565057412193427526125139943062483732353369059902143126883533765447106960495457745029916237144786993152
>
>  *r6^4
>
>
>  -635169949538565681576982326160549436240883584439338432673525994990913080807950260860564946300674411801050636541037207842492256551427372153039955618506945366919221198007650439524793012000679388450770855239677149593727590431005134028800
>
>   *r6^3
>
>
>  +3773582944487787455830551492803259388211349008151340804131552138188093815385623045926247772506066538790464420557034914421623994643207226503772289573557474677051463750616572107720270280803074163957471223808122510464828502175400748646400000
>
>   *r6^2
>
>
>  -9948279623627893869859289209594080416041808413076560304612414814587118716991627738173731600116720301970551077968613578942034278006729503574412645544866229928697333169298505177614238627894769858092583862465518854305173865700801694924800000000
>
>   *r6
>
>
>  +9822922763110022709311739852374932544701082989677604695371182275825523040365218991617759102209164976234430910614344753842949200199390345412495844846726154719177727206292919047920074582348967155833787042629327031586343766578291539968000000000000)
>
>  /582076609134674072265625$
>
>
>
> trace(bfloat)$
>
>
>
> solve(eq, r6);
>
> 1 Enter bfloat [0.0 + 1..#INF00e+000E+2678440]
>
>  2 Enter bfloat [0.0]
>
>  2 Exit  bfloat 0.0b0
>
>  2 Enter bfloat [1..#INF00e+000E+69]
>
> bfloat: attempted conversion of floating-point infinity.
>
>  -- an error. To debug this try: debugmode(true);
>
>
>
>
>
>
> ------------------------------------------------------------------------------
> Slashdot TV.
> Video for Nerds.  Stuff that matters.
> http://tv.slashdot.org/
> _______________________________________________
> Maxima-discuss mailing list
> Max...@li...
> https://lists.sourceforge.net/lists/listinfo/maxima-discuss
>
>

Re: [Maxima-discuss] Problem with bfloat

[Barton W. <wi...@un...>]

Maxima compiled with CCL gives--due to subtractive cancellation, the floating point values might be

 totaly bogus, but no errors:


(%i7) rectform(solve(eq,r6))$

(%i8) bfloat(%);
(%o8) [r6 = 3.78775493814778b3 - 3.156703346873123b2 %i,
r6 = 3.156703346873123b2 %i + 3.78775493814778b3, r6 = 3.673147182356886b3,
r6 = 4.62662337510712b3]
(%i9) build_info();
(%o9)
Maxima version: "branch_5_33_base_182_gc59006e_dirty"
Maxima build date: "2014-08-08 06:41:23"
Host type: "i686-pc-cygwin"
Lisp implementation type: "Clozure Common Lisp"
Lisp implementation version: "Version 1.9-r15764  (WindowsX8632)"


--Barton

Re: [Maxima-discuss] Problem with bfloat

[Robert D. <rob...@gm...>]

On 2014-08-31, Raymond Toy <toy...@gm...> wrote:

> I didn't have any bfloat issues when trying to solve this with cmucl and a
> very recent git version of maxima. I don't see why solve should use bfloats
> since it's supposed to be doing a symbolic solution.

The error is in SIGN1, not solve itself -- there is some code to
attempt a floating or bigfloat evaluation under some circumstances.
I haven't tried to understand the rationale. Maybe you can take a
look at it.

best

Robert Dodier

Re: [Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[Robert D. <rob...@gm...>]

On 2014-09-03, Raymond Toy <toy...@gm...> wrote:

> On the other hand, didn't you also write some kind of float-infinity-p
> function already? If not, we should probably add a float-infinity-p
> and float-nan-p functions just for cases like this.

Hrm, so I did ... FLOAT-INF-P and FLOAT-NAN-P in src/float.lisp.
It's been a few years. I'll revise the patch accordingly.

best

Robert Dodier

Re: [Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[Richard F. <fa...@be...>]

On 9/3/2014 12:28 PM, Raymond Toy wrote:
>>>>>> "Richard" == Richard Fateman <fa...@be...> writes:
>      Richard> A solution to this problem of inf-inf (and a number of others)
>      Richard> is to tag each inf uniquely.  So then inf-inf is zero  if those infs are
>      Richard> identically tagged.
>
> I think this "inf" is the IEEE754 infinity value, not maxima's inf.
>
> --
> Ray
>
>
Seems to me that if we had a theory of how to deal with infinities, we 
should convert the IEEE
infinities to "our" infinities.  I'm thinking of a hack that will do 
most of what I proposed (in conjunction
with intervals), that can be read in as a patch.  If there are any 
people eager to look at this,
send me a note.

Re: [Maxima-discuss] SIGN1 strangeness, was: Problem with bfloat

[Henry B. <hb...@pi...>]

Shouldn't 'inf be converted to IEEE inf only when "numer" is active?

Shouldn't IEEE inf be converted back to 'inf by rat, just like it converts floats to rational numbers?

At 08:46 PM 9/4/2014, Richard Fateman wrote:
>On 9/3/2014 12:28 PM, Raymond Toy wrote:
>>>>>>> "Richard" == Richard Fateman <fa...@be...> writes:
>>      Richard> A solution to this problem of inf-inf (and a number of others)
>>      Richard> is to tag each inf uniquely.  So then inf-inf is zero  if those infs are
>>      Richard> identically tagged.
>>
>> I think this "inf" is the IEEE754 infinity value, not maxima's inf.
>>
>> --
>> Ray
>
>Seems to me that if we had a theory of how to deal with infinities, we should convert the IEEE
>infinities to "our" infinities.  I'm thinking of a hack that will do most of what I proposed (in conjunction
>with intervals), that can be read in as a patch.  If there are any people eager to look at this,
>send me a note.