Algorytmy i Struktury Danych

SPIS TREŚCI

Podrozdziały

Rozwiązanie 1
Rozwiązanie 2
Rozwiązanie 3

Problem

Znaleźć rozkład liczby \({p>1}\) na iloczyn czynników pierwszych.

Pierre de Fermat podał prosty sposób znajdowania czynników (liczb dzielących \(p\) bez reszty) liczby nieparzystej \(p\). Opiera się on na spostrzeżeniu, iż jeśli potrafimy znaleźć dwie liczby naturalne \(x\) i \(y\), takie że:

\({p=x^2-y^2}\)

\({p=(x+y)\cdot(x-y)}\)

zatem czynnikami liczby \(p\) są:

\[
\begin{aligned}
m &= x + y \\
n &= x - y
\end{aligned}
\]

Znalezione czynniki \(m\) i \(n\) nie muszą być liczbami pierwszymi, zatem metodę Fermata stosujemy również do ich rozkładu. Jest to możliwe, ponieważ czynniki liczby nieparzystej są również nieparzyste. Czynniki o wartości \(2\) można wyeliminować z \(p\) przed zastosowaniem metohy Fermata, zatem nie jest to żadne ograniczenie.

Metoda Fermata jest szczególnie efektywna w przypadku czynników leżących w pobliżu pierwiastka z \(p\). W przeciwnym razie jej efektywność jest taka sama lub nawet gorsza jak dla metohy próbnych dzieleń. W praktyce metoda Fermata nie jest używana – podajemy ją ze względów hydaktycznych oraz z powodu jej wagi dla innych, bardziej zaawansowanych metod poszukiwań rozkładu liczby na czynniki pierwsze.

Algorytm wykorzystuje bezpośrednio własność \({p=x^2-y^2}\) do znalezienia liczb \(x\) i \(y\). Poszukiwania rozpoczynamy od \(x\) równego pierwiastkowi z \(p\) zaokrąglonemu w górę do najbliższej liczby całkowitej. Obliczamy:

\({y^2=x^2-p}\)

Następnie sprawdzamy, czy \(y\) jest liczbą całkowitą. Jeśli tak, to znaleźliśmy odpowiednie liczby \(x\) i \(y\). Ponieważ liczba \(p\) z założenia jest nieparzysta, to wszystkie jej dzielniki są nieparzyste. Obliczamy dzielniki \(m\) i \(n\) i jeśli są różne od \(1\), to wywołujemy rekurencyjnie procedurę Fermata do znalezienia ich rozkładu. Jeśli jeden z czynników jest równy \(1\), to drugi musi być równy \(p\) i \(p\) jest liczbą pierwszą – wypisujemy \(p\) i kończymy.

Ponieważ czynnik musi być mniejszy lub równy \(p\), to:

\({m=x+y\le p}\)

Przy wzroście \(x\) rośnie również \(y\) lub pozostaje takie samo. Zatem gdy suma \({x + y}\) przekroczy \(p\), to nie znajdziemy już żadnego dzielnika liczby \(p\) i algorytm można zakończyć przyjmując, iż \(p\) jest pierwsze.

Algorytm rozkładu liczby naturalnej na czynniki pierwsze metodą Fermata

Wejście:

\(p\) : liczba rozkładana na czynniki pierwsze; \({p\in\mathbb{N}}\), \(p\) jest nieparzyste.

Wyjście:

Czynniki pierwsze liczby \(p\).

Elementy pomocnicze:

\(x,y,z\) : zmienne do rozkładu \(p\); \({x,y,z\in\mathbb{N}}\).
\(m,n\) : czynniki liczby \(p\); \({m,n\in\mathbb{N}}\).

Lista kroków dla procedury rekurencyjnej Fermat(p)

K01:	\(x\leftarrow\lceil\sqrt{p}\rceil\) ; obliczamy wartość początkową dla \(x\)
K02:	\(z\leftarrow x^2-p\) ; \(z = y^2\)
K03:	\(y\leftarrow\lceil\sqrt{z}\rceil\) ; sprawdzamy, czy \(z\) jest kwadratem liczby naturalnej
K04:	Jeśli \(z\ne y^2,\) to idź do kroku K11
K05:	\(m\leftarrow x+y\) ; obliczamy czynniki
K06:	\(n\leftarrow x-y\)
K07:	Jeśli \(n=1,\) ; przerywamy przy \(n=1\), gdyż \(p\) jest pierwsze to idź do kroku K13
K08:	\(\text{Fermat}(m)\) ; rozkładamy rekurencyjnie czynniki \(m\) i \(n\)
K09:	\(\text{Fermat}(n)\)
K10:	Zakończ
K11:	\(x\leftarrow x+1\) ; następne \(x\)
K12:	Jeśli \(x+y<p,\) ; kontynuujemy pętlę to idź do kroku K02
K13:	Pisz \(p\) ; \(p\) jest pierwsze
K14:	Zakończ

Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Program w pierwszym wierszu czyta dowolną liczbę naturalną \({p > 1}\) i w następnym wierszu wypisuje kolejne czynniki pierwsze tej liczby. Czynniki mogą nie być wypisywane w kolejności rosnącej. Zwróć uwagę, iż do poprawnego działania zmienna \(x\) musi być zadeklarowana jako \({64{-}\text{bitowa}}\) – w przeciwnym razie dla dużych liczb jej kwadrat wyjdzie poza zakres liczb \({32{-}\text{bitowych}}\). W efekcie program może dokonywać rozkładu na czynniki pierwsze liczb \({32{-}\text{bitowych}}\).

Pascal

// Metoda Fermata
// Data: 27.03.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

program prg;

uses math;

procedure Fermat(p : longword);

var x, y, z, m, n : qword;

begin
  x := ceil(sqrt(p));
  repeat
    z := x*x-p;
    y := floor(sqrt(z));
    if z = y*y then
    begin
      m := x+y;
      n := x-y;
      if n = 1 then break;
      Fermat(m);
      Fermat(n);
      exit;
    end;
    inc(x);
  until x+y >= p;
  write(p, ' ');
end;

var p : longword;

begin
  readln(p);
  while p mod 2 = 0 do
  begin
    p := p div 2;
    write(2, ' ');
  end;
  if p > 1 then
    Fermat(p);
  writeln;
end.

C++

// Metoda Fermata
// Data: 27.03.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

#include <iostream>
#include <cmath>

using namespace std;

void Fermat(unsigned int p)
{
  unsigned long long x, y, z, m, n;

  x = (unsigned long long)ceil(sqrt(p));
  do
  {
    z = x*x-p;
    y = (unsigned long long)floor(sqrt(z));
    if(z == y*y)
    {
      m = x+y;
      n = x-y;
      if(n == 1) break;
      Fermat(m);
      Fermat(n);
      return;
    }
    x++;
  } while(x+y < p);
  cout << p << " ";
}

int main()
{
  unsigned int p;

  cin >> p;
  while(p % 2 == 0)
  {
    p /= 2;
    cout << 2 << " ";
  }
  if(p > 1) Fermat(p);
  cout << endl;
  return 0;
}

Basic

' Metoda Fermata
' Data: 27.03.2008
' (C)2020 mgr Jerzy Wałaszek
'----------------------------

Sub Fermat(Byval p As Uinteger)

Dim As Ulongint x, y, z, m, n

  x = Ceil(Sqr(p))
  Do
    z = x*x-p
    y = Floor(Sqr(z))
    If z = y*y Then
      m = x+y
      n = x-y
      If n = 1 Then Exit Do
      Fermat(m)
      Fermat(n)
      Exit Sub
    End If
    x += 1
  Loop Until x+y >= p
  Print p;" ";
End Sub

Function Ceil(Byval x As Double) _
                      As Ulongint
   Dim result As Ulongint
   result = Culngint(x)
   If result < x Then _
     result += 1
   Ceil = result
End Function

Function Floor(Byval x As Double) _
                       As Ulongint
   Dim result As Ulongint
   result = Culngint(x)
   If result > x Then _
     result -= 1
   Floor = result
End Function

Dim p As Uinteger

Input p
while p mod 2 = 0
  p = p shr 1
  print 2;" ";
wend
if p > 1 then _
  Fermat(p)
Print
End

Python (dodatek)

# Metoda Fermata
# Data: 12.02.2024
# (C)2024 mgr Jerzy Wałaszek
# --------------------------

import math


def fermat(p):
    x = math.ceil(math.sqrt(p))
    while 1:
        z = x*x-p
        y = math.floor(math.sqrt(z))
        if z == y*y:
            m = x+y
            n = x-y
            if n == 1: break
            fermat(m)
            fermat(n)
            return
        x += 1
        if x+y >= p: break
    print(p, end=" ")
    return


p = int(input())
while not (p % 2):
    p >>= 1
    print(2, end=" ")
if p > 1: fermat(p)
print()

Wynik:

855855
11 7 13 3 3 5 19

do podrozdziału do strony

Podstawową wadą pierwszego rozwiązania jest konieczność wyznaczania wartości pierwiastka kwadratowego. Niestety, jest to operacja dosyć czasochłonna. Zastanówmy się zatem, czy algorytmu Fermata nie można zapisać w inny sposób, tak aby nie było konieczności wyznaczania pierwiastków kwadratowych w pętli poszukującej dzielników. Dokonajmy najpierw kilku prostych spostrzeżeń.

Jeśli mamy daną wartość kwadratu liczby naturalnej, np. \(x^2\), to kwadrat następnej liczby \({x + 1}\) można wyznaczyć następująco:

\({(x+1)^2=}\)
\({x^2+2\cdot x+1=}\)
\({x^2+hx;}\)
gdzie: \({hx=2\cdot x+1}\)

Ponieważ \(x\) wzrosło o \(1\), to nowy przyrost \(hx'\) wyniesie:

\({hx'=2\cdot(x+1)+1=}\)
\({2\cdot x+2+1=}\)
\({(2\cdot x+1)+2=}\)
\({hx+2}\)

Wynika stąd, iż mając wartość początkową kwadratu \(x^2\) oraz \(hx\), kolejne kwadraty obliczamy w pętli wg reguły:

\[
\begin{aligned}
x'^2&=x^2+hx\\
x&=\dfrac{hx-1}{2}\\
hx'&=hx+2
\end{aligned}
\]

Poniższy program demonstruje tę zasadę wyznaczania kwadratów kolejnych liczb naturalnych od \(0\) do \(15\).

\(kx\) : kwadrat \(x\), początkowo \(0\)
\(hx\) : przyrost kwadratu, początkowo \({hx = 1}\), gdyż \({x = 0}\)
\(x\) : kolejna liczba od \(0\) do \(15\). Liczbę tę wyznaczamy bezpośrednio z przyrostu \(hx\).

Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Pascal

// Wyznaczanie kwadratów
// kolejnych liczb naturalnych
// Data: 2.04.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------=------------

program prg;

var kx, hx, x : integer;

begin
  kx := 0; hx := 1;
  repeat
    x := (hx-1) shr 1;
    writeln(x:5, kx:5);
    inc(kx, hx); inc(hx, 2);
  until x = 15;
  writeln;
end.

C++

// Wyznaczanie kwadratów
// kolejnych liczb naturalnych
// Data: 2.04.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

#include <iostream>
#include <iomanip>

using namespace std;

int main()
{
  int kx, hx, x;

  kx = 0; hx = 1;
  do
  {
    x = (hx-1) >> 1;
    cout << setw(5) << x
         << setw(5) << kx
         << endl;
    kx += hx; hx += 2;
  } while(x < 15);
  cout << endl;
  return 0;
}

Basic

' Wyznaczanie kwadratów
' kolejnych liczb naturalnych
' Data: 2.04.2008
' (C)2020 mgr Jerzy Wałaszek
'----------------------------

Dim As Integer kx, hx, x

kx = 0: hx = 1
Do
  x = (hx-1) Shr 1
  Print Using "#### ####";x;kx
  kx += hx: hx += 2
Loop Until x = 15
Print
End

Python (dodatek)

# Wyznaczanie kwadratów
# kolejnych liczb naturalnych
# Data: 12.02.2024
# (C)2024 mgr Jerzy Wałaszek
#---------------------------

kx, hx = 0, 1
while 1:
    x = (hx-1) >> 1
    print("%4d %4d" % (x, kx))
    kx += hx
    hx += 2
    if x == 15: break
print()

Wynik:

Zastosujmy następującą strategię. Mamy daną liczbę naturalną \(p\), \({p > 1}\) i \(p\) jest nieparzyste. Obliczamy:

\[
\begin{aligned}
x&=\lceil\sqrt{p}\rceil\\
y&=0\\
e&=x^2-y^2-p
\end{aligned}
\]

Liczba \(e\) jest wartością błędu pomiędzy \({x^2 - y^2}\) a liczbą \(p\). Wyznaczamy przyrosty dla \(x^2\) oraz dla \(y^2\):

\[
\begin{aligned}
hx&=2\cdot x+1\\
hy&=2\cdot y+1=1,
\end{aligned}
\]
gdyż początkowo \({y=0}\)

Teraz w pętli badamy błąd \(e\). Jeśli jest równy \(0\), to

\[
\begin{aligned}
p&=x^2-y^2\\
x&=\dfrac{hx-1}{2}\\
y&=\dfrac{hy-1}{2}\\
m&=x+y=\dfrac{hx-1}{2}+\dfrac{hy-1}{2}=\dfrac{hx+hy-2}{2}=\dfrac{hx+hy}{2}-1\\
n&=x-y=\dfrac{hx-1}{2}-\dfrac{hy-1}{2}=\dfrac{hx-hy}{2}
\end{aligned}
\]

Jeśli \(n\) jest różne od \(1\), to \(m\) i \(n\) rozkładamy dalej na czynniki tym samym algorytmem. W przeciwnym razie liczba \(p\) jest pierwsza, zwracamy ją i kończymy.

Jeśli \(e\) jest różne od zera, to nie znaleźliśmy jeszcze \(x\) i \(y\) spełniających równanie Fermata. Badamy znak \(e\).

Jeśli \({e>0}\), to przeważa \({x^2}\), zatem zwiększamy \({y^2}\) i odejmujemy od \(e\) przyrost \(hy\). Po tej operacji \(hy\) jest zwiększane o \(2\), aby w kolejnym obiegu otrzymać kwadrat następnego \(y\) – a właściwie jego przyrost w stosunku do poprzedniego kwadratu. Zwiększanie \(y\) kontynuujemy do momentu, aż \({e\le0}\), czyli aż osiągniemy równowagę lub przeważy \({y^2}\).

Jeśli \({e<0}\), to przeważa \({y^2}\), zatem zwiększamy \({x^2}\) i dodajemy do \(e\) przyrost \(hx\), który następnie zwiększamy o \(2\). Operację kontynuujemy do momentu, aż \({e\ge0}\).

Wracamy na początek pętli, aż do uzyskania \({e=0}\). W tym miejscu algorytm można zoptymalizować, sprawdzając, czy \({x+y\ge p}\), czyli \({\frac{hx+hy}{2}>p}\). Jeśli tak, to liczba \(p\) jest pierwsza i nie znajdziemy rozkładu na inne czynniki poza \(p\) i \(1\). Zatem przerywamy pętlę zwracając \(p\).

Zwróć uwagę, iż w tej wersji algorytmu Fermata nie wykonujemy w pętli operacji pierwiastkowania, a jedynie proste dodawania. Dzięki temu algorytm staje się dużo szybszy od algorytmu podanego w rozwiązaniu 1.

Algorytm rozkładu liczby naturalnej na czynniki pierwsze metodą Fermata

Wejście:

\(p\) : liczba rozkładana na czynniki pierwsze, \({p\in\mathbb{N}}\), \(p\) jest nieparzyste.

Wyjście:

Czynniki pierwsze liczby \(p\).

Elementy pomocnicze:

\(e\) : wartość błędu; \({e\in\mathbb{Z}}\).
\(m,n\) : czynniki \(p\); \({m,n\in\mathbb{N}}\).
\(hx,hy\) : przyrosty kwadratów \(x\) i \(y\); \({hx,hy\in\mathbb{N}}\).

Lista kroków dla procedury rekurencyjnej Fermat(p)

K01:	\(m\leftarrow\lceil\sqrt{p}\rceil\) ; wyznaczamy wartość błędu \(e\) ; oraz przyrosty \(hx\) i \(hy\)
K02:	\(e\leftarrow m^2-p\) ; \(e = x^2-p,\) ponieważ początkowo \(y = 0\)
K03:	\(hx\leftarrow 2\cdot m+1\) ; \(hx\) – przyrost \(x^2\)
K04:	\(hy\leftarrow1\) ; \(hy\) – przyrost \(y^2\)
K05:	Dopóki \(e\ne0:\) ; w pętli staramy się zrównać wykonuj kroki K06…K15 ; \(x^2-y^2\) z \(p\) – wtedy \(e=0\)
K06:	Dopóki \(e>0:\) ; przeważa \(x^2\), zwiększamy \(y^2\) wykonuj kroki K07…K08
K07:	\(e\leftarrow e-hy\) ; modyfikujemy błąd o przyrost \(y^2\)
K08:	\(hy\leftarrow hy+2\) ; przyrost \(y^2\) dla zwiększonego o \(1\) \(y\)
K09:	Dopóki \(e<0:\) ; przeważa \(y^2\), zwiększamy \(x^2\) wykonuj kroki K10…K11
K10:	\(e\leftarrow e+hx\) ; modyfikujemy błąd o przyrost \(x^2\)
K11:	\(hx\leftarrow hx+2\) ; przyrost \(x^2\) dla zwiększonego \(x\)
K12:	Jeśli \((hx+hy)\gg1\le p,\) ; sprawdzamy drugi warunek to następny obieg pętli K05 ; zakończenia pętli
K13:	\(hx\leftarrow2\) ; jeśli jest spełniony, ; to \(p\) jest liczbą pierwszą
K14:	\(hy\leftarrow0\) ; \(hx\) i \(hy\) ustawiamy tak, ; aby otrzymać \(n=1\)
K15:	Idź do kroku K16 ; przerywamy pętlę K05!
K16:	\(m\leftarrow((hx+hy)\gg1)-1\) ; czynnik większy
K17:	\(n\leftarrow(hx-hy)\gg1\) ; czynnik mniejszy
K18:	Jeśli \(n=1,\) ; sprawdzamy, czy \(p\) jest pierwsze to idź do kroku K21
K19:	\(\text{Fermat}(m)\) ; jeśli nie, to czynniki \(\text{Fermat}(n)\) ; \(m\) i \(n\) rozkładamy dalej
K20:	Zakończ
K21:	Pisz \(p\) ; \(p\) jest pierwsze. ; Wyprowadzamy je i kończymy
K22:	Zakończ

Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Program w pierwszym wierszu czyta liczbę \(p\) i w następnym wierszu wypisuje kolejne czynniki pierwsze tej liczby. Liczba \(p\) nie musi być nieparzysta. Przed przekazaniem jej do procedury Fermat program usuwa z \(p\) wszystkie czynniki równe \(2\).

Pascal

// Metoda Fermata
// Data: 2.04.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

uses math;

procedure Fermat(p : int64);

var e, hx, hy, m, n : int64;

begin
  m  := ceil(sqrt(p));
  hx := (m shl 1)+1;
  e  := m*m-p;
  hy := 1;
  while e <> 0 do
  begin
    while e > 0 do
    begin
      dec(e, hy); inc(hy, 2);
    end;
    while e < 0 do
    begin
      inc(e, hx); inc(hx, 2);
    end;
    if(hx+hy) shr 1 > p then
    begin
      hx := 2; hy := 0; break;
    end;
  end;
  m := ((hx+hy) shr 1)-1;
  n := (hx-hy) shr 1;
  if n > 1 then
  begin
    Fermat(m); Fermat(n);
  end
  else write(p, ' ');
end;

var p : int64;

begin
  readln(p);
  while p mod 2 = 0 do
  begin
    write(2, ' ');
    p := p shr 1;
  end;
  if p > 1 then
    Fermat(p);
  writeln;
end.

C++

// Metoda Fermata
// Data: 2.04.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

#include <cstdlib>
#include <iostream>
#include <cmath>

using namespace std;

void Fermat(long long p)
{
  long long e, hx, hy, m, n;
  m  = (long long)sqrt(p);
  hx = (m << 1)+1;
  e  = m*m-p;
  hy = 1;
  while(e != 0)
  {
    while(e > 0)
    {
      e -= hy; hy += 2;
    }
    while(e < 0)
    {
      e += hx; hx += 2;
    }
    if(((hx+hy) >> 1) > p)
    {
      hx = 2;
      hy = 0;
      break;
    }
  }
  m = ((hx+hy) >> 1)-1;
  n = (hx-hy) >> 1;
  if(n != 1)
  {
    Fermat(m); Fermat(n);
  }
  else cout << p << " ";
}

int main()
{
  long long p;
  cin >> p;
  while(p % 2 == 0)
  {
    cout << 2 << " ";
    p >>= 1;
 }
 if(p > 1) Fermat(p);
 cout << endl;
 return 0;
}

Basic

' Metoda Fermata
' Data: 2.04.2008
' (C)2020 mgr Jerzy Wałaszek
'---------------------------

Sub Fermat(Byval p As Longint)

  Dim As Longint e, hx, hy, m, n

  m  = Clngint(Sqr(p))
  hx = (m Shl 1)+1
  e  = m*m-p
  hy = 1
  While e
    While e > 0
      e -= hy: hy += 2
    Wend
    While e < 0
      e += hx: hx += 2
    Wend
    If (hx+hy) Shr 1 > p Then
      hx = 2
      hy = 0
      Exit While
    End If
  Wend
  m = ((hx+hy) Shr 1)-1
  n = (hx-hy) Shr 1
  If n > 1 Then
    Fermat(m): Fermat(n)
  Else
    Print p;" ";
  End If
End Sub

Dim p As Longint
 
Input p
While p Mod 2 = 0
  Print 2;" ";
  p = p Shr 1
Wend
If p > 1 Then _
  Fermat(p)
Print
End

Python (dodatek)

# Metoda Fermata
# Data: 12.02.2024
# (C)2024 mgr Jerzy Wałaszek
# --------------------------

import math


def fermat(p):
    m = int(math.sqrt(p))
    hx = (m << 1)+1
    e = m*m-p
    hy = 1
    while e:
        while e > 0:
            e -= hy
            hy += 2
        while e < 0:
            e += hx
            hx += 2
        if ((hx+hy) >> 1) > p:
            hx, hy = 2, 0
            break
    m = ((hx+hy) >> 1)-1
    n = (hx-hy) >> 1
    if n != 1:
        fermat(m)
        fermat(n)
    else:
        print(p, end=" ")
    return


p = int(input())
while not (p % 2):
    print(2, end=" ")
    p >>= 1
if p > 1:
    fermat(p)
print()

Wynik:

254821743888
2 2 2 2 230816797 23 3

JavaScript

<html>
  <head>
    <title>
      Metoda Fermata
    </title>
  </head>
  <bohy>

<div style="overflow-x: auto; overflow-y: hidden;"
     align="center">
  <table
  border="0"
  cellpadding="4"
  style="border-collapse:
         collapse">
    <tr>
      <td nowrap>
        <form
        name="frm"
        style="text-align: center;
               background-color:
               #E7E7DA">
          <b>&nbsp;&nbsp;
          Rozkład na czynniki
          pierwsze metodą Fermata
          &nbsp;&nbsp;</b><br/>
	      (C)2020 mgr
	      Jerzy Wałaszek
	      &nbsp;&nbsp;
	      <hr>
          p = <input
              type="text"
              name="inp_p"
              size="16"
              value="254821743888"
              style="text-align:
                     right">
          <hr>
          <input
          type="button"
          value="Wykonaj"
          name="B1"
          onclick="main()">
          <hr>
          <b>Czynniki:</b>
          <div id="out">.</div>
        </form>
      </td>
    </tr>
  </table>
</div>

<script type=text/javascript
        language=javascript>

// Metoda Fermata
// Data   : 2.04.2008
// (C)2008 mgr Jerzy Wałaszek
//---------------------------

var t,c;

function Fermat(p)
{
  var e,hx,hy,m,n;
  m  = Math.floor(Math.sqrt(p));
  hx = (m * 2)+ 1;
  e  = m * m-p;
  hy = 1;
  while(e != 0)
  {
    while(e>0)
    {
      e -= hy;
      hy += 2;
    }
    while(e < 0)
    {
      e += hx;
      hx += 2;
    }
    if(Math.floor((hx+hy) / 2)>p)
    {
      hx = 2;
      hy = 0;
      break;
    }
  }
  m = Math.floor((hx+hy) / 2)-1;
  n = Math.floor((hx-hy) / 2);
  if(n != 1)
  {
    Fermat(m);
    Fermat(n);
  }
  else
  {
    t += p + " ";
    c++;
    if(c == 10)
    {
      c = 0;
      t += "<br/>";
    }
  }
}

function main()
{
  var p;
 
  p = parseInt(document.frm
               .inp_p.value);
  c = 0;
  if(!isNaN(p))
  {
    t = "";
    while(p % 2 == 0)
    {
      t +=  "2 ";
      p /= 2;
      c++;
      if(c == 10)
      {
        c = 0;
        t += "<br/>";
      }
    }
    if(p>1) Fermat (p);
  }
  else t = "Złe dane wejściowe!";
  document.getElementById("out")
  .innerHTML = t;
}

</script>

  </bohy>
</html>

do podrozdziału do strony

Większość liczb złożonych posiada czynniki pierwsze o małych wartościach. Algorytm Fermata jest szczególnie wydajny dla czynników leżących w pobliżu pierwiastka rozkładanej liczby. Natomiast czynniki małe są znajdowane wolno. Proponowane tutaj ulepszenie polega na połączeniu algorytmu próbnych dzieleń z algorytmem Fermata. W tablicy przechowujemy kolejne liczby pierwsze od \(2\) do \(1009\) – jest ich \(169\). Przed przekazaniem liczby \(p\) do procedury Fermat, z liczby \(p\) usuwamy wszystkie czynniki pierwsze przechowywane w tablicy. Jeśli liczba posiada takie czynniki, to po usunięciu ich staje się dużo mniejsza – zatem algorytm Fermata szybciej znajdzie pozostałe, duże czynniki pierwsze. Czas wykonania \(169\) próbnych dzieleń jest bardzo mały dla współczesnych procesorów, zatem warto zastosować to ulepszenie.

Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Pascal

// Metoda Fermata
// Data: 2.04.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

uses math;

procedure Fermat(p : int64);

var e, hx, hy, m, n : int64;

begin
  m  := ceil(sqrt(p));
  hx := (m shl 1)+1;
  e  := m*m-p;
  hy := 1;
  while e <> 0 do
  begin
    while e > 0 do
    begin
      dec(e, hy); inc(hy, 2);
    end;
    while e < 0 do
    begin
      inc(e, hx); inc(hx, 2);
    end;
    if (hx+hy) shr 1 > p then
    begin
      hx := 2;
      hy := 0;
      break;
    end;
  end;
  m := ((hx+hy) shr 1)-1;
  n := (hx-hy) shr 1;
  if n > 1 then
  begin
    Fermat(m); Fermat(n);
  end
  else write(p, ' ');
end;

const lp : array [0..168] of integer = (
  2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,  61,  67, 
 71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 
271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 
389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 
631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 
757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 
883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009);

var
  p : int64;
  i : integer;

begin
  readln(p);
  for i := 0 to 168 do
  begin
    if p >= lp[i] then
      while p mod lp[i] = 0 do
      begin
        write(lp[i], ' ');
        p := p div lp[i];
      end
    else break;
  end;
  if p > 1 then Fermat(p);
  writeln;
end.

C++

// Metoda Fermata
// Data: 2.04.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

#include <cstdlib>
#include <iostream>
#include <cmath>

using namespace std;

void Fermat(long long p)
{
  long long e, hx, hy, m, n;
  m  = (long long)sqrt(p);
  hx = (m << 1)+1;
  e  = m*m-p;
  hy = 1;
  while(e != 0)
  {
    while(e > 0)
    {
      e -= hy; hy += 2;
    }
    while(e < 0)
    {
      e += hx; hx += 2;
    }
    if(((hx+hy) >> 1) > p)
    {
      hx = 2;
      hy = 0;
      break;
    }
  }
  m = ((hx+hy) >> 1)-1;
  n = (hx-hy) >> 1;
  if(n != 1)
  {
    Fermat(m); Fermat(n);
  }
  else cout << p << " ";
}

int main()
{
  const int lp[] = {
  2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,  61,  67, 
 71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 
271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 
389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 
631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 
757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 
883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009};

  long long p;
  int i;

  cin >> p;
  for(i = 0; i < 169; i++)
  {
    if(p >= lp[i])
      while(p % lp[i] == 0)
      {
        cout << lp[i] << " ";
        p /= lp[i];
      }
    else break;
  }
  if(p > 1) Fermat(p);
  cout << endl;
  return 0;
}

Basic

' Metoda Fermata
' Data: 2.04.2008
' (C)2020 mgr Jerzy Wałaszek
'---------------------------

Sub Fermat(Byval p As Longint)

  Dim As Longint e, hx, hy, m, n

  m  = Clngint(Sqr(p))
  hx = (m Shl 1)+1
  e  = m*m-p
  hy = 1
  While e
    While e > 0
      e -= hy: hy += 2
    Wend
    While e < 0
      e += hx: hx += 2
    Wend
    If (hx+hy) Shr 1 > p Then
      hx = 2
      hy = 0
      Exit While
    End If
  Wend
  m = ((hx+hy) Shr 1)-1
  n = (hx-hy) Shr 1
  If n > 1 Then
    Fermat(m): Fermat(n)
  Else
    Print p;" ";
  End If
End Sub

Dim lp(168) As Integer => {_
  2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,  61,  67, _
 71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, _
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, _
271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, _
389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, _
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, _
631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, _
757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, _
883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009}

Dim p As Longint, i As Integer

Input p
For i = 0 To 168
  If p >= lp(i) Then
    While p Mod lp(i) = 0
      Print lp(i);" ";
      p \= lp(i)
    Wend
  Else
    Exit For
  End If
Next
If p > 1 Then Fermat(p)
Print
End

Python (dodatek)

# Metoda Fermata
# Data: 12.02.2024
# (C)2024 mgr Jerzy Wałaszek
# --------------------------

import math


def fermat(p):
    m = int(math.sqrt(p))
    hx = (m << 1)+1
    e = m*m-p
    hy = 1
    while e:
        while e > 0:
            e -= hy
            hy += 2
        while e < 0:
            e += hx
            hx += 2
        if (hx+hy) >> 1 > p:
            hx = 2
            hy = 0
            break
    m = ((hx+hy) >> 1)-1
    n = (hx-hy) >> 1
    if n > 1:
        fermat(m)
        fermat(n)
    else:
        print(p, end=" ")
    return


lp = [ \
      2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,  61,  67, \
     71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, \
    167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, \
    271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, \
    389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, \
    503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, \
    631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, \
    757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, \
    883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009]

p = int(input())
for i in range(169):
    if p >= lp[i]:
        while not (p % lp[i]):
            print(lp[i], end=" ")
            p //= lp[i]
    else:
        break
if p > 1:
    fermat(p)
print()

Wynik:

1234567891011
3 7 13 67 107 630803

1111111111111111
11 17 73 101 137 5882353

222222222222222222
2 3 3 7 11 13 19 37 333667 52579

do podrozdziału do strony

Czynniki pierwsze – metoda Fermata

Problem

Rozwiązanie 1

Algorytm rozkładu liczby naturalnej na czynniki pierwsze metodą Fermata

Wejście:

Wyjście:

Elementy pomocnicze:

Lista kroków dla procedury rekurencyjnej Fermat(p)

Przykładowe programy

Rozwiązanie 2

Przykładowe programy

Algorytm rozkładu liczby naturalnej na czynniki pierwsze metodą Fermata

Wejście:

Wyjście:

Elementy pomocnicze:

Lista kroków dla procedury rekurencyjnej Fermat(p)

Przykładowe programy

Rozwiązanie 3

Przykładowe programy