Sum of the Euler Totient function

May 10, 2018 Leave a comment

Given a positive integer {n}, the Euler totient function {\varphi(n)} is defined as the number of positive integers less than {n} which are co-prime with {n} (i.e. they have no common factors with {n}). There are formulas for computing {\varphi(n)} starting from the factorization of {n}. One such formula is

\displaystyle \varphi(n) = n \prod_{p|n} \left(1-\frac{1}{p}\right),

where the product is made over all primes dividing {n}.

If you have to compute {\varphi(n)} for all numbers less than a threshold then another property could be useful: {\varphi} is arithmetic, that is, {\varphi(mn) = \varphi(m)\varphi(n)} whenever {\gcd(m,n)=1}. Therefore you could store all values computed until {k} and for computing the value {\varphi(k+1)} there are two possibilities: {k+1=p^\alpha} is a prime power and then {\varphi(k+1) = p^{\alpha}-p^{\alpha-1}} or {k+1} is composite and {k+1 = mn} with {m,n\leq k,\ \gcd(m,n)=1}. Then use the stored values to compute {\varphi(k+1)=\varphi(m)\varphi(n)}.

I now come to the main point of this post: computing the sum of all values of the totient function up to a certain {N}:

\displaystyle \text{Compute } S(N) = \sum_{i=1}^N \varphi(i).

One approach is to compute each {\varphi(i)} and sum them. I will call this the brute-force approach. For all numerical purposes I will use Pari-GP in this post. On my computer it takes less than a second to compute {S(10^6)} and about {12} seconds to compute {S(10^7)}. This is super linear in time, since the algorithm computes the factorization for each {n} and then sums the values. Using the sieve approach could improve the timing a bit, but the algorithm is still super linear.

In some Project Euler problems it is not uncommon to have to compute something like {S(10^{11})} or even larger. Therefore, there must be more efficient ways to compute {S(N)} out there, so let’s study some of the properties of {S(N)}. In another post I dealt with the acceleration of the computation of the sum of the divisor function.

We have {S(N) = \sum_{i=1}^N \varphi(i)} which is the number of pairs {(a,b)} with {1\leq a\leq b \leq N} such that {\gcd(a,b) =1}. It is not difficult to see that the total number of such pairs is {n(n+1)/2}. Moreover, the possible values of {\gcd(a,b)} are {1,2,...,N}. Now, if for {m \leq N} we search instead for pairs satisfying {\gcd(a,b)=m} then we have {a = ma',\ b = mb'} with {\gcd(a',b')=1} and we get

\displaystyle 1 \leq a' \leq b' \leq N/m,\ \gcd(a',b')=1.

There fore the number of pairs with gcd equal to {m} is {S(\lfloor N/m\rfloor )}. Now we arrive at an interesting recursive formula:

\displaystyle S(N) = \frac{n(n+1)}{2} - \sum_{m=2}^N S(\lfloor N/m \rfloor ).

At a first sight this looks more complicated, but there is a trick to keep in mind whenever you see a summation over {m} of terms of the form {\lfloor N/m \rfloor}: these quantities are constant on large intervals. Indeed,

\displaystyle \lfloor N/m \rfloor = d \Leftrightarrow md \leq N < m(d+1)\Leftrightarrow N/(d+1)<m \leq N/d.

Therefore we can change the index of summation from {m} to {d=\lfloor N/m \rfloor}. The range of {d} for which the interval {I_d = [N/(d+1),N/d]} contains more than one integer is of order {\sqrt(N)}. Indeed, {N/d-N/(d+1) = N/(d(d+1))}. Therefore for {d\leq \sqrt(N)} we should have at least one integer in the interval {I_d}. The part where {d} is larger than {\sqrt{N}} corresponds to {m} smaller than {\sqrt{N}}. Therefore, we can split {S(N)} into two sums, each of order {\sqrt{N}}. and get that

\displaystyle S(N) = \frac{n(n+1)}{2}- \sum_{m=2}^{\sqrt{N}} S(\lfloor N/m \rfloor)-\sum_{d=1}^{\sqrt{N}}\left(\lfloor N/d \rfloor-\lfloor N/(d+1)\rfloor\right)S(d),

where in the last sum we must make sure that {d \neq \lfloor N/d \rfloor} in order to avoid duplicating terms in the sum.

Therefore we replaced a sum until {N} to two sums with upper bound {\sqrt{N}}. The complexity is not {\sqrt{N}}, but something like {N^{2/3}} since we have a recursive computation. Nevertheless, with this new formula and using memoization, to keep track of the values of {S} already computed, we can compute {S(N)} very fast:

{S(10^6) = 303963552392} is computed instantly (vs {1} second with brute force)

{S(10^7) = 30396356427242} takes {1} second (vs {12} seconds with brute force)

{S(10^8)} takes {5} seconds (vs over {3} minutes with brute force)

{S(10^9)} takes {30} seconds

{S(10^{11})} takes about {12} minutes

etc. Recall that these computations are done in Pari GP, which is not too fast. If you use C++ you can compute {S(10^8)} in {0.2} seconds, {S(10^9)} in {1} second and {S(10^{10})} in {6} seconds and {S(10^{11})} in under a minute, if you manage to get past overflow errors.

Advertisements

FreeFem++ Tutorial – Part 2

Click here for the first part. Some other posts related to FreeFem: link 1, link 2.

Here are a few tricks, once you know the basics of FreeFem. If your plan is straightforward: define the domain, build the mesh, define the problem, solve the problem, plot the result… then things are rather easy. If you want to do stuff in a loop, like an optimization problem, things may get more complicated, but FreeFem still has lots of tricks up its sleeves. I’ll go through some of them.

  1. Defining the geometry of the domain using border might be tricky at the beginning. Keep in mind that the domain should be on the left side of the curves definining the boundaries. This could be acheived in the parametrization chosen or you could reverse the parametrization of a particular part of the domain in the following way. Suppose you have the pieces of boundary C1,C2,C3,C4. You wish to build a mesh with the command mesh Th = buildmesh(C1(100)+C2(100)+C3(100)+C4(100));but you get an error concerning a bad sens on one of the boundaries. If you identify, for example, that you need to change the orientation of C3, you can acheive this by changing the sign of the integer defining the number of points on that part of the boundary: mesh Th = buildmesh(C1(100)+C2(100)+C3(-100)+C4(100));You could make sure that you have the right orientations and good connectivities for all the boundaries by running a plot command before trying mesh: something like

    plot(C1(100)+C2(100)+C3(-100)+C4(100));

    should produce a graph of all your boundaries with arrows showing the orientations. This is good as a debug tool when you don’t know where the error comes from when you define the domain.

  2. Keep in mind that there are also other ways to build a mesh, like square, which meshes a rectangular domain with quite a few options and trunc which truncates or modifies a mesh following some criteria. Take a look in the documentation to see all the options.
  3. Let’s say that you need to build a complex boundary, but with many components with similar properties. Two examples come to mind: polygons and domains with many circular holes. In order to do this keep in mind that it is possible to define a some kind of “vectorial boundary”. Let’s say that you want to mesh a polygon and you have the coordinates stored in the arrays xs,ys. Furthermore, you have another array of integers ind which point to the index of the next vertex. Then the boundary of the polygon could be defined with the following syntax:
    border poly(t=0,1; i){
    x=(1-t)*xx[i]+t*xx[ind(i)];
    y=(1-t)*yy[i]+t*yy[ind(i)];
    label=i;
    }
    Now the mesh could be constructed using a vector of integers NC containing the number of desired points on each of the sides of the polygon:

    mesh Th = buildmesh (poly(NC));
  4. You can change a 2D mesh using the command adaptmesh. There are various options and you’ll need to search in the documentation for a complete list. I use it in order to improve a mesh build with buildmesh. An example of command which gives good results in some of my codes is:

    Th = adaptmesh(Th,0.02,IsMetric=1,nbvx=30000);

    The parameters are related to the size of the triangles, the geometric properties of the griangles and the maximal number of vertices you want in your mesh. Experimenting a bit might give you a better idea of how this command works in practice.
  5. There are two ways of defining the problems in FreeFem. One is with solve and the other one is with problem. If you use solve then FreeFem solves the problem where it is defined. If you use problem then FreeFem remembers the problem as a variable and will solve it whenever you call this variable. This is useful when solving the same problem multiple times. Note that you can modify the coefficients of the PDE and FreeFem will build the new problem with the updated coefficients.
  6. In a following post I’ll talk about simplifying your code using macros and functions. What is good to keep in mind is that macros are verbatim code replacements, which are quite useful when dealing with complex formulas in your problem definition or elsewhere. Functions allow you to run a part of the code with various parameters (just like functions in other languages like Matlab).

I’ll finish with a code which computes the eigenvalue of the Laplace operator on a square domain with multiple holes. Try and figure out what the commands and parameters do.

int N=5; //number of holes
int k=1; //number of eigenvalue

int nb = 100;  // parameter for mesh size
verbosity = 10; // parameter for the infos FreeFem gives back

real delta = 0.1;
int nbd = floor(1.0*nb/N*delta*2*pi);

int bsquare = 0;
int bdisk   = 1;

// vertices of the squares
real[int] xs(N^2),ys(N^2);

real[int] initx = 0:(N^2-1);
for(int i=0;i&lt;N^2;i++){
  xs[i] = floor(i/N)+0.5;
  ys[i] = (i%N)+0.5;
}

cout &lt;&lt; xs &lt;&lt; endl;
cout &lt;&lt; ys &lt;&lt; endl;

int[int] Nd(N^2);

Nd = 1;
Nd = -nbd*Nd;

cout &lt;&lt; Nd &lt;&lt; endl;

// sides of the square
border Cd(t=0,N){x = t;y=0;label=bsquare;}
border Cr(t=0,N){x = N;y=t;label=bsquare;}
border Cu(t=0,N){x = N-t;y=N;label=bsquare;}
border Cl(t=0,N){x = 0;y=N-t;label=bsquare;}

border disks(t=0,2*pi; i){
    x=xs[i]+delta*cos(t);
    y=ys[i]+delta*sin(t);
    label=bdisk;
}

plot(Cd(nb)+Cr(nb)+Cu(nb)+Cl(nb)+disks(Nd));

mesh Th = buildmesh(Cd(nb)+Cr(nb)+Cu(nb)+Cl(nb)+disks(Nd));

plot(Th);

int[int] bc = [0,1]; // Dirichlet boundary conditions

load &quot;Element_P3&quot;;
fespace Vh(Th,P2);
// variables on the mesh
Vh u1,u2;

// Define the problem in weak form
varf a(u1,u2) = int2d(Th)
(dx(u1)*dx(u2) + dy(u1)*dy(u2))+on(1,u1=0)//on(C1,C2,C3,C4,u1=1)
                               +on(bc,u1=0);
varf b([u1],[u2]) = int2d(Th)(  u1*u2 ) ;
// define matrices for the eigenvalue problem
matrix A= a(Vh,Vh,solver=Crout,factorize=1);
matrix B= b(Vh,Vh,solver=CG,eps=1e-20); 

// we are interested only in the first eigenvalue
int eigCount = k;
real[int] ev(eigCount); // Holds eigenvalues
Vh[int] eV(eigCount);   // holds eigenfunctions
// Solve Ax=lBx
int numEigs = EigenValue(A,B,sym=true,sigma=0,value=ev,vector=eV);

cout &lt;&lt; ev[k-1] &lt;<span id="mce_SELREST_start" style="overflow:hidden;line-height:0;"></span>&lt; endl;

plot(eV[k-1],fill=1,nbiso=50,value=1);

Here is the mesh and the solution given by the above program:

MeshEig5

Spherical triangles of area Pi

April 6, 2018 Leave a comment

Recently I stumbled upon this page and found out a very nice result:

If a spherical triangle has area \pi then four copies of it can tile the sphere.

Here we are talking about triangles on the unit sphere whose edges are geodesics. The above result is a simple consequence of the following facts:

  1. If a spherical triangle has angles \alpha,\beta,\gamma then its area is \alpha+\beta+\gamma-\pi. Therefore if a triangle has area \pi, then \alpha+\beta+\gamma = 2\pi.
  2. If ABC is a triangle of area \pi and D is obtained by symmetrizing A with respect to the midpoint of BC on the sphere, then DCB and ABC are congruent triangles.
  3. Using angles and the fact that on the sphere similar triangles are congruent, we obtain that the triangles ABC,BCD,CDA,DAC are all congruent.

Here are a few examples of such partitions:

Mesh a hollow sphere

April 5, 2018 Leave a comment

This is an interesting experiment I’ve done today and I’d like to share it. There are nice software out there which allow you to mesh regular surfaces, like the sphere, torus or other examples.

I needed something different: I wanted to mesh a sphere with a hole and an inner surface. This kind of surfaces may be useful for those who want to design objects for 3D printing. When using 3D printers, the cost is usually proportional to the volume of material used, since this also consumes impression time. If you manage to make a hollow object, the cost goes down immediately (as the resistance of the object…). However, to make a hole, you need to consider an inner parallel surface. I’ll not handle this here, since it’s more complicated. I’ll just mesh a hollow sphere.

I did this in Matlab using the nice library Distmesh, which is free and really easy to use. The whole difficulty was making a function whose zero level set is the hollow sphere. Recall that it is possible to find the union of two shapes using the minimum of two level set function and the intersection using the maximum. This is all that you need to understand what I did below. You can see the details in the Matlab code below.

function [p,t] = Hole_Spher_Mesh(mh)

  R = 0.5;
  r = 0.4;
  rh = 0.1;

  fh=@(p) sqrt(min(p(:,1).^2+p(:,2).^2,ones(size(p(:,1)))));
  fh=@(p) sqrt((p(:,1).^2+p(:,2).^2));
  [p,t]=distmeshsurface(@(x)fd(x,R,r,rh),@huniform,mh,[-1.1,-1.1,-1.1;1.1,1.1,1.1]);
  points = p';
  npt = prod(size(points))/3

clf
patch('Faces',t,'Vertices',p,'FaceColor',[0,0.7,0.9],'EdgeColor','none','FaceAlpha',0.5);

function res = fd(p,R,r,rh)

res1 = sqrt(sum(p.^2,2))-R;
res2 = (rh-sqrt(p(:,1).^2+p(:,2).^2));
res3 = p(:,3);
res2 = max(-res2,res3);
res4 = sqrt(sum(p.^2,2))-r;
res2 = min(res2,res4);
res = max(res1,-res2);

And here is the resulting surface for the parameter mh = 0.01 (if you put a higher precision computations will take longer).

HollowSphere

When is arccos a rational multiple of pi?

March 31, 2018 Leave a comment

Recently I had to test if the {\arccos} of some algebraic number is a rational multiple of {\pi} or not. I found this solution online and since it is quite nice, I write it here in a more generalized form.

Suppose that {r} is an algebraic number. Under which circumstances we have

\displaystyle \arccos(r) \in \pi\Bbb{Q}?

Suppose that {\arccos(r) = \frac{p}{q} \pi} with {p,q \in \Bbb{Z}, q > 0}. Then {q\arccos(r) = p\pi} and applying {\cos} we see that

\displaystyle \cos(q\arccos(r)) -(-1)^p = 0.

It is well known that the function {P_q\ :\ x \mapsto \cos(q\arccos(x))} is polynomial for {|x|\leq 1}. If we denote {f(x) = P_q(x)-(-1)^p} then {f} is a polynomial of degree {q} with leading coefficient {2^{q-1}}. Since we have {f(r) = 0}, the minimal polynomial {Q_r} associated to {r} must divide {f(r)}. This gives us the following necessary condition:

If {\arccos(r) \in \pi\Bbb{Q}} then there exist {p,q \in \Bbb{Z}}, {q > 0} such that {Q_r} (the minimal polynomial of {r}) divides {P_q\pm 1} (where {P_q} is the Chebyshev polynomial of degree {q}).

It is not difficult to see that this condition is also sufficient. Now, using this we have the following application:

Consequence. If {r} is an algebraic number whose minimal polynomial has a leading coefficient which is not a power of {2}, then

\displaystyle \arccos(r) \notin \pi \Bbb{Q}.

Since the minimal polynomial {Q_r} has a leading coefficient which is not a power of {2}, it cannot divide {P_q\pm 1} for some positive integer {q}. Therefore, following the arguments stated above, {\arccos(r) \notin \pi\Bbb{Q}}.

Applications.

  • {\arccos(\sqrt{n}/3)\notin \pi\Bbb{Q}} if {\gcd(n,3) \leq 3}.
  • {\arccos(\sqrt{m/n}) \notin \pi\Bbb{Q}} for {m,n \in \Bbb{N}^*} if {n>1},\ \gcd(m,n)=1 and {n} is not a power of {2}.

Romanian Masters in Mathematics contest – 2018

March 9, 2018 Leave a comment

Problem 1. Let {ABCD} be a cyclic quadrilateral an let {P} be a point on the side {AB.} The diagonals {AC} meets the segments {DP} at {Q.} The line through {P} parallel to {CD} mmets the extension of the side {CB} beyond {B} at {K.} The line through {Q} parallel to {BD} meets the extension of the side {CB} beyond {B} at {L.} Prove that the circumcircles of the triangles {BKP} and {CLQ} are tangent .

Problem 2. Determine whether there exist non-constant polynomials {P(x)} and {Q(x)} with real coefficients satisfying

\displaystyle P(x)^{10}+P(x)^9 = Q(x)^{21}+Q(x)^{20}.

Problem 3. Ann and Bob play a game on the edges of an infinite square grid, playing in turns. Ann plays the first move. A move consists of orienting any edge that has not yet been given an orientation. Bob wins if at any point a cycle has been created. Does Bob have a winning strategy?

Problem 4. Let {a,b,c,d} be positive integers such that {ad \neq bc} and {gcd(a,b,c,d)=1}. Let {S} be the set of values attained by {\gcd(an+b,cn+d)} as {n} runs through the positive integers. Show that {S} is the set of all positive divisors of some positive integer.

Problem 5. Let {n} be positive integer and fix {2n} distinct points on a circle. Determine the number of ways to connect the points with {n} arrows (oriented line segments) such that all of the following conditions hold:

  • each of the {2n} points is a startpoint or endpoint of an arrow;
  • no two arrows intersect;
  • there are no two arrows {\overrightarrow{AB}} and {\overrightarrow{CD}} such that {A}, {B}, {C} and {D} appear in clockwise order around the circle (not necessarily consecutively).

Problem 6. Fix a circle {\Gamma}, a line {\ell} to tangent {\Gamma}, and another circle {\Omega} disjoint from {\ell} such that {\Gamma} and {\Omega} lie on opposite sides of {\ell}. The tangents to {\Gamma} from a variable point {X} on {\Omega} meet {\ell} at {Y} and {Z}. Prove that, as {X} varies over {\Omega}, the circumcircle of {XYZ} is tangent to two fixed circles.

Source: Art of Problem Solving forums

Some quick ideas: For Problem 1 just consider the intersection of the circle {(BKP)} with the circle {(ABCD)}. You’ll notice immediately that this point belongs to the circle {(CLQ)}. Furthermore, there is a common tangent to the two circles at this point.

For Problem 2 we have {10\deg P = 21 \deg Q}. Eliminate the highest order term from both sides and look at the next one to get a contradiction.

Problem 4 becomes easy after noticing that if {q} divides {an+b} and {cn+d} then {q} divides {ad-bc}.

In Problem 5 try to prove that the choice of start points determines that of the endpoints. Then you have a simple combinatorial proof.

Problem 6 is interesting and official solutions use inversions. Those are quite nice, but it may be worthwhile to understand what happens in the non-inverted configuration.

I will come back to some of these problems in some future posts.

SEEMOUS 2018 – Problems

March 1, 2018 Leave a comment

Problem 1. Let {f:[0,1] \rightarrow (0,1)} be a Riemann integrable function. Show that

\displaystyle \frac{\displaystyle 2\int_0^1 xf^2(x) dx }{\displaystyle \int_0^1 (f^2(x)+1)dx }< \frac{\displaystyle \int_0^1 f^2(x) dx}{\displaystyle \int_0^1 f(x)dx}.

Problem 2. Let {m,n,p,q \geq 1} and let the matrices {A \in \mathcal M_{m,n}(\Bbb{R})}, {B \in \mathcal M_{n,p}(\Bbb{R})}, {C \in \mathcal M_{p,q}(\Bbb{R})}, {D \in \mathcal M_{q,m}(\Bbb{R})} be such that

\displaystyle A^t = BCD,\ B^t = CDA,\ C^t = DAB,\ D^t = ABC.

Prove that {(ABCD)^2 = ABCD}.

Problem 3. Let {A,B \in \mathcal M_{2018}(\Bbb{R})} such that {AB = BA} and {A^{2018} = B^{2018} = I}, where {I} is the identity matrix. Prove that if {\text{tr}(AB) = 2018} then {\text{tr}(A) = \text{tr}(B)}.

Problem 4. (a) Let {f: \Bbb{R} \rightarrow \Bbb{R}} be a polynomial function. Prove that

\displaystyle \int_0^\infty e^{-x} f(x) dx = f(0)+f'(0)+f''(0)+...

(b) Let {f} be a function which has a Taylor series expansion at {0} with radius of convergence {R=\infty}. Prove that if {\displaystyle \sum_{n=0}^\infty f^{(n)}(0)} converges absolutely then {\displaystyle \int_0^{\infty} e^{-x} f(x)dx} converges and

\displaystyle \sum_{n=0}^\infty f^{(n)}(0) = \int_0^\infty e^{-x} f(x).

Source: official site of SEEMOUS 2018 

Hints: 1. Just use 2f(x) \leq f^2(x)+1  and xf^2(x) < f^2(x). The strict inequality comes from the fact that the Riemann integral of strictly positive function cannot be equal to zero. This problem was too simple…

2. Use the fact that ABCD = AA^t, therefore ABCD is symmetric and positive definite. Next, notice that (ABCD)^3 = ABCDABCDABCD = D^tC^tB^tA^t = (ABCD)^t = ABCD. Notice that ABCD  is diagonalizable and has eigenvalues among -1,0,1. Since it is also positive definite, -1 cannot be an eigenvalue. This allows to conclude.

3. First note that the commutativity allows us to diagonalize A,B  using the same basis. Next, note that A,B both have eigenvalues of modulus one. Then the trace of AB is simply the sum \sum \lambda_i\mu_i where \lambda_i,\mu_i are eigenvalues of A and B, respectively. The fact that the trace equals 2018  and the triangle inequality shows that eigenvalues of A are a multiple of eigenvalues of B. Finish by observing that they have the same eigenvalues.

4. (a) Integrate by parts and use a recurrence. (b) Use (a) and the approximation of a continuous function by polynomials on compacts to conclude.

I’m not sure about what others think, but the problems of this year seemed a bit too straightforward.

%d bloggers like this: