2

As the title says we are given an odd number $n$ and wish to find the number of bases $b$ such that $n$ is an Euler pseudoprime; That is, $\gcd(b,n)=1$ and $b^{(n-1)/2} \equiv \left( \frac{b}{n} \right) \mod n$. I wrote a simple function in Matlab to compute the Jacobi symbol and used this in a simple for-loop over the numbers that are coprime to $n$. I tested it for some random numbers and although the numbers seem to check out I feel that it somehow gives me the wrong answer. Here is my function:

function [count] = EulerPP(n)
count=0;
for i=2:(n-1)
    if gcd(i,n)==1
        if mod(modexp(i,(n-1)/2,n)-JacSym(i,n),n)==0
            count=count+1;
        end
    end
end
end

Here is the input/output for some numbers:

$\begin{align*} 1387 &\mapsto 126 \\ 1189 &\mapsto 6 \\ 1537 &\mapsto 7 \\ 105 &\mapsto 7 \\ 781 &\mapsto 55 \\ \end{align*} $

Also, is it possible to calculate this in a simpler way (i.e. by not checking all numbers)?

Edit: I have also included the code for the functions modexp and JacSym:

function [s] = JacSym(m,n)
%Jacobi symbol of (m/n) where n is always odd.

m=mod(m,n);

if mod(n,2)==0
    s=0;
elseif m==1;
    s=1;
elseif m==0
    s=0;

elseif mod(m,2)==0
    if abs(mod(n,8))==1
        s=JacSym(m/2,n);
    else
        s=-JacSym(m/2,n);
    end
else
    if mod(n,4)==3 && mod(m,4)==3
        s=-JacSym(n,m);
    else
        s=JacSym(n,m);
    end
end    
end


function result = modexp (x, y, n)
    %anything raised to 0th power = 1 so return 1
    if (y == 0)
        result = 1;
        return;
    end

    %recurse
    z = modexp(x, floor(y/2), n);

    %if even square the result
    if (mod(y, 2) == 0)
        result = mod(z*z, n);
        return;
    else
        %odd so square the result & multiply by itself
        result = mod(x*z*z, n);
        return;
    end
end

Found the latter one here: https://gist.github.com/ttezel/4635562

KryptoJon
  • 53
  • With Pari/GP, for $n=1387$ I get $161$, for $n=1189$ I get $7$, for $n=781$ I get $49$. Check your functions. Generally, for a given $n$, such $b$ including $b=1$ form a subgroup of $(\mathbb{Z}/n\mathbb{Z})^\times$. To take advantage of this, you will probably need to factor some numbers however. – ccorn Dec 12 '13 at 17:11
  • Your code has Matlab syntax. Maybe it has issues with a floating-point mod operation applied to integers? – ccorn Dec 12 '13 at 17:19
  • Is the base 1 usually included when counting the bases? I updated the first post to include the code for the functions called in the EulerPP function. It might be the case that there is a floating-point mod operation somewhere, but I cannot spot it ... – KryptoJon Dec 12 '13 at 17:47
  • Your JacSym(2,7) yields $-1$ but should yield $1$. – ccorn Dec 12 '13 at 17:52
  • Whether the base $1$ is usually included? I do not know. Certainly, it is not used in pseudoprimality tests. But you include $n-1$ although one may argue that $-1\pmod{n}$ is useless in pseudoprimality tests as well. I'd include $1$ in the base list because it fulfills the criteria and completes the subgroup. Besides, as we are talking about convention, I remember that only non-primes are called pseudoprimes although they typically pass the pseudoprime test. – ccorn Dec 12 '13 at 18:04
  • When testing larger $n$, keep in mind that the x*z*z in your modexp function typically uses IEEE double precision ($53$ bits) which may lose accuracy for x or z above $17$ bits. – ccorn Dec 12 '13 at 18:09
  • Thanks! I fixed the JacSym function and now I get the same results as those you posted. I forgot that Matlab does not realize that $7 \equiv -1 \mod 8$ which caused some problems. Also, this is just a simple test for "small $n$'s", so I have not optimized the code at all. I will definitively take your advices for the modexp into consideration if I ever am to make it better. Again, thanks for your time! – KryptoJon Dec 12 '13 at 18:54

0 Answers0