Fermat numbers

You are encouraged to solve this task according to the task description, using any language you may know.
In mathematics, a Fermat number, named after Pierre de Fermat who first studied them, is a positive integer of the form Fn = 22n + 1 where n is a non-negative integer.
Despite the simplicity of generating Fermat numbers, they have some powerful mathematical properties and are extensively used in cryptography & pseudo-random number generation, and are often linked to other number theoric fields.
As of this writing, (mid 2019), there are only five known prime Fermat numbers, the first five (F0 through F4). Only the first twelve Fermat numbers have been completely factored, though many have been partially factored.
- Task
- Write a routine (function, procedure, whatever) to generate Fermat numbers.
- Use the routine to find and display here, on this page, the first 10 Fermat numbers - F0 through F9.
- Find and display here, on this page, the prime factors of as many Fermat numbers as you have patience for. (Or as many as can be found in five minutes or less of processing time). Note: if you make it past F11, there may be money, and certainly will be acclaim in it for you.
- See also
ALGOL 68
Uses Algol 68G's LONG LONG INT which has programmer specifiable precision.
The source of primes.incl.a68 is on another page on Rosetta Code, see the above link.
Note that ALGOL 68 Genie version 3 issues a warning: "261 digits precision impacts performance".
BEGIN # find and factorise some Fermat numbers: F(n) = 2^(2^n) + 1 #
PR read "primes.incl.a68" PR # include prime utilities #
PR precision 256 PR # set the precision of LONG LONG INT #
PROC gcd = ( LONG LONG INT x, y )LONG LONG INT: # iterative gcd #
BEGIN
LONG LONG INT a := x, b := y;
WHILE b /= 0 DO
LONG LONG INT next a = b;
b := a MOD b;
a := next a
OD;
ABS a
END # gcd # ;
# returns a prime factor (if possible) of n, if n is prime, n is returned #
PROC pollard rho = ( LONG LONG INT n )LONG LONG INT:
IF is probably prime( n )
THEN n
ELIF LONG LONG INT x := 2, y := 2, d := 1;
PROC g = ( LONG LONG INT gx )LONG LONG INT: ( ( gx * gx ) + 1 ) MOD n;
WHILE d = 1 DO
x := g( x );
y := g( g( y ) );
d := gcd( ABS( x - y ), n )
OD;
d = n
THEN print( ( "pollard rho found non probably prime n for: ", n, newline ) );
n
ELIF LONG LONG INT other d = n OVER d;
d > other d
THEN other d
ELSE d
FI # pollard rho # ;
# returns the lowest prime factor of n, or n if n is prime #
PROC prime factor = ( LONG LONG INT n )LONG LONG INT:
IF LONG LONG INT d := pollard rho( n );
d = n
THEN d
ELSE # check for a lower factor #
LONG LONG INT other d := n OVER d;
LONG LONG INT d1 := pollard rho( other d );
WHILE d1 < d DO
d := d1;
other d := other d OVER d;
d1 := pollard rho( other d )
OD;
IF d1 < d THEN d1 ELSE d FI
FI # prime factor # ;
# task #
INT p2 := 1;
FOR i FROM 0 TO 9 DO
LONG LONG INT fn = 1 + ( LONG LONG 2 )^p2;
print( ( "F(", whole( i, 0 ), "): ", whole( fn, 0 ) ) );
IF i < 7 THEN
print( ( ", " ) );
LONG LONG INT pf = prime factor( fn );
IF pf = fn THEN
print( ( "prime" ) )
ELSE
print( ( whole( pf, 0 ), " x ", whole( fn OVER pf, 0 ) ) )
FI
FI;
print( ( newline ) );
p2 *:= 2
OD
END
- Output:
F(0): 3, prime F(1): 5, prime F(2): 17, prime F(3): 257, prime F(4): 65537, prime F(5): 4294967297, 641 x 6700417 F(6): 18446744073709551617, 274177 x 67280421310721 F(7): 340282366920938463463374607431768211457 F(8): 115792089237316195423570985008687907853269984665640564039457584007913129639937 F(9): 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097
Arturo
nPowers: [1 2 4 8 16 32 64 128 256 512]
fermatSet: map 0..9 'x -> 1 + 2 ^ nPowers\[x]
loop 0..9 'i ->
print ["F(" i ") =" fermatSet\[i]]
print ""
loop 0..9 'i ->
print ["Prime factors of F(" i ") =" factors.prime fermatSet\[i]]
C
Compile with :
gcc -o fermat fermat.c -lgmp
#include <stdlib.h>
#include <stdio.h>
#include <gmp.h>
void mpz_factors(mpz_t n) {
int factors = 0;
mpz_t s, m, p;
mpz_init(s), mpz_init(m), mpz_init(p);
mpz_set_ui(m, 3);
mpz_set(p, n);
mpz_sqrt(s, p);
while (mpz_cmp(m, s) < 0) {
if (mpz_divisible_p(p, m)) {
gmp_printf("%Zd ", m);
mpz_fdiv_q(p, p, m);
mpz_sqrt(s, p);
factors ++;
}
mpz_add_ui(m, m, 2);
}
if (factors == 0) printf("PRIME\n");
else gmp_printf("%Zd\n", p);
}
int main(int argc, char const *argv[]) {
mpz_t fermat;
mpz_init_set_ui(fermat, 3);
printf("F(0) = 3 -> PRIME\n");
for (unsigned i = 1; i < 10; i ++) {
mpz_sub_ui(fermat, fermat, 1);
mpz_mul(fermat, fermat, fermat);
mpz_add_ui(fermat, fermat, 1);
gmp_printf("F(%d) = %Zd -> ", i, fermat);
mpz_factors(fermat);
}
return 0;
}
F(0) = 3 -> PRIME F(1) = 5 -> PRIME F(2) = 17 -> PRIME F(3) = 257 -> PRIME F(4) = 65537 -> PRIME F(5) = 4294967297 -> 641 6700417 F(6) = 18446744073709551617 -> 274177 67280421310721 F(7) = 340282366920938463463374607431768211457 -> 59649589127497217 5704689200685129054721 F(8) = 115792089237316195423570985008687907853269984665640564039457584007913129639937 -> 1238926361552897 93461639715357977769163558199606896584051237541638188580280321 ......
C++
Built and tested on macOS 10.15, CPU: 3.2 GHz Intel Core i5. Execution time is about 12 minutes.
#include <iostream>
#include <vector>
#include <boost/integer/common_factor.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/miller_rabin.hpp>
typedef boost::multiprecision::cpp_int integer;
integer fermat(unsigned int n) {
unsigned int p = 1;
for (unsigned int i = 0; i < n; ++i)
p *= 2;
return 1 + pow(integer(2), p);
}
inline void g(integer& x, const integer& n) {
x *= x;
x += 1;
x %= n;
}
integer pollard_rho(const integer& n) {
integer x = 2, y = 2, d = 1, z = 1;
int count = 0;
for (;;) {
g(x, n);
g(y, n);
g(y, n);
d = abs(x - y);
z = (z * d) % n;
++count;
if (count == 100) {
d = gcd(z, n);
if (d != 1)
break;
z = 1;
count = 0;
}
}
if (d == n)
return 0;
return d;
}
std::vector<integer> get_prime_factors(integer n) {
std::vector<integer> factors;
for (;;) {
if (miller_rabin_test(n, 25)) {
factors.push_back(n);
break;
}
integer f = pollard_rho(n);
if (f == 0) {
factors.push_back(n);
break;
}
factors.push_back(f);
n /= f;
}
return factors;
}
void print_vector(const std::vector<integer>& factors) {
if (factors.empty())
return;
auto i = factors.begin();
std::cout << *i++;
for (; i != factors.end(); ++i)
std::cout << ", " << *i;
std::cout << '\n';
}
int main() {
std::cout << "First 10 Fermat numbers:\n";
for (unsigned int i = 0; i < 10; ++i)
std::cout << "F(" << i << ") = " << fermat(i) << '\n';
std::cout << "\nPrime factors:\n";
for (unsigned int i = 0; i < 9; ++i) {
std::cout << "F(" << i << "): ";
print_vector(get_prime_factors(fermat(i)));
}
return 0;
}
- Output:
First 10 Fermat numbers: F(0) = 3 F(1) = 5 F(2) = 17 F(3) = 257 F(4) = 65537 F(5) = 4294967297 F(6) = 18446744073709551617 F(7) = 340282366920938463463374607431768211457 F(8) = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F(9) = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Prime factors: F(0): 3 F(1): 5 F(2): 17 F(3): 257 F(4): 65537 F(5): 641, 6700417 F(6): 274177, 67280421310721 F(7): 59649589127497217, 5704689200685129054721 F(8): 1238926361552897, 93461639715357977769163558199606896584051237541638188580280321
Common Lisp
This uses the 'factor' function defined in the page Prime Decomposition http://rosettacode.org/wiki/Prime_decomposition#Common_Lisp
(defun fermat-number (n)
"Return the n-th Fermat number"
(1+ (expt 2 (expt 2 n))) )
(defun factor (n &optional (acc '()))
"Return the list of factors of n"
(when (> n 1) (loop with max-d = (isqrt n)
for d = 2 then (if (evenp d) (1+ d) (+ d 2)) do
(cond ((> d max-d) (return (cons (list n 1) acc)))
((zerop (rem n d))
(return (factor (truncate n d) (if (eq d (caar acc))
(cons
(list (caar acc) (1+ (cadar acc)))
(cdr acc))
(cons (list d 1) acc)))))))))
- Output:
(dotimes (i 8) (format t "~d: ~d = ~d~%" i (fermat-number i) (factors (fermat-number i)))) 0: 3 = (3) 1: 5 = (5) 2: 17 = (17) 3: 257 = (257) 4: 65537 = (65537) 5: 4294967297 = (641 6700417) 6: 18446744073709551617 = (274177 (67280421310721)) 7: 340282366920938463463374607431768211457 = ((340282366920938463463374607431768211457)) 8: 115792089237316195423570985008687907853269984665640564039457584007913129639937 = ((115792089237316195423570985008687907853269984665640564039457584007913129639937)) 9: 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 = ((13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097)) 10: 179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137217 = ((179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137217)) 11: 32317006071311007300714876688669951960444102669715484032130345427524655138867890893197201411522913463688717960921898019494119559150490921095088152386448283120630877367300996091750197750389652106796057638384067568276792218642619756161838094338476170470581645852036305042887575891541065808607552399123930385521914333389668342420684974786564569494856176035326322058077805659331026192708460314150258592864177116725943603718461857357598351152301645904403697613233287231227125684710820209725157101726931323469678542580656697935045997268352998638215525166389437335543602135433229604645318478604952148193555853611059596230657 = (319489 (101152171346465785365739905563790778275446424351747584524444802254614885454171789617787158279386499891040749324458425859713854183244152133860909616251101863039512713637405344446131784663602352840930541077733717180487566766438342966931062084574042206368862921265008513729385286790910065162204496552694867070609361616173540692858549041708993328392702647150062512506151403207406283761595736673720405375033810606080158013948717662760215065784116654734290374983906448207065425365852408720566771005345821995341556493590254118091846659097349200248570452085641250044738949182704974520704370036542579394575574913724915713)) 12: 1044388881413152506691752710716624382579964249047383780384233483283953907971557456848826811934997558340890106714439262837987573438185793607263236087851365277945956976543709998340361590134383718314428070011855946226376318839397712745672334684344586617496807908705803704071284048740118609114467977783598029006686938976881787785946905630190260940599579453432823469303026696443059025015972399867714215541693835559885291486318237914434496734087811872639496475100189041349008417061675093668333850551032972088269550769983616369411933015213796825837188091833656751221318492846368125550225998300412344784862595674492194617023806505913245610825731835380087608622102834270197698202313169017678006675195485079921636419370285375124784014907159135459982790513399611551794271106831134090584272884279791554849782954323534517065223269061394905987693002122963395687782878948440616007412945674919823050571642377154816321380631045902916136926708342856440730447899971901781465763473223850267253059899795996090799469201774624817718449867455659250178329070473119433165550807568221846571746373296884912819520317457002440926616910874148385078411929804522981857338977648103126085903001302413467189726673216491511131602920781738033436090243804708340403154190337 = (114689 (9106268965752186405773463110818163752233991481723476361152625650968740750826648212547208641935996986118024454955917854702609434541985662158212523327009262247869952450049350838706079834460006786304075107567909269645531121898331250125751682239313156601738683820643686003638396435055834553570682260579462973839574318172464558815116581626749391315641251152532705571615644886981829338611134458123396450764186936496833100701185274214915961723337127995182593580031119299575446791424418154036863609858251201843852076223383379133238000289598800458955855329052103961332983048473420515918928565951506637819342706575976725030506905683310915700945442329953941604008255667676914945655757474715779252371155778495946746587469464160684843488975918662295274965457887082037460184558511575570318625886351712499453155527762335682281851520733417380809781252979478377941937578568481859702438295520231435016188495646093490407803983345420364088331996467459309353537828143080691834120737157445502646809195267166779721413577366833939771467773331873590129210913628329073978766992198221682739812652450408607796042492802295258713711959073218748776359806123717024800451461326745599716651128725627280643537507664130920416107218492950792456907321580171946770433))
Crystal
This uses the `factor` function from the `coreutils` library that comes standard with most GNU/Linux, BSD, and Unix systems. https://www.gnu.org/software/coreutils/ https://en.wikipedia.org/wiki/GNU_Core_Utilities
require "big"
def factors(n)
factors = `factor #{n}`.split(' ')[1..-1].map(&.to_big_i)
factors.group_by(&.itself).map { |prime, exp| [prime, exp.size] }
end
def fermat(n); (1.to_big_i << (1 << n)) | 1 end
puts "Value for each Fermat Number F0 .. F9."
(0..9).each { |n| puts "F#{n} = #{fermat(n)}" }
puts
puts "Factors for each Fermat Number F0 .. F8."
(0..8).each { |n| puts "F#{n} = #{factors fermat(n)}" }
System: Lenovo V570 (2011), I5-2410M, 2.9 GHz, Crystal 0.34 Run as: $ time crystal fermat.cr --release
- Output:
Value for each Fermat Number F0 .. F9. F0 = 3 F1 = 5 F2 = 17 F3 = 257 F4 = 65537 F5 = 4294967297 F6 = 18446744073709551617 F7 = 340282366920938463463374607431768211457 F8 = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F9 = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors for each Fermat Number F0 .. F8. F0 = [[3, 1]] F1 = [[5, 1]] F2 = [[17, 1]] F3 = [[257, 1]] F4 = [[65537, 1]] F5 = [[641, 1], [6700417, 1]] F6 = [[274177, 1], [67280421310721, 1]] F7 = [[59649589127497217, 1], [5704689200685129054721, 1]] F8 = [[1238926361552897, 1], [93461639715357977769163558199606896584051237541638188580280321, 1]] crystal fermat.cr --release 174.19s user 0.19s system 100% cpu 2:54.08 total
Factor
USING: formatting io kernel lists lists.lazy math math.functions
math.primes.factors sequences ;
: lfermats ( -- list )
0 lfrom [ [ 1 2 2 ] dip ^ ^ + ] lmap-lazy ;
CHAR: ₀ 10 lfermats ltake list>array [
"First 10 Fermat numbers:" print
[ dupd "F%c = %d\n" printf 1 + ] each drop nl
] [
"Factors of first few Fermat numbers:" print [
dupd factors dup length 1 = " (prime)" "" ?
"Factors of F%c: %[%d, %]%s\n" printf 1 +
] each drop
] 2bi
- Output:
First 10 Fermat numbers: F₀ = 3 F₁ = 5 F₂ = 17 F₃ = 257 F₄ = 65537 F₅ = 4294967297 F₆ = 18446744073709551617 F₇ = 340282366920938463463374607431768211457 F₈ = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F₉ = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors of first few Fermat numbers: Factors of F₀: { 3 } (prime) Factors of F₁: { 5 } (prime) Factors of F₂: { 17 } (prime) Factors of F₃: { 257 } (prime) Factors of F₄: { 65537 } (prime) Factors of F₅: { 641, 6700417 } Factors of F₆: { 274177, 67280421310721 } ^D
Go
The first seven Fermat numbers are factorized almost instantly by the Pollard's rho algorithm and, on switching to Lenstra's elliptical curve method, F₇ only takes a couple of seconds (compared to over 12 minutes when using Pollard's rho).
It's back to Pollard rho for F₈ and the first prime factor of F₉ which takes a further 40 seconds or so. ECM doesn't seem to be effective for numbers as big as these.
As the second and third prime factors of F₉ are respectively 49 and 99 digits long there would be no chance of finding these any time soon so I haven't bothered.
The ECM implementation is based on the Python code here.
The timings are for my Intel Core i7-8565U laptop using Go 1.14.1 on Ubuntu 18.04.
package main
import (
"fmt"
"github.com/jbarham/primegen"
"math"
"math/big"
"math/rand"
"sort"
"time"
)
const (
maxCurves = 10000
maxRnd = 1 << 31
maxB1 = uint64(43 * 1e7)
maxB2 = uint64(2 * 1e10)
)
var (
zero = big.NewInt(0)
one = big.NewInt(1)
two = big.NewInt(2)
three = big.NewInt(3)
four = big.NewInt(4)
five = big.NewInt(5)
)
// Uses algorithm in Wikipedia article, including speed-up.
func pollardRho(n *big.Int) (*big.Int, error) {
// g(x) = (x^2 + 1) mod n
g := func(x, n *big.Int) *big.Int {
x2 := new(big.Int)
x2.Mul(x, x)
x2.Add(x2, one)
return x2.Mod(x2, n)
}
x, y, d := new(big.Int).Set(two), new(big.Int).Set(two), new(big.Int).Set(one)
t, z := new(big.Int), new(big.Int).Set(one)
count := 0
for {
x = g(x, n)
y = g(g(y, n), n)
t.Sub(x, y)
t.Abs(t)
t.Mod(t, n)
z.Mul(z, t)
count++
if count == 100 {
d.GCD(nil, nil, z, n)
if d.Cmp(one) != 0 {
break
}
z.Set(one)
count = 0
}
}
if d.Cmp(n) == 0 {
return nil, fmt.Errorf("Pollard's rho failure")
}
return d, nil
}
// Gets all primes under 'n' - uses a Sieve of Atkin under the hood.
func getPrimes(n uint64) []uint64 {
pg := primegen.New()
var primes []uint64
for {
prime := pg.Next()
if prime < n {
primes = append(primes, prime)
} else {
break
}
}
return primes
}
// Computes Stage 1 and Stage 2 bounds.
func computeBounds(n *big.Int) (uint64, uint64) {
le := len(n.String())
var b1, b2 uint64
switch {
case le <= 30:
b1, b2 = 2000, 147396
case le <= 40:
b1, b2 = 11000, 1873422
case le <= 50:
b1, b2 = 50000, 12746592
case le <= 60:
b1, b2 = 250000, 128992510
case le <= 70:
b1, b2 = 1000000, 1045563762
case le <= 80:
b1, b2 = 3000000, 5706890290
default:
b1, b2 = maxB1, maxB2
}
return b1, b2
}
// Adds two specified P and Q points (in Montgomery form). Assumes R = P - Q.
func pointAdd(px, pz, qx, qz, rx, rz, n *big.Int) (*big.Int, *big.Int) {
t := new(big.Int).Sub(px, pz)
u := new(big.Int).Add(qx, qz)
u.Mul(t, u)
t.Add(px, pz)
v := new(big.Int).Sub(qx, qz)
v.Mul(t, v)
upv := new(big.Int).Add(u, v)
umv := new(big.Int).Sub(u, v)
x := new(big.Int).Mul(upv, upv)
x.Mul(x, rz)
if x.Cmp(n) >= 0 {
x.Mod(x, n)
}
z := new(big.Int).Mul(umv, umv)
z.Mul(z, rx)
if z.Cmp(n) >= 0 {
z.Mod(z, n)
}
return x, z
}
// Doubles a point P (in Montgomery form).
func pointDouble(px, pz, n, a24 *big.Int) (*big.Int, *big.Int) {
u2 := new(big.Int).Add(px, pz)
u2.Mul(u2, u2)
v2 := new(big.Int).Sub(px, pz)
v2.Mul(v2, v2)
t := new(big.Int).Sub(u2, v2)
x := new(big.Int).Mul(u2, v2)
if x.Cmp(n) >= 0 {
x.Mod(x, n)
}
z := new(big.Int).Mul(a24, t)
z.Add(v2, z)
z.Mul(t, z)
if z.Cmp(n) >= 0 {
z.Mod(z, n)
}
return x, z
}
// Multiplies a specified point P (in Montgomery form) by a specified scalar.
func scalarMultiply(k, px, pz, n, a24 *big.Int) (*big.Int, *big.Int) {
sk := fmt.Sprintf("%b", k)
lk := len(sk)
qx := new(big.Int).Set(px)
qz := new(big.Int).Set(pz)
rx, rz := pointDouble(px, pz, n, a24)
for i := 1; i < lk; i++ {
if sk[i] == '1' {
qx, qz = pointAdd(rx, rz, qx, qz, px, pz, n)
rx, rz = pointDouble(rx, rz, n, a24)
} else {
rx, rz = pointAdd(qx, qz, rx, rz, px, pz, n)
qx, qz = pointDouble(qx, qz, n, a24)
}
}
return qx, qz
}
// Lenstra's two-stage ECM algorithm.
func ecm(n *big.Int) (*big.Int, error) {
if n.Cmp(one) == 0 || n.ProbablyPrime(10) {
return n, nil
}
b1, b2 := computeBounds(n)
dd := uint64(math.Sqrt(float64(b2)))
beta := make([]*big.Int, dd+1)
for i := 0; i < len(beta); i++ {
beta[i] = new(big.Int)
}
s := make([]*big.Int, 2*dd+2)
for i := 0; i < len(s); i++ {
s[i] = new(big.Int)
}
// stage 1 and stage 2 precomputations
curves := 0
logB1 := math.Log(float64(b1))
primes := getPrimes(b2)
numPrimes := len(primes)
idxB1 := sort.Search(len(primes), func(i int) bool { return primes[i] >= b1 })
// compute a B1-powersmooth integer 'k'
k := big.NewInt(1)
for i := 0; i < idxB1; i++ {
p := primes[i]
bp := new(big.Int).SetUint64(p)
t := uint64(logB1 / math.Log(float64(p)))
bt := new(big.Int).SetUint64(t)
bt.Exp(bp, bt, nil)
k.Mul(k, bt)
}
g := big.NewInt(1)
for (g.Cmp(one) == 0 || g.Cmp(n) == 0) && curves <= maxCurves {
curves++
st := int64(6 + rand.Intn(maxRnd-5))
sigma := big.NewInt(st)
// generate a new random curve in Montgomery form with Suyama's parameterization
u := new(big.Int).Mul(sigma, sigma)
u.Sub(u, five)
u.Mod(u, n)
v := new(big.Int).Mul(four, sigma)
v.Mod(v, n)
vmu := new(big.Int).Sub(v, u)
a := new(big.Int).Mul(vmu, vmu)
a.Mul(a, vmu)
t := new(big.Int).Mul(three, u)
t.Add(t, v)
a.Mul(a, t)
t.Mul(four, u)
t.Mul(t, u)
t.Mul(t, u)
t.Mul(t, v)
a.Quo(a, t)
a.Sub(a, two)
a.Mod(a, n)
a24 := new(big.Int).Add(a, two)
a24.Quo(a24, four)
// stage 1
px := new(big.Int).Mul(u, u)
px.Mul(px, u)
t.Mul(v, v)
t.Mul(t, v)
px.Quo(px, t)
px.Mod(px, n)
pz := big.NewInt(1)
qx, qz := scalarMultiply(k, px, pz, n, a24)
g.GCD(nil, nil, n, qz)
// if stage 1 is successful, return a non-trivial factor else
// move on to stage 2
if g.Cmp(one) != 0 && g.Cmp(n) != 0 {
return g, nil
}
// stage 2
s[1], s[2] = pointDouble(qx, qz, n, a24)
s[3], s[4] = pointDouble(s[1], s[2], n, a24)
beta[1].Mul(s[1], s[2])
beta[1].Mod(beta[1], n)
beta[2].Mul(s[3], s[4])
beta[2].Mod(beta[2], n)
for d := uint64(3); d <= dd; d++ {
d2 := 2 * d
s[d2-1], s[d2] = pointAdd(s[d2-3], s[d2-2], s[1], s[2], s[d2-5], s[d2-4], n)
beta[d].Mul(s[d2-1], s[d2])
beta[d].Mod(beta[d], n)
}
g.SetUint64(1)
b := new(big.Int).SetUint64(b1 - 1)
rx, rz := scalarMultiply(b, qx, qz, n, a24)
t.Mul(two, new(big.Int).SetUint64(dd))
t.Sub(b, t)
tx, tz := scalarMultiply(t, qx, qz, n, a24)
q, step := idxB1, 2*dd
for r := b1 - 1; r < b2; r += step {
alpha := new(big.Int).Mul(rx, rz)
alpha.Mod(alpha, n)
limit := r + step
for q < numPrimes && primes[q] <= limit {
d := (primes[q] - r) / 2
t := new(big.Int).Sub(rx, s[2*d-1])
f := new(big.Int).Add(rz, s[2*d])
f.Mul(t, f)
f.Sub(f, alpha)
f.Add(f, beta[d])
g.Mul(g, f)
g.Mod(g, n)
q++
}
trx := new(big.Int).Set(rx)
trz := new(big.Int).Set(rz)
rx, rz = pointAdd(rx, rz, s[2*dd-1], s[2*dd], tx, tz, n)
tx.Set(trx)
tz.Set(trz)
}
g.GCD(nil, nil, n, g)
}
// no non-trivial factor found, return an error
if curves > maxCurves {
return zero, fmt.Errorf("maximum curves exceeded before a factor was found")
}
return g, nil
}
// find prime factors of 'n' using an appropriate method.
func primeFactors(n *big.Int) ([]*big.Int, error) {
var res []*big.Int
if n.ProbablyPrime(10) {
return append(res, n), nil
}
le := len(n.String())
var factor1 *big.Int
var err error
if le > 20 && le <= 60 {
factor1, err = ecm(n)
} else {
factor1, err = pollardRho(n)
}
if err != nil {
return nil, err
}
if !factor1.ProbablyPrime(10) {
return nil, fmt.Errorf("first factor is not prime")
}
factor2 := new(big.Int)
factor2.Quo(n, factor1)
if !factor2.ProbablyPrime(10) {
return nil, fmt.Errorf("%d (second factor is not prime)", factor1)
}
return append(res, factor1, factor2), nil
}
func fermatNumbers(n int) (res []*big.Int) {
f := new(big.Int).SetUint64(3) // 2^1 + 1
for i := 0; i < n; i++ {
t := new(big.Int).Set(f)
res = append(res, t)
f.Sub(f, one)
f.Mul(f, f)
f.Add(f, one)
}
return res
}
func main() {
start := time.Now()
rand.Seed(time.Now().UnixNano())
fns := fermatNumbers(10)
fmt.Println("First 10 Fermat numbers:")
for i, f := range fns {
fmt.Printf("F%c = %d\n", 0x2080+i, f)
}
fmt.Println("\nFactors of first 10 Fermat numbers:")
for i, f := range fns {
fmt.Printf("F%c = ", 0x2080+i)
factors, err := primeFactors(f)
if err != nil {
fmt.Println(err)
continue
}
for _, factor := range factors {
fmt.Printf("%d ", factor)
}
if len(factors) == 1 {
fmt.Println("- prime")
} else {
fmt.Println()
}
}
fmt.Printf("\nTook %s\n", time.Since(start))
}
- Output:
First 10 Fermat numbers: F₀ = 3 F₁ = 5 F₂ = 17 F₃ = 257 F₄ = 65537 F₅ = 4294967297 F₆ = 18446744073709551617 F₇ = 340282366920938463463374607431768211457 F₈ = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F₉ = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors of first 10 Fermat numbers: F₀ = 3 - prime F₁ = 5 - prime F₂ = 17 - prime F₃ = 257 - prime F₄ = 65537 - prime F₅ = 641 6700417 F₆ = 274177 67280421310721 F₇ = 59649589127497217 5704689200685129054721 F₈ = 1238926361552897 93461639715357977769163558199606896584051237541638188580280321 F₉ = 2424833 (second factor is not prime) Took 41.683532956s
Haskell
import Data.Numbers.Primes (primeFactors)
import Data.Bool (bool)
fermat :: Integer -> Integer
fermat = succ . (2 ^) . (2 ^)
fermats :: [Integer]
fermats = fermat <$> [0 ..]
--------------------------- TEST ---------------------------
main :: IO ()
main =
mapM_
putStrLn
[ fTable "First 10 Fermats:" show show fermat [0 .. 9]
, fTable
"Factors of first 7:"
show
showFactors
primeFactors
(take 7 fermats)
]
------------------------- DISPLAY --------------------------
fTable :: String -> (a -> String) -> (b -> String) -> (a -> b) -> [a] -> String
fTable s xShow fxShow f xs =
unlines $
s : fmap (((++) . rjust w ' ' . xShow) <*> ((" -> " ++) . fxShow . f)) xs
where
rjust n c = drop . length <*> (replicate n c ++)
w = maximum (length . xShow <$> xs)
showFactors :: [Integer] -> String
showFactors x
| 1 < length x = show x
| otherwise = "(prime)"
- Output:
First 10 Fermats: 0 -> 3 1 -> 5 2 -> 17 3 -> 257 4 -> 65537 5 -> 4294967297 6 -> 18446744073709551617 7 -> 340282366920938463463374607431768211457 8 -> 115792089237316195423570985008687907853269984665640564039457584007913129639937 9 -> 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors of first 7: 3 -> (prime) 5 -> (prime) 17 -> (prime) 257 -> (prime) 65537 -> (prime) 4294967297 -> [641,6700417] 18446744073709551617 -> [274177,67280421310721]
J
fermat =: 1 1 p. 2 ^ 2 ^ x: (,. fermat)i.10 0 3 1 5 2 17 3 257 4 65537 5 4294967297 6 18446744073709551617 7 340282366920938463463374607431768211457 8 115792089237316195423570985008687907853269984665640564039457584007913129639937 9 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 (; q:@:fermat)&>i.7 +-+---------------------+ |0|3 | +-+---------------------+ |1|5 | +-+---------------------+ |2|17 | +-+---------------------+ |3|257 | +-+---------------------+ |4|65537 | +-+---------------------+ |5|641 6700417 | +-+---------------------+ |6|274177 67280421310721| +-+---------------------+
Java
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.stream.Collectors;
public class FermatNumbers {
public static void main(String[] args) {
System.out.println("First 10 Fermat numbers:");
for ( int i = 0 ; i < 10 ; i++ ) {
System.out.printf("F[%d] = %s\n", i, fermat(i));
}
System.out.printf("%nFirst 12 Fermat numbers factored:%n");
for ( int i = 0 ; i < 13 ; i++ ) {
System.out.printf("F[%d] = %s\n", i, getString(getFactors(i, fermat(i))));
}
}
private static String getString(List<BigInteger> factors) {
if ( factors.size() == 1 ) {
return factors.get(0) + " (PRIME)";
}
return factors.stream().map(v -> v.toString()).map(v -> v.startsWith("-") ? "(C" + v.replace("-", "") + ")" : v).collect(Collectors.joining(" * "));
}
private static Map<Integer, String> COMPOSITE = new HashMap<>();
static {
COMPOSITE.put(9, "5529");
COMPOSITE.put(10, "6078");
COMPOSITE.put(11, "1037");
COMPOSITE.put(12, "5488");
COMPOSITE.put(13, "2884");
}
private static List<BigInteger> getFactors(int fermatIndex, BigInteger n) {
List<BigInteger> factors = new ArrayList<>();
BigInteger factor = BigInteger.ONE;
while ( true ) {
if ( n.isProbablePrime(100) ) {
factors.add(n);
break;
}
else {
if ( COMPOSITE.containsKey(fermatIndex) ) {
String stop = COMPOSITE.get(fermatIndex);
if ( n.toString().startsWith(stop) ) {
factors.add(new BigInteger("-" + n.toString().length()));
break;
}
}
factor = pollardRhoFast(n);
if ( factor.compareTo(BigInteger.ZERO) == 0 ) {
factors.add(n);
break;
}
else {
factors.add(factor);
n = n.divide(factor);
}
}
}
return factors;
}
private static final BigInteger TWO = BigInteger.valueOf(2);
private static BigInteger fermat(int n) {
return TWO.pow((int)Math.pow(2, n)).add(BigInteger.ONE);
}
// See: https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
@SuppressWarnings("unused")
private static BigInteger pollardRho(BigInteger n) {
BigInteger x = BigInteger.valueOf(2);
BigInteger y = BigInteger.valueOf(2);
BigInteger d = BigInteger.ONE;
while ( d.compareTo(BigInteger.ONE) == 0 ) {
x = pollardRhoG(x, n);
y = pollardRhoG(pollardRhoG(y, n), n);
d = x.subtract(y).abs().gcd(n);
}
if ( d.compareTo(n) == 0 ) {
return BigInteger.ZERO;
}
return d;
}
// Includes Speed Up of 100 multiples and 1 GCD, instead of 100 multiples and 100 GCDs.
// See Variants section of Wikipedia article.
// Testing F[8] = 1238926361552897 * Prime
// This variant = 32 sec.
// Standard algorithm = 107 sec.
private static BigInteger pollardRhoFast(BigInteger n) {
long start = System.currentTimeMillis();
BigInteger x = BigInteger.valueOf(2);
BigInteger y = BigInteger.valueOf(2);
BigInteger d = BigInteger.ONE;
int count = 0;
BigInteger z = BigInteger.ONE;
while ( true ) {
x = pollardRhoG(x, n);
y = pollardRhoG(pollardRhoG(y, n), n);
d = x.subtract(y).abs();
z = z.multiply(d).mod(n);
count++;
if ( count == 100 ) {
d = z.gcd(n);
if ( d.compareTo(BigInteger.ONE) != 0 ) {
break;
}
z = BigInteger.ONE;
count = 0;
}
}
long end = System.currentTimeMillis();
System.out.printf(" Pollard rho try factor %s elapsed time = %d ms (factor = %s).%n", n, (end-start), d);
if ( d.compareTo(n) == 0 ) {
return BigInteger.ZERO;
}
return d;
}
private static BigInteger pollardRhoG(BigInteger x, BigInteger n) {
return x.multiply(x).add(BigInteger.ONE).mod(n);
}
}
Output includes composite numbers attempted to factor with Pollard rho.
- Output:
First 10 Fermat numbers: F[0] = 3 F[1] = 5 F[2] = 17 F[3] = 257 F[4] = 65537 F[5] = 4294967297 F[6] = 18446744073709551617 F[7] = 340282366920938463463374607431768211457 F[8] = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F[9] = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 First 12 Fermat numbers factored: F[0] = 3 (PRIME) F[1] = 5 (PRIME) F[2] = 17 (PRIME) F[3] = 257 (PRIME) F[4] = 65537 (PRIME) Pollard rho try factor 4294967297 elapsed time = 2 ms (factor = 641). F[5] = 641 * 6700417 Pollard rho try factor 18446744073709551617 elapsed time = 6 ms (factor = 274177). F[6] = 274177 * 67280421310721 Pollard rho try factor 340282366920938463463374607431768211457 elapsed time = 574251 ms (factor = 59649589127497217). F[7] = 59649589127497217 * 5704689200685129054721 Pollard rho try factor 115792089237316195423570985008687907853269984665640564039457584007913129639937 elapsed time = 31640 ms (factor = 1238926361552897). F[8] = 1238926361552897 * 93461639715357977769163558199606896584051237541638188580280321 Pollard rho try factor 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 elapsed time = 30 ms (factor = 2424833). F[9] = 2424833 * (C148) Pollard rho try factor 179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137217 elapsed time = 120 ms (factor = 45592577). Pollard rho try factor 3942951359960012586542991835686376608231592127249807732373409846031135195659174148737161255930050543559319182152642816343958573976075461198274610155058226350701077796608546283231637018483208223116080561800334422176622099740983337736621316898600121619871377542107047343253864459964167331555646795960321 elapsed time = 5026 ms (factor = 6487031809). F[10] = 45592577 * 6487031809 * (C291) Pollard rho try factor 32317006071311007300714876688669951960444102669715484032130345427524655138867890893197201411522913463688717960921898019494119559150490921095088152386448283120630877367300996091750197750389652106796057638384067568276792218642619756161838094338476170470581645852036305042887575891541065808607552399123930385521914333389668342420684974786564569494856176035326322058077805659331026192708460314150258592864177116725943603718461857357598351152301645904403697613233287231227125684710820209725157101726931323469678542580656697935045997268352998638215525166389437335543602135433229604645318478604952148193555853611059596230657 elapsed time = 8 ms (factor = 974849). Pollard rho try factor 33150781373639412155846573868024639672856106606987835072026893834352453701925006737655987144186344206834820532125383540932102878651453631377873037143648178457002958783669056532601662155256508553423204658756451069116132055982639112479817996775373591674794399801442382402697828988429044712163168243619196804348072710121945390948428910309765481110260333687910970886853046635254307274981520537180895290310783635953818082306553996934497908037349010876970379631341148456045116407475229712217130141926525362871253437794629422541384355185626695660779797862427347553871011957167960991543632376506466281643163047416635393 elapsed time = 98 ms (factor = 319489). F[11] = 974849 * 319489 * (C606) Pollard rho try factor 1044388881413152506691752710716624382579964249047383780384233483283953907971557456848826811934997558340890106714439262837987573438185793607263236087851365277945956976543709998340361590134383718314428070011855946226376318839397712745672334684344586617496807908705803704071284048740118609114467977783598029006686938976881787785946905630190260940599579453432823469303026696443059025015972399867714215541693835559885291486318237914434496734087811872639496475100189041349008417061675093668333850551032972088269550769983616369411933015213796825837188091833656751221318492846368125550225998300412344784862595674492194617023806505913245610825731835380087608622102834270197698202313169017678006675195485079921636419370285375124784014907159135459982790513399611551794271106831134090584272884279791554849782954323534517065223269061394905987693002122963395687782878948440616007412945674919823050571642377154816321380631045902916136926708342856440730447899971901781465763473223850267253059899795996090799469201774624817718449867455659250178329070473119433165550807568221846571746373296884912819520317457002440926616910874148385078411929804522981857338977648103126085903001302413467189726673216491511131602920781738033436090243804708340403154190337 elapsed time = 75 ms (factor = 114689). Pollard rho try factor 9106268965752186405773463110818163752233991481723476361152625650968740750826648212547208641935996986118024454955917854702609434541985662158212523327009262247869952450049350838706079834460006786304075107567909269645531121898331250125751682239313156601738683820643686003638396435055834553570682260579462973839574318172464558815116581626749391315641251152532705571615644886981829338611134458123396450764186936496833100701185274214915961723337127995182593580031119299575446791424418154036863609858251201843852076223383379133238000289598800458955855329052103961332983048473420515918928565951506637819342706575976725030506905683310915700945442329953941604008255667676914945655757474715779252371155778495946746587469464160684843488975918662295274965457887082037460184558511575570318625886351712499453155527762335682281851520733417380809781252979478377941937578568481859702438295520231435016188495646093490407803983345420364088331996467459309353537828143080691834120737157445502646809195267166779721413577366833939771467773331873590129210913628329073978766992198221682739812652450408607796042492802295258713711959073218748776359806123717024800451461326745599716651128725627280643537507664130920416107218492950792456907321580171946770433 elapsed time = 301 ms (factor = 26017793). Pollard rho try factor 350001591824186871106763863899530746217943677302816436473017740242946077356624684213115564488348300185108877411543625345263121839042445381828217455916005721464151569276047005167043946981206545317048033534893189043572263100806868980998952610596646556521480658280419327021257968923645235918768446677220584153297488306270426504473941414890547838382804073832577020334339845312084285496895699728389585187497914849919557000902623608963141559444997044721763816786117073787751589515083702681348245404913906488680729999788351730419178916889637812821243344085799435845038164784900107242721493170135785069061880328434342030106354742821995535937481028077744396075726164309052460585559946407282864168038994640934638329525854255227752926198464207255983432778402344965903148839661825873175277621985527846249416909718758069997783882773355041329208102046913755441975327368023946523920699020098723785533557579080342841062805878477869513695185309048285123705067072486920463781103076554014502567884803571416673251784936825115787932810954867447447568320403976197134736485611912650805539603318790667901618038578533362100071745480995207732506742832634459994375828162163700807237997808869771569154136465922798310222055287047244647419069003284481 elapsed time = 1616 ms (factor = 63766529). F[12] = 114689 * 26017793 * 63766529 * (C1213)
jq
Works with gojq, the Go implementation of jq
Preliminaries
# To take advantage of gojq's arbitrary-precision integer arithmetic:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);
def gcd(a; b):
# subfunction expects [a,b] as input
# i.e. a ~ .[0] and b ~ .[1]
def rgcd: if .[1] == 0 then .[0]
else [.[1], .[0] % .[1]] | rgcd
end;
[a,b] | rgcd;
# This is fast because the state of `until` is just a number
def is_prime:
. as $n
| if ($n < 2) then false
elif ($n % 2 == 0) then $n == 2
elif ($n % 3 == 0) then $n == 3
elif ($n % 5 == 0) then $n == 5
elif ($n % 7 == 0) then $n == 7
elif ($n % 11 == 0) then $n == 11
elif ($n % 13 == 0) then $n == 13
elif ($n % 17 == 0) then $n == 17
elif ($n % 19 == 0) then $n == 19
elif ($n % 23 == 0) then $n == 23
elif ($n % 29 == 0) then $n == 29
elif ($n % 31 == 0) then $n == 31
elif ($n % 37 == 0) then $n == 37
elif ($n % 41 == 0) then $n == 41
else 43
| until( (. * .) > $n or ($n % . == 0); . + 2)
| . * . > $n
end;
Fermat Numbers
def fermat:
. as $n
| (2 | power( 2 | power($n))) + 1;
# https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
def pollardRho($x):
. as $n
| def g: (.*. + 1) % $n ;
{x:$x, y:$x, d:1}
| until(.d != 1;
.x |= g
| .y |= (g|g)
| .d = gcd((.x - .y)|length; $n) )
| if .d == $n then 0
else .d
end ;
def rhoPrimeFactors:
. as $n
| pollardRho(2)
| if . == 0
then [$n, 1]
else [., ($n / .)]
end ;
"The first 10 Fermat numbers are:",
[ range(0;10) | fermat ] as $fns
| (range(0;10) | "F\(.) is \($fns[.])"),
("\nFactors of the first 7 Fermat numbers:",
range(0;7) as $i
| $fns[$i]
| rhoPrimeFactors as $factors
| if $factors[1] == 1
then "F\($i) : rho-prime", " ... => \(if is_prime then "prime" else "not" end)"
else "F\($i) => \($factors)"
end )
- Output:
The first 10 Fermat numbers are: F0 is 3 F1 is 5 F2 is 17 F3 is 257 F4 is 65537 F5 is 4294967297 F6 is 18446744073709551617 F7 is 340282366920938463463374607431768211457 F8 is 115792089237316195423570985008687907853269984665640564039457584007913129639937 F9 is 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors of the first 7 Fermat numbers: F0 : rho-prime ... => prime F1 : rho-prime ... => prime F2 : rho-prime ... => prime F3 : rho-prime ... => prime F4 : rho-prime ... => prime F5 => [641,6700417] F6 => [274177,67280421310721]
Julia
using Primes
fermat(n) = BigInt(2)^(BigInt(2)^n) + 1
prettyprint(fdict) = replace(replace(string(fdict), r".+\(([^)]+)\)" => s"\1"), r"\=\>" => "^")
function factorfermats(max, nofactor=false)
for n in 0:max
fm = fermat(n)
if nofactor
println("Fermat number F($n) is $fm.")
continue
end
factors = factor(fm)
println("Fermat number F($n), $fm, ",
length(factors) < 2 ? "is prime." : "factors to $(prettyprint(factors)).")
end
end
factorfermats(9, true)
factorfermats(10)
- Output:
Fermat number F(0) is 3. Fermat number F(1) is 5. Fermat number F(2) is 17. Fermat number F(3) is 257. Fermat number F(4) is 65537. Fermat number F(5) is 4294967297. Fermat number F(6) is 18446744073709551617. Fermat number F(7) is 340282366920938463463374607431768211457. Fermat number F(8) is 115792089237316195423570985008687907853269984665640564039457584007913129639937. Fermat number F(9) is 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097. Fermat number F(0), 3, is prime. Fermat number F(1), 5, is prime. Fermat number F(2), 17, is prime. Fermat number F(3), 257, is prime. Fermat number F(4), 65537, is prime. Fermat number F(5), 4294967297, factors to 641^1,6700417^1. Fermat number F(6), 18446744073709551617, factors to 274177^1,67280421310721^1. ...waited >5 minutes
Kotlin
import java.math.BigInteger
import kotlin.math.pow
fun main() {
println("First 10 Fermat numbers:")
for (i in 0..9) {
println("F[$i] = ${fermat(i)}")
}
println()
println("First 12 Fermat numbers factored:")
for (i in 0..12) {
println("F[$i] = ${getString(getFactors(i, fermat(i)))}")
}
}
private fun getString(factors: List<BigInteger>): String {
return if (factors.size == 1) {
"${factors[0]} (PRIME)"
} else factors.map { it.toString() }
.joinToString(" * ") {
if (it.startsWith("-"))
"(C" + it.replace("-", "") + ")"
else it
}
}
private val COMPOSITE = mutableMapOf(
9 to "5529",
10 to "6078",
11 to "1037",
12 to "5488",
13 to "2884"
)
private fun getFactors(fermatIndex: Int, n: BigInteger): List<BigInteger> {
var n2 = n
val factors: MutableList<BigInteger> = ArrayList()
var factor: BigInteger
while (true) {
if (n2.isProbablePrime(100)) {
factors.add(n2)
break
} else {
if (COMPOSITE.containsKey(fermatIndex)) {
val stop = COMPOSITE[fermatIndex]
if (n2.toString().startsWith(stop!!)) {
factors.add(BigInteger("-" + n2.toString().length))
break
}
}
//factor = pollardRho(n)
factor = pollardRhoFast(n)
n2 = if (factor.compareTo(BigInteger.ZERO) == 0) {
factors.add(n2)
break
} else {
factors.add(factor)
n2.divide(factor)
}
}
}
return factors
}
private val TWO = BigInteger.valueOf(2)
private fun fermat(n: Int): BigInteger {
return TWO.pow(2.0.pow(n.toDouble()).toInt()).add(BigInteger.ONE)
}
// See: https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
@Suppress("unused")
private fun pollardRho(n: BigInteger): BigInteger {
var x = BigInteger.valueOf(2)
var y = BigInteger.valueOf(2)
var d = BigInteger.ONE
while (d.compareTo(BigInteger.ONE) == 0) {
x = pollardRhoG(x, n)
y = pollardRhoG(pollardRhoG(y, n), n)
d = (x - y).abs().gcd(n)
}
return if (d.compareTo(n) == 0) {
BigInteger.ZERO
} else d
}
// Includes Speed Up of 100 multiples and 1 GCD, instead of 100 multiples and 100 GCDs.
// See Variants section of Wikipedia article.
// Testing F[8] = 1238926361552897 * Prime
// This variant = 32 sec.
// Standard algorithm = 107 sec.
private fun pollardRhoFast(n: BigInteger): BigInteger {
val start = System.currentTimeMillis()
var x = BigInteger.valueOf(2)
var y = BigInteger.valueOf(2)
var d: BigInteger
var count = 0
var z = BigInteger.ONE
while (true) {
x = pollardRhoG(x, n)
y = pollardRhoG(pollardRhoG(y, n), n)
d = (x - y).abs()
z = (z * d).mod(n)
count++
if (count == 100) {
d = z.gcd(n)
if (d.compareTo(BigInteger.ONE) != 0) {
break
}
z = BigInteger.ONE
count = 0
}
}
val end = System.currentTimeMillis()
println(" Pollard rho try factor $n elapsed time = ${end - start} ms (factor = $d).")
return if (d.compareTo(n) == 0) {
BigInteger.ZERO
} else d
}
private fun pollardRhoG(x: BigInteger, n: BigInteger): BigInteger {
return (x * x + BigInteger.ONE).mod(n)
}
langur
val fermat = fn i: 2 ^ 2 ^ i + 1
val sz = '\u2080' # small zero code point
val factors = fn(var x) {
for[f=[]] i, s = 2, trunc(x ^/ 2); i < s; i += 1 {
if x div i {
f ~= [i]
x \= i
s = trunc(x ^/ 2)
}
} ~ [x]
}
writeln "first 10 Fermat numbers"
for i in 0..9 {
writeln "F{{i + sz:cp}} = {{fermat(i)}}"
}
writeln()
writeln "factors of first few Fermat numbers"
for i in 0..9 {
val ferm = fermat(i)
val fac = factors(ferm)
if len(fac) == 1 {
writeln "F{{i + sz:cp}} is prime"
} else {
writeln "F{{i + sz:cp}} factors: {{fac}}"
}
}
- Output:
first 10 Fermat numbers F₀ = 3 F₁ = 5 F₂ = 17 F₃ = 257 F₄ = 65537 F₅ = 4294967297 F₆ = 18446744073709551617 F₇ = 340282366920938463463374607431768211457 F₈ = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F₉ = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 factors of first few Fermat numbers F₀ is prime F₁ is prime F₂ is prime F₃ is prime F₄ is prime F₅ factors: [641, 6700417] F₆ factors: [274177, 67280421310721]
Mathematica / Wolfram Language
ClearAll[Fermat]
Fermat[n_] := 2^(2^n) + 1
Fermat /@ Range[0, 9]
Scan[FactorInteger /* Print, %]
- Output:
{3,5,17,257,65537,4294967297,18446744073709551617,340282366920938463463374607431768211457,115792089237316195423570985008687907853269984665640564039457584007913129639937,13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097} {{3,1}} {{5,1}} {{17,1}} {{257,1}} {{65537,1}} {{641,1},{6700417,1}} {{274177,1},{67280421310721,1}} {{59649589127497217,1},{5704689200685129054721,1}} {{1238926361552897,1},{93461639715357977769163558199606896584051237541638188580280321,1}} $Aborted
Nim
import math
import bignum
import strformat
import strutils
import tables
import times
const Composite = {9: "5529", 10: "6078", 11: "1037", 12: "5488", 13: "2884"}.toTable
const Subscripts = ["₀", "₁", "₂", "₃", "₄", "₅", "₆", "₇", "₈", "₉"]
let One = newInt(1)
#---------------------------------------------------------------------------------------------------
func fermat(n: int): Int {.inline.} = 2^(culong(2^n)) + 1
#---------------------------------------------------------------------------------------------------
template isProbablyPrime(n: Int): bool = n.probablyPrime(25) != 0
#---------------------------------------------------------------------------------------------------
func pollardRhoG(x, n: Int): Int {.inline.} = (x * x + 1) mod n
#---------------------------------------------------------------------------------------------------
proc pollardRhoFast(n: Int): Int =
let start = getTime()
var
x = newInt(2)
y = newInt(2)
count = 0
z = One
while true:
x = pollardRhoG(x, n)
y = pollardRhoG(pollardRhoG(y, n), n)
result = abs(x - y)
z = z * result mod n
inc count
if count == 100:
result = gcd(z, n)
if result != One: break
z = One
count = 0
let duration = (getTime() - start).inMilliseconds
echo fmt" Pollard rho try factor {n} elapsed time = {duration} ms (factor = {result})."
if result == n:
result = newInt(0)
#---------------------------------------------------------------------------------------------------
proc factors(fermatIndex: int; n: Int): seq[Int] =
var n = n
var factor: Int
while true:
if n.isProbablyPrime():
result.add(n)
break
if fermatIndex in Composite:
let stop = Composite[fermatIndex]
let s = $n
if s.startsWith(stop):
result.add(newInt(-s.len))
break
factor = pollardRhoFast(n)
if factor.isZero():
result.add(n)
break
result.add(factor)
n = n div factor
#---------------------------------------------------------------------------------------------------
func `$`(factors: seq[Int]): string =
if factors.len == 1:
result = fmt"{factors[0]} (PRIME)"
else:
result = $factors[0]
let start = result.high
for factor in factors[1..^1]:
result.addSep(" * ", start)
result.add(if factor < 0: fmt"(C{-factor})" else: $factor)
#---------------------------------------------------------------------------------------------------
func subscript(n: Natural): string =
var n = n
while true:
result.insert(Subscripts[n mod 10], 0)
n = n div 10
if n == 0: break
#———————————————————————————————————————————————————————————————————————————————————————————————————
echo "First 10 Fermat numbers:"
for i in 0..9:
echo fmt"F{subscript(i)} = {fermat(i)}"
echo ""
echo "First 12 Fermat numbers factored:"
for i in 0..12:
echo fmt"F{subscript(i)} = {factors(i, fermat(i))}"
- Output:
First 10 Fermat numbers: F₀ = 3 F₁ = 5 F₂ = 17 F₃ = 257 F₄ = 65537 F₅ = 4294967297 F₆ = 18446744073709551617 F₇ = 340282366920938463463374607431768211457 F₈ = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F₉ = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 First 12 Fermat numbers factored: F₀ = 3 (PRIME) F₁ = 5 (PRIME) F₂ = 17 (PRIME) F₃ = 257 (PRIME) F₄ = 65537 (PRIME) Pollard rho try factor 4294967297 elapsed time = 0 ms (factor = 641). F₅ = 641 * 6700417 Pollard rho try factor 18446744073709551617 elapsed time = 1 ms (factor = 274177). F₆ = 274177 * 67280421310721 Pollard rho try factor 340282366920938463463374607431768211457 elapsed time = 516705 ms (factor = 59649589127497217). F₇ = 59649589127497217 * 5704689200685129054721 Pollard rho try factor 115792089237316195423570985008687907853269984665640564039457584007913129639937 elapsed time = 20765 ms (factor = 1238926361552897). F₈ = 1238926361552897 * 93461639715357977769163558199606896584051237541638188580280321 Pollard rho try factor 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 elapsed time = 3 ms (factor = 2424833). F₉ = 2424833 * (C148) Pollard rho try factor 179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137217 elapsed time = 46 ms (factor = 45592577). Pollard rho try factor 3942951359960012586542991835686376608231592127249807732373409846031135195659174148737161255930050543559319182152642816343958573976075461198274610155058226350701077796608546283231637018483208223116080561800334422176622099740983337736621316898600121619871377542107047343253864459964167331555646795960321 elapsed time = 595 ms (factor = 6487031809). F₁₀ = 45592577 * 6487031809 * (C291) Pollard rho try factor 32317006071311007300714876688669951960444102669715484032130345427524655138867890893197201411522913463688717960921898019494119559150490921095088152386448283120630877367300996091750197750389652106796057638384067568276792218642619756161838094338476170470581645852036305042887575891541065808607552399123930385521914333389668342420684974786564569494856176035326322058077805659331026192708460314150258592864177116725943603718461857357598351152301645904403697613233287231227125684710820209725157101726931323469678542580656697935045997268352998638215525166389437335543602135433229604645318478604952148193555853611059596230657 elapsed time = 1 ms (factor = 974849). Pollard rho try factor 33150781373639412155846573868024639672856106606987835072026893834352453701925006737655987144186344206834820532125383540932102878651453631377873037143648178457002958783669056532601662155256508553423204658756451069116132055982639112479817996775373591674794399801442382402697828988429044712163168243619196804348072710121945390948428910309765481110260333687910970886853046635254307274981520537180895290310783635953818082306553996934497908037349010876970379631341148456045116407475229712217130141926525362871253437794629422541384355185626695660779797862427347553871011957167960991543632376506466281643163047416635393 elapsed time = 11 ms (factor = 319489). F₁₁ = 974849 * 319489 * (C606) Pollard rho try factor 1044388881413152506691752710716624382579964249047383780384233483283953907971557456848826811934997558340890106714439262837987573438185793607263236087851365277945956976543709998340361590134383718314428070011855946226376318839397712745672334684344586617496807908705803704071284048740118609114467977783598029006686938976881787785946905630190260940599579453432823469303026696443059025015972399867714215541693835559885291486318237914434496734087811872639496475100189041349008417061675093668333850551032972088269550769983616369411933015213796825837188091833656751221318492846368125550225998300412344784862595674492194617023806505913245610825731835380087608622102834270197698202313169017678006675195485079921636419370285375124784014907159135459982790513399611551794271106831134090584272884279791554849782954323534517065223269061394905987693002122963395687782878948440616007412945674919823050571642377154816321380631045902916136926708342856440730447899971901781465763473223850267253059899795996090799469201774624817718449867455659250178329070473119433165550807568221846571746373296884912819520317457002440926616910874148385078411929804522981857338977648103126085903001302413467189726673216491511131602920781738033436090243804708340403154190337 elapsed time = 11 ms (factor = 114689). Pollard rho try factor 9106268965752186405773463110818163752233991481723476361152625650968740750826648212547208641935996986118024454955917854702609434541985662158212523327009262247869952450049350838706079834460006786304075107567909269645531121898331250125751682239313156601738683820643686003638396435055834553570682260579462973839574318172464558815116581626749391315641251152532705571615644886981829338611134458123396450764186936496833100701185274214915961723337127995182593580031119299575446791424418154036863609858251201843852076223383379133238000289598800458955855329052103961332983048473420515918928565951506637819342706575976725030506905683310915700945442329953941604008255667676914945655757474715779252371155778495946746587469464160684843488975918662295274965457887082037460184558511575570318625886351712499453155527762335682281851520733417380809781252979478377941937578568481859702438295520231435016188495646093490407803983345420364088331996467459309353537828143080691834120737157445502646809195267166779721413577366833939771467773331873590129210913628329073978766992198221682739812652450408607796042492802295258713711959073218748776359806123717024800451461326745599716651128725627280643537507664130920416107218492950792456907321580171946770433 elapsed time = 18 ms (factor = 26017793). Pollard rho try factor 350001591824186871106763863899530746217943677302816436473017740242946077356624684213115564488348300185108877411543625345263121839042445381828217455916005721464151569276047005167043946981206545317048033534893189043572263100806868980998952610596646556521480658280419327021257968923645235918768446677220584153297488306270426504473941414890547838382804073832577020334339845312084285496895699728389585187497914849919557000902623608963141559444997044721763816786117073787751589515083702681348245404913906488680729999788351730419178916889637812821243344085799435845038164784900107242721493170135785069061880328434342030106354742821995535937481028077744396075726164309052460585559946407282864168038994640934638329525854255227752926198464207255983432778402344965903148839661825873175277621985527846249416909718758069997783882773355041329208102046913755441975327368023946523920699020098723785533557579080342841062805878477869513695185309048285123705067072486920463781103076554014502567884803571416673251784936825115787932810954867447447568320403976197134736485611912650805539603318790667901618038578533362100071745480995207732506742832634459994375828162163700807237997808869771569154136465922798310222055287047244647419069003284481 elapsed time = 114 ms (factor = 63766529). F₁₂ = 114689 * 26017793 * 63766529 * (C1213)
PARI/GP
\\ Define a function to calculate Fermat numbers
fermat(n) = 2^(2^n) + 1;
\\ Define a function to factor Fermat numbers up to a given maximum
\\ and optionally just print them without factoring
factorfermats(mymax, nofactor=0) =
{
local(fm, factors, n);
for (n = 0, mymax,
fm = fermat(n);
if (nofactor,
print("Fermat number F(", n, ") is ", fm, ".");
next;
);
factors = factorint(fm);
if (matsize(factors)[1] == 1 && factors[1,2] == 1, \\ Check if it has only one factor with exponent 1 (i.e., prime)
print("Fermat number F(", n, "), ", fm, ", is prime.");
,
print("Fermat number F(", n, "), ", fm, ", factors to ", (factors), ".");
);
);
}
{
\\ Example usage
factorfermats(9, 1); \\ Print Fermat numbers without factoring
print(""); \\ Just to add a visual separator in the output
factorfermats(10); \\ Factor Fermat numbers
}
- Output:
Fermat number F(0) is 3. Fermat number F(1) is 5. Fermat number F(2) is 17. Fermat number F(3) is 257. Fermat number F(4) is 65537. Fermat number F(5) is 4294967297. Fermat number F(6) is 18446744073709551617. Fermat number F(7) is 340282366920938463463374607431768211457. Fermat number F(8) is 115792089237316195423570985008687907853269984665640564039457584007913129639937. Fermat number F(9) is 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097. Fermat number F(0), 3, is prime. Fermat number F(1), 5, is prime. Fermat number F(2), 17, is prime. Fermat number F(3), 257, is prime. Fermat number F(4), 65537, is prime. Fermat number F(5), 4294967297, factors to [641, 1; 6700417, 1]. Fermat number F(6), 18446744073709551617, factors to [274177, 1; 67280421310721, 1]. Fermat number F(7), 340282366920938463463374607431768211457, factors to [59649589127497217, 1; 5704689200685129054721, 1]. Fermat number F(8), 115792089237316195423570985008687907853269984665640564039457584007913129639937, factors to [1238926361552897, 1; 93461639715357977769163558199606896584051237541638188580280321, 1].
PascalABC.NET
##
function IsPrime(n: biginteger): boolean;
begin
if (n = 2) or (n = 3) then Result := true
else if (n <= 1) or ((n mod 2) = 0) or ((n mod 3) = 0) then Result := false
else
begin
var i := 5bi;
Result := False;
while i * i < n do
begin
if (n mod i) = 0 then
begin
println(i, n div i);
exit;
end
else if (n mod (i + 2)) = 0 then
begin
println(i + 2, n div (i + 2));
exit;
end;
i += 6;
end;
Result := True;
end;
end;
function fermat(n: integer) := power(2bi, int64(power(2, n))) + 1;
for var n := 0 to 9 do writeln('F', n, ' = ', fermat(n));
for var n := 0 to 6 do
begin
write('F', n, ' = ');
if isprime(fermat(n)) then println(fermat(n));
end;
- Output:
F0 = 3 F1 = 5 F2 = 17 F3 = 257 F4 = 65537 F5 = 4294967297 F6 = 18446744073709551617 F7 = 340282366920938463463374607431768211457 F8 = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F9 = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 F0 = 3 F1 = 5 F2 = 17 F3 = 257 F4 = 65537 F5 = 641 6700417 F6 = 274177 67280421310721
Perl
use strict;
use warnings;
use feature 'say';
use bigint try=>"GMP";
use ntheory qw<factor>;
my @Fermats = map { 2**(2**$_) + 1 } 0..9;
my $sub = 0;
say 'First 10 Fermat numbers:';
printf "F%s = %s\n", $sub++, $_ for @Fermats;
$sub = 0;
say "\nFactors of first few Fermat numbers:";
for my $f (map { [factor($_)] } @Fermats[0..8]) {
printf "Factors of F%s: %s\n", $sub++, @$f == 1 ? 'prime' : join ' ', @$f
}
- Output:
First 10 Fermat numbers: F0 = 3 F1 = 5 F2 = 17 F3 = 257 F4 = 65537 F5 = 4294967297 F6 = 18446744073709551617 F7 = 340282366920938463463374607431768211457 F8 = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F9 = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors of first few Fermat numbers: Factors of F0: prime Factors of F1: prime Factors of F2: prime Factors of F3: prime Factors of F4: prime Factors of F5: 641 6700417 Factors of F6: 274177 67280421310721 Factors of F7: 59649589127497217 5704689200685129054721 Factors of F8: 1238926361552897 93461639715357977769163558199606896584051237541638188580280321
Phix
with javascript_semantics -- demo\rosetta\Fermat.exw include mpfr.e procedure fermat(mpz res, integer n) integer pn = power(2,n) mpz_ui_pow_ui(res,2,pn) mpz_add_si(res,res,1) end procedure mpz fn = mpz_init() constant lim = iff(platform()=JS?18:29), -- (see note) print_lim = iff(platform()=JS?16:20) for i=0 to lim do fermat(fn,i) if i<=print_lim then printf(1,"F%d = %s\n",{i,shorten(mpz_get_str(fn))}) else -- (since printing it takes too long...) printf(1,"F%d has %,d digits\n",{i,mpz_sizeinbase(fn,10)}) end if end for printf(1,"\n") constant flimit = iff(platform()=JS?11:13) for i=0 to flimit do atom t = time() fermat(fn,i) sequence f = mpz_prime_factors(fn, 200000) t = time()-t string fs = "", ts = elapsed(t) if length(f[$])=1 then -- (as per docs) mpz_set_str(fn,f[$][1]) if not mpz_prime(fn) then if length(f)=1 then fs = " (not prime)" else fs = " (last factor is not prime)" end if end if f = deep_copy(f) f[$][1] = shorten(f[$][1]) elsif length(f)=1 and mpz_prime(fn) then fs = " (prime)" end if fs = mpz_factorstring(f)&fs printf(1,"Factors of F%d: %s [%s]\n",{i,fs,ts}) end for
- Output:
Note that mpz_prime_factors(), a phix-specific extension to gmp, is designed to find small factors quickly and give up early, however it works by maintaining a table of primes, so any prime factor over 10 digits or so is beyond reach. You could increase the maxprime parameter, here set at 200,000, which guarantees all factors up to 2,750,159 (obviously 7 digits), but it will just get exponentially slower without getting close to finding anything more, specifically in this case 1,238,926,361,552,897 (16 digits) or 59,649,589,127,497,217 (17 digits).
Calculating F0..F29 is pretty quick, but F30 and above hit integer limits on 32 bit, F32 and above exceed my physical memory on 64 bit.
As noted above, there is not really much point, and it just takes far too long to bother printing out any numbers with more than 500,000 digits.
Attempting to factor F14 and above gets nowhere, with each attempt some 5-10 times slower than the previous, until F18 which eventually crashes.
F0 = 3 F1 = 5 F2 = 17 F3 = 257 F4 = 65537 F5 = 4294967297 F6 = 18446744073709551617 F7 = 340282366920938463463374607431768211457 F8 = 115792089...<78 digits>...129639937 F9 = 134078079...<155 digits>...006084097 F10 = 179769313...<309 digits>...224137217 F11 = 323170060...<617 digits>...596230657 F12 = 104438888...<1,234 digits>...154190337 F13 = 109074813...<2,467 digits>...715792897 F14 = 118973149...<4,933 digits>...964066817 F15 = 141546103...<9,865 digits>...712377857 F16 = 200352993...<19,729 digits>...719156737 F17 = 401413218...<39,457 digits>...934173697 F18 = 161132571...<78,914 digits>...298300417 F19 = 259637056...<157,827 digits>...185773057 F20 = 674114012...<315,653 digits>...335579137 F21 has 631,306 digits F22 has 1,262,612 digits F23 has 2,525,223 digits F24 has 5,050,446 digits F25 has 10,100,891 digits F26 has 20,201,782 digits F27 has 40,403,563 digits F28 has 80,807,125 digits F29 has 161,614,249 digits Factors of F0: 3 (prime) [0.0s] Factors of F1: 5 (prime) [0s] Factors of F2: 17 (prime) [0s] Factors of F3: 257 (prime) [0s] Factors of F4: 65537 (prime) [0s] Factors of F5: 641*6700417 [0s] Factors of F6: 274177*67280421310721 [0.0s] Factors of F7: 340282366920938463463374607431768211457 (not prime) [0.2s] Factors of F8: 115792089...<78 digits>...129639937 (not prime) [0.2s] Factors of F9: 2424833*552937374...<148 digits>...393118209 (last factor is not prime) [0.2s] Factors of F10: 179769313...<309 digits>...224137217 (not prime) [0.2s] Factors of F11: 319489*974849*103761886...<606 digits>...591348737 (last factor is not prime) [0.3s] Factors of F12: 114689*910626896...<1,228 digits>...946770433 (last factor is not prime) [0.6s] Factors of F13: 109074813...<2,467 digits>...715792897 (not prime) [1.3s]
PicoLisp
(seed (in "/dev/urandom" (rd 8)))
(de **Mod (X Y N)
(let M 1
(loop
(when (bit? 1 Y)
(setq M (% (* M X) N)) )
(T (=0 (setq Y (>> 1 Y)))
M )
(setq X (% (* X X) N)) ) ) )
(de isprime (N)
(cache '(NIL) N
(if (== N 2)
T
(and
(> N 1)
(bit? 1 N)
(let (Q (dec N) N1 (dec N) K 0 X)
(until (bit? 1 Q)
(setq
Q (>> 1 Q)
K (inc K) ) )
(catch 'composite
(do 16
(loop
(setq X
(**Mod
(rand 2 (min (dec N) 1000000000000))
Q
N ) )
(T (or (=1 X) (= X N1)))
(T
(do K
(setq X (**Mod X 2 N))
(when (=1 X) (throw 'composite))
(T (= X N1) T) ) )
(throw 'composite) ) )
(throw 'composite T) ) ) ) ) ) )
(de gcd (A B)
(until (=0 B)
(let M (% A B)
(setq A B B M) ) )
(abs A) )
(de g (A)
(% (+ (% (* A A) N) C) N) )
(de pollard-brent (N)
(let
(A (dec N)
Y (rand 1 (min A 1000000000000000000))
C (rand 1 (min A 1000000000000000000))
M (rand 1 (min A 1000000000000000000))
G 1
R 1
Q 1 )
(ifn (bit? 1 N)
2
(loop
(NIL (=1 G))
(setq X Y)
(do R
(setq Y (g Y)) )
(zero K)
(loop
(NIL (and (> R K) (=1 G)))
(setq YS Y)
(do (min M (- R K))
(setq
Y (g Y)
Q (% (* Q (abs (- X Y))) N) ) )
(setq
G (gcd Q N)
K (+ K M) )
)
(setq R (* R 2)) )
(when (== G N)
(loop
(NIL (> G 1))
(setq
YS (g YS)
G (gcd (abs (- X YS)) N) ) ) )
(if (== G N)
NIL
G ) ) ) )
(de factors (N)
(sort
(make
(loop
(setq N (/ N (link (pollard-brent N))))
(T (isprime N)) )
(link N) ) ) )
(de fermat (N)
(inc (** 2 (** 2 N))) )
(for (N 0 (>= 8 N) (inc N))
(println N ': (fermat N)) )
(prinl)
(for (N 0 (>= 8 N) (inc N))
(let N (fermat N)
(println
N
':
(if (isprime N) 'PRIME (factors N)) ) ) )
- Output:
$ pil VU.l 0 : 3 1 : 5 2 : 17 3 : 257 4 : 65537 5 : 4294967297 6 : 18446744073709551617 7 : 340282366920938463463374607431768211457 8 : 115792089237316195423570985008687907853269984665640564039457584007913129639937 3 : PRIME 5 : PRIME 17 : PRIME 257 : PRIME 65537 : PRIME 4294967297 : (641 6700417) 18446744073709551617 : (274177 67280421310721) 340282366920938463463374607431768211457 : (59649589127497217 5704689200685129054721) 115792089237316195423570985008687907853269984665640564039457584007913129639937 : (1238926361552897 93461639715357977769163558199606896584051237541638188580280321)
Python
Procedural
def factors(x):
factors = []
i = 2
s = int(x ** 0.5)
while i < s:
if x % i == 0:
factors.append(i)
x = int(x / i)
s = int(x ** 0.5)
i += 1
factors.append(x)
return factors
print("First 10 Fermat numbers:")
for i in range(10):
fermat = 2 ** 2 ** i + 1
print("F{} = {}".format(chr(i + 0x2080) , fermat))
print("\nFactors of first few Fermat numbers:")
for i in range(10):
fermat = 2 ** 2 ** i + 1
fac = factors(fermat)
if len(fac) == 1:
print("F{} -> IS PRIME".format(chr(i + 0x2080)))
else:
print("F{} -> FACTORS: {}".format(chr(i + 0x2080), fac))
- Output:
First 10 Fermat numbers: F₀ = 3 F₁ = 5 F₂ = 17 F₃ = 257 F₄ = 65537 F₅ = 4294967297 F₆ = 18446744073709551617 F₇ = 340282366920938463463374607431768211457 F₈ = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F₉ = 1340780792994259709957402499820584612747936582059239337772356144372176403007354697680187429816690342769003185818648605085375388281194656994643364 9006084097 Factors of first few Fermat numbers: F₀ IS PRIME F₁ IS PRIME F₂ IS PRIME F₃ IS PRIME F₄ IS PRIME F₅ FACTORS: [641, 6700417] F₆ FACTORS: [274177, 67280421310721] F₇ FACTORS: [59649589127497217, 5704689200685129054721] F₈ FACTORS: [1238926361552897, 93461639715357977769163558199606896584051237541638188580280321]
Functional
'''Fermat numbers'''
from itertools import count, islice
from math import floor, sqrt
# fermat :: Int -> Int
def fermat(n):
'''Nth Fermat number.
Nth term of OEIS A000215.
'''
return 1 + (2 ** (2 ** n))
# fermats :: () -> [Int]
def fermats():
'''Non-finite series of Fermat numbers.
OEIS A000215.
'''
return (fermat(x) for x in enumFrom(0))
# --------------------------TEST---------------------------
# main :: IO ()
def main():
'''First 10 Fermats, and factors of first 7.'''
print(
fTable('First ten Fermat numbers:')(str)(str)(
fermat
)(enumFromTo(0)(9))
)
print(
fTable('\n\nFactors of first seven:')(str)(
lambda xs: repr(xs) if 1 < len(xs) else '(prime)'
)(primeFactors)(
take(7)(fermats())
)
)
# -------------------------DISPLAY-------------------------
# fTable :: String -> (a -> String) ->
# (b -> String) -> (a -> b) -> [a] -> String
def fTable(s):
'''Heading -> x display function -> fx display function ->
f -> xs -> tabular string.
'''
def go(xShow, fxShow, f, xs):
ys = [xShow(x) for x in xs]
w = max(map(len, ys))
return s + '\n' + '\n'.join(map(
lambda x, y: y.rjust(w, ' ') + ' -> ' + fxShow(f(x)),
xs, ys
))
return lambda xShow: lambda fxShow: lambda f: lambda xs: go(
xShow, fxShow, f, xs
)
# -------------------------GENERIC-------------------------
# enumFrom :: Enum a => a -> [a]
def enumFrom(x):
'''A non-finite stream of enumerable values,
starting from the given value.
'''
return count(x) if isinstance(x, int) else (
map(chr, count(ord(x)))
)
# enumFromTo :: Int -> Int -> [Int]
def enumFromTo(m):
'''Enumeration of integer values [m..n]'''
def go(n):
return list(range(m, 1 + n))
return lambda n: go(n)
# primeFactors :: Int -> [Int]
def primeFactors(n):
'''A list of the prime factors of n.
'''
def f(qr):
r = qr[1]
return step(r), 1 + r
def step(x):
return 1 + (x << 2) - ((x >> 1) << 1)
def go(x):
root = floor(sqrt(x))
def p(qr):
q = qr[0]
return root < q or 0 == (x % q)
q = until(p)(f)(
(2 if 0 == x % 2 else 3, 1)
)[0]
return [x] if q > root else [q] + go(x // q)
return go(n)
# take :: Int -> [a] -> [a]
# take :: Int -> String -> String
def take(n):
'''The prefix of xs of length n,
or xs itself if n > length xs.
'''
return lambda xs: (
xs[0:n]
if isinstance(xs, (list, tuple))
else list(islice(xs, n))
)
# until :: (a -> Bool) -> (a -> a) -> a -> a
def until(p):
'''The result of repeatedly applying f until p holds.
The initial seed value is x.
'''
def go(f, x):
v = x
while not p(v):
v = f(v)
return v
return lambda f: lambda x: go(f, x)
# MAIN ---
if __name__ == '__main__':
main()
- Output:
First ten Fermat numbers: 0 -> 3 1 -> 5 2 -> 17 3 -> 257 4 -> 65537 5 -> 4294967297 6 -> 18446744073709551617 7 -> 340282366920938463463374607431768211457 8 -> 115792089237316195423570985008687907853269984665640564039457584007913129639937 9 -> 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors of first seven: 3 -> (prime) 5 -> (prime) 17 -> (prime) 257 -> (prime) 65537 -> (prime) 4294967297 -> [641, 6700417] 18446744073709551617 -> [274177, 67280421310721]
Raku
(formerly Perl 6)
I gave up on factoring F₉ after about 20 minutes.
use ntheory:from<Perl5> <factor>;
my @Fermats = (^Inf).map: 2 ** 2 ** * + 1;
my $sub = '₀';
say "First 10 Fermat numbers:";
printf "F%s = %s\n", $sub++, $_ for @Fermats[^10];
$sub = '₀';
say "\nFactors of first few Fermat numbers:";
for @Fermats[^9].map( {"$_".&factor} ) -> $f {
printf "Factors of F%s: %s %s\n", $sub++, $f.join(' '), $f.elems == 1 ?? '- prime' !! ''
}
- Output:
First 10 Fermat numbers: F₀ = 3 F₁ = 5 F₂ = 17 F₃ = 257 F₄ = 65537 F₅ = 4294967297 F₆ = 18446744073709551617 F₇ = 340282366920938463463374607431768211457 F₈ = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F₉ = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors of first few Fermat numbers: Factors of F₀: 3 - prime Factors of F₁: 5 - prime Factors of F₂: 17 - prime Factors of F₃: 257 - prime Factors of F₄: 65537 - prime Factors of F₅: 641 6700417 Factors of F₆: 274177 67280421310721 Factors of F₇: 59649589127497217 5704689200685129054721 Factors of F₈: 1238926361552897 93461639715357977769163558199606896584051237541638188580280321
REXX
factoring by trial division
/*REXX program to find and display Fermat numbers, and show factors of Fermat numbers.*/
parse arg n . /*obtain optional argument from the CL.*/
if n=='' | n=="," then n= 9 /*Not specified? Then use the default.*/
numeric digits 20 /*ensure enough decimal digits, for n=9*/
do j=0 to n; f= 2** (2**j) + 1 /*calculate a series of Fermat numbers.*/
say right('F'j, length(n) + 1)': ' f /*display a particular " " */
end /*j*/
say
do k=0 to n; f= 2** (2**k) + 1; say /*calculate a series of Fermat numbers.*/
say center(' F'k": " f' ', 79, "═") /*display a particular " " */
p= factr(f) /*factor a Fermat number, given time. */
if words(p)==1 then say f ' is prime.'
else say 'factors: ' p
end /*k*/
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
factr: procedure; parse arg x 1 z,,?
do k=1 to 11 by 2; j= k; if j==1 then j= 2; if j==9 then iterate
call build /*add J to the factors list. */
end /*k*/ /* [↑] factor X with some low primes*/
do y=0 by 2; j= j + 2 + y // 4 /*ensure not ÷ by three. */
parse var j '' -1 _; if _==5 then iterate /*last digit a "5"? Skip it.*/
if j*j>x | j>z then leave
call build /*add Y to the factors list. */
end /*y*/ /* [↑] factor X with other higher #s*/
j= z
if z\==1 then ?= build()
if ?='' then do; @.1= x; ?= x; #= 1; end
return ?
/*──────────────────────────────────────────────────────────────────────────────────────*/
build: do while z//j==0; z= z % j; ?= ? j; end; return strip(?)
- output when using the default input:
F0: 3 F1: 5 F2: 17 F3: 257 F4: 65537 F5: 4294967297 F6: 18446744073709551617 F7: 340282366920938463463374607431768211457 F8: 115792089237316195423570985008687907853269984665640564039457584007913129639937 F9: 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 ═══════════════════════════════════ F0: 3 ════════════════════════════════════ 3 is prime. ═══════════════════════════════════ F1: 5 ════════════════════════════════════ 5 is prime. ═══════════════════════════════════ F2: 17 ═══════════════════════════════════ 17 is prime. ══════════════════════════════════ F3: 257 ═══════════════════════════════════ 257 is prime. ═════════════════════════════════ F4: 65537 ══════════════════════════════════ 65537 is prime. ═══════════════════════════════ F5: 4294967297 ═══════════════════════════════ factors: 641 6700417 ══════════════════════════ F6: 18446744073709551617 ══════════════════════════ ■ ■ ■ (the REXX program stopped via Ctrl─Alt─Break) ■ ■ ■
factoring via Pollard's rho algorithm
/*REXX program to find and display Fermat numbers, and show factors of Fermat numbers.*/
parse arg n . /*obtain optional argument from the CL.*/
if n=='' | n=="," then n= 9 /*Not specified? Then use the default.*/
numeric digits 200 /*ensure enough decimal digits, for n=9*/
do j=0 to n; f= 2** (2**j) + 1 /*calculate a series of Fermat numbers.*/
say right('F'j, length(n) + 1)': ' f /*display a particular " " */
end /*j*/
say
do k=5 to n; f= 2** (2**k) + 1; say /*calculate a series of Fermat numbers.*/
say center(' F'k": " f' ', 79, "═") /*display a particular " " */
a= rho(f) /*factor a Fermat number, given time. */
b= f % a
if a==b then say f ' is prime.'
else say 'factors: ' commas(a) " " commas(b)
end /*k*/
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _; do ?=length(_)-3 to 1 by -3; _=insert(',', _, ?); end; return _
/*──────────────────────────────────────────────────────────────────────────────────────*/
rho: procedure; parse arg n; y= 2; d= 1 /*initialize X, Y, and D variables.*/
do x=2 until d==n /*try rho method with X=2 for 1st time.*/
do while d==1
x= (x*x + 1) // n
v= (y*y + 1) // n
y= (v*v + 1) // n
parse value x-y with xy 1 sig 2 /*obtain sign of the x-y difference. */
if sig=='-' then parse var xy 2 xy /*Negative? Then use absolute value. */
nn= n
do until nn==0
parse value xy//nn nn with nn xy /*assign two variables: NN and XY */
end /*until*/ /*this is an in-line GCD function. */
d= xy /*assign variable D with a new XY */
end /*while*/
end /*x*/
return d /*found a factor of N. Return it.*/
- output when using the default input:
F0: 3 F1: 5 F2: 17 F3: 257 F4: 65537 F5: 4294967297 F6: 18446744073709551617 F7: 340282366920938463463374607431768211457 F8: 115792089237316195423570985008687907853269984665640564039457584007913129639937 F9: 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 ═══════════════════════════════ F5: 4294967297 ═══════════════════════════════ factors: 641 6,700,417 ══════════════════════════ F6: 18446744073709551617 ══════════════════════════ factors: 274,177 67,280,421,310,721 ════════════════ F7: 340282366920938463463374607431768211457 ═════════════════ ■ ■ ■ (the REXX program stopped via Ctrl─Alt─Break) ■ ■ ■
Using REXX libraries
Libraries: How to use
Libraries: Source code
include Settings
numeric digits 200
say 'Fermat numbers - Using REXX libraries'
parse version version; say version; say
do i = 0 to 9
p = 2**(2**i)+1; l = Length(p)
say 'Fermat number' i',' p',' l 'digits'
if i < 7 then do
f = Factors(p)
if f = 1 then
say 'is prime'
else do
call Charout ,'has factors '
do j = 1 to f
call Charout, fact.factor.j' '
end
say
end
end
else do
if IsPrime(p) then
say 'is prime'
else
say 'is composite, but could not be factorized'
end
say
end
say Format(Time('e'),,3) 'seconds'; say
exit
include Functions
include Numbers
include Abend
Procedure Factors() is in Sequences and Prime() in Numbers.
- Output:
REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022 Fermat numbers Fermat number 0, 3, 1 digits is prime Fermat number 1, 5, 1 digits is prime Fermat number 2, 17, 2 digits is prime Fermat number 3, 257, 3 digits is prime Fermat number 4, 65537, 5 digits is prime Fermat number 5, 4294967297, 10 digits has factors 641 6700417 Fermat number 6, 18446744073709551617, 20 digits has factors 274177 67280421310721 Fermat number 7, 340282366920938463463374607431768211457, 39 digits is composite, but could not be factorized Fermat number 8, 115792089237316195423570985008687907853269984665640564039457584007913129639937, 78 digits is composite, but could not be factorized Fermat number 9, 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097, 155 digits is composite, but could not be factorized
Because Fermat number 7 has a smallest factor of 59649589127497217, this number and the higher ones take far too long to factorize in REXX, even when using the PollardRho method.
Ring
decimals(0)
load "stdlib.ring"
see "working..." + nl
see "The first 10 Fermat numbers are:" + nl
num = 0
limit = 9
for n = 0 to limit
fermat = pow(2,pow(2,n)) + 1
mod = fermat%2
if n > 5
ferm = string(fermat)
tmp = number(right(ferm,1))+1
fermat = left(ferm,len(ferm)-1) + string(tmp)
ok
see "F(" + n + ") = " + fermat + nl
next
see "done..." + nl
Output:
working... The first 10 Fermat numbers are: F(0) = 3 F(1) = 5 F(2) = 17 F(3) = 257 F(4) = 65537 F(5) = 4294967297 F(6) = 18446744073709551617 F(7) = 340282366920938463463374607431768211457 F(8) = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F(9) = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097
RPL
« 2 DUP ROT ^ ^ 1 + » 'FERMAT' STO « n FERMAT » 'n' 0 9 1 SEQ « { } 0 9 FOR n n FERMAT IF DUP ISPRIME? THEN + ELSE DROP END NEXT » TEVAL
- Output:
3: { 3 5 17 257 65537 4294967297 18446744073709551617 340282366920938463463374607431768211457 115792089237316195423570985008687907853269984665640564039457584007913129639937 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 } 2: { 3 5 17 257 65537 } 1: s:273.0092
It takes 4 min 33 s to test primality of the first 10 Fermat numbers on a HP-50g.
Ruby
This uses the `factor` function from the `coreutils` library that comes standard with most GNU/Linux, BSD, and Unix systems. https://www.gnu.org/software/coreutils/ https://en.wikipedia.org/wiki/GNU_Core_Utilities
def factors(n)
factors = `factor #{n}`.split(' ')[1..-1].map(&:to_i)
factors.group_by { _1 }.map { |prime, exp| [prime, exp.size] } # Ruby 2.7 or later
#factors.group_by { |prime| prime }.map { |prime, exp| [prime, exp.size] } # for all versions
end
def fermat(n); (1 << (1 << n)) | 1 end
puts "Value for each Fermat Number F0 .. F9."
(0..9).each { |n| puts "F#{n} = #{fermat(n)}" }
puts
puts "Factors for each Fermat Number F0 .. F8."
(0..8).each { |n| puts "F#{n} = #{factors fermat(n)}" }
System: Lenovo V570 (2011), I5-2410M, 2.9 GHz, Ruby 2.7.1 Run as: $ time ruby fermat.rb
- Output:
Value for each Fermat Number F0 .. F9. F0 = 3 F1 = 5 F2 = 17 F3 = 257 F4 = 65537 F5 = 4294967297 F6 = 18446744073709551617 F7 = 340282366920938463463374607431768211457 F8 = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F9 = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors for each Fermat Number F0 .. F8. F0 = [[3, 1]] F1 = [[5, 1]] F2 = [[17, 1]] F3 = [[257, 1]] F4 = [[65537, 1]] F5 = [[641, 1], [6700417, 1]] F6 = [[274177, 1], [67280421310721, 1]] F7 = [[59649589127497217, 1], [5704689200685129054721, 1]] F8 = [[1238926361552897, 1], [93461639715357977769163558199606896584051237541638188580280321, 1]] ruby fermat.rb 175.26s user 0.01s system 99% cpu 2:55.27 total
Rust
struct DivisorGen {
curr: u64,
last: u64,
}
impl Iterator for DivisorGen {
type Item = u64;
fn next(&mut self) -> Option<u64> {
self.curr += 2u64;
if self.curr < self.last{
None
} else {
Some(self.curr)
}
}
}
fn divisor_gen(num : u64) -> DivisorGen {
DivisorGen { curr: 0u64, last: (num / 2u64) + 1u64 }
}
fn is_prime(num : u64) -> bool{
if num == 2 || num == 3 {
return true;
} else if num % 2 == 0 || num % 3 == 0 || num <= 1{
return false;
}else{
for i in divisor_gen(num){
if num % i == 0{
return false;
}
}
}
return true;
}
fn main() {
let fermat_closure = |i : u32| -> u64 {2u64.pow(2u32.pow(i + 1u32))};
let mut f_numbers : Vec<u64> = Vec::new();
println!("First 4 Fermat numbers:");
for i in 0..4 {
let f = fermat_closure(i) + 1u64;
f_numbers.push(f);
println!("F{}: {}", i, f);
}
println!("Factor of the first four numbers:");
for f in f_numbers.iter(){
let is_prime : bool = f % 4 == 1 && is_prime(*f);
let not_or_not = if is_prime {" "} else {" not "};
println!("{} is{}prime", f, not_or_not);
}
}
- Output:
First 4 Fermat numbers: F0: 5 F1: 17 F2: 257 F3: 65537 Factor of the first four numbers: 5 is prime 17 is prime 257 is prime 65537 is prime
Alternative using rug crate
Based on the C++ code, which was based on the Java solution.
// [dependencies]
// rug = "1.9"
use rug::Integer;
fn fermat(n: u32) -> Integer {
Integer::from(Integer::u_pow_u(2, 2u32.pow(n))) + 1
}
fn g(x: Integer, n: &Integer) -> Integer {
(Integer::from(&x * &x) + 1) % n
}
fn pollard_rho(n: &Integer) -> Integer {
use rug::Assign;
let mut x = Integer::from(2);
let mut y = Integer::from(2);
let mut d = Integer::from(1);
let mut z = Integer::from(1);
let mut count = 0;
loop {
x = g(x, n);
y = g(g(y, n), n);
d.assign(&x - &y);
d = d.abs();
z *= &d;
z %= n;
count += 1;
if count == 100 {
d.assign(z.gcd_ref(n));
if d != 1 {
break;
}
z.assign(1);
count = 0;
}
}
if d == *n {
return Integer::from(0);
}
d
}
fn get_prime_factors(n: &Integer) -> Vec<Integer> {
use rug::integer::IsPrime;
let mut factors = Vec::new();
let mut m = Integer::from(n);
loop {
if m.is_probably_prime(25) != IsPrime::No {
factors.push(m);
break;
}
let f = pollard_rho(&m);
if f == 0 {
factors.push(m);
break;
}
factors.push(Integer::from(&f));
m = m / f;
}
factors
}
fn main() {
for i in 0..10 {
println!("F({}) = {}", i, fermat(i));
}
println!("\nPrime factors:");
for i in 0..9 {
let f = get_prime_factors(&fermat(i));
print!("F({}): {}", i, f[0]);
for j in 1..f.len() {
print!(", {}", f[j]);
}
println!();
}
}
- Output:
Execution time is about 8.5 minutes on my system (macOS 10.15, 3.2 GHz Quad-Core Intel Core i5).
F(0) = 3 F(1) = 5 F(2) = 17 F(3) = 257 F(4) = 65537 F(5) = 4294967297 F(6) = 18446744073709551617 F(7) = 340282366920938463463374607431768211457 F(8) = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F(9) = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Prime factors: F(0): 3 F(1): 5 F(2): 17 F(3): 257 F(4): 65537 F(5): 641, 6700417 F(6): 274177, 67280421310721 F(7): 59649589127497217, 5704689200685129054721 F(8): 1238926361552897, 93461639715357977769163558199606896584051237541638188580280321
Scala
import scala.collection.mutable
import scala.collection.mutable.ListBuffer
object FermatNumbers {
def main(args: Array[String]): Unit = {
println("First 10 Fermat numbers:")
for (i <- 0 to 9) {
println(f"F[$i] = ${fermat(i)}")
}
println()
println("First 12 Fermat numbers factored:")
for (i <- 0 to 12) {
println(f"F[$i] = ${getString(getFactors(i, fermat(i)))}")
}
}
private val TWO = BigInt(2)
def fermat(n: Int): BigInt = {
TWO.pow(math.pow(2.0, n).intValue()) + 1
}
def getString(factors: List[BigInt]): String = {
if (factors.size == 1) {
return s"${factors.head} (PRIME)"
}
factors.map(a => a.toString)
.map(a => if (a.startsWith("-")) "(C" + a.replace("-", "") + ")" else a)
.reduce((a, b) => a + " * " + b)
}
val COMPOSITE: mutable.Map[Int, String] = scala.collection.mutable.Map(
9 -> "5529",
10 -> "6078",
11 -> "1037",
12 -> "5488",
13 -> "2884"
)
def getFactors(fermatIndex: Int, n: BigInt): List[BigInt] = {
var n2 = n
var factors = new ListBuffer[BigInt]
var loop = true
while (loop) {
if (n2.isProbablePrime(100)) {
factors += n2
loop = false
} else {
if (COMPOSITE.contains(fermatIndex)) {
val stop = COMPOSITE(fermatIndex)
if (n2.toString.startsWith(stop)) {
factors += -n2.toString().length
loop = false
}
}
if (loop) {
val factor = pollardRhoFast(n2)
if (factor == 0) {
factors += n2
loop = false
} else {
factors += factor
n2 = n2 / factor
}
}
}
}
factors.toList
}
def pollardRhoFast(n: BigInt): BigInt = {
var x = BigInt(2)
var y = BigInt(2)
var z = BigInt(1)
var d = BigInt(1)
var count = 0
var loop = true
while (loop) {
x = pollardRhoG(x, n)
y = pollardRhoG(pollardRhoG(y, n), n)
d = (x - y).abs
z = (z * d) % n
count += 1
if (count == 100) {
d = z.gcd(n)
if (d != 1) {
loop = false
} else {
z = BigInt(1)
count = 0
}
}
}
println(s" Pollard rho try factor $n")
if (d == n) {
return 0
}
d
}
def pollardRhoG(x: BigInt, n: BigInt): BigInt = ((x * x) + 1) % n
}
- Output:
First 10 Fermat numbers: F[0] = 3 F[1] = 5 F[2] = 17 F[3] = 257 F[4] = 65537 F[5] = 4294967297 F[6] = 18446744073709551617 F[7] = 340282366920938463463374607431768211457 F[8] = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F[9] = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 First 12 Fermat numbers factored: F[0] = 3 (PRIME) F[1] = 5 (PRIME) F[2] = 17 (PRIME) F[3] = 257 (PRIME) F[4] = 65537 (PRIME) Pollard rho try factor 4294967297 F[5] = 641 * 6700417 Pollard rho try factor 18446744073709551617 F[6] = 274177 * 67280421310721 Pollard rho try factor 340282366920938463463374607431768211457 F[7] = 59649589127497217 * 5704689200685129054721 Pollard rho try factor 115792089237316195423570985008687907853269984665640564039457584007913129639937 F[8] = 1238926361552897 * 93461639715357977769163558199606896584051237541638188580280321 Pollard rho try factor 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 F[9] = 2424833 * (C148) Pollard rho try factor 179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137217 Pollard rho try factor 3942951359960012586542991835686376608231592127249807732373409846031135195659174148737161255930050543559319182152642816343958573976075461198274610155058226350701077796608546283231637018483208223116080561800334422176622099740983337736621316898600121619871377542107047343253864459964167331555646795960321 F[10] = 45592577 * 6487031809 * (C291) Pollard rho try factor 32317006071311007300714876688669951960444102669715484032130345427524655138867890893197201411522913463688717960921898019494119559150490921095088152386448283120630877367300996091750197750389652106796057638384067568276792218642619756161838094338476170470581645852036305042887575891541065808607552399123930385521914333389668342420684974786564569494856176035326322058077805659331026192708460314150258592864177116725943603718461857357598351152301645904403697613233287231227125684710820209725157101726931323469678542580656697935045997268352998638215525166389437335543602135433229604645318478604952148193555853611059596230657 Pollard rho try factor 33150781373639412155846573868024639672856106606987835072026893834352453701925006737655987144186344206834820532125383540932102878651453631377873037143648178457002958783669056532601662155256508553423204658756451069116132055982639112479817996775373591674794399801442382402697828988429044712163168243619196804348072710121945390948428910309765481110260333687910970886853046635254307274981520537180895290310783635953818082306553996934497908037349010876970379631341148456045116407475229712217130141926525362871253437794629422541384355185626695660779797862427347553871011957167960991543632376506466281643163047416635393 F[11] = 974849 * 319489 * (C606) Pollard rho try factor 1044388881413152506691752710716624382579964249047383780384233483283953907971557456848826811934997558340890106714439262837987573438185793607263236087851365277945956976543709998340361590134383718314428070011855946226376318839397712745672334684344586617496807908705803704071284048740118609114467977783598029006686938976881787785946905630190260940599579453432823469303026696443059025015972399867714215541693835559885291486318237914434496734087811872639496475100189041349008417061675093668333850551032972088269550769983616369411933015213796825837188091833656751221318492846368125550225998300412344784862595674492194617023806505913245610825731835380087608622102834270197698202313169017678006675195485079921636419370285375124784014907159135459982790513399611551794271106831134090584272884279791554849782954323534517065223269061394905987693002122963395687782878948440616007412945674919823050571642377154816321380631045902916136926708342856440730447899971901781465763473223850267253059899795996090799469201774624817718449867455659250178329070473119433165550807568221846571746373296884912819520317457002440926616910874148385078411929804522981857338977648103126085903001302413467189726673216491511131602920781738033436090243804708340403154190337 Pollard rho try factor 9106268965752186405773463110818163752233991481723476361152625650968740750826648212547208641935996986118024454955917854702609434541985662158212523327009262247869952450049350838706079834460006786304075107567909269645531121898331250125751682239313156601738683820643686003638396435055834553570682260579462973839574318172464558815116581626749391315641251152532705571615644886981829338611134458123396450764186936496833100701185274214915961723337127995182593580031119299575446791424418154036863609858251201843852076223383379133238000289598800458955855329052103961332983048473420515918928565951506637819342706575976725030506905683310915700945442329953941604008255667676914945655757474715779252371155778495946746587469464160684843488975918662295274965457887082037460184558511575570318625886351712499453155527762335682281851520733417380809781252979478377941937578568481859702438295520231435016188495646093490407803983345420364088331996467459309353537828143080691834120737157445502646809195267166779721413577366833939771467773331873590129210913628329073978766992198221682739812652450408607796042492802295258713711959073218748776359806123717024800451461326745599716651128725627280643537507664130920416107218492950792456907321580171946770433 Pollard rho try factor 350001591824186871106763863899530746217943677302816436473017740242946077356624684213115564488348300185108877411543625345263121839042445381828217455916005721464151569276047005167043946981206545317048033534893189043572263100806868980998952610596646556521480658280419327021257968923645235918768446677220584153297488306270426504473941414890547838382804073832577020334339845312084285496895699728389585187497914849919557000902623608963141559444997044721763816786117073787751589515083702681348245404913906488680729999788351730419178916889637812821243344085799435845038164784900107242721493170135785069061880328434342030106354742821995535937481028077744396075726164309052460585559946407282864168038994640934638329525854255227752926198464207255983432778402344965903148839661825873175277621985527846249416909718758069997783882773355041329208102046913755441975327368023946523920699020098723785533557579080342841062805878477869513695185309048285123705067072486920463781103076554014502567884803571416673251784936825115787932810954867447447568320403976197134736485611912650805539603318790667901618038578533362100071745480995207732506742832634459994375828162163700807237997808869771569154136465922798310222055287047244647419069003284481 F[12] = 114689 * 26017793 * 63766529 * (C1213)
Sidef
func fermat_number(n) {
2**(2**n) + 1
}
func fermat_one_factor(n) {
fermat_number(n).ecm_factor
}
for n in (0..9) {
say "F_#{n} = #{fermat_number(n)}"
}
say ''
for n in (0..13) {
var f = fermat_one_factor(n)
say ("F_#{n} = ", join(' * ', f.shift,
f.map { <C P>[.is_prime] + .len }...))
}
- Output:
F_0 = 3 F_1 = 5 F_2 = 17 F_3 = 257 F_4 = 65537 F_5 = 4294967297 F_6 = 18446744073709551617 F_7 = 340282366920938463463374607431768211457 F_8 = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F_9 = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 F_0 = 3 F_1 = 5 F_2 = 17 F_3 = 257 F_4 = 65537 F_5 = 641 * P7 F_6 = 274177 * P14 F_7 = 59649589127497217 * P22 F_8 = 1238926361552897 * P62 F_9 = 2424833 * C148 F_10 = 45592577 * C301 F_11 = 319489 * C612 F_12 = 114689 * C1228 F_13 = 2710954639361 * C2454
Tcl
namespace import ::tcl::mathop::*
package require math::numtheory 1.1.1; # Buggy before tcllib-1.20
proc fermat n {
+ [** 2 [** 2 $n]] 1
}
for {set i 0} {$i < 10} {incr i} {
puts "F$i = [fermat $i]"
}
for {set i 1} {1} {incr i} {
puts -nonewline "F$i... "
flush stdout
set F [fermat $i]
set factors [math::numtheory::primeFactors $F]
if {[llength $factors] == 1} {
puts "is prime"
} else {
puts "factors: $factors"
}
}
- Output:
$ tclsh fermat.tcl F0 = 3 F1 = 5 F2 = 17 F3 = 257 F4 = 65537 F5 = 4294967297 F6 = 18446744073709551617 F7 = 340282366920938463463374607431768211457 F8 = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F9 = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 F1... is prime F2... is prime F3... is prime F4... is prime F5... factors: 641 6700417 F6...
Wren
Wren-cli
The first seven Fermat numbers are factorized in about 0.75 seconds but as it would take 'hours' to get any further I've stopped there.
import "./big" for BigInt
var fermat = Fn.new { |n| BigInt.two.pow(2.pow(n)) + 1 }
var fns = List.filled(10, null)
System.print("The first 10 Fermat numbers are:")
for (i in 0..9) {
fns[i] = fermat.call(i)
System.print("F%(String.fromCodePoint(0x2080+i)) = %(fns[i])")
}
System.print("\nFactors of the first 7 Fermat numbers:")
for (i in 0..6) {
System.write("F%(String.fromCodePoint(0x2080+i)) = ")
var factors = BigInt.primeFactors(fns[i])
System.write("%(factors)")
if (factors.count == 1) {
System.print(" (prime)")
} else if (!factors[1].isProbablePrime(5)) {
System.print(" (second factor is composite)")
} else {
System.print()
}
}
- Output:
The first 10 Fermat numbers are: F₀ = 3 F₁ = 5 F₂ = 17 F₃ = 257 F₄ = 65537 F₅ = 4294967297 F₆ = 18446744073709551617 F₇ = 340282366920938463463374607431768211457 F₈ = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F₉ = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors of the first 7 Fermat numbers: F₀ = [3] (prime) F₁ = [5] (prime) F₂ = [17] (prime) F₃ = [257] (prime) F₄ = [65537] (prime) F₅ = [641, 6700417] F₆ = [274177, 67280421310721]
Embedded
This uses Lenstra's elliptical curve method (ECM) to find F₇ and Pollard's Rho to find the rest, including the first factor of F₉.
The overall time for this particular run was about 86 seconds with F₇ being found in only a couple of them. However, being non-deterministic, ECM is quite erratic time-wise and could take up to 3 minutes longer.
/* fermat_numbers_2.wren */
import "./gmp" for Mpz
import "./ecm" for Ecm
import "random" for Random
var fermat = Fn.new { |n| Mpz.two.pow(2.pow(n)) + 1 }
var fns = List.filled(10, null)
System.print("The first 10 Fermat numbers are:")
for (i in 0..9) {
fns[i] = fermat.call(i)
System.print("F%(String.fromCodePoint(0x2080+i)) = %(fns[i])")
}
System.print("\nFactors of the first 8 Fermat numbers:")
for (i in 0..8) {
System.write("F%(String.fromCodePoint(0x2080+i)) = ")
var factors = (i != 7) ? Mpz.primeFactors(fns[i]) : Ecm.primeFactors(fns[i])
System.write("%(factors)")
if (factors.count == 1) {
System.print(" (prime)")
} else if (!factors[1].probPrime(15)) {
System.print(" (second factor is composite)")
} else {
System.print()
}
}
System.print("\nThe first factor of F₉ is %(Mpz.pollardRho(fns[9])).")
- Output:
The first 10 Fermat numbers are: F₀ = 3 F₁ = 5 F₂ = 17 F₃ = 257 F₄ = 65537 F₅ = 4294967297 F₆ = 18446744073709551617 F₇ = 340282366920938463463374607431768211457 F₈ = 115792089237316195423570985008687907853269984665640564039457584007913129639937 F₉ = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097 Factors of the first 8 Fermat numbers: F₀ = [3] (prime) F₁ = [5] (prime) F₂ = [17] (prime) F₃ = [257] (prime) F₄ = [65537] (prime) F₅ = [641, 6700417] F₆ = [274177, 67280421310721] F₇ = [59649589127497217, 5704689200685129054721] F₈ = [1238926361552897, 93461639715357977769163558199606896584051237541638188580280321] The first factor of F₉ is 2424833.
zkl
GNU Multiple Precision Arithmetic Library
for big ints and primes
fermatsW:=[0..].tweak(fcn(n){ BI(2).pow(BI(2).pow(n)) + 1 });
println("First 10 Fermat numbers:");
foreach n in (10){ println("F",n,": ",fermatsW.next()) }
- Output:
First 10 Fermat numbers: F0: 3 F1: 5 F2: 17 F3: 257 F4: 65537 F5: 4294967297 F6: 18446744073709551617 F7: 340282366920938463463374607431768211457 F8: 115792089237316195423570985008687907853269984665640564039457584007913129639937 F9: 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097
fcn primeFactorsBI(n){ // Return a list of the prime factors of n
acc:=fcn(n,k,acc,maxD){ // k is primes
if(n==1 or k>maxD) acc.close();
else{
q,r:=n.div2(k); // divr-->(quotient,remainder)
if(r==0) return(self.fcn(q,k,acc.write(k.copy()),q.root(2)));
return(self.fcn(n, k.nextPrime(), acc,maxD)) # both are tail recursion
}
}(n,BI(2),Sink(List),n.root(2));
m:=acc.reduce('*,BI(1)); // mulitply factors
if(n!=m) acc.append(n/m); // opps, missed last factor
else acc;
}
fermatsW:=[0..].tweak(fcn(n){ BI(2).pow(BI(2).pow(n)) + 1 });
println("Factors of first few Fermat numbers:");
foreach n in (7){
println("Factors of F",n,": ",factorsBI(fermatsW.next()).concat(" "));
}
- Output:
Factors of first few Fermat numbers: Factors of F0: 3 Factors of F1: 5 Factors of F2: 17 Factors of F3: 257 Factors of F4: 65537 Factors of F5: 641 6700417 Factors of F6: 274177 67280421310721