\documentstyle[12pt,kol,emlines2,twoside]{iarticle}

\begin{document}


\abovedisplayskip=4pt plus 2pt minus 2pt 
\belowdisplayskip=4pt plus 2pt minus 3pt

\font\smr=plr7 scaled\magstep0
\font\hsmr=plr9 scaled\magstep0
\font\srm=plr7 scaled\magstep0
\font\sc=plcsc10 scaled\magstep1
\font\ssrm=plr5 scaled\magstep0
\font\lasy=lasy10 scaled\magstep1
\def \Box {{\lasy 2}}
\def\aaa{\LARGE\bf}



{\centerline{\Large \bf
PSPrimeQ - version 1, 10.07.2009
}}
{\centerline{\Large \bf
by Jaros/law Wr/oblewski
}}
{\centerline{\tt jwr@math.uni.wroc.pl}}
{\centerline{\tt http:////www.math.uni.wroc.pl//$\sim$jwr//AP26//GPU}}



\bigskip

This document is an attempt to produce pseudoprimality testing algorithm suitable for implementing on GPU, i.e. using only 32-bit operations.

\bigskip


$p$ - unsigned integer with slightly less than 64 significant bits. This number is being tested for base 2 pseudoprimality:

{\Large $$2^p\equiv2({\hbox{mod}}\;p)$$}


\bigskip

$z$ - number of leading zeros of $p$ when represented as 64-bit unsigned int. Hence $p$ has $64-z$ significant bits:

{\Large $$2^{64-z-1}\le p<2^{64-z}$$}

\bigskip

$R$ - reciprocal of $p$ with 128 significant bits, rounded down:

{\Large $$R=\left[{2^{191-z}\over p}\right]$$}

$R$ is represented as 128-bit unsigned int.

$r=1//p$ can be computed using recursive formula:

$$r={1-e\over p}$$

$$r_{new}=2r-p r^2$$

Then

$$r_{new}={1-e^2\over p}$$

Note that the above recursion assumes $r$ to be a real number close to $1//p$.

\bigskip

$A$ - a number with at most $127-z$ significant bits. Our short term goal is to compute $A({\hbox{mod}}\;p)$.

$As$ - number $A$ shifted left by $z+1$ bits.

\bigskip


$f$ - 64-bit product of two 128-bit integers $As$ and $R$. We do not care about 128 highest bits of the 256-bit product, so in principle we could multiply those as 128-bit integers, where the overflow bits are being discarded. Also we do not care about the 64 lowest bit, but we do care about some precision of the lowest bit we keep. So perhaps one could compute 80 bits (ignore computations on the lowest 48 bits) and discard the 16 lowest bits from the result.

\bigskip

$out$ - 64-bit product of two 64-bit integers $f$ and $p$, where we want to keep only the 64 highest bits. It should be enought to perform computations on the highest 80 bits and then discard 16 lowest bits from the result.

We now have

{\Large $$A({\hbox{mod}}\;p) = out+1$$}

provided errors aren't too large - this has to be traced carefully.

Also we assume that $A$ is not divisible by $p$, but that is always the case if $A$ arises from pseudoprimality testing and $p$ is not a power of 2.

\end{document}
