Listable compiled function crashes the kernel with "no more memory available"

26

2

Bug introduced in 9.0.0 and fixed in 10.3


I am trying to speed up a function that involves continued fractions. Since ContinuedFraction cannot be compiled, and the alternative way with Nestlist of FractionalPart does not produce the same results due to ContinuedFraction's internal algorithm for the choice of the number of terms, I think the best way is to compute the continued fractions and then pass the integer list to the compiled function. Here's what I did:

fc = Compile[{{z, _Real}, {xcf, _Integer, 1}}, 
  Fold[z/(#1 + #2) &, 0, 1/xcf]/z, RuntimeAttributes -> Listable]
z = {-.5, .5}
xcf = Reverse /@ Rest /@ ContinuedFraction[{.123455, .546452}]
fc[z, xcf]

This approach should work, however on my computer (Mathematica 10.0.2, Xubuntu 15.04), evaluating the compiled function over a list of arguments always crashes the kernel and sometimes displays this message:

No more memory available.
Mathematica kernel has shut down.
Try quitting other applications and then retry.

Why doesn't it work and how to fix this?

By the way, I am trying to recreate a higher quality version of the "Road to the Horizon" transformation of continued fractions from Linas Vepstas's math art gallery.

Update: Setting Parallelization -> False in the compiled function as per @Joel Klein's answer in 17971 works, and it seems this is a confirmed bug. However, in my case it does not seem to be fixed as reported in that thread.

shrx

Posted 2015-08-07T08:12:32.950

Reputation: 7 577

4FWIW, it also crashes on Mathematica 9.0.1, but not on 8.0.4. – Oleksandr R. – 2015-08-07T10:53:12.513

1I can also confirm it crashes with Mathematica 10.2. – shrx – 2015-08-07T11:59:55.260

1As a workaround, it does not crash for me if you change the type of xcf to Real and use Parallelization -> False. – rcollyer – 2015-08-07T13:08:30.870

1@rcollyer is changing the type of xcf even necessary? See the update in my question, disabling parallelization alone prevents the crash (but slows down the execution). – shrx – 2015-08-07T13:47:16.400

Nope, not necessary. It appears I was changing two things at once, and did not isolate which actually caused it not to crash. – rcollyer – 2015-08-07T13:53:15.480

I'd be interested if the kernel crashes with Mathematica 8.04 if Parallelization is explicitly set to True. Can anyone check? – shrx – 2015-08-07T13:56:03.857

1@shrx no, it does not. – Oleksandr R. – 2015-08-07T14:08:55.637

I should say I didn't try this in 9.0.0 (I don't have it installed). So we do not know for sure that this appeared only in 9.0.1. – Oleksandr R. – 2015-08-07T14:12:07.527

6Bug reported internally. Thank you! – ilian – 2015-08-07T15:22:39.900

@OleksandrR. you are right, I edited the notice. – shrx – 2015-08-07T16:29:06.860

3@shrx 9.0.0 is correct. I checked before making my edit. – ilian – 2015-08-07T16:33:45.017

3Fixed in the development version. – ilian – 2015-08-11T15:16:49.380

Answers

4

This bug was speedily confirmed...

Bug reported internally. Thank you! – ilian Aug 7 at 15:22

and fixed...internally at least.

Fixed in the development version. – ilian Aug 11 at 15:16

dr.blochwave

Posted 2015-08-07T08:12:32.950

Reputation: 8 483

2

This answer is not intended to address the bug described in the OP. Rather, my goal is to address

Since ContinuedFraction cannot be compiled...

I don't know the exact specifics of what goes on inside ContinuedFraction[], so the implementation I am about to present will slightly differ in the results. Nevertheless, I think it would be useful for me to present the algorithm I use on non-Mathematica systems.

The following algorithm is due to Pyzalski and Vala. Their algorithm is based on the elementary algorithm for generating the terms of a simple continued fraction; the wrinkle is in their choice of termination criterion for stopping the iteration. You can read their paper for more details. Here is a compiled routine implementing their algorithm:

ccf = With[{m = 1*^6}, Compile[{{x, _Real}},
           Module[{e = 1/m/(2 m - 1), d, f, h, n, nb, r, s, v, w},
                  r = IntegerPart[x]; s = Sign[x]; w = v = Abs[x - r];
                  f = {{1, 0}, {0, 1}};
                  nb = Internal`Bag[{r}];
                  While[True, h = 1/v; r = Floor[h];
                        Internal`StuffBag[nb, s r];
                        f = f.{{0, 1}, {1, r}};
                        {n, d} = f[[All, 2]];
                        If[Abs[w - n/d] >= e, v = h - r, Break[]]];
                  Internal`BagPart[nb, All]], RuntimeAttributes -> {Listable}]];

m is an adjustable parameter that controls the termination criterion; crudely stated, ccf stops when the denominator of the CF convergent being considered has $\approx\log_{10}(\mathtt m)$ digits.

Test:

ccf[{0.123455, 0.546452}]
   {{0, 8, 9, 1, 84, 4, 7}, {0, 1, 1, 4, 1, 7, 2, 7, 6, 5, 1, 3, 1}}

ContinuedFraction[{0.123455, 0.546452}]
   {{0, 8, 9, 1, 84, 4}, {0, 1, 1, 4, 1, 7, 2, 7, 6, 5, 1, 3, 1}}

ContinuedFraction[N[π]]
   {3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14}

ccf[N[π]]
   {3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3}

Even with the slight differences, I think ccf[] does a good job.

(Note that with some work, one can modify ccf[] to be a compiled substitute for Convergents[].)

J. M.'s ennui

Posted 2015-08-07T08:12:32.950

Reputation: 115 520