% This is a LaTeX input file for p146.ps and p146.pdf

% To produce p146.ps and p146.pdf:
%   latex p146.tex && dvips -o p146.ps p146.dvi && ps2pdf p146.ps p146.pdf

% document class options
%   <drivers>              I always use only dvips
%   10pt|11pt|12pt         Choose the standard font size (default: 10pt)
%   onecolum|twocolumn     Number of column (default: onecolumn)
%   oneside|twoside        The output is printed on one side or on two side.
%                          (defaults: oneside for article, letter and report,
%                          twoside for book)
%   notitlepage|titlepage  Choose if the document has a title page.
%                          (default: titlepage for book and report, notitlepage
%                          for the others)
%   final|draft            final or draft (default: final)
%   leqno                  left equation number
%                          Number of the equation is left instead of rigth.
%   openbib                not standar form for the bibliography
%   <paper_format>         a4paper        (297 x 210 mm)
%                          a5paper        (210 x 148 mm)
%                          b5paper        (250 x 176 mm)
%                          letterpaper    (11 x 8.5 in)
%                          legalpaper     (14 x 8.5 in)
%                          executivepaper (10.5 x 7.25 in)
%   openright|openany      with openright a new chapter is always on a 
%                          right page.
%
\documentclass[dvips,10pt,onecolumn,twoside,a4paper]{article}

% \overfullrule=4pt % Comment this line (or use 0pt) to get rid of the overfull box rules.

%%------------------------------------------------------------------------------
%% Enable latin1 accents.
%%------------------------------------------------------------------------------
\usepackage[latin1]{inputenc}
\usepackage[T1]{fontenc}

% Load the fancyvrb package.
%   You can use the \DefineVerbatimEnvironment to create custom verbatim
%   environment. See lextensions for an ex.: the 'scrivimi' verbatim.
\usepackage{fancyvrb}

%%------------------------------------------------------------------------------
%% Document specific macros.
%%------------------------------------------------------------------------------
%
% Define VerbatimFORTH based on Verbatim (to be used for FORTH verbatim).
%\DefineVerbatimEnvironment{VerbatimFORTH}{Verbatim}{
%  xleftmargin=\parindent,   % indentazione di \parindent a sinistra
%  xrightmargin=\parindent,  % indentazione di \parindent a destra
% commandchars=\\\{\},      % escape and command charachters
% commentchar=\%,           % comment character
%  fontfamily=tt,            % courier-like font
%  fontsize=\footnotesize    % ..footnotesize
% frame=lines,              % there is a top line and a bottom line
% framesep=1em,             % ..separated by 1 em
% label={\bfseries\ Label}, % and the label
% labelposition=topline     % is positioned at the center of the topline
%}
%
% \url{<url>}
% ex. \url{http://www.gnu.org}
%
%   TeXbook pag. 93: 'Sometimes you want to permit a linebreak after a '/' just
%   as if it werea hypen. For this purpose plain TeX allows you to say '\slash';
%   for example, '\input\slash output' produces 'input/output' with an optional break.
%
%   The text of a URL cannot be hypenated so I suggest TeX a linebreak before it.
%\newcommand{\url}[1]{\linebreak[3]\texttt{\textsl{\small#1}}}   % \linebreak[0=consiglio..4=ordine}
%
\begin{document}%                     % End of preamble and beginning of text.
%
\title{Projekt Euler — Problem p146}
\author{Luca Masini
\thanks{\protect\email{luca.masini@member.fsf.org}}}
\date{\today}
\maketitle
\label{p146}

%
\begin{multicols}{2}
\section{Einleitung zum Projekt Euler}
%
Die Projekt--Euler\footnote{\url{http://projecteuler.net/index.php?section=about}}--Internet--Seite bietet eine Reihe von mathematischen
Programmierproblemen an.

Neben der Mathematik selbst, die hilft, elegante und effi\-ziente
Lösungsmethoden zu finden, sind für die vorgestellten Probleme ein
Rechner und Programmierfähigkeiten wesentliche Zutaten, um überhaupt zu
einer Lösung zu kommen.

Wer ein Interesse an der Mathematik hat, kann hier sehr interessante Probleme finden.

Die Probleme haben verschiedene Schwierigkeitsgrade, aber mit steigender
Erfahrung wird man lernen, auch die schwierigen Probleme zu lösen.

Am Ende ist es auch eine gute Gelegenheit, um ein bisschen Erfahrung mit
Forth zu sammeln!
%
\section{Ein Beispiel: Problem 146}
%
Als Beispiel nehme ich das Problem 146, hier zu finden:
\href{http://projecteuler.net/index.php?section=problems\&id=146}{\tt
  http://projecteuler.net/index.php?section=prob lems\&id=146} und mit
folgendem Satz beschrieben:

\begin{quote}
The smallest positive integer $n$ for which the numbers $n^2+1$, $n^2+3$, $n^2+7$, $n^2+9$, $n^2+13$, and $n^2+27$ are consecutive primes is 10.
The sum of all such integers $n$ below 1 million is 1242490.
What is the sum of all such integers $n$ below 150 million?
\end{quote}

Hier gilt es also, sechs Primzahlen zu erkennen, von denen die ersten vier Primzahlen sog.\ \emph{Primzahlvierlinge} sein sollen.

\emph{Primzahlvierlinge} sind vier Primzahlen der Form $(p, p+2, p+6, p+8)$.
Wenn also $p$ die kleinste Primzahl von Primzahlvierlingen ist, dann soll man jeden einzelnen Primzahlvierling so durch die natürliche Zahl $n$ beschreiben können:
$$
\begin{array}{r@{\quad=\quad}l}
n^2 + 1 & p     \\
n^2 + 3 & p + 2 \\
n^2 + 7 & p + 6 \\
n^2 + 9 & p + 8
\end{array}
$$

Eine kurze Suche im Internet findet auf der Seite
\href{http://mathworld.wolfram.com/PrimeQuadruplet.html}{\tt
  http://mathworld.wolfram.com/PrimeQuadruplet. html} eine 
Erklärung mit weiteren Eigenschaften von Primzahlvierlingen.

Insbesondere gilt folgendes:

\begin{quote}
With the exception of (5, 7, 11, 13), a prime quadruple must be of the form ( 30 m + 11, 30 m + 13, 30 m + 17, 30 m + 19 )
\end{quote}

Aus $$p = n^2 + 1 = ( 30 m + 10 ) + 1$$ erhalten wir $$n^2 = ( 30 m +
10 ) = 10 ( 3 m + 1 )$$

Deshalb ist $n$ immer ein Vielfaches von 10, aber nie ein Vielfaches von
3.

Die Anzahl der Tests wird damit von 150 Millionen auf 10 Millionen
verringert, aber es sind immer noch für jedes $n$ teure Primzahltests zu machen.

Mit einigen mathematischen Schritten, die relativ lange und
langweilig sind, ist es auch möglich zu zeigen, dass
$$( n\ \mbox{mod}\ 7 = 3 )\ or\ ( n\ \mbox{mod}\ 7 = 4 )$$ ist.
Diese Bedingung erlaubt, die Zahl der Tests weiter zu vermindern und
wird im Code verwendet.
%
\subsection{Algorithmus}
%
Die äußere Schleife des Algorithmus erhöht den Index um 10. Außerdem
kann sie alle durch drei teilbaren Zahlen weglassen, und $n\ {\tt mod}\ 7$
muss gleich $3$ oder $4$ sein.

Die Indizes, die diese Bedingungen erfüllen, werden an
\texttt{(main-ba)} übergeben, das eine Reihe von Tests über die
Restwerte ausführt. Dann macht es einen Primzahltest mit
\texttt{miller\_rabin} und gibt schließlich TOS zurück, wenn der
Parameter die Bedingungen des Problems erfüllt, sonst 0.

Dieser Wert wird verwendet, um die lokale Variable \texttt{sum} zu
vergrößern, die am Ende der Schleife die Lösung des Problems
enthalten wird.

\begin{verbatim}
: main { W: limit -- d }
    0. { D: sum }
    limit 1+ 10 do
        \ i must be divisible by 10,
        \ loop only through multiples of 10
        i 3 mod 0<> if
            \ i can't be divisible by 3
            i 4 + 7 mod 1 <= if
                \ i % 7 must be either 3 or 4
                i (main-ba) s>d sum d+  to sum
            then
        then
    10 +loop
    ." sum=" sum ud. cr
;
\end{verbatim}
%
\subsection{Der Miller--Rabin--Primzahltest}
%
Der Primzahltest ist durch das Wort \texttt{miller\_rabin} (siehe
\href{http://en.wikipedia.org/wiki/Miller-Rabin\_primality\_test}{\tt
  http://en.wikipedia.org/wiki/Miller-Rabin\_ primality\_test})
realisiert.

Es handelt sich um einen probabilistischen Primzahltest.  Erhält er
eine natürliche Zahl $n$ als Eingabe, so gibt er aus, ob diese Zahl
eine Primzahl ist oder nicht.  Gibt er dabei \emph{FALSE ($n$ ist
  keine Primzahl)} aus, so ist das immer richtig.  Gibt er allerdings
\emph{TRUE ($n$ ist wahrscheinlich eine Primzahl)} aus, so ist dies
mit geringer Wahrscheinlichkeit falsch.

Für die Implementierung dieses Tests braucht man einen
Zufallszahlengenerator. Hier wurde einer der \emph{Scientific Forth
  Library\footnote{\url{http://www.taygeta.com/fsl/sciforth.html}}}
benutzt.

%
\subsection{Logik-- und Shift--Wörter für Doppelwortoperanden}
%
Der Zufallszahlengenerator \texttt{ran4} und der Primzahltest
\texttt{miller\_rabin} brauchen logische Worte, die mit Doppelzellen
arbeiten.  Diese Worte wurden wie folgt implementiert:

{\small
\begin{verbatim}
: dand    ( d d -- d    ) rot and >r and r>       ;
: dinvert (   d -- d    ) swap invert swap invert ;
: dlshift ( d u -- d    ) 0 ?do d2* loop          ;
: dor     ( d d -- d    ) rot or  >r or  r>       ;
: dnot    (   d -- 0.|1.) d0= abs s>d             ;
: drshift ( d u -- d    ) 0 ?do d2/ loop dabs     ;
: dxor    ( d d -- d    ) rot xor >r xor r>       ;
\end{verbatim}
}
%
\subsection{Arithmetische und modular--arithmetische Wörter für Doppelwortoperanden}
%
Der Teil, der den größten Aufwand verursachte, war die Implementierung
der modulo--arithmetischen Worte für Doppelzellen.  Folgende Worte sind im Quellcode zu finden:

\begin{verbatim}
  d**mod  { D: base D: power D: m -- d  }
  \ Modular exponentiation
  d*mod   { D: a    D: b     D: m -- d  }
  \ Modular multiplication
  d+mod   { D: a    D: b     D: m -- d  }
  \ Modular addition
  d*      ( d1 d2                 -- d  )
  \ Multiplication for double cells
  d/      { D: u    D: v          -- d  }
  \ Division for double cells
\end{verbatim}

Einige Worte wurden auch aus \emph{Long Divisors and Short Fractions (Nathaniel Grossman\footnote{http://www.forth.org/fd/FD-V06N3.pdf})} übernommen.

Versuchsweise wurde das Wort \texttt{d/} auch mit \emph{float} implementiert, und
war dann auch deutlich schneller, aber das hatte eine nicht korrekte
Lösung zur Folge:

\begin{Verbatim}
Found n=10
Found n=315410
Found n=927070
Found n=2525870
Found n=8146100
Found n=16755190
Found n=39313460
Found n=97387280
Found n=119571820
Found n=121288430
Found n=130116970
Found n=139985660
Found n=144774340    <---- WRONG
sum=821107610
\end{Verbatim}
%
\section{Zusammenfassung}
%
Der präsentierte Code löste das Problem in 6 Minuten und einer Sekunde
auf meinem Laptop (Intel Core 2 Duo Prozessor T7200, Merom, 2.00 GHz,
667 MHz FSB, 4). Ich bin immer noch nicht glücklich über die
Performance, und denke, es sollte möglich sein, unter einer Minute bleiben
zu können.
%
\end{multicols}
\section{Source Code}
%
\begin{quote}
{\small\listinginput[1]{1}{2008-01/p146.fs}}
\end{quote}
%
% \end{document}                 % End of document.
