2
$\begingroup$

so I am extremely new to Mathematica and currently doing a project that involves calculating a specific function of the Farey Fractions. That is given an interval $I: (0,\alpha)$, I want to calculate the number of distinct pairs $i$,$j$ s.t $i-j$ $\in I$ in a normalized list of farey fractions. The normalization factor here is $\frac{1}{|F_d-1|}$, where $|F_d|$ is the size of the list of farey fractions with denominator d. So, for example with the denominator being 5, the normal list of Farey fractions is $FareySequence[5]$, while the normalized list is $|F_d-1|(FareySequence[5])$.

Below is my attempt to count the number of distinct pairs $i,j$ which are part of a normalized Farey Sequence such that $i-j<1$, where $i>j$. Here, I am counting it based on fractions with denominator 1000. There are 304193 such fractions. I used a counter, where a count is added everytime there is a pair whose difference is less than 1. However, everytime I run it, the program never evaluates. Any help would be appreciated.

count = 0

For[i = 1, i < 304193, i++,
 
For[j = 1, j < i, j++,
 
 if[(304193 FareySequence[1000, i]) - (304193 FareySequence[1000, 
        j]) < 1, count ++, count = count ]
  
]
 ]

Print[count]
$\endgroup$
3
  • $\begingroup$ I tried it, it still doesn't work $\endgroup$ – NotSoTrivial Oct 12 '20 at 16:08
  • 1
    $\begingroup$ I think it's working fine it's just going to take forever. First pull the calculation of (304193 FareySequence[1000, i]) out of the j loop because it's just wasting time recalculating the same number every time you increment j, then try a smaller number for the outer loop. Even setting i<500 takes 2 minutes on my machine, and it will only take longer as i increases beyond that. You'll likely need a cleverer calculation if you want this to happen in any reasonable time. $\endgroup$ – N.J.Evans Oct 12 '20 at 16:56
  • $\begingroup$ Looks like @N.J.Evans has the right analysis on this. See my extended comment below. $\endgroup$ – Jagra Oct 12 '20 at 17:06
4
$\begingroup$

Here's solution that works quickly:

FareySequence[1000] is only 304193 elements long, plenty short enough to store in memory, so just calculate the sequence once then access the elements for a speed up.

fs = 304193*FareySequence[1000];

Then notice that the elements of FareySequence are monotonically increasing, so if you take element fs[[i]] and the preceding element, fs[[i-1]] has a distance father than 1 from f[[i]] you don't need to test any of the elements earlier in the sequence. So instead of testing j from 0 to i, test j from i-1 down and break every time your condition isn't met.

For[
 i = 1, i < 304193, i++,
 For[j = (i - 1), j > 0, j--,
  If[(fs[[i]] - fs[[j]]) < 1, count++, Break[]]
  ];
 ]

Print[count]

232074

This completes in just under 2 seconds.

Here's a more mathematica-like ('functional') solution that doesn't work very quickly:

First write a function that takes a sequence then uses TakeWhile to traverse all of the sequence except for the last element using the last element for the comparison and breaking if the difference between the two elements is larger than lim, Length counts the length of this list.

countBelowWithin[seq_, lim_:1] := 
  Length@TakeWhile[Reverse[Most[seq]], Last[seq] - # < lim &];

Then use Table to take subsequences of the full sequence and pass them to countBelowWithin. The resulting table will be the number of elements that met the criteria for that subsequence, so Total adds those for a total number of elements meeting the criteria.

With[
 {a = fs[[;; 10000]]},
 Total[
  Table[
   countBelowWithin[a[[;; i]]]
   , {i, 2, Length@a}
   ]
  ]
 ]

I use the With to isolate a portion of the sequence fs for testing purposes.

Very often the functional solution can be quicker - in this case it's not. It is often less error prone as you can think in bigger concepts than indices. In this case you'd probably need to think a little more if you wanted a functional solution that worked in a reasonable time.

If you interested in learning more, there's a lot of fundamentals in here Part i.e. ...[[;;]], Table, TakeWhile. I'd start just reading the docs on this stuff, as well as the links other posters provided.

$\endgroup$
1
  • 1
    $\begingroup$ Very nicely conceived and done ;-) $\endgroup$ – Jagra Oct 12 '20 at 19:21
2
$\begingroup$

More an extended comment than an answer...

If you reduce the limit of your i iterator to something small, e.g., 10 rather than 304193, the corrected code ("IF" rather than "if") seems to work fine.

Try the following:

count = 0;
For[i = 1, i < 10, i++,
 For[j = 1, j < i, j++,
  If[
   (304193 *FareySequence[1000, i]) - (304193 * FareySequence[1000, j]) < 1,
   count++, count = count]]]
Print[count]

which outputs: 18

It appears that setting the iterative limit of i at 304193 just produces an enormous number of evaluations in the nested For loops.

I think you can produce your result in far less time by "functional programing" means (the real strength of Mathematica) rather than iterative looping.

Simply Mapping the IF statement across the appropriate range(s) will likely pick up a great deal of speed. Then do a separate "count" of the instances you want to identify.

I'll give some thought to this later in the day and come back if no one else as supplied an answer yet.

$\endgroup$
3
  • $\begingroup$ I see now, I guess the time complexity is complicating things. Could you explain what you mean by function programming, and how I can implement this for my problem? I only have been using Mathematica for 3 days so I'm very new to this. Thanks. $\endgroup$ – NotSoTrivial Oct 12 '20 at 17:11
  • 2
    $\begingroup$ I suggest starting with: mathematica.stackexchange.com/questions/7924/… $\endgroup$ – Jagra Oct 12 '20 at 17:19
  • 2
    $\begingroup$ Also: Leonid's ebook: Mathematica Programming: An Advanced Introduction -- by Leonid Shifrin, free download at wolfram.com/books/profile.cgi?id=8101 $\endgroup$ – Jagra Oct 12 '20 at 17:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.