Wolfram Mathematica

Oleksandr R.

12:18 AM

@ℛ.ℳ Block[{Internal`$EqualTolerance = -Infinity}, (1.1 + 1.1 + 1.1) > 3.3]

Interesting. I did not know this tolerance was applied for relational operators other than == and ===

Though I suppose it makes sense. Having it both equal and greater would not be good.

Mr.Wizard

My "reputation" is approaching a nice round number. I expect you all to throw me a party. ;^)

Brett Champion

12:34 AM

@Mr.Wizard As far as I can tell you're nowhere near 2^2^2^2, but if you come out for the conference we can go to Jarlings.

You now need to have at least 4K rep to be on the first page of users.

J. M.

1:00 AM

@halirutan Here, the space in between + and 1.1 matters. ;)

@OleksandrR. There's also Internal`$SameQTolerance...

Oleksandr R.

@J.M. yup. But no Internal`$LessTolerance...

J. M.

Yeah, I suppose it's all a matter of checking "are they the same? if not, what's the sign of their difference?"

Oleksandr R.

1:12 AM

@J.M. do you know any convenient identity for Log[Cosh[x]] that doesn't give numerical problems for very large x? I was thinking of this, specifically.

To tell the truth I am getting discouraged over it. He started from completely the wrong place in asking about parallelization whereas checking the function is well enough behaved would have IMO been more preferable.

J. M.

@OleksandrR. Oh joy, hyperbolics within logs are a pain... let me think.

On the other hand, since $\cosh\,x \asymp \frac12\exp\,x$, there's a convenient approximation in the large...

Oleksandr R.

@J.M. well, I figured that already, and a Padé approximant is okay for small arguments. The problem is joining the two up.

I was hoping you might have an idea for an identity that doesn't involve piecewise approximations. :)

J. M.

@OleksandrR. ...and you've hit on the usual problem of people who write routines for special functions. :)

halirutan

@OleksandrR. What do you think they are doing internally, when they don't use the machine > for comparison?

J. M.

In any event, the problem could be reduced to implementing $\log(\exp\,x\pm 1)$ sans overflow and/or cancellation.

halirutan

1:28 AM

How do they manage to make 1.1+1.1+1.1>3.3 return false?

J. M.

@halirutan I'd guess the same machinery behind N[]...

Oleksandr R.

@halirutan an unfortunate consequence of having 1.1 + 1.1 + 1.1 == 3.3 return True.

J. M.

"adaptive precision" and such, if we are to believe the docs on this.

Well, (1.1 + 1.1 + 1.1) - 3.3 is "sufficiently tiny", whatever that means...

Oleksandr R.

@J.M. no, adaptive precision and significance arithmetic are not that bad. It is a quasi-analytic first order approximation of the error, which works surprisingly well provided you add a few extra digits of working precision for luck.

halirutan

@OleksandrR. Ok, but in compiled functions they use double

J. M.

1:31 AM

@OleksandrR. I didn't say it was bad; halirutan was asking how Less[] might work, and if my reading of the docs is correct, the relational operators use the same machinery as N[].

Oleksandr R.

@J.M. yes, the actual function has Log[Exp[-x] + Exp[x]] but I felt Cosh was easier to deal with.

@J.M. okay. Well, AIUI, that's not correct. I think the documentation is very unclear on how numerics work in general.

J. M.

@OleksandrR. Not that much simpler, except maybe the writing is shorter. :)

Oleksandr R.

@J.M. easier to look up identities on Wolfram Functions, anyway. :) (Didn't find any.)

J. M.

@OleksandrR. Yes, I wish the precision handling could be made less mysterious...

@OleksandrR. Heh. Yeah, I think patching up the asymptotics and the Padé approximant is the best you can do at the moment... I'm drawing a blank on the sets of identities I'm aware of.

Oleksandr R.

@J.M. what really irks me is that they removed useful information in more recent versions in favor of a more black-box approach.

J. M.

1:36 AM

@OleksandrR. Ah, I miss that. IIRC Paul Abbott wrote something nice as well. Too bad it was no longer applicable starting with version 5. Architecture changes, I guess.

@halirutan: still running my stuff; I'll get back to you later on the FTP bit.

halirutan

@OleksandrR. @J.M. Let me clear my point. Maybe I have too less experience in this: When you compile a function down to "C" the data-type they are using is a simple double or float. What they do is to call a special comparison function. Lets say I compute 1.1+1.1+1.1 and look at the bit pattern and I want to compare it with the bitpattern of 3.3: it is obviously not the same, because the calculation and the roundoff error happened before. So how can they return False?

This here:

g = Compile[{{a, _Real, 0}, {b, _Real, 0}},
 (a + a + a) > b,
  CompilationTarget -> "C"]

<< CCodeGenerator`

CCodeStringGenerate[g, "fun"]

Oleksandr R.

@halirutan I don't think it's any different than just calling Less or Equal at the top level, is it?

halirutan

@OleksandrR. No, I don't think so but in "C" we really see that there happens no special thing except of the comparison itself

mbool B0_ 0;
mreal R0_ 0;
mreal R0_ 1;
mreal R0_ 2;
R0_ 0 = A1;
R0_ 1 = A2;
R0_ 2 = R0_ 0 + R0_ 0 + R0_ 0;
{
mreal S0[2];
S0[0] = R0_ 2;
S0[1] = R0_ 1;
B0_ 0 = funStructCompile->Compare_R(7, R0_ 3, 2, S0);
}

Oleksandr R.

@halirutan yes, sure. Well, I think that makes sense. Not sure what the 7 and the 2 mean.

J. M.

Yes, I'd say Mathematica and your favorite compiler would treat relational operators differently...

halirutan

1:44 AM

@J.M. Yes, because, as expected, if you substitute a normal > there, you get True ;-)

As I said, at the point of the comparison they have only the different binary representations of R0_ and S0_.

Oleksandr R.

@halirutan personally I'm not keen on this behavior being carried over to compiled code. But, it is what it is, I suppose.

J. M.

Compare_R() is in the WolframCompileLibrary header file, right?

halirutan

So how they can possibly decide wether R0 and So should be equal?

Oleksandr R.

@halirutan mask off the last 7 bits of the significand, see if they match... if they do, it's good enough for government work!

I suppose the 7 means 7 bits? Though not necessarily, of course.

halirutan

@OleksandrR. Cool, I was afraid of that.

Yes, maybe..

Maybe it just means, return false no matter what.

Oleksandr R.

1:49 AM

Anyway, the tolerance is 7 bits. Log[10] Internal`$EqualTolerance/Log[2] == 7

halirutan

@J.M. No, to the functions they use when you compile something down you don't have direct access.

You don't know what the parameters mean.

J. M.

Tsk, tsk. Black box.

Oleksandr R.

Okay, Looking at CCodeGenerator.m, the 7 is the tolerance.

ℛ.ℳ

@OleksandrR. Thanks, that helps. I've always assumed that machine precision numbers would behave exactly as a double in other systems, but never came across a case where the output differed... it's all the more confusing that the tolerance used is that of a single

@halirutan I tried compiling it to C, but still gave me false... I don't know C, so the output from CompiledPrint didn't make much sense :P

(I only use it to check if there's a MainEvaluate hanging around :D)

Oleksandr R.

@ℛ.ℳ hm? Doubles have a 53-bit significand; minus 7 gives 46 bits. A single has a 23-bit significand.

halirutan

2:01 AM

@ℛ.ℳ Yes, because they use there special comparison even in compiled C code.

Try the following:

ℛ.ℳ

@OleksandrR. Oh, 7 bits... sorry, was reading quickly. I thought 7 decimal digits

halirutan

src = "
  #include \"WolframLibrary.h\"

  DLLEXPORT mint WolframLibrary_getVersion(){
    return WolframLibraryVersion;
  }

  DLLEXPORT int WolframLibrary_initialize( WolframLibraryData
libData) {
  	return 0;
  }

  DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData
libData) {
  	return;
  }

  DLLEXPORT int cmp(WolframLibraryData libData,
  			mint Argc, MArgument *Args, MArgument Res) {
  	mreal I0;
  	mreal I1;
  	I0 = MArgument_getReal(Args[0]);
  	I1 = MArgument_getReal(Args[1]);

@ℛ.ℳ Note the (I0+I0+I0)>I1 where I use the machine comparison.

ℛ.ℳ

@halirutan Neat! Thanks for the help @halirutan @OleksandrR. @J.M. :)

Oleksandr R.

Error in Sin@Log[x^2]. Points: Monte Carlo. Yellow: analytic 1st-order. Green: analytic 2nd-order. Red: significance arithmetic.

halirutan

@OleksandrR. It's the error compared to the true Sin[Log[x^2]]?

Oleksandr R.

2:16 AM

@halirutan the error assuming an uncertain x, sorry.

Sorry again... I do significance arithmetic a disservice here by using a WorkingPrecision of 1. It should look like this:

halirutan

@OleksandrR. I'm definitely too less experienced in this area; meaning I have really no plan about error calculation of approximations. Maybe you should write a blog post one day.

Oleksandr R.

@halirutan yeah, may be a good idea. It is really easy. Just Taylor-expand the covariance matrix.

Actually, it's a very nice demonstration of Mathematica's support for tensors of any order as well. Otherwise it would be more tricky.

Okay, time for me to go. Goodnight!

J. M.

@OleksandrR. Ah, better. I was scratching my head on why the red curve was so rough...

halirutan

@OleksandrR. Good night.

J. M.

Good night!

Hmm, I had no idea there was a lot of research on reassembling shredded documents with image processing. I wonder if anybody's ever used Mathematica for this...

halirutan

2:42 AM

@J.M. If you have shredded paper.. who turns, de-folds and scans the small parts? I mean, that must be a lot of work anyway. You could just start puzzling by yourself then.

Ok, 4:49:36 am. Time for me too to collect my wife from the couch and carry her to bed. Good night all.

J. M.

3:00 AM

@halirutan I don't really have shredded docs; I was reading an old article that brought up the idea of reassembling shredded documents. I suppose the tedium of scanning in the shreds would be better spent actually doing the jigsaw puzzle, though. :)

Transcript for