Message Boards Message Boards

[Numberphile] - Amazing Graphs II

This is part of a series where I explore some of the videos of Numberphile, see also the other ones:

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]

enter image description here

The second can also be easily made using the following code:

n = 1000;
nums = # - (Times @@ DeleteCases[IntegerDigits[#], 0]) & /@ Range[n];
ListPlot[nums]

giving:

enter image description here

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:

enter image description here

Hope you enjoyed these codes, perhaps you can modify them and make them more intricate/faster/better!

POSTED BY: Sander Huisman
5 Replies

enter image description here - 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!

POSTED BY: EDITORIAL BOARD
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

Neat! Thanks for sharing!

POSTED BY: Sander Huisman

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

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
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract