Rosetta Code

Fermat numbers

Revision as of 13:16, 29 May 2020 by PureFox (talk | contribs) (→‎{{header|Wren}}: Now uses Wren-math module.)

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.

Task
Fermat numbers
You are encouraged to solve this task according to the task description, using any language you may know.

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


Arturo

<lang arturo>nPowers: #(1 2 4 8 16 32 64 128 256 512) fermatSet: map 0..9 -> 2^nPowers.[&]+1

loop 0..9 { print "F(" + & + ") = " + fermatSet.[&] }

print ""

loop 0..9 { print "prime factors of F(" + & + ") = " + [primeFactors fermatSet.[&]] }</lang>

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

prime factors of F(0) = #(3)
prime factors of F(1) = #(5)
prime factors of F(2) = #(17)
prime factors of F(3) = #(257)
prime factors of F(4) = #(65537)
prime factors of F(5) = #(641 6700417)
prime factors of F(6) = #(274177 67280421310721)
prime factors of F(7) = #(59649589127497217 5704689200685129054721)
prime factors of F(8) = #(1238926361552897 93461639715357977769163558199606896584051237541638188580280321)
prime factors of F(9) = #(2424833)

C

Compile with : 

gcc -o fermat fermat.c -lgmp

<lang c>#include <stdlib.h>

  1. include <stdio.h>
  2. 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;

}</lang> 

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++

Translation of: Java
Library: Boost

Built and tested on macOS 10.15.2, CPU: 3.2 GHz Intel Core i5. Execution time is about 18 minutes. <lang cpp>#include <iostream>

  1. include <vector>
  2. include <boost/integer/common_factor.hpp>
  3. include <boost/multiprecision/cpp_int.hpp>
  4. 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 integer g(integer x, integer n) { 

   return (x * x + 1) % n;

}

integer pollard_rho(integer n) { 

   integer x = 2, y = 2, d = 1, z = 1;
   int count = 0;
   for (;;)
   {
       x = g(x, n);
       y = g(g(y, n), 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;

}</lang>

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

Crystal

Translation of: 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

<lang ruby>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)}"</lang> 

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

<lang 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</lang>

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. <lang go>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))

}</lang>

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

<lang 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
       ((bool "(prime)" . show) <*> ((1 <) . length))
       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)</lang>
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

<lang 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);
   }

} </lang> 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)

Julia

<lang 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)

</lang>

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

langur

Translation of: Python
Works with: langur version 0.10

<lang langur>val .fermat = f 2 ^ 2 ^ .n + 1

val .factors = f(var .x) { 

   for[.f=[]] .i, .s = 2, truncate .x ^/ 2; .i < .s; .i += 1 {
       if .x div .i {
           .f ~= [.i]
           .x \= .i
           .s = truncate .x ^/ 2
       }
   } ~ [.x]

}

writeln "first 10 Fermat numbers" for .i in 0..9 { 

   writeln $"F\(.i + 16x2080: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 + 16x2080:cp) is prime"
   } else {
       writeln $"F\(.i + 16x2080:cp) factors: ", .fac
   }

} </lang>

Output:

I just ran an initial test. Maybe I'll take the time to calculate more factors later. 

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]

Perl

Library: ntheory
Translation of: Raku

<lang 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

}</lang>

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

Library: Phix/mpfr

<lang Phix>-- 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() for i=0 to 29 do -- (see note) 

   fermat(fn,i)
   if i<=20 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") randstate state = gmp_randinit_mt() for i=0 to 13 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_probable_prime_p(fn, state) then
           if length(f)=1 then
               fs = " (not prime)"
           else
               fs = " (last factor is not prime)"
           end if
       end if
       f[$][1] = shorten(f[$][1])
   elsif length(f)=1
     and mpz_probable_prime_p(fn, state) then
       fs = " (prime)"
   end if
   fs = mpz_factorstring(f)&fs
   printf(1,"Factors of F%d: %s [%s]\n",{i,fs,ts})

end for</lang>

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]

Python

Procedural

<lang python>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))</lang>
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

Works with: Python version 3.7

<lang python>Fermat numbers

from itertools import count, islice from math import floor, sqrt


  1. fermat :: Int -> Int

def fermat(n): 

   Nth Fermat number.
      Nth term of OEIS A000215.
   
   return 1 + (2 ** (2 ** n))


  1. fermats :: () -> [Int]

def fermats(): 

   Non-finite series of Fermat numbers.
      OEIS A000215.
   
   return (fermat(x) for x in enumFrom(0))


  1. --------------------------TEST---------------------------
  2. 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())
       )
   )


  1. -------------------------DISPLAY-------------------------
  1. fTable :: String -> (a -> String) ->
  2. (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
   )


  1. -------------------------GENERIC-------------------------
  1. 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)))
   )


  1. 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)


  1. 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)


  1. take :: Int -> [a] -> [a]
  2. 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))
   )


  1. 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)


  1. MAIN ---

if __name__ == '__main__': 

   main()</lang>
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)

Works with: Rakudo version 2019.07.1

I gave up on factoring F₉ after about 20 minutes.

<lang perl6>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' !!

}</lang>

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

<lang rexx>/*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=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 /*──────────────────────────────────────────────────────────────────────────────────────*/ 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   /*forever*/
       return strip(?)</lang>
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

<lang rexx>/*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 /*──────────────────────────────────────────────────────────────────────────────────────*/ 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                                   /*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*/
       if d==n  then iterate                    /*Current X failure; bump X; try again.*/
       return d                                 /*found a factor of  N.      Return it.*/
       end     /*x*/</lang>
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)   ■  ■  ■

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

<lang ruby>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)}" }</lang> 

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

<lang 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);
   }

}</lang>

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

Sidef

<lang ruby>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 }...))

}</lang>

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

Wren

Library: Wren-math

In the absence of 'big integer' support, limited to computing and factorizing the first six Fermat numbers. <lang ecmascript>import "/math" for Int

var fermat = Fn.new { |n| 2.pow(2.pow(n)) + 1 }

for (i in 0..5) { 

   var n = fermat.call(i)
   var pf = Int.primeFactors(n)
   var kind = (pf.count == 1) ? "prime" : "composite"
   System.print("F%(String.fromCodePoint(0x2080+i)): %(n) -> factors = %(pf) -> %(kind)")

}</lang>

Output:
F₀: 3 -> factors = [3] -> prime
F₁: 5 -> factors = [5] -> prime
F₂: 17 -> factors = [17] -> prime
F₃: 257 -> factors = [257] -> prime
F₄: 65537 -> factors = [65537] -> prime
F₅: 4294967297 -> factors = [641, 6700417] -> composite

zkl

Library: GMP

GNU Multiple Precision Arithmetic Library

for big ints and primes <lang zkl>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()) }</lang>

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

<lang zkl>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;

}</lang> <lang zkl>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(" "));

}</lang>

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
Retrieved from "https://rosettacode.org/wiki/Fermat_numbers?oldid=292925"