# How do I plot Thomae's function in Mathematica?

I wanted to plot this function

$$f(x) =\begin{cases} 1 & \text{if } x= 0 \\ \tfrac1{q} & \text{if } x = \tfrac{p}{q}\\ 0 & \text{if } x \in \mathbb{R}-\mathbb{Q} \end{cases}$$

so I wrote

FuncThomae[x_] := If[ExactNumberQ[Rationalize[x]], If[x == 0, a = 1,
L = #^-1 & /@ Divisors[Numerator[Rationalize[x]]]], a = 0]


and

ListDomain[xmin_, xmax_] := Table[Outer[List, {x}, FuncThomae[x]],
{x,xmin,xmax,0.001}] // Flatten[#, 1] &


My result doesn't take all the real numbers (nor negatives) in its domain, and for $-1$ to $1$, it should have looked like so:

but my function does not cater to negatives, nor does it look like the above plot. It looks like this from 2 to 5:

Nearly similar, but not quite. Can someone help me to perfect the function?

Wanted to duplicate this

-
Hi there, is the code you supplied complete, I didn't seem to be able to get it to run successfully ? –  image_doctor May 13 '12 at 11:38
@image_doctor it does run. what was the problem? BTW you need to plot the generated list using the ListDomain function i coded as argument for the ListPlot Function..... –  The-Ever-Kid May 13 '12 at 12:42
Outer::ipnfm: "Positive machine-sized integer or Infinity expected at position 3 in Outer[List,{2.},0]." Is the error I get. a and L seem not to be defined in the code segment you have posted. But as you have an answer, perhaps it isn't important :) –  image_doctor May 13 '12 at 17:31
A tiny reminder: in floating-point arithmetic (which is what Mathematica internally uses when plotting), all numbers are rational. –  Ｊ. Ｍ. May 14 '12 at 11:15

I'd suggest producing a list of rational numbers and then plot the function there, like so:

maxq = 100;
fracs = Table[p/q, {q, 2, maxq}, {p, 2, q}] // Flatten // DeleteDuplicates;
pq = {#, 1/Denominator @ #} & /@ fracs;

ListPlot[pq, PlotRange -> {0, 1}]


-
If you define pq as pq = {#, 1/Denominator[#]}& /@ fracs you could ListPlot pq directly without having to transform it first. –  Heike May 13 '12 at 12:39
@acl thanks for the lovely answer but the plot is missing the first bit it isnt plotting 1 when x is 0 but anyway thanks.... –  The-Ever-Kid May 13 '12 at 13:19
btw could you do this for -1 to 1 –  The-Ever-Kid May 13 '12 at 13:20
@The-Ever-Kid yes that bit is straightforward, so I didn't bother implementing it –  acl May 13 '12 at 13:20
yeah cool thanks... –  The-Ever-Kid May 13 '12 at 13:23

Another possibility that avoids the generation of fractions not in lowest terms (and thus the use of Union[] or DeleteDuplicates[]) rests on generating a Farey sequence, and then applying the Dirichlet-Thomae function to that:

farey[n_Integer?Positive] := Module[{v = 0, w = 1, p = 1, q = n, t},
Join[{0, 1/n}, Flatten[Last[
Reap[While[p < q,
t = Quotient[n + w, q];
{{p, v}, {q, w}} = {{t p - v, p}, {t q - w, q}};
Sow[p/q]]]
]]]]

ListPlot[{#, 1/Denominator[#]} & /@ farey[100], Frame -> True, PlotRange -> {0, 1}]


-
Why not use DiscretePlot directly?
DiscretePlot[1/Denominator[Rationalize[x]], {x, -1, 1, 1/(2*3*4*5*6*7)},

The RegionFunction throws out the cases where $\frac{p}{q}=1$.