This is really a response to Dror Speiser's comment (but it is too long to give as a comment), and gives an analytic $O(n^{1/3+\epsilon})$ time and $O(n^{1/6+\epsilon})$ space algorithm to count the number of square free numbers less than $x$ (thus improving on the complexity of the algorithm given in the Jakub Pawlewicz paper cited in the question). If we do exactly in the same way as in my other answer for
$$\sum_{n\leq x} |\mu(n)|$$
instead of
$$\sum_{n\leq x} \mu(n)$$
by using the fact that $$\frac{\zeta(s)}{\zeta(2s)} = \sum_{n=1}^\infty |\mu(n)| n^{-s}$$
it is true that the time complexity will be $O(x^{1/2+\epsilon})$. There is however a way to combine the analytic method (which is essentially the same as in Lagarias-Odlyzkos 1987 paper "Computing $\pi(x)$ an analytic method" http://www.dtc.umn.edu/~odlyzko/doc/arch/analytic.pi.of.x.pdf) with the elementary identity
$$
S(n)=\sum_{d=1}^{\lfloor \sqrt n\rfloor} \mu(d) \lfloor \frac n {d^2} \rfloor
$$
that is Pawleviec's starting point in obtaining his algorithm. Let $\Phi$ be defined as in my other answer, i.e. let $\Phi(x)$ be a smooth test function such that $\Phi(x)=1$ for $x<0$ and $\Phi(x)=0$ for $x>1$. Consider the identity
$$ \sum_{n\leq x} |\mu(n)| =\sum_{n=1}^\infty |\mu(n)| \Phi \left( \frac {n-x} {x^{2/3}} \right )
- \sum_{x< n < x + x^{2/3} } |\mu(n)| \Phi \left(\frac {n-x} {x^{2/3}} \right)$$
The first sum can be estimated by a complex integral of length $x^{1/3+\epsilon}$ and by the Odlyzko-Schönhage algorithm it can be calculated in $O(x^{1/3+\epsilon})$ time (and $O(x^{1/6+\epsilon})$ space). The remaining sum will have length $O(x^{2/3})$ so might at first glance not be so simple to treat. I claim however that in fact it can be calculated in $O(x^{1/3+\epsilon})$ time. This is where using the elementary identity is handy. Let $\Psi(x)=0$ for $x<0$ and $\Psi(x)=\Phi(x)$ for $x \geq 0$ (this will have a discontinuity at $x=0$) By a similar elementary identity we obtain
$$\sum_{x< n < x + x^{2/3} } |\mu(n)| \Phi \left(\frac {n-x} {x^{2/3}} \right) = \sum_{d=1}^{\lfloor \sqrt {x+x^{2/3}}\rfloor} \mu(d) \sum_{k=1}^\infty \Psi \left(\frac {d^2k-x} {x^{2/3}} \right).
$$
For any $d$ the inner sum can be calculated fast (certainly in $O(x^\epsilon)$ time), by either noticing that the sum is empty, just contains one element or the Poisson summation formula. The trick now is to see that there will only be $O(x^{1/3})$ integers $d$ where the inner sum (in $k$) is non empty, namely $d$ must either be $O(x^{1/3})$ or belong to an interval $(\sqrt{ x/k},\sqrt{(x+x^{2/3})/k})$ for $k \leq x^{1/3}$ and there are only $O(x^{1/3})$ such $d$. Now if we can calculate the values of $\mu(d)$ fast (simple sieving seems more difficult in this case) we obtain the desired algorithm. Determining $\mu(d)$ basically comes down to factoring the number. However it is sufficient to determine the factors less than $d^{1/3}$ since if we know that all prime factors are greater than $d^{1/3}$, then either it has one or two prime factors. We can determine if it is a prime fast by the Agrawal-Kayal-Saxena algorithm, then of course $\mu(d)=-1$. If it has two prime factors either $\mu(d)=0$, but then it has to be a square, and it can be easily checked, or $\mu(d)=1$. Now factoring the $O(x^{1/3})$ numbers of size $O(x^{1/2})$ into all prime factors less than $x^{1/6}$ will take $O(x^{1/3+\epsilon})$ time and $O(x^{1/6+\epsilon})$ space by Bernstein's algorithm http://cr.yp.to/papers/sf.pdf. I might possibly write up this answer as a paper proper. If anyone has seen these arguments anywhere else or have references I would be interested.
It depends on what you mean by "calculate".
If you want a "closed-form" expression for $k(m)$, the period of the Fibonacci sequence modulo $m$: One should expect that there is no such thing, because finding the period of the Fibonacci sequence modulo $m$ is a very similar problem to that of finding the period of the sequence $1,a,a^2,\ldots$ modulo $m$, a.k.a. the discrete logarithm problem, which is widely believed to be "hard". (The Fibonacci sequence is a 2nd-order linear recurrence sequence, while the powers of $a$ are a 1st-order linear recurrence sequence.)
As it sounds like you know already, if $m=p_1^{a_1}\cdots p_k^{a_k}$, then the period of a sequence modulo $m$ is the $\mathrm{lcm}$ of the periods modulo $p_1^{a_1}$, modulo $p_2^{a_2}$, etc. (This is a trivial consequence of the Chinese remainder theorem.)
The next natural step would be to reduce the problem from computing the Pisano periods for prime powers to computing them only for primes. As far as we know, $k(p^n)=p^{n-1}\cdot k(p)$ for any $n$ and any prime $p$, but it's still conjectured that there are infinitely many primes $p$ for which this fails. I don't think there's been any significant progress other than higher and higher computer checks, but then again, I really haven't been keeping up with this problem. This formula fails for a prime $p$ if and only if the prime is a Wall-Sun-Sun prime. The analog of such primes for the 1st-order linear recurrence sequences are Wieferich primes and their generalizations.
As far as finding a formula for $k(p)$ in terms of $p$, that should be hopeless, for the reason I mentioned at the start of my post.
If you want an algorithm, i.e. a program your computer can run to find out the period experimentally, then I'm afraid that's something I know nothing about (beyond the obvious naive one, of course); it would probably be helpful to you to look at the research papers of those mathematicians who have carried out the computer searches.
Best Answer
This is not directly an answer to your question (graph theoretic relation, rather than number theoretic), but it's too long for a comment.
Euler's formula is a special case of a disharmonicity function
$$D(x)=\sum_i |e_i| g(p_i)$$
with $g(p_i)>0$ a function of the prime factors $p_i$ of a rational number $x=p_1^{e_1} p_2^{e_2}\cdots $ (each with positive or negative multiplicity $e_i$). Euler's disharmonicity has $g(p_i)=p_i-1$ (and adds unity to the sum over $i$). An alternative disharmonicity, due to Barlow, takes $g(p_i)=(2/p_i)(p_i-1)^2$.
Disharmonicity functions can be used to define the notion of a "harmonic distance" of two rational numbers, and to formulate the problem of the harmonization of a musical scale as a problem in graph theory.
See Musical scale rationalization – a graph-theoretic approach, by Albert Gräf.