Message Boards

WOLFRAM COMMUNITY

11101 Views

5 Replies

16 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Staff Picks Mathematics Recreation Mathematica Compilation Discrete Mathematics Graphics and Visualization Wolfram Language Computer-Based Maths

[Numberphile] - Amazing Graphs II

Sander Huisman

Sander Huisman, University of Twente

Posted 5 years ago

This is part of a series where I explore some of the videos of Numberphile, see also the other ones: [Numberphile] - The Illumination Problem [Numberphile] - Frog Jumping - Solving the puzzle [Numberphile] - Sandpiles - Done in the Wolfram Language [Numberphile] - The Square-Sum Problem [Numberphile] - How to make a binary prime logo Today we are gonna look at the graphs from https://www.youtube.com/watch?v=o8c4uYnnNnc The first can be really easily made: n = 3^5 - 1; nums = FromDigits[IntegerDigits[#, 3] /. 2 -> -1, 3] & /@ Range[n]; ListPlot[nums] The second can also be easily made using the following code: n = 1000; nums = # - (Times @@ DeleteCases[IntegerDigits[#], 0]) & /@ Range[n]; ListPlot[nums] giving: The last sequence is harder to program, especially if one wants a fast solution Here is the code: x = ConstantArray[1, 10^5]; x[[;; 2]] = {1, 1}; Dynamic[i] AbsoluteTiming@Do[ max = Floor[(i - 1)/2]; invalid = 2 x[[i - max ;; i - 1]] - x[[i - 2 max ;; i - 2 ;; 2]]; invalid = Pick[invalid, 1 - UnitStep[-invalid], 1]; (* remove nonpositive *) invalid = Sort[DeleteDuplicates[invalid]]; new = -1; If[Last[invalid] == Length[invalid], new = Length[invalid] + 1; , Do[ If[invalid[[k]] =!= k, new = k; Break[]; ] , {k, 1, Length[invalid]} ]; If[new == -1, new = Last[invalid] + 1]; ]; x[[i]] = new; , {i, 3, Length[x]} ] ListPlot[x] giving: Hope you enjoyed these codes, perhaps you can modify them and make them more intricate/faster/better!

POSTED BY: Sander Huisman

5 Replies

Sort By:

Shenghui Yang

Shenghui Yang, WOLFRAM

Posted 5 years ago

I feel all plots in both episodes are related to something A New Kind of Science. Meanwhile I put your code directly within a `Compile` function and the speedup is definitely noticeable with minimum modification: cf = Compile[{}, Module[{x, max, invalid, new}, x = ConstantArray[1, 10^5]; Do[ max = Floor[(i - 1)/2]; invalid = 2 x[[i - max ;; i - 1]] - x[[i - 2 max ;; i - 2 ;; 2]]; invalid = Pick[invalid, 1 - UnitStep[-invalid], 1];(remove nonpositive) invalid = Sort[DeleteDuplicates[invalid]]; new = -1; If[Last[invalid] == Length[invalid], new = Length[invalid] + 1;, Do[If[invalid[[k]] =!= k, new = k; Break[];], {k, 1, Length[invalid]}]; If[new == -1, new = Last[invalid] + 1];]; x[[i]] = new;, {i, 3, Length[x]}]; x ] , CompilationTarget -> "C"]

I feel all plots in both episodes are related to something A New Kind of Science. Meanwhile I put your code directly within a Compile function and the speedup is definitely noticeable with minimum modification:

cf = Compile[{},
  Module[{x, max, invalid, new},
   x = ConstantArray[1, 10^5];
   Do[
    max = Floor[(i - 1)/2];
    invalid = 2 x[[i - max ;; i - 1]] - x[[i - 2 max ;; i - 2 ;; 2]];
    invalid = 
     Pick[invalid, 1 - UnitStep[-invalid], 1];(*remove nonpositive*)
    invalid = Sort[DeleteDuplicates[invalid]];
    new = -1;
    If[Last[invalid] == Length[invalid], new = Length[invalid] + 1;, 
     Do[If[invalid[[k]] =!= k, new = k;
       Break[];], {k, 1, Length[invalid]}];
     If[new == -1, new = Last[invalid] + 1];];
    x[[i]] = new;, {i, 3, Length[x]}];
   x
   ]
  , CompilationTarget -> "C"]

speedup

POSTED BY: Shenghui Yang

Sander Huisman

Sander Huisman, University of Twente

Posted 5 years ago

Hi Shenghui! Thanks for checking the code out. Yes I did not try any compilation. I think it can be further optimized, because Pick, UnitStep and DeleteDuplicates are handled by the Wolfram Language again. They are not compiled to C But anyhow, it is a nice speed up. Note that I made it speed-up for Mathematica, I think a more 'raw' approach would be (a lot) faster in compiled code. Cheers! Sander

POSTED BY: Sander Huisman

Paul Cleary

Posted 5 years ago

Here is a snippet of code I wrote to test the Fly straight damnit from the Amazing graphs 1 episode. an = 1; lst = Reap[Do[gf = GCD[an, n]; If[gf == 1, an = an + n + 1, an = an/gf]; Sow[{n, an}], {n, 2, 1000}]]; lst = Flatten[lst[[2]], 1]; ListPlot[lst]

POSTED BY: Paul Cleary

EDITORIAL BOARD

EDITORIAL BOARD, WOLFRAM

Posted 5 years ago

- Congratulations! This post is now featured in our Staff Pick column as distinguished by a badge on your profile of a Featured Contributor! Thank you, keep it coming, and consider contributing your work to the The Notebook Archive!